Article Text
PDF
Abstract
Background The existence of immunologically ‘cold tumors’ frequently found across a wide spectrum of tumor types represents a significant challenge for cancer immunotherapy. Cold tumors have poor baseline pan-leukocyte infiltration, including a low prevalence of cytotoxic lymphocytes, and not surprisingly respond unfavorably to immune checkpoint (IC) inhibitors. We hypothesized that cold tumors harbor a mechanism of immune escape upstream and independent of ICs that may be driven by tumor biology rather than differences in mutational neoantigen burden.
Methods Using a bioinformatic approach to analyze TCGA (The Cancer Genome Atlas) RNA sequencing data we identified genes upregulated in cold versus hot tumors across four different smoking-related cancers, including squamous carcinomas from the oral cavity (OCSCC) and lung (LUSC), and adenocarcinomas of the bladder (BLCA) and lung (LUAD). Biological significance of the gene most robustly associated with a cold tumor phenotype across all four tumor types, glutathione peroxidase 2 (GPX2), was further evaluated using a combination of in silico analyses and functional genomic experiments performed both in vitro and in in vivo with preclinical models of oral cancer.
Results Elevated RNA expression of five metabolic enzymes including GPX2, aldo-keto reductase family 1 members AKR1C1, AKR1C3, and cytochrome monoxygenases (CP4F11 and CYP4F3) co-occurred in cold tumors across all four smoking-related cancers. These genes have all been linked to negative regulation of arachidonic acid metabolism—a well-established inflammatory pathway—and are also known downstream targets of the redox sensitive Nrf2 transcription factor pathway. In OCSCC, LUSC, and LUAD, GPX2 expression was highly correlated with Nrf2 activation signatures, also elevated in cold tumors. In BLCA, however, GPX2 correlated more strongly than Nrf2 signatures with decreased infiltration of multiple leukocyte subtypes. GPX2 inversely correlated with expression of multiple pro- inflammatory cytokines/chemokines and NF-kB activation in cell lines and knockdown of GPX2 led to increased secretion of prostaglandin E2 (PGE2) and interleukin-6. Conversely, GPX2 overexpression led to reduced PGE2 production in a murine OCSCC model (MOC1). GPX2 overexpressing MOC1 tumors had a more suppressive tumor immune microenvironment and responded less favorably to anti-cytotoxic T-lymphocytes-associated protein 4 IC therapy in mice.
Conclusion GPX2 overexpression represents a novel potentially targetable effector of immune escape in cold tumors.
- tumor escape
- immune tolerance
- immune evation
- immunotherapy
- head and neck neoplasms
Data availability statement
Data are available in a public, open access repository. Files containing recalculated RNA-seq FPKM-UQ values derived from publically availabe Toil datasets corresponding to the TCGA cohorts used in this study will be made available upon individual request.
This is an open access article distributed in accordance with the Creative Commons Attribution Non Commercial (CC BY-NC 4.0) license, which permits others to distribute, remix, adapt, build upon this work non-commercially, and license their derivative works on different terms, provided the original work is properly cited, appropriate credit is given, any changes made indicated, and the use is non-commercial. See http://creativecommons.org/licenses/by-nc/4.0/.
Statistics from Altmetric.com
WHAT IS ALREADY KNOWN ON THIS TOPIC
Immune escape is a hallmark of cancer with multiple mechanisms identified, but the biology behind existence of immunologically cold tumors with little baseline infiltration of leukocytes is poorly understood.
WHAT THIS STUDY ADDS
We have identified a common set of metabolic genes, including glutathione peroxidase 2, that are at the nexus of arachidonic acid metabolism and neutralization of reactive oxygen species that likely quench intracellular inflammation, shutting down NF-κB signaling and extracellular triggering of host immunity to create immunologically cold tumors despite high mutational neoantigen burdens.
HOW THIS STUDY MIGHT AFFECT RESEARCH, PRACTICE OR POLICY
We have identified novel potentially targetable enzymes that are likely responsible for a poor host antitumor response, which can be used theranostically to both predict which patient tumors may be poor candidates for immune therapy and to develop new approaches for overcoming treatment failure.
Background
Antitumor immunity is a critical regulator of tumorigenesis, metastasis and response to conventional chemoradiation regimens across a wide variety of solid tumor histologies. Large scale characterizations of the tumor immune microenvironment (TIME) across multiple cancer types have consistently revealed at least two kinds of immunosuppressive phenotypes.1–3 In the first category are tumors characterized by abundant infiltration of myeloid derived suppressor cells and regulatory T-cells (T-regs) often present simultaneously with cytotoxic CD8 +lymphocytes, sometimes referred to as immunologically ‘hot’ tumors. In the second category are immunologically ‘cold’ tumors that have minimal baseline leukocyte infiltration. While both scenarios likely contribute to immune tolerance, it is generally believed that hot tumors may respond better to immune checkpoint inhibitors (ICIs) than cold tumors because the former have pre-existing cytotoxic cells potentially capable of recognizing tumor derived neoantigens.4–7 For cold tumors, on the other hand, antitumor immunity is likely abrogated at a much earlier step—making this subtype potentially more challenging to treat with ICIs.
Suppression of the TIME is driven by multiple mechanisms including increased expression of immune checkpoint molecules like programmed death-ligand 1 (PD-L1) on tumors and their cognate programmed cell death protein-1 (PD-1) receptors on lymphocytes, or decreased generation/presentation of mutation-based or expression-based neoantigens by cancer cells.8 9 However, there is now increasing evidence that altered metabolic activity inside tumor cells plays a critical role in the recruitment of suppressive cells and the exhaustion of effector cells including CD8 +tumor infiltrating lymphocytes (TILs).10–13 A reduction of oxidative phosphorylation and increased production of lactate has broad based suppressive effects on functional immunocytes and has been shown by us and others to generate ‘cold’ tumors depleted of effector cells.12 14 Modulation of the kynurenine pathway through differential activation of indoleamine 2,3 dioxygenase has also been shown to cause T-cell anergy, increase proliferation of T-regs and alter the overall balance of the Th1/Th2 response.15 Just how broadly relevant these phenomena are remains unclear.
We hypothesized that there may be tumor-derived factors other than neoantigen burden driving an immunologically cold phenotype. Using an unbiased bioinformatics approach, we sought to identify candidate genes or pathways that may explain why a substantial subset of human tumors exhibit relatively minimal baseline leukocyte infiltration. Our analysis identified overexpression of glutathione peroxidase 2 (GPX2) to be robustly associated with cold tumors in four different tobacco associated cancers, including squamous cell carcinomas of oral cavity (OCSCC), lung (LUSC), lung adenocarcinomas (LUAD), and bladder cancer (BLCA). Forced overexpression of GPX2 in preclinical models of OCSCC reduced production of inflammatory mediators by tumor cells, generated an aggressive in vivo tumorigenic profile, reduced infiltration of effector immunocytes and reduced effectiveness of cytotoxic T-lymphocytes-associated protein 4 (CTLA-4) inhibition. Collectively, our data indicate that activation of GPX2 represents a previously understudied and potentially targetable potent driver of a suppressive TIME which could be responsible for limited ICI effectiveness in a significant subset of tumors.
Methods
Genomic data from TCGA cohorts
Files containing raw RNA-seq by Expectation Maximization (RSEM) count data from the The Cancer Genome Atlas (TCGA) RNA sequencing (RNA-seq) data sets that were previously harmonized computationally and mapped through the Toil open source calling pipeline16 were downloaded from Amazon web services S3 as previously described.17 RSEM counts were normalized to get upper quartile Fragments per Kilobase of transcript per Million mapped reads (FPKM) values (FPKM-UQ) with an in-house script in R studio. Essentially, non-protein coding RNAs were removed using annotations from Ensembl, which left 19,345 genes. The 75th percentile RSEM value (excluding zeros) for each TCGA sample was calculated and used to derive FPKM-UQ values. For most statistical analyses, the log2 (FPKM-UQ +0.01) transformation was applied. TCGA clinical and somatic mutation data were downloaded from cbioportal.18
Bioinformatic analysis of TIME in TCGA cohorts
Single sample gene set enrichment analysis (ssGSEA) using published gene lists specific for 13 different leukocyte subpopulations was performed using GenePattern software as we described previously,12 except FPKM normalized RNA-seq data from 20,000 different TCGA samples representing >30 cancer types was used to exclude genes with poor cross-correlation of expression from leukocyte gene lists. The final vetted immune subtype gene lists used for ssGSEA appear in online supplemental table 1. Because only one gene, FOXP3, represented T-regs, the log2 expression value was used in place of a ssGSEA score in downstream analyses for this leukocyte subtype. To derive immune phenotypes, two-way consensus hierarchical ward’s clustering based on random resampling was performed and optimal cluster numbers chosen using normalized Euclidian distances as we previously described,12 except a new in-house MATLAB script was employed (see online supplemental methods). To generate final heatmaps, ward’s clustering of the consensus clustering similarity matrices was performed with JMP V.13 software (SAS, Cary, North Carolina, USA).
Supplemental material
Supplemental material
Cells
Murine oral cancer (MOC) cell line MOC1 was obtained from our collaborator and maintained as previously described.19 Human head and neck squamous cell (HNSCC) lines, their maintenance, and genomic characterization have also been previously described.20 The viral packaging cell line 293FT was purchased from Thermo Fisher Scientific (Waltham, Massachusetts, USA) and maintained in Dulbecco’ s Modified Eagle Medium high glucose containing 10% heat-inactivated fetal bovine serum, non-essential amino acids (0.1 mM), sodium pyruvate (1 mM), 2 mM L-glutamine, penicillin (100 units/mL) and streptomycin (50 mg/mL).
GPX2 overexpression and knockdown
A mouse complementary DNA (cDNA) encoding GPX2 and the 3’ regulatory selenocysteine insertion element was cloned (online supplemental methods) into pLVX-IRES-Puro (online supplemental figure 1), which lacks exogenous tags. For in vitro GPX2 knockdown, human cells were infected with lentiviral pGIPZ-shRNA (Horizon, Waltham, Massachusetts, USA) targeting human GPX2 (V3LHS_368762 with antisense TCCAGGCCACATCTGAGCG) or an empty vector (EV) control. Details regarding infection with lentiviral constructs and isolation of cells with genetically altered GPX2 expression are provided in online supplemental methods.
Supplemental material
Western blotting
Cells were lysed in standard radioimmunoprecipitation assay buffer with detergent, and equal amounts of protein diluted in Laemmli SDS loading buffer plus reducing agent were resolved on SDS gels (Invitrogen, Carlsbad, California, USA) before transferring onto polyvinylidene difluoride (PVDF) membrane (Millipore, Temecula, California, USA). Blots were blocked with 5% non-fat milk, washed in TBS-tween, probed with primary antibody, and further washed before developing with horseradish peroxidase (HRPO)-conjugated secondary antibody and the ECL Western detection reagent (Thermo Scientific). Primary antibodies used in this study were as follows: from R&D systems, GPX2 (clone 496010); from Millipore Sigma, β-actin (clone AC-15).
Flow cytometry assessment of TIME in MOC1 tumors overexpressing GPX2
A total of 4×106 MOC1/wt, MOC1/pLVX-IRES-puro (EV), or MOC1/pLVX-IRES-puro-GPX2 cells were implanted subcutaneously in the left flank of female C57BL/6J mice (8–10 weeks). Tumors of equal size (minimum five per group) were dissociated in RPMI 1640 (Sigma-Aldrich) containing DNase I (20 U/mL; Sigma-Aldrich), Collagenase I (1 mg/mL; EMD Millipore) and Collagenase IV (250 U/mL; Worthington Biochemical Corporation) using a gentleMACS Dissociator (Miltenyi Biotec, California, USA) and then incubated at 37°C for 30 min to complete digestion. Tumor infiltrating leukocytes were enriched from single cell suspensions using Lymphoprep (STEMCELL Technologies) and analyzed in unblinded fashion using a panel of subset specific markers (online supplemental methods) as we previously described11 21 using the gating strategy illustrated in online supplemental figure 2.
In vivo response to immune checkpoint inhibitor anti-cytotoxic T-lymphocytes-associated protein 4
A total of 5×106 MOC1/wt, MOC1/pLVX-IRES-puro (EV), or MOC1/pLVX-IRES-puro-GPX2 cells were implanted subcutaneously in the left flank of female C57BL/6J mice (8–10 weeks). When tumors reached palpable size around day 7 post tumor inoculation, 10 mice were randomized into control or treatment groups. Treatment regimens included a total of three doses each of intraperitoneal injections of either InVivoMAb anti-mouse cytotoxic T-lymphocytes-associated protein 4 (CTLA-4) (clone 9H10; BioXCell; 100 µg per dose) or InVivoMAb polyclonal Syrian hamster IgG (Isotype Control; BioXCell; 100 µg per dose) once every 5 days. The mice were monitored, in unblinded fashion, two times per week and the tumor growth was measured using calipers. Tumor area (mm2) was calculated as A=L×W, where L is length and W is width. Mice were sacrificed when their tumors reached maximum allowable burden (above 200 mm2), if tumors became ulcerated, or animals lost >20% body weight.
Metabolomics analysis
Cells were plated at low density in biologic triplicate and grown exponentially (50%–70% confluence) for 48 hours. Samples were normalized to an internal standard depending on the method used during processing: L-creatine was used for the Luna method, with a coefficient of variation (CV) of 0.12; L-zeatine (CV of 0.11) was used in the negative method; and L-tryptophan (CV of 0.11) was used for the positive method. After normalization, data were log2 transformed and then differential expression was determined. Individual metabolite shifts were quantified using linear fold changes compared with the parental cell line. Only metabolites which generated a significant shift (false discovery rate (FDR)<0.1) in two GPX2 overexpression (OE) clones and did not generate a significant shift in the EV expressing cells were further analyzed. Individual metabolite functional analysis and pathway assigned was manually curated using the Human Metabolome Database (V.4.0). Secondary pathway enrichment analysis was performed using MetaboAnalyst (V.5.0).
Other statistical analyses
Differential analysis of gene expression (DEG) between hot and cold tumors was analyzed individually for each TCGA cohort after removing low expression genes (ie, max average log2 transformed FPKM-UQ for hot or cold groups <2) and performing multiple (two-tailed) t-tests with a Benjamini-Hochberg correction (FDR=0.1) using the Response Screening module in JMP. To find the intersection of up or downregulated genes across all four cohorts a minimum threshold of 1.5-fold change and FDR<0.1 was applied. Multiple t-testing with Benjamini-Hochberg correction or analysis of variance with post-hoc Tukey multiple comparison testing for other data sets were performed with Prism V.8 software (GraphPad, San Diego, California, USA). The Pearson product-moment correlation (ie, Spearman’s rho) for gene correlations, along with their FDR-adjusted p values were calculated in JMP. For some analyses, tumor purity values for the TCGA data sets were obtained through the ‘TCGA_mastercalls.abs_tables’ file available from the National Cancer Institute (NCI) genomic data commons website.
Results
Identification of immunologically hot and cold tumor samples
TCGA RNA-seq data were used to computationally quantify the relative tumor infiltration levels of leukocytes within solid tumors derived from four different smoking-related cancers, including OCSCC), LUSC, LUAD, and BLCA. Vetted gene lists (online supplemental table 1) defining different leukocyte subpopulations were used to derive ssGSEA scores for each leukocyte subtype in individual tumor samples from all four cancers (online supplemental table 2). For each TCGA cohort, the ssGSEA scores plus RNA values from FOXP1 (ie, T-regs) were then used to perform consensus hierarchical clustering with optimal selection of cluster numbers (ie, online supplemental figure 3) using approaches we previously described.12 Heatmaps from hierarchical clustering revealed easily identifiable and highly similar patterns of leukocyte infiltration among all four cancer types (figure 1). In each cohort we identified a subset of tumors with high levels of infiltration from nearly every immune subset examined (hot tumors) and at the opposite extreme a substantial number of samples with relatively little infiltration by any of the different leukocytes, which we refer to as ‘cold tumors’ (figure 1, online supplemental table 2). There were also varying groups of samples with an intermediate phenotype that sometimes separated based on whether cells mediating adaptive immunity or non-adaptive immunity predominated. Tumors with a higher influx of CD8 +T cells or cytotoxic cells in their TIME have frequently been associated with increased survival. Consistent with this observation, we found an increased tendency for improved overall median survival time in patients whose tumors were from cluster 1 (heavy leukocyte infiltration) compared with patients whose tumors had far reduced leukocyte infiltration, across all four tumor types (online supplemental figure 4) although statistical significance was not always reached.
Identifying correlates of tumor immunity in tobacco-related malignancies. (A) Clustering of OCSCC, LUSC, LUAD and BLCA based on single sample gene set enrichment analysis score for 13 leukocyte subtypes (and log2 FOXP3 RNA for regulatory T-cells) identifies subsets of immunologically hot and cold tumors. OCSCC samples included tumors from the oral cavity, oral tongue (excluding base of tongue), floor of mouth, buccal mucosa, alveolar ridge, hard palate, and lip. (B) Differential analysis of genes upregulated in hot and cold tumors across all four tumor histologies. AKR1, aldo-keto reductase family 1; BLCA, bladder cancer; CYP, cytochrome P; GCLC, glutamine-cysteine ligase catalytic subunit; GPX2, glutathione peroxidase 2; LUAD, lung adenocarcinomas; LUSC, lung squamous cell carcinomas; OCSCC, squamous cell carcinomas of oral cavity; NQO1, NAD(P)H quinone dehydrogenase 1.
For subsequent analysis, tumor clusters were grouped into the categories of hot, intermediate, or cold tumors based on robustness of pan-leukocyte infiltration (figure 1, online supplemental table 2). As a further validation we performed unbiased DEG expression to confirm higher expression of immune-related genes in hot tumors versus cold tumors. The number of significantly upregulated genes in hot tumors with at least a 1.5-fold average increase ranged from 1710 to 2285 among the four cancers (figure 1, online supplemental tables S3–S6), with 846 commonly upregulated genes (online supplemental table 7). Numerous immune-related genes and known regulators of inflammation, including many not part of gene lists used for ssGSEA, were commonly upregulated across hot tumors. Among these were 20 chemokines/chemokine receptors, 19 interleukins, and 17 human leukocyte antigen (HLA) genes (figure 1, online supplemental table 7). A more intense immune signature has sometimes been linked to a higher mutational tumor burden (ie, increased potential for neoantigens), while other investigators have found the exact opposite relationship depending on the tumor type or study.22 We found no significant differences in the total number of non-synonymous or truncating mutations between hot and cold tumors across all four smoking-related cancers (online supplemental figure 5) that could explain differences in leukocyte infiltration.
Overexpression of enzymes negatively regulating arachidonic acid metabolism dominates the profile of immunologically cold tumors
Next, we focused on genes overexpressed in the cold tumors. Fewer genes were significantly upregulated in the cold tumors (figure 1 and online supplemental tables S3–S6), with LUAD having the least (N=1049) and BLCA the most (N=1628). Only 117 genes were consistently upregulated across cold tumors from all four smoking-related cancers (figure 1, online supplemental table 8) and these were significantly enriched for multiple Gene Ontology (GO) pathways involving metabolism (online supplemental table 9), including oxidoreductase activity (FDR=0.0297). GPX2 showed the highest fold increase in cold tumors when averaged across all four cancer types (9.4-fold), ranging from 6.3-fold to 17.8-fold (figure 1, online supplemental table 8), followed by aldo-keto reductase family 1 member C2 (AKR1C2) which had the second highest average increase in expression (5.5-fold). Other AKR1 members (AKR1C1 and AKR1C3) as well as two cytochrome P450 genes (CYP4F11 and CYP4F3) that encode monooxygenases were also among the 117 common genes upregulated in cold tumors. An anti-inflammatory role of GPX2 is suggested by previous work showing that knockout of the gene in mice increases inflammation and tumorigenesis in an inflammatory model of colon carcinogenesis.23 Furthermore, depletion of GPX2 expression through siRNA and shRNA reportedly causes increased interleukin (IL)-1β-mediated activation of NF-kB, which is associated with increased expression of cyclooxygenase-2 (COX-2) and downstream production of pro-inflammatory prostaglandin E2 (PGE2) in human colon cancer cell lines.24 25 All three of the AKR1C family members elevated in cold tumors are capable of reducing PGE2 production by conversion to PGF2a,26 while CYP4F3 and CYP4F11 inactivate leukotriene 4B (a potent inducer of neutrophil chemotaxis) as well as other and inflammatory eicosanoids derived from arachidonic acid metabolism.27 28 Although all of these metabolic genes may contribute to the cold tumor phenotype, we focused on GPX2 for further analysis because of its strong link to inflammatory models of colon cancer in mice and previous studies correlating its expression to tumor progression and poor prognosis in multiple human cancers.29–31 Levels of GPX2 messenger RNA broadly correlated with the degree of pan-leukocyte infiltration among all four smoking-related cancers (online supplemental figure 6A). Because GPX2 has a restricted expression pattern in normal tissues with little RNA found in fibroblasts or leukocytes (online supplemental figure 6B), we examined whether increased tumor purity in cold tumors could account for large increases in GPX2 expression. Applying a conservative approach, linear GPX2 RNA data were normalized to tumor purity values publicly available for TCGA samples and still found to be significantly elevated in cold tumors (ie, threefold to ninefold increase, FDR<0.01) across all four cancer types (online supplemental figure 6C). Analysis of single cell RNA-seq data (GSE103322) for OCSCC specimens available through the Gene Expression Omnibus (GEO) data portal confirmed that GPX2 is predominately expressed by tumor cells and not leukocytes or other stromal compartments (online supplemental figure 6D). As further validation, we compared GPX2 expression in hot and cold tumors using data deposited in GEO (ie, GSE81089) from a TCGA-independent Swedish cohort of lung tumors. Using our same ssGSEA and clustering approach for the Swedish cohort (online supplemental figure 7), GPX2 RNA expression was found to be 13.6-fold higher in cold tumors (p<0.002).
GPX2—Nrf2 correlations with TIME
Eleven out of the 117 genes commonly upregulated in cold tumors are known downstream targets of Nrf2 activation (online supplemental table 8),32 33 including GPX2, AKR1C family members, CYP4F genes, glutamine-cysteine ligase catalytic subunit (GCLC), and NAD(P)H quinone dehydrogenase 1. Since mutations in the KEAP1/NRF2 pathway have been recently linked to poor leukocyte infiltration,34 we sought to determine whether the correlation between GPX2 and a cold TIME were specific to GPX2 or a byproduct of a more general Nrf2 pathway activation. We first used unbiased gene correlation methods to define a 138 gene signature of Nrf2 activation across all four tumor sites (online supplemental methods and figure 8A). GO analysis (online supplemental table 10) confirmed that this 138 gene signature was enriched for genes involved in processing of oxidative, metabolic, and xenobiotic stress consistent with the function of Nrf2 as a master oxidative stress sensor. The Nrf2-gene set (online supplemental table 11) and ssGSEA scores defined by this signature correlated strongly with Nrf2 pathway mutational status across all four tumor types (online supplemental figure 8B,C). GPX2 expression levels also correlated strongly with Nrf2 ssGSEA scores for OCSCC, LUSC, LUAD and less so for BLCA (online supplemental figure 9). Both GPX2 expression and Nrf2 activation (ie, ssGSEA scores) were significantly correlated with specific immune subsets when examined individually in OCSCC and BLCA (online supplemental figure 10A). However, the significance of correlations was greater for GPX2 than Nrf2 score in OCSCC and both the correlations and the significance levels were stronger for GPX2 than Nrf2 score in BLCA. Collectively, the data suggest that GPX2 was a stronger driver of poor leukocyte infiltration than Nrf2 activation, particularly for BLCA (online supplemental figure 10B).
To determine whether the relationship between GPX2 and Nrf2 is maintained in preclinical models of OCSCC we re-analyzed RNA-seq data from a panel of 65 HNSCC cell lines that we previously published.20 A strong correlation between GPX2 expression and Nrf2 ssGSEA score was found in the HNSCC cell lines (figure 2A, r=0.64; p<0.0001) regardless of HPV and NRF2 mutational status. We also found a strong correlation between GPX2 levels in the cell lines and expression of TP63 (figure 2B), which has previously been described as a transcriptional activator of GPX2.35 To better identify possible regulatory mechanisms in the cell lines we clustered the expression of GPX2, TP63, Nrf2 score, and GCLC—another downstream target of Nrf2 that strongly correlated with the Nrf2 ssGSEA scores both in cell lines and TCGA samples. Cell lines with the strongest GPX2 levels were found mainly in two clusters (figure 2C), including one with high Nrf2 scores and elevated expression of the Nrf2 target GCLC that was enriched for Nrf2 pathway mutations and another cluster with much lower Nrf2 scores but higher TP63 levels. Additional evidence for Nrf2-independent mechanisms also regulating GPX2 is apparent from clustering TCGA OCSCC specimens by the 117 genes upregulated in cold tumors and their Nrf2 activation scores, which revealed a group of cold tumors with relatively high GPX2 but low Nrf2 activation (online supplemental figure 11). Collectively, the data suggest a dual regulatory mechanism of GPX2 in HNSCC that may involve both stress-induced signaling by Nrf2 and endogenous oncogenic signaling through TP63.
Nrf2-dependent and Nrf2-indepdent GPX2 regulation in HNSCC tumor lines. (A) Correlation between GPX2 RNA expression and Nrf2 activation (ie, Nrf2 single sample gene set enrichment analysis ssGSEA score) in a large panel of HNSCC cell lines annotated by KEAP1/NRF2 mutational status or relative amount of Nrf2 activation (U.Q., upper quartile). Remaining WT samples with Nrf2 scores below the UQ are represented by black circles. (B) Correlation between GPX2 RNA and TP63 expression in the same panel of HNSCC cell lines. (C) Clustering HNSCC cell lines by Nrf2 activation and expression of GPX2, TP63, or a second downstream Nrf2 target (GCLC) revealed five clusters or patterns of expression. GCLC, glutamine-cysteine ligase catalytic subunit; GPX2, glutathione peroxidase 2; HNSCC, human head and neck squamous cell; ssGSEA, single sample gene set enrichment analysis.
GPX2 reduces inflammatory mediator production in tumors
We hypothesized that GPX2 could negatively regulate intracellular inflammation in tumor cells causing reduced secretion of pro-inflammatory mediators into the surrounding TIME. Conceptually, this would be difficult to validate from in vivo RNA-seq data sets because the mere presence of leukocytes in hot tumors is both a source and cause of inflammatory gene expression. Therefore, we leveraged the RNA data from our panel of HNSCC cell lines,20 representing pure tumor cells, to examine correlations between GPX2 expression and 69 different cytokines/chemokines. Remarkably, 15 different cytokines/chemokines were significantly associated (FDR<0.1) with GPX2 expression in the cell lines (online supplemental table 12 and figure 12A,B). The majority (ie, 14/15) represented pro-inflammatory genes or cytokines that had moderate but significant anti-correlations, supporting our hypothesis that GPX2 negatively regulates inflammation. GPX2 expression in the cell lines was also significantly anti-correlated with NF-kB activation assessed by ssGSEA (online supplemental figure 12C), consistent with an earlier report that GPX2 reduces NF-kB activity in colon cancer cell lines. Collectively, the complimentary associations of GPX2 and cytokine/chemokine expression with NF-kB activation in HNSCC cell lines (online supplemental figure 12D,E). suggest the former may broadly impact inflammatory mediators by attenuating NF-kB activity
Among chemokines negatively correlated with GPX2 were molecules strongly chemotactic for neutrophils (CXCL1,2,3) or monocytes, activated CD8 +T cells, NK cells, macrophages and dendritic cells (ie, CCL2 and CCL5). The pro-inflammatory cytokines IL-1A, IL-1B, IL-6, IL-8, IL-32, IL-20, IL-24, IL-23A were also inversely correlated with GPX2 expression. Interestingly, the immune checkpoint PD-L1 was also negatively correlated while IL-17D known to suppress CD8 +T cells was the only cytokine positively correlated with GPX2 expression. To understand better whether GPX2 expression or Nrf2 activation was a stronger driver of the anti-inflammatory genotype we performed the same analysis using Nrf2 ssGSEA scores derived from the cell lines and compared their anti-correlations and FDR adjusted p-values (online supplemental table 12 and figure 12A). Pro-inflammatory genes were more strongly anti-correlated with GPX2, than with Nrf2 activation scores. Finally, we analyzed correlations with GPX2 expression in an independent RNA-seq data set from a different smoking-related tumor type consisting of 24 LUSC tumor lines, which are part of the Cell Line Encyclopedia Project. GPX2 expression was similarly inversely correlated with multiple pro-inflammatory genes in LUSC cell lines (online supplemental table 13). IL-6, IL-32, along with PD-L1, were all strongly anti-correlated with GPX2, along with additional pro-inflammatory genes uniquely significant for the LUSC tumors. The suppressive cytokine IL-17D was also positively correlated with GPX2 in LUSC cell lines as in HNSCC.
To test the functional link between GPX2 and downstream proinflammatory mediators, we knocked down expression of the former using shRNA in two different human HNSCC cell lines (figure 3A), which expressed high (FaDu) or low (UMSC22A) baseline levels of GPX2 protein. Suppression of GPX2 with targeting shRNA but not an EV led to a twofold to threefold increase in levels of secreted IL-6 detectable in tumor supernatants collected at 48 hours and 72 hours following knockdown in both FaDu and UMSCC22A (adj p<0.01, figure 3B). GPX2 knockdown also resulted in increased secretion (1.6–1.8-fold) of proinflammatory PGE2 from both cell lines at each time point (adj p<0.001, figure 3C). Conversely, we overexpressed mouse GPX2 cDNA in MOC1, a preclinical murine model of oral carcinoma previously shown to have relatively high leukocyte infiltration when grown in vivo.19 Increased GPX2 protein expression by western blotting and a doubling of intracellular peroxidase activity in overexpressing clones were confirmed (figure 4A,B). Consistent with the human cell line data, overexpressing GPX2 drove levels of PGE2 in the opposite direction resulting in more than a twofold drop (figure 4C). IL-6 levels in parental cells were not detectable in either control MOC1 or after GPX2 OE.
GPX2 suppression increases inflammatory mediator production by HNSCC cell lines. FaDu and UM22A cells were infected with lentiviral constructs containing either empty vector (EV) or shRNA targeting GPX2. (A) Western blot validation of GPX2 protein knockdown 72 hours post infection. Insert demonstrates suppression of GPX2 protein levels at 48 hours post infection. Conditioned media was analyzed at 48 and 72 hours post-infection using ELISA to quantitate levels of secreted IL-6 (B) or PGE2 (C). Soluble mediator levels are normalized to total number of cells. Data are presented as means, with error bars denoting SD. * denotes p<0.05; ** denotes p<0.01; *** denotes p<0.001, **** denotes p<0.0001. IL, interleukin; GPX2, glutathione peroxidase 2; HNSCC, human head and neck squamous cell; PGE2, prostaglandin E2.
GPX2 overexpression suppresses PGE2 production in murine HNSCC cell lines. (A) Western blot validation of murine GPX2 OE in MOC1 cells. Clones 1 and 2 (C1 and C2) were derived from polyclonal stable cell lines selected with puromycin following infection with either an empty vector or one with the mouse GPX2 complementary DNA insert. (B) GPX2 OE doubled the intracellular GPX2 enzymatic activity measured in MOC1 cell lines. (C) PGE2 secretion was reduced by twofold in MOC1 clones with GPX2 OE. Stable overexpression of either empty vector or GPX2 containing constructs were accomplished in the MOC1 murine cell line (insert denotes Western blot). Individual clones (C) were established for both the empty vector control (EV) and for the GPX2 overexpressing constructs. and reduced PGE2 secretion. Data are presented as means, with error bars denoting SD. *** denotes p<0.001; **** denotes p<0.0001. GPX2, glutathione peroxidase 2; HNSCC, human head and neck squamous cell; PGE2, MOC, murine oral cancer; OE, overexpression; PGE2, prostaglandin E2.
GPX2 regulates cellular metabolism
To test whether GPX2 OE shifted cellular metabolism in a manner consistent with its previously described function of a stress response enzyme we performed steady state metabolomics analysis. GPX2 OE cells demonstrated a dramatic decrease in carnitine derivatives (online supplemental table 14) consistent with suppression of fatty acid beta oxidation, decreased levels of oxidized glutathione, glucose and disruption of multiple metabolic regulators of oxidative stress including taurine (−7.7, –4.5 linear decreases, respectively for each of 2 GPX2 OE clones, compared with MOC parental line), betaine (−3.7, –3.4, respectively), octulose-8-phosphate (1.8, 2.1, respectively) and sialic acid (3.6, 3.5, respectively). Several metabolites with known immunomodulatory activity were significantly increased in GPX2 OE cells including 3- hydroxykynurenine (2.1, 2.5, respectively), phenylalanine (2.3, 2.2, respectively) and phosphocholine (2.7, 2.1, respectively). These metabolic shifts were paralleled by changes in expression of metabolic genes measured through RNA-seq. GO enrichment analysis of genes significantly upregulated in both GPX2 OE clones compared with the EV and parental cell lines mapped to multiple metabolic pathways (online supplemental table 15) including the following pathways with enrichment ≥2 fold and FDR<0.0001: peptide metabolic process, translation and cellular nitrogen compound biosynthetic process. Genes downregulated in GPX2 OE cells mapped to pathways (enrichment ≥2 fold and FDR<0.05) including: double-stranded break repair via non-homologous end joining, inter-strand cross-link repair, glycoprotein biosynthetic process, phospholipid biosynthetic process, cellular response to stress and cellular response to DNA damage stimulus. The following genes mapped (Reactome) to selenocysteine synthesis were significantly upregulated in both OE clones: RPL13, FAU, RPS29, RPL18A, RPL12, RPL5, RPL9, RPL29 and UBA52.
GPX2 drives shifts in TIME and impacts ICI response
In three independent experiments (figure 5, online supplemental figure 13), tumors grown subcutaneously from MOC1 overexpressing GPX2 in immunocompetent C57BL/6J mice had significantly fewer overall percentage of T cells and decreased percentages of M1-macrophages, with reduced levels of CD8 +T cells significant in 2/3 experiments and trending towards significance in a third replicate experiment. Conversely, MOC1 GPX2 OE tumors had significantly increased levels (ie, almost double) of infiltrating granulocyte-myeloid derived suppressor cells (G-MDSC) compared with MOC1 EV or parental MOC1 in both experiments. For most leukocyte subsets, no statistically significant differences were found between the MOC1 parental and MOC1 EV tumors, but the average percentages were sometimes intermediate in value between MOC1 and MOC1 GPX2 OE, suggesting that the presence of viral vector sequences could have weakly impacted the TIME.
GPX2 drives tumor immune microenvironment shifts. MOC1 tumors (N=6) generated from parental (wild-type (WT)), empty vector (EV) or GPX2 containing constructs were harvested at the same size and analyzed via flow cytometry. All data are presented as means, with error bars denoting SD and individual tumor volumes showed using individual circles. * denotes p<0.05; ** denotes p<0.01; *** denotes p<0.001, **** denotes p<0.0001, after performing a Tukey multiple comparison test and further adjusting p values with a Benjamini-Hochberg correction (FDR<0.1) to control the family wise error within a cell line group due to multiple testing across different leukocyte subpopulations. The arcsine transformation was applied to percentages prior to statistical analysis. EV, empty vector; G-MDSC, granulocyte-myeloid derived suppressor cells; GPX2, glutathione peroxidase 2; MOC, murine oral cancer.
We tested whether the altered TIME in MOC1 GPX OE tumors, characterized by reduced T-cell infiltration and increased presence of G-MDSC, associated with GPX2 OE would impair antitumor efficacy of immune checkpoint therapy. In two independent experiments, single agent blockade of CTLA-4 was less effective at blocking tumor progression in MOC1 GPX2 OE tumors compared with MOC1 EV control. In a pooled analysis from the two experiments (figure 6), just 11/20 (ie, 55%) mice with MOC1 GPX2 OE tumors responded compared with 19/20 (ie, 95%) of mice harboring MOC1 EV tumors (p<0.009). This was accompanied by a highly significant reduction of survival time (p<0.005) following anti-CTLA-4 treatment in mice harboring tumors overexpressing GPX2. Interestingly, the anti-CTLA-4 response observed for MOC1 EV tumors was more robust than previously reported for MOC1 parental tumors,36 raising the possibility that the presence of viral vector sequences does impart some added immunogenicity that was reversed by GPX2 OE. Nevertheless, MOC1 EV tumors grew in 100% of control (isoform) treated mice, although the GPX2 OE tumors were comparatively more aggressive in the absence of anti-CTLA-4 (p<0.0001, figure 6C, D and E). Our preclinical data suggest that patients with higher GPX2 expression may be less responsive to ICIs. To that end, we examined the relationship between GPX2 expression in TCGA samples (OCSCC and BLCA) and their tumor inflammatory score (TIS) derived using an 18 gene pan-cancer immune signature that robustly predicts ICI response across multiple cancers.37 We found highly significant inverse correlations (p<5×10–11) between levels of GPX2 RNA and the TIS (online supplemental figure 14). In BLCA the anti-correlation with TIS was much greater for GPX2 than Nrf2, supporting that the immunophenotype was more strongly driven by GPX2.
GPX2 reduces ICI effectiveness. Pooled analysis of two independent experiments in which MOC1 tumors with EV (A,B) or GPX2 OE (C,D) were established and allowed to grow in the flank model. Established tumors underwent treatment with three total injections of either isotype (iso) control (A,C) or anti-CTLA4 antibody (B,D). Tumor measurements are presented as individual symbols (square=experiment 1; circle=experiment 2) and each tumor is represented as an individual curve for the growth panels (red=EV; blue=OE). Each treatment group had 20 tumors (exp I=12; exp II=8) (E) Survival is denoted using Kaplan-Meier curves as an aggregate of the treatment groups. (F) Table of responders and non-responders following anti-CTLA4 antibody treatment. CTCLA-4, cytotoxic T-lymphocytes-associated protein 4; EV, empty vector; GPX2, glutathione peroxidase 2; ICI, immune checkpoint inhibitor; MOC, murine oral cancer; OE, overexpression.
GPX2 expression is a poor prognostic marker in OCSCC
Consistent with an immunologically cold phenotype, high levels of GPX2 RNA were strongly correlated with worse survival in the OCSCC TCGA cohort. Patients with OCSCC with high GPX2 expression in their tumors had a median survival time of just 32.4 months compared with 71.2 months for patients with low GPX2 expression (p=0.0092, figure 7), which represents a greater than twofold difference. High Nrf2 activation was similarly associated with a substantial reduction in survival time as well (p=0.039), although it was statistically not as significant as GPX2 expression.
Both GPX2 and Nrf2 activation correlate with poor survival in OSCC. Recursive partitioning was used to dichotomize patient tumors as either high or low for GPX2 RNA expression (left) or Nrf2 activation (right). Kaplan-Meier curves of overall survival were analyzed for statistical significance with a log-rank test. HRs for tumors expressing high GPX2 and Nrf2 activation are 1.593 and 1.452, respectively. GPX2, glutathione peroxidase 2; OSCC, squamous cell carcinomas of oral cavity.
Pan-cancer analysis of GPX2 expression in cold tumors
The link between GPX2 expression and cold tumors was examined in 23 other cancer types and in the two other main subsites of head and neck cancer: squamous cell carcinomas of the oropharynx (OPSCC) and larynx/hypopharynx (LHSCC). Although differences in GPX2 expression did not correlate with the TIME in OPSCC (online supplemental figure 15A,B), a majority of these tumors harbored high-risk oncogenic human papilloma virus-16 which was itself strongly associated with hot tumors (p=0.02). Analysis of LHSCC identified significant OE of arginase 2 in cold tumors (online supplemental figure 15C), rather than GPX2, suggesting these tumors utilize an alternative metabolic mechanism for immunosuppression. We examined the correlation between GPX2 RNA expression and leukocyte infiltration (based on ssGSEA scores) across an additional 23 cancer types (online supplemental figure 16A) to identify tumors that behaved similarly to smoking-related cancers. Levels of GPX2 RNA (after adjusting for tumor purity) were significantly elevated in cold tumors compared with hot tumors in colon and rectal adenocarcinomas (online supplemental figure 16B,C), but the opposite relationship was found for a few cancers where GPX2 expression levels were comparatively lower, such as prostate adenocarcinoma (online supplemental figure 16A,D). In other cancers, such as cholangiocarcinoma, moderate inverse correlations with GPX2 and leukocyte infiltration did not persist after adjusting for tumor purity (data not shown).
Discussion
For nearly a decade we have been both dazzled by the amazing effectiveness of ICIs in the context of previously untreatable solid tumors (eg, metastatic melanoma) and disappointed when these effects fail to materialize in less immunogenic tumors such as OCSCC.38 We have come to realize that tumor cells are able to modulate their immunological phenotype much like they can internal signaling networks (eg, protein phosphorylation cascades). ICI effectiveness is critically linked to the presence and activation of TILs.4 Similarly, CAR-T therapies require that lymphocytes infiltrate and maintain activation inside of solid tumors. In both cases, if tumor cells can either prevent infiltration or suppress activation once TILs have reached the tumor, the therapy will be ineffective. Our data indicate that GPX2 activation likely does both by preventing secretion of inflammatory mediators. In particular, the potential of GPX2 to blunt very early stages of tumor inflammation by suppressing pro-inflammatory chemokines and cytokines sets it apart mechanistically from the more widely studied processes of immune escape tied to expression of immune checkpoint molecules like PD-L1.
GPXs are ubiquitously expressed selenium-dependent enzymes that protect cells against oxidative damage by reducing hydrogen peroxide and a wide range of organic peroxides using glutathione as a primary source of reducing potential. GPX2, one of the eight members of the GPX family (GPX1–8), is also called gastrointestinal GPX because of its relatively high expression in the gastrointestinal tract. Because of their antioxidative activity, several isoforms of GPXs have been investigated in association with inflammation and cancer.39 GPX2 in particular has been shown to play an important role in the progression of malignant tumors with GPX2 OE detected in liver,40 colorectal,41 42 breast,43 lung44 and castration-resistant prostate tumors.31 In the context of tumorigenesis, a preventive role of GPX2 has been proposed because mice, in which both GPX1 and GPX2 had been knocked out, progressively developed ileocolitis and subsequently intestinal cancer, especially when not raised under specific pathogen-free conditions.45
As with all other regulators of oxidative stress, expression of GPX2 likely generates a push-stop phenomenon whereby sufficient oxidative stress is allowed for tumor cell proliferation, growth and metastasis but excess levels are kept in check to prevent premature tumor cell death. It has been described as a tobacco response gene in both lung and oral cavity epithelium.44 46 GPX2 and AKR1Cs both have both been linked to development of resistance to chemotherapy and radiation in multiple tumor models including HNSCC and LUAD.47 48 We previously showed that HNSCC cell lines adapted to oxidative stress rewire metabolism49 in a manner similar to that observed when GPX2 is exogenously overexpressed in the current manuscript.
Adaptation to oxidative stress, whether induced by genotoxic agents (eg, cisplatin or radiation), inflammatory cells (eg, neutrophils) or tobacco exposure likely triggers fundamental metabolic reprogramming which involves, at least in part hyperactivation of GPX2. What we have now shown mechanistically is that when this happens, there are important immunologic consequences. OE of GPX2 in an immune competent HNSCC model not only reduced effector TIL levels, but also increased levels of suppressive immunocytes including G-MDSC. Future experiments with additional mouse tumor lines where GPX2 is manipulated are needed to dissect how existing genetic background and intracellular tumor signaling may be influenced by GPX2 to alter the inflammatory secretome. One possible way GPX2 may shift the TIME is by attenuating NF-kB mediated inflammation as suggested by prior reports in which GPX2 knockdown in human colon cancer cells reduced NF-kB activation, downstream expression of COX-2, and increased PGE2 synthesis25 NF-kB is stimulated by intracellular reactive oxygen species,50 and GPX2 is thought to inhibit NF-kB indirectly by reducing ROS.25 In support of this model, we found a plethora of pro-inflammatory cytokines and chemokines whose expression inversely correlated with GPX2 levels in a large panel of established human HNSCC lines. Many of these negatively correlating cytokines, including IL-1, IL-6, IL-8, CCL2, CCL5, CXCL1, and CXCL2 are all downstream targets of NF-kB activation We have confirmed in vitro that GPX2 knockdown reduced IL-6 production in two different human HNSCC cell lines tested. Our model suggests that GPX2 OE in tumors limits inflammation by dampening initial recruitment and subsequent activation of effectors from both innate and adaptive immunity. This model could account for why GPX2 is linked to the cold TIME in a subset of tumor types. Possibly, NF-kB -induced inflammation may be more relevant to smoking-related cancers and colorectal tumors. In cancer types with a comparatively lower range of GPX2 expression, such as prostate adenocarcinoma, neutralization of intracellular ROS may be incomplete and thus elevated GPX2 would become a surrogate for unresolved intracellular inflammation that leads to increased leukocyte recruitment.
While this explanation is inherently simple, it may in fact be incomplete as several of the metabolites which increase following constitutive GPX2 OE have been shown to directly regulate immunocyte activity (eg, taurine, kynurenine derivatives). Basic shifts in energy and biomass balance have been shown by us and others to correlate with shifts in TIME at least in part to differential balance of oxidative phosphorylation and lactate producing glycolytic activity.12 Data generated here suggest that additional metabolic shifts may further complicate this metabolic link between tumor cells and immunocytes. Furthermore, a number of other enzymes including the AKR1 family members and cytochrome monoxygenases were upregulated in parallel with GPX2 in cold tumors. These enzymes can directly regulate the arachidonic acid cascade, altering production of pro-inflammatory prostaglandins and leukotrienes, complementing the ability of GPX2 to inhibit PGE2 expression via COX2 downregulation.
Because GPX2 has been functionally liked to chemotherapy resistance47 we hypothesize that treatment failure to first line genotoxic therapies will select for patient tumors with a metabolic phenotype less likely to respond to subsequent ICIs. Two future directions remain our current focus. First, development of small molecule GPX2 inhibitors is critical to further preclinical validation of this immunomodulatory model in combination with clinically available ICIs. Second, expansion of functional genomic experiments to include other Nrf2-regulated metabolic enzymes upregulated in cold tumors will allow us to dissect their individual and collective contributions to modulating the TIME in preclinical models of HNSCC, BLCA, and LUAD. Such experiments will help clarify how best to clinically target these novel putative metabolic drivers of immunological escape in smoking-related cancers.
Data availability statement
Data are available in a public, open access repository. Files containing recalculated RNA-seq FPKM-UQ values derived from publically availabe Toil datasets corresponding to the TCGA cohorts used in this study will be made available upon individual request.
Ethics statements
Patient consent for publication
Ethics approval
Human subjects research restricted to public databases and therefore exempt; all mouse studies performed with institutionally approved IACUC protocols.
References
Supplementary materials
Supplementary Data
This web only file has been produced by the BMJ Publishing Group from an electronic file supplied by the author(s) and has not been edited for content.
Footnotes
Contributors MAK: Performed experiments and analyzed data. RV: Performed experiments and analyzed data. DD: Performed experiments and analyzed data. NRP: Performed experiments and analyzed data. VP: Performed experiments and analyzed data. MFC: Wrote code, downloaded and analyzed data. DAW: Provided intellectual input and edited manuscript. WD: Provided intellectual input and edited manuscript. AIF: Developed and algorithms and analyzed data. SK: Performed experiments and analyzed data. AS: Provided intellectual input and edited manuscript. VS: Performed experiments, analyzed data, intellectual input, and writing. MJF: Project conception, data analysis, writing, intellectual input, supervision. MJF as the guarantor for the current manuscript accepts full responsibility for the finished work and/or the conduct of the study, had access to the entire dataset, and controlled the decision to publish.
Funding VCS is supported by a Career Development Award from the Veterans Administration Clinical Science Research and Development division (1IK2CX001953) and an American Cancer Society Research Scholar Grant. This work was supported by the National Institute of Dental and Craniofacial Research through R03DE028858 (VCS) and the Sid W. Richardson Foundation (MJF&VCS). This project was also supported by the Cytometry and Cell Sorting Core at Baylor College of Medicine with funding from the CPRIT Core Facility Support Award (CPRIT-RP180672), the NIH (P30 CA125123 and S10 RR024574) and the expert assistance of Joel M. Sederstrom.
Disclaimer The contents of this manuscript do not represent the views of the U.S. Department of Veterans Affairs or the United States Government.
Competing interests None declared.
Provenance and peer review Not commissioned; externally peer reviewed.
Supplemental material This content has been supplied by the author(s). It has not been vetted by BMJ Publishing Group Limited (BMJ) and may not have been peer-reviewed. Any opinions or recommendations discussed are solely those of the author(s) and are not endorsed by BMJ. BMJ disclaims all liability and responsibility arising from any reliance placed on the content. Where the content includes any translated material, BMJ does not warrant the accuracy and reliability of the translations (including but not limited to local regulations, clinical guidelines, terminology, drug names and drug dosages), and is not responsible for any error and/or omissions arising from translation and adaptation or otherwise.