Redox^high phenotype mediated by KEAP1/STK11/SMARCA4/NRF2 mutations diminishes tissue-resident memory CD8+ T cells and attenuates the efficacy of immunotherapy in lung adenocarcinoma

ABSTRACT

Metabolism reprogramming within the tumor microenvironment (TME) can have a profound impact on immune cells. Identifying the association between metabolic phenotypes and immune cells in lung adenocarcinoma (LUAD) may reveal mechanisms of resistance to immune checkpoint inhibitors (ICIs). Metabolic phenotypes were classified by expression of metabolic genes. Somatic mutations and transcriptomic features were compared across the different metabolic phenotypes. The metabolic phenotype of LUAD is predominantly determined by reductase-oxidative activity and is divided into two categories: redox^high LUAD and redox^low LUAD. Genetically, redox^high LUAD is mainly driven by mutations in KEAP1, STK11, NRF2, or SMARCA4. These mutations are more prevalent in redox^high LUAD (72.5%) compared to redox^low LUAD (17.4%), whereas EGFR mutations are more common in redox^low LUAD (19.0% vs. 0.7%). Single-cell RNA profiling of pre-treatment and post-treatment samples from patients receiving neoadjuvant chemoimmunotherapy revealed that tissue-resident memory CD8+ T cells are responders to ICIs. However, these cells are significantly reduced in redox^high LUAD. The redox^high phenotype is primarily attributed to tumor cells and is positively associated with mTORC1 signaling. LUAD with the redox^high phenotype demonstrates a lower response rate (39.1% vs. 70.8%, p = 0.001), shorter progression-free survival (3.3 vs. 14.6 months, p = 0.004), and overall survival (12.1 vs. 31.2 months, p = 0.022) when treated with ICIs. The redox^high phenotype in LUAD is predominantly driven by mutations in KEAP1, STK11, NRF2, and SMARCA4. This phenotype diminishes the number of tissue-resident memory CD8+ T cells and attenuates the efficacy of ICIs.

KEYWORDS:

• Immune checkpoint inhibitors
• KEAP1/STK11/SMARCA4/NRF2 mutations
• lung adenocarcinoma
• metabolic phenotypes
• tissue-resident memory CD8+ T cells

Introduction

Lung cancer is the leading cause of cancer-related mortality worldwide.¹ In recent years, Immune checkpoint inhibitors (ICIs) targeting PD-1 or PD-L1 have increased the five-year overall survival rate from 8.5% to 21.9% in patients with advanced non-small cell lung cancer (NSCLC).² However, approximately half of NSCLC do not respond to ICIs and exhibit a poor prognosis, highlighting an unmet clinical need.³ On the one hand, a robust predictive biomarker is needed to identify patients who will respond to ICIs. On the other hand, the mechanism underlying primary resistance to ICIs is still not fully understood.

Metabolic reprogramming is a recognized hallmark of cancer.⁴ Cancer cells genetically and epigenetically alter their metabolism to meet their bioenergetic and biosynthetic demands. Tumor metabolism demonstrates individual heterogeneity, affected by genetic and environmental factors.⁵ Recently, an association between diverse bioenergetic phenotypes and metabolic dependencies has been identified in NSCLC, suggesting that metabolic dependency of cancer can result in specific metabolic phenotypes.⁶ Emerging evidence suggests that competition for nutrients between tumor and immune cells, as well as by-products of tumor cells within the tumor microenvironment (TME), substantially affects anti-tumor immunity.⁷ Upregulation of glutamine metabolism in cancer cells can lead to glutamine depletion, suppressing T cell proliferation and differentiation.^8,9 Tumor-derived metabolites such as fumarate and kynurenine can weaken antitumor immunity by reducing the anti-tumor capacity of CD8+ T cells¹⁰ and promoting regulatory T cell expansion.^11,12 This metabolic interplay helps the tumor to evade immune surveillance.¹³

Given that the complex interplay of metabolism within the TME can have a profound impact on immune cells, it is important to advance our understanding of these diverse metabolic phenotypes in the TME. Identifying the association between metabolic phenotypes and immune cells, especially CD8+ T cells, may reveal mechanisms of primary resistance to ICIs and help discover predictive biomarkers. To address the questions above, we classified the metabolic phenotypes of lung adenocarcinoma (LUAD), the most common pathological subtype of lung cancer. The mechanisms shaping these metabolic phenotypes, as well as their impact on CD8+ T cells are investigated.

Material and methods

Study design

A study flowchart was illustrated in Figure 1. Metabolic phenotypes in The Cancer Genome Atlas (TCGA) LUAD cohort were classified using unsupervised consensus clustering based on the RNA expression profiles of metabolic genes. The dataset was split into a training set and a validation set at a 7:3 ratio to establish a metabolic scoring model. Somatic mutations and transcriptomic features were compared across the different metabolic phenotypes. The mechanisms shaping these metabolic phenotypes and their impact on CD8+ T cells were investigated by single cell profiling analysis of treatment-naïve or neoadjuvant chemoimmunotherapy samples. Furthermore, the impact of metabolic phenotypes on the efficacy of ICIs was investigated in two independent immunotherapy cohorts.

Figure 1. Overview of study flow chart.

Patient cohorts

Discovery cohorts

LUAD patients with paired RNA sequencing and whole-exome sequencing data (n = 471) were included from TCGA. An additional LUAD cohort (n = 204) from the Gene Expression Omnibus (GSE31210)¹⁴ and a pan-cancer cohort from TCGA (n = 10,103) were used to validate the prognostic stratification capability of the metabolic scoring model.

Immunotherapy cohorts

Two independent immunotherapy cohorts (the Ravi cohort and the OAK cohort) were included to compare the efficacy of ICIs between metabolic phenotypes.^15,16 The inclusion criteria were LUAD histology, the availability of RNA-seq data, confirmed response, and availability of progression-free survival (PFS) and overall survival (OS) data. Responders were defined as patients achieving a partial or complete response, or stable disease with PFS of at least six months. Non-responders were those with disease progression or stable disease with PFS of less than six months. The response rate was defined as the ratio of responders to the sum of responders and non-responders.

Single-cell RNA sequencing cohorts

Eleven fresh tumor tissues from ten NSCLC patients undergoing lung surgery were obtained at the Guangdong Lung Cancer Institute, Guangdong Provincial People’s Hospital (GLCI cohort). Following surgical resection, samples were immediately processed for tissue preparation and single-cell RNA sequencing (scRNA-seq). The metabolic phenotype shaping mechanism and impact on CD8+T cells were explored in this scRNA-seq cohort. Additionally, another scRNA-seq cohort, which included three pre-treatment and twelve post-treatment samples from fifteen NSCLC patients who had received neoadjuvant chemoimmunotherapy, was included as a validation cohort (Hu cohort).¹⁷ This study was approved by the Ethics Committee of Guangdong Provincial People’s Hospital (approval No. GDRECKY2020-199-01). The study adhered to the principles of the Declaration of Helsinki.

Classifying metabolic phenotypes

A total of 1,540 metabolic genes were gathered from 84 pathways involving human metabolism in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Supplemental Table S1). Metabolic phenotypes were classified by unsupervised consensus clustering algorithm based on RNA expression of the metabolic genes. Key metabolic modules were identified by the weighted gene co-expression network analysis (WGCNA) in each metabolic phenotype. Metabolic modules were annotated by the Gene Ontology (GO) and KEGG pathway using clusterProfiler (version 4.8.2). Metabolic flux balance was assessed using METAFlux to substantiate the metabolic phenotype.¹⁸

Analysis of molecular features and tumor ecosystem

Somatic variants identified by the MuTect2 pipeline were analyzed and visualized using maftools (version 2.16.0). We compared frequently mutated genes (≥5%) between metabolic phenotypes. Copy number aberrations at the gene level were also analyzed across metabolic phenotypes. The signature scores of the hallmark gene sets from the Molecular Signatures Database, and those of the tissue-resident memory CD8+ T cells (CD8+ Trm) (CD8A, CD8B, CD69, CD103), were calculated using GSVA (version 1.48.3). Cell communities and states were identified using EcoTyper.¹⁹ Cell composition and cell abundance were compared between metabolic phenotypes. Abundance of twenty-four immune cell was also estimated by ImmuCellAI.²⁰

Establishment of the redox scoring model

The least absolute shrinkage and selection operator (LASSO) binomial analysis was performed to reduce redundant metabolic genes. Ten-fold cross-validation was used and lambda.min was chosen. Logistic regression model was used to construct the redox scoring model. Expression of selected genes was normalized as follows: normalizedGene_{i =} (gene_i *10^5)/[ACTB_i*(gene₁+gene₂+ … +gene_i)]. gene_i and ACTB_i indicates expression level of selected gene and actin beta (ACTB). The standard deviation of our normalization approach was the smallest when compared with those of other gene expressions, suggesting its stability across tumor samples (Supplemental Fig S1A). Redox score was calculated as: redox score =${\sum }_{i=1}^n\mathit{Co}{e}_i\ast \mathit{normalizedGen}{e}_i$. Coe_i is the corresponding model coefficient for the gene. Cut-off point to define low or high redox score was derived in corresponding datasets using maximally selected rank statistics method.

Immunohistochemistry validation

Twenty-eight patients diagnosed with LUAD from the Human Protein Atlas (https://www.proteinatlas.org/) were included. The expression of genes in the redox scoring model was evaluated by immunohistochemistry in the tumor tissue. Staining intensity was evaluated as negative, weak, moderate, or strong; the fraction of stained cells was evaluated as < 25%, 25–75%, or > 75%.

Single cell RNA sequencing

Single-cell RNA-seq experiment was performed by NovelBio Bio-Pharm Technology Co., Ltd. The detailed method of scRNA-seq pipeline were the same as previous described.²¹ The Seurat package (version 4.3.0.1) was employed to create the object, filtering out cells of poor quality. The library size was normalized using the Seurat function NormalizeData. The top 2000 variable genes, identified by the FindVariableFeatures function, were used for principal components analysis. The top 10 principal components were retained for uniform manifold approximation and projection visualization and clustering. Batch effects were corrected using the fastMNN function from the batchelor package (version 1.16.0), employing the mutual nearest-neighbor method. Cell clustering was conducted using the FindClusters function with a graph-based clustering method. Marker genes were identified using the FindAllMarkers function with the Wilcoxon rank-sum algorithm, under the criteria: lnFC > 0.25, p < 0.05, and min.pct >0.1. Cell types were manually annotated based on representative cell markers. Malignant cells were identified by inferCNV analysis. Redox signature and mTORC1 signaling were quantified using the AddModuleScore function.

Cell-cell communication analysis

We performed Nichenet analysis²² to infer the interactions between tumor cells and CD8+ T cells. Considering genes expressed in more than 10% of cells in each cluster, we extracted the top 20 ligands and top 30 receptor targets based on differential expression between ‘sender cells’ (tumor cells) and ‘receiver cells’ (CD8+T cells) for paired ligand-receptor activity analysis. Tumor cells with a low redox signature served as reference cells, and feature genes of CD8+Trm cells were identified as differentiating genes.

Statical analysis

Continuous variables were compared using the Wilcoxon rank-sum test. Categorical data were analyzed using the chi-square test or Fisher’s exact test. Kaplan-Meier curve was used to estimate median OS and PFS. Survival significance between groups was assessed using the log-rank test. Odds ratios (ORs) and 95% confidence intervals (CIs) for confounding factors were calculated via binary logistic regression analysis. Hazard ratios (HRs) and 95% CIs were computed using Cox proportional hazards survival analysis. The Spearman test was used to quantify the strength of the association between variables. A two-sided p-value of less than 0.05 was considered statistically significant unless if indicated. All statistical analyzes were conducted using R software (version 4.3.0).

Results

Metabolic clustering distinguishes a metabolic phenotype of high reductase-oxidative activity in LUAD

The TCGA LUAD cohort was clustered into two metabolic phenotypes (k = 2, Cluster A and Cluster B) (Figure 2(a)). Cluster B demonstrated decreased OS compared to Cluster A (Figure 2(b), p < 0.001). Multivariate analysis identified Cluster B (HR: 1.51, p = 0.016) and stage (HR: 2.33, p < 0.001) as independent risk factors for OS (Supplemental Fig S1B), with stratified survival analysis indicating a decrease in OS for both early and advanced stages in Cluster B (Supplemental Fig S1C).

Figure 2. Identification of metabolic phenotypes in lung adenocarcinoma. (a) Two metabolic phenotypes of LUAD (Cluster a and (b) identified by unsupervised consensus clustering. (b) Kaplan–Meier analysis reveal poor prognosis in Cluster B. (c) a dendrogram of clustered genes in WGCNA, with each color representing a co-expressed gene module and each branch representing a gene. (d) Correlation heatmap showing the relationship between metabolic phenotypes and metabolic gene modules. Numbers indicate the correlation coefficient, followed by a statistical p-value in parentheses. (e) Differential expression analysis of three gene modules (MEbrown, MEyellow, and MEblue) between two metabolic phenotypes. (f) Cluster B identified as LUAD with high reductase-oxidative activity (redox^high LUAD). (g) Output of superoxide anion (O₂-) and pathway score of ROS detoxification analyzed by metabolic flux balance. LUAD, lung adenocarcinoma.

Five metabolic gene modules were identified by WGCNA between two metabolic phenotypes (Figure 2(c)). MEyellow module was the most significantly associated with Cluster B (r = 0.53, p < 0.001), followed by MEblue (r = 0.42, p < 0.001), while MEbrown module was more associated with Cluster A (r = 0.28, p < 0.001) (Figure 2(d)). Expression of MEbrown module genes tended to be higher in Cluster A, while MEyellow or MEblue module genes were more expressed in Cluster B (Figure 2(e), Supplemental Table S2). Compared with Cluster A, oxidative phosphorylation, reactive oxygen species (ROS) production and reductive stress pathways (glutathione metabolism, biosynthesis of cofactors, pentose and guluronate interconversions) were significantly enriched in Cluster B (Figure 2(f)). Metabolic flux balance analysis revealed a significant increase in superoxide anion (O2-) production in Cluster B (0.0056 vs. 0.0045, p = 0.003), with no marked difference observed for H₂O₂ or hydroxide (Figure 2(g)). Furthermore, the pathway scoring for ROS detoxification was reduced in Cluster B (0.0258 vs. 0.00269, p = 0.036), suggesting ROS accumulation (Figure 2(g)). These results indicate that Cluster B represents a metabolic phenotype characterized by high reductase-oxidative activity, termed redox^high LUAD, as opposed to Cluster A (redox^low LUAD).

Redox^high LUAD was characterized as highly proliferative tumor driven by KEAP1/STK11/SMARCA4/NRF2 mutations

Analysis revealed no significant differences in age, smoking history, or clinical stage distribution between metabolic phenotypes, with the exception of a higher male prevalence in the redox^high LUAD (Supplemental Fig S2A). For somatic mutations, KEAP1 mutation (55.6% vs. 2.4%, p < 0.001), STK11 mutation (22.2% vs. 9.8%, p = 0.001), SMARCA4 mutation (16.3% vs. 4.6%, p < 0.001), and LRRIQ3 mutation (11.1% vs. 3.1%, p = 0.001) were significantly higher in redox^high LUAD, while EGFR mutation (0.7% vs. 19.0%, p < 0.001) was more common in redox^low LUAD (Figure 3(a-b)). Despite NRF2 activating mutations occurring at a frequency (2.4%) below our initial threshold, a comparative analysis between phenotypes revealed a higher occurrence of NRF2 activating mutations within the redoxhigh phenotype (5.2% vs. 1.2%, p = 0.02). No difference was found in the mutation frequency of TP53 (52% vs. 44%) and KRAS (28% vs. 30%) between two metabolic phenotypes. EGFR mutation did not co-mutate with KEAP1, STK11, or SMARCA4 mutation (Supplemental Fig S2B). The incidence of AKR1C2 and AKR1C3 amplification were higher in the redox^high phenotype (AKR1C2: 10.5% vs. 0.6%, p = 0.03; AKR1C3: 9.7% vs. 0.6%, p = 0.04). It was noted that 75.0% (12/16) of patients with AKR1C2 or AKR1C3 amplification coexisted with mutations in STK11, SMARCA4, KEAP1, or NRF2. Copy number aberrations of STK11 (3.0% vs. 0.9%), KEAP1 (1.5% vs. 0.6%), SMARCA4 (4.5% vs. 0.6%), or NRF1 (3.0% vs. 0.9%) were numerically higher in the redox^high phenotype, though no statistical difference was observed.

Figure 3. Comparison of Molecular features and Enrichment Scores of Hallmark Gene Sets between Metabolic Phenotypes. (a) Volcano plots comparing somatic mutations across metabolic phenotypes, with the horizontal axis representing the odds ratio and the vertical axis representing the statistical significance; (b) Bar plot shows mutation frequencies of representative genes in redox^high and redox^low LUAD; (c) Horizontal axis representing pathway signatures and the vertical axis representing statistical significance. Dots above the horizontal axis indicate pathways upregulated in redox^high LUAD, below indicate pathways upregulated in redox^low LUAD.

Cell proliferation (G2M checkpoint, unfolded protein response, DNA repair, and mitotic spindle) was more active in redox^high LUAD, while apoptosis and immune response pathways (inflammatory response, complement, interferon-γ and interferon-α response) were more active in redox^low LUAD (Figure 3(c)). In cell signaling pathways, mTORC1 signaling was the most significantly elevated pathway in redox^high LUAD. In contrast, KRAS signaling, NOTCH signaling, IL2-STAT5, and IL6-JAK-STAT3 signaling were more prominent in redox^low LUAD (Figure 3(c)).

Redox^high LUAD displays as an ecosystem lacking lymphocytes infiltration

The TME of LUAD was classified into ten carcinoma ecosystems (CE1-CE10) based on cell states and cellular communities (Figure 4(a), Supplemental Table S3). CE2 was enriched in redox^high LUAD, while CE6 and CE10 were dominantly observed in redox^low LUAD (Figure 4(b)). CE2 was characterized by an ecosystem rich in basal-like epithelial cells, pro-migratory-like fibroblasts and M2-like macrophages, but deficient in lymphocytes (Supplemental Table S3). CE6 was enriched for immune cells (normal-enriched B cells, resting CD4+ T cells, mature DCs, normal-enriched mast cells, M2-like macrophages, and activated neutrophils), stromal cells (normal-enriched endothelial cells and fibroblasts), and normal-enriched epithelial cells. CE10 was characterized as a proinflammatory tumor with abundant leukocytes, including naïve B cells, plasma cells, naïve/central memory CD4+ T cells, naïve/central memory CD8+ T cells, mast cells, monocytes, and cDC1s.

Figure 4. Comparative analysis of cell communities and states across metabolic Phenotypes. (a) Distribution of cell states and carcinoma ecotypes (CEs) within the TCGA LUAD cohort. The horizontal axis denotes tumor samples, while the vertical axis categorizes cell states (S01-S09) and their respective cell types (e.g., CD8+T cells) found in each tumor sample. The heatmap illustrates the proportion of each cell state within its corresponding cell type. CEs were identified by clustering the predominant cell state of each cell type within individual samples. Details on cell state annotation and marker genes are provided in Supplemental Table S3. (b) CE comparison across metabolic phenotypes. (c) Comparison of predominant cell states by cell type across metabolic phenotypes.

Distinct differences in cell state compositions were noted between the two metabolic phenotypes (Figure 4(c)). Metabolic epithelial cells were enriched in redox^high LUAD, while normal-enriched epithelial cells were prevalent in redox^low LUAD. B cells, naïve/central memory CD4+ T cells, and naïve/central memory CD8+ T cells were more abundant in redox^low LUAD, while the abundance of migratory and mature immunogenic DCs, tumor-associated endothelial cells, basal-like epithelial cells, pro-migratory fibroblasts, and M2-like proliferative macrophages was higher in redox^high LUAD (Supplemental Fig S3A). Estimations by ImmuCellAI also indicated lymphocytes enrichment in redox^low LUAD, whereas myeloid cells were more prevalent in redox^high LUAD (Supplemental Fig S3B). Cytokines involved in T cell co-stimulation and the formation of tertiary lymphoid structures were significantly decreased in redox^high LUAD (Supplemental Fig S3C).

Redox^high phenotype is associated with poor prognosis in LUAD

Thirty-nine hub genes from the MEyellow module were initially selected to construct the redox scoring model. Twenty-two genes were retained after LASSO regression (Figure 5(a), Supplemental Table S4). The redox scoring model demonstrated good predictive ability for metabolic phenotypes with an area under the curve (AUC) of 0.89 (sensitivity: 0.84, specificity: 0.80) in the validation set (Figure 5(b)). High redox score was indicative of redox^high phenotype. Immunohistochemistry of tumor tissue validated the expression of genes in redox scoring model (Figure 5(c), Supplemental Table S5). In the GSE31210 LUAD cohort, patients with high redox score showed significantly decreased OS (Figure 5(d), 5-years OS rate: 71.7% vs. 91.1%, p = 0.00033). Similarly, pan-cancer patients with high redox score exhibited decreased OS (Figure 5(e), median OS: 5.1 years vs. 9.5 years, p < 0.0001).

Figure 5. Construction and validation of redox scoring model in LUAD. (a) Tenfold cross-validation of the LASSO Cox regression model for tuning parameter selection, with the horizontal axis representing log(lambda) and the vertical axis representing misclassification error. (b) ROC curve analysis of redox scoring model for predicting metabolic phenotypes in the validation cohort. (c) the immunohistochemistry-based expression of genes from the redox scoring model in LUAD tumor tissues, utilizing data sourced from the Human Protein Atlas, not conducted in-house. (d) Kaplan–Meier analysis of OS between high and low redox score groups in the GSE31210 LUAD cohort. (e) Kaplan–Meier analysis of OS between high and low redox score groups in the TCGA pan-cancer cohort.

ScRNA-seq reveals the Redox^high phenotype is predominantly shaped by malignant cells

The clinical information of the GLCI cohort is listed in Supplemental Table S6. After quality control, transcriptomes from 44,352 cells were analyzed. Cell clusters were annotated as CD4+ T cells, NK/CD8+ T cells, myeloid cells, epithelial cells, B cells, neutrophils, plasma cells, fibroblasts, endothelial cells, and mast cells, based on the expression of canonical marker genes (Figure 6(a-b)). Notably, most metabolic genes were highly expressed in epithelial cells (Figure 6(c)), which also demonstrated the highest redox signature (Figure 6(d)). Epithelial cells were re-clustered into ten subclusters (epi 0 to epi 9), and all subclusters, except for epi 2, were annotated as malignant cells by inferCNV analysis (Figure 6(e-f)). The redox signature was significantly higher in malignant cells, as well as in EGFR wild-type tumors (p < 2.2e-16, Figure 6(g)). The redox signature was positively associated with mTORC1 signaling in tumor cells (Figure 6(h); rho = 0.43, p < 2.2e-16). These findings were also observed in the Hu cohort (Supplemental Fig S4).

Figure 6. ScRNA-seq analysis on the metabolic phenotype shaping mechanism. (a) UMAP plot of the 44,352 cells, visualizing the distribution and clustering of cell types. (b) Dot plot showcasing the expression of marker genes across all identified cell types. (c) Expression of twenty-two genes of redox scoring model in each cell type at the single-cell level. (d) Redox signature in each cell type. (e) Re-clustering of epithelial cells. (f) Results of identifying malignant cells by Infer-CNV analysis. (g) Comparison of redox signature between malignant cells and normal cells or between EGFR wild-type epithelial cells and EGFR-mutant epithelial cells. (h) the correlation between mTORC1 signaling and redox signature.

Redox^high phenotype leads to a diminished presence of tissue-resident memory CD8+T cells

We next focused on the impact of redox activity on CD8+ T cells, the direct anti-tumor lymphocytes. A total of 5,743 CD8+ T cells in the GLCI cohort were re-clustered into nine subclusters, including four central memory CD8+ T cells (Tcm1–4), two effector memory CD8+ T cells (Tem1–2), two tissue-resident memory CD8+ T cells (Trm1–2), and one cytotoxic CD8+ T cell (Tcyt) (Figure 7(a)). Trm1 and Trm2 highly expressed the tissue-resident marker CD103 (ITGAE) and immune checkpoints including LAG3, TIGIT, PDCD1, HAVCR2, and CTLA-4 (Figure 7(b)). No difference was observed in the composition of CD8+ T cells between patients with a high and low redox signature (p > 0.05, Figure 7(c)). However, the abundance of CD8+ Trm, CD8+ Tcm, CD8+ Tem, and CD8+ Tcyt was negatively correlated with the redox signature (Figure 7(d)). An evaluation of the CD8+ Trm signature alongside specific Trm markers (CD69 and CD103) through RNA expression profiling by tissue microarray from the GSE31210 cohort revealed that both the Trm marker CD69 (p = 0.004) and the CD8+ Trm signature (p = 0.031) were significantly diminished in redox^high LUAD (Figure 7(e)). The Nichenet analysis revealed that high redox-signature tumor cells exhibited increased ligand expression of FGA, TIMP1, FURIN, while showing decreased expression of HLA-E, CCL4, CCL2, S100A9, and TYROBP (Figure 7(f)).

Figure 7. The association between redox signature and CD8+T cells. (a) Re-clustering of CD8+T cells. (b) Normalized Z-score expression of marker genes in subtypes of CD8+T cells. (c) Proportions of subtypes of CD8+T cells in individual samples (up) and across redox^high and redox^low (bottom). (d) Correlation between redox signature and abundance of CD8+T cells. (e) Analysis of CD8+Trm signature, expression of specific Trm markers (CD69 and CD103) using tissue RNA microarray of GSE31210GSE cohort without EGFR mutation. (f) Predicted modulation of tumor cells to CD8+T cells by NicheNet. Relative expression of predicted ligands in tumor cells of redox^high vs. redox^low (the left panel) and heatmap showing the interaction potential of top-ranked ligands in tumor cells with corresponding receptors in CD8+T cells (the right panel). (g) Expression of cytotoxic chemokines, MKI67, and CXCL13 in CD8+T cells at pre-treatment and post-treatment of neoadjuvant chemoimmunotherapy. (h) Expression of cytotoxic chemokines, MKI67, and CXCL13 in CD8+T cells between MPR and NMPR after neoadjuvant chemoimmunotherapy. (i) Abundance of CD8+T cells at pre-treatment and post-treatment of neoadjuvant chemoimmunotherapy in Hu cohort. (j) Abundance of CD8+T cells in major pathological response (MPR) and non-major pathological response (NMPR) after neoadjuvant chemoimmunotherapy (post-treatment) in Hu cohort. Trm, tissue-resident memory T cell; Tcyt, cytotoxic T cell; Tem, effector memory T cell; Tcm, central memory T cell.

We further compared the abundance of CD8+ T cells before and after neoadjuvant chemoimmunotherapy in the Hu cohort. We utilized the ProjecTILs²³ to annotate cell states of CD8+T cells in the Hu cohort, with the GLCI cohort’s CD8+T cell states serving as a reference (Supplemental Fig S4E). The expression of cytotoxic chemokines (GZMA, GNLY, PRF1, GZMB) in Trm was found to increase after neoadjuvant chemoimmunotherapy, as shown in Figure 7(g). Interestingly, there was no change observed between patients exhibiting a major pathological response (MPR) and those with a non-major pathological response (NMPR), according to Figure 7(h). CXCL13, a chemokine crucial for the recruitment and activation of CXCR5+ CD8+ T cells, was predominantly expressed in Trm cells and saw enhanced expression after treatment. Moreover, the abundance of Trm was numerically greater in patients post-neoadjuvant chemoimmunotherapy, including in those with MPR, with post-treatment versus pre-treatment levels at 0.051 versus 0.010 (Figure 7(i)) and MPR versus NMPR at 0.070 versus 0.016 (Figure 7(j)).

Redox^high phenotype attenuates the efficacy of immune checkpoint inhibitors

Redox phenotypes were stratified by the redox scoring model: a high redox score indicated a redox^high phenotype while a low redox score indicated a redox^low phenotype. Consistent with previous findings, KEAP1 and STK11 mutations were significantly higher in patients with a high redox score (KEAP1, high vs. low: 50.9% vs. 9.0%, p < 0.001; STK11, high vs. low: 39.6% vs. 20.5%, p = 0.025), whereas EGFR mutations were more prevalent in patients with a low redox score (16.7% vs. 1.9%, p = 0.016) (Figure 8(a)). PD-L1 expression (p = 0.245) and tissue tumor mutation burden (tTMB, p = 0.779) were similar between patients with low and high redox score.

Figure 8. Impact of metabolic phenotypes on response to ICIs in lung adenocarcinoma. (a) Distribution of tTMB, PD-L1 expression, KEAP1/STK11/EGFR/KRAS status in high and low redox score groups. (b) Comparison of response rate, PFS, and OS between high and low redox score groups in patients receiving immunotherapy in Ravi cohort. (c) Comparison of response rate, PFS, and OS between high and low redox score groups in OAK cohort patients receiving atezolizumab (upper) or docetaxel (lower). (d) ROC curve analysis comparing redox score (low vs. high), PD-L1 TPS (<1% vs. >1%), and tTMB (<16 vs. ≥16) in predicting responders to atezolizumab. (e) Changes of redox signature between pre-treatment and post-treatment or between patients achieved NMPR and MPR in neoadjuvant immunochemotherapy cohort. (f) Changes of mTORC1 signaling between pre-treatment and post-treatment or between patients achieved NMPR and MPR in neoadjuvant immunochemotherapy cohort. tTMB, tissue tumor mutation burden; PFS, progression-free survival; OS, overall survival; NMPR, non-major pathological response; MPR, major pathological response.

In the Ravi cohort, LUAD patients with a high redox score demonstrated a significantly decreased response rate (39.1% vs. 70.8%, p = 0.001), shorter PFS (3.3 vs. 14.6 months, p = 0.004), and OS (12.1 vs. 31.2 months, p = 0.022) compared to those with a low redox score (Figure 8(b)). In patients receiving atezolizumab in the OAK cohort, those with a high redox score also showed a decreased response rate (16.5% vs. 36.8%, p = 0.001), shorter PFS (1.6 vs. 3.0 months, p = 0.011), and OS (11.1 vs. 16.8 months, p = 0.095) (Figure 8(c)). For patients treated with docetaxel, no significant differences were observed in response rate (26.4% vs. 37.8%, p = 0.101), PFS (3.8 vs. 4.2 months, p = 0.18), and OS (10.5 vs. 14.3 months, p = 0.24) between the high and low redox score groups, indicating that the redox score may specifically predict immunotherapy outcomes (Figure 8(c)). In a multivariate analysis, a high redox score was an independent risk factor for response rate, PFS, and OS (Supplemental Fig S5). The AUC of the redox score, tTMB, and PD-L1 for predicting responders to atezolizumab was 0.72, 0.60, and 0.53 (Figure 8(d)).

The redox signature was higher in patients who did not achieve MPR (NMPR vs. MPR: 0.509 vs. 0.267, p < 0.001, Figure 8(e)). Moreover, activity of mTORC1 signaling was significantly lower in patients who achieved an MPR (NMPR vs. MPR: 0.321 vs. 0.178, p < 0.001, Figure 8(f)).

Discussion

Metabolic reprogramming within the TME profoundly impacts immune cells. Clarifying the interplay between genetic alterations and metabolic changes, and their collective role in immune evasion, is pivotal for the development of effective therapies to bypass immunotherapy resistance. Previous studies has found a link between high redox signatures,^24–27 mutations in genes associated with redox homeostasis (such as STK11 or KEAP1),^28–30 and adverse outcomes in lung cancer. However, the dynamic between metabolic states, specific genetic alterations, and TME composition in LUAD has not been fully elucidated. Our investigation highlights that a redox^high phenotype, predominantly induced by mutations in KEAP1, STK11, SMARCA4, and NRF2, results in a depletion of tissue-resident memory CD8+ T cells within the TME, undermining ICI efficacy. Hence, the metabolic phenotypes identified in our study offer valuable insights for patient stratification potentially responsive to ICIs or trials targeting metabolic pathways to alleviate oxidative stress.

The metabolic microenvironment of tumor is a focal point in current cancer research, but few studies have assessed whether specific metabolic phenotypes can serve as criteria for tumor classification or ICIs treatment before.³¹ Our study revealed that metabolic phenotype of LUAD is dominantly decided by reductase-oxidative activity, divided into redox^high or redox^low phenotype. The redox^high phenotype accounted for about 30% of LUAD patients. The reductase-oxidative activity involves oxidative stress and reduction reactions.³² The imbalance of reductase-oxidative activity can lead to damage to DNA, proteins, and lipids within immune cells, contributing to immune evasion and tumor progression.³³ In our study, the mutation frequency of KEAP1, STK11, NRF1, and SMARCA4 was up to 72.5% in redox^high LUAD, compared to the mutation frequency of 17.4% in redox^low LUAD. Copy number aberrations of KEAP1, STK11, NRF1, and SMARCA4 were low in this population (4.1% in total population). It suggested that the mutation of KEAP1, STK11, NRF1, and SMARCA4 are predominantly responsible for the redox^high phenotype. But not all KEAP1 mutations were associated with activation of the NRF2 pathway.³⁴ Recent two studies reported the utility of transcriptomic phenotypes to define KEAP1/NRF2-activated or STK11-deficient lung cancers, respectively.^35,36 It is interesting that the gene sets applying for identifying the transcriptomic phenotypes were similar between two studies. It suggests that mutations of different genes can caused similar functional downstream effects.³⁷ It is supported by Arolt et al.‘s finding that SMARCA4 and CUL3 mutations act as mimics of KEAP1/NRF2 mutations.³⁶ In addition, our study found LRRIQ3 mutation was also significantly higher in redox^high LUAD. Taken together, redox^high phenotype represented similar functional downstream effects caused by aberrations of different genes.

Cellular communities and cell states within the TME were identified and validated in previous pan-cancer study,¹⁹ but not illustrated and characterized based on metabolic phenotypes. Here, the TME of redox^high LUAD was displayed as an ecosystem with enriched metabolic epithelial cells and deficient lymphocytes, suggesting a prominent impact of reductase-oxidative activity on local immune activity. Deficiencies in STK11 or activation of the KEAP1/NRF2 pathway have been shown to promote ROS accumulation.^38,39 This accumulation can lead macrophages to adopt an immunosuppressive phenotype, achieved through the activation of the NF-κB signaling pathway, which promotes PD-L1 transcription and the release of immunosuppressive chemokines.⁴⁰ Excess ROS can also selectively upregulate PAC1 in CD8+ T cells via the activation of EGR1.⁴¹ PAC1 recruits the Mi-2β nucleosome-remodeling and histone-deacetylase complex, leading to chromatin remodeling of effector T cells, which ultimately results in the loss of their proliferative and effector capacities.⁴¹ CD8+ Trm has been associated with the response to ICIs.^21,42 In our study, the abundance and cytotoxic chemokines of CD8+ Trm increased in responders to ICIs after treatment. However, a high redox signature was associated with a decrease in CD8+ Trm within the TME, suggesting that redox activity reprogramming caused by genetic aberrations might impair CD8+ Trm. Our cell-cell communication analysis highlighted that FURIN was upregulated in high redox-signature tumor cells. Known for its overexpression in various malignancies and potential as an anti-tumor target, FURIN deficiency in T cells has been linked to impaired immune tolerance and accelerated immune responses.^43,44 In mouse models, FURIN is instrumental in regulating the effector functions of cytotoxic T cells and in the development of CD103+ tissue-resident memory T cells.⁴⁵ Moreover, we noted a reduction in CCL4 expression, critical for recruiting CCR5+ CD8+ T cells,⁴⁶ in high redox-signature tumor cells. These observations suggest that redox-high tumor cells may employ multiple mechanisms to evade immunosurveillance, potentially contributing to resistance to ICIs. Further research is warranted to validate these findings and explore combining ICIs with anti-redox activity therapies as a strategy to counter immunotherapy resistance.

Several limitations exist in this study. Firstly, metabolic phenotyping is currently only applied to LUAD, and applicability to other tumors requires further study. Secondly, nearly 30% of redox^high LUAD could not be explained by the mutation and copy number aberration of KEAP1, STK11, NRF2, or SMARCA4, emphasizing that an elucidation of the impact of other genetic aberrations such as epigenetic modification might be needed in future studies. Thirdly, although the performance of redox scoring model to predict the efficacy of ICIs was validated in two independent cohorts, it still needs validation in a prospective cohort.

In summary, we discover a redox^high phenotype of LUAD that arises from tumor metabolic reprogramming. This redox^high phenotype is predominantly driven by mutations in KEAP1, STK11, and SMARCA4. The phenotype decreased number of tissue-resident memory CD8+ T cells within the TME, and it attenuates the efficacy of ICIs. Our metabolic phenotyping study provides new insights into the resistance mechanisms against ICIs and proposes the redox phenotype as a predictive biomarker for ICI responsiveness.

Acknowledgments

We extend our gratitude to the Genentech Data Access Committee for granting access to RNA-seq and corresponding clinical data from the OAK trial.

Disclosure statement

Prof. Yi-Long Wu reports personal financial interests: consulting and advisory services, speaking engagements of Roche, AstraZeneca, Eli Lilly, Boehringer Ingelheim, Sanofi, MSM, and BMS. Prof. Qing Zhou reports lecture and presentations fees from AstraZeneca, Boehringer Ingelheim, BMS, Eli Lilly, MSM, Pfizer, Roche, and Sanofi outside the submitted work. The remaining authors declare no competing interests that may have influenced the reported work.

Data availability statement

TCGA datasets were available in https://portal.gdc.cancer.gov/. GEO dataset was available at the following URL: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc= GSE31210. Single cell sequencing data of GLCI cohort and Hu cohort was available at the following URL: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc= GSE171145 and https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE207422. Bulk RNA-seq and clinical data of OAK cohort was available upon request at the following URL: https://ega-archive.org/studies/EGAS00001005013. Bulk RNA-seq and clinical data of Ravi cohort was available at the following URL: https://www.nature.com/articles/s41588-023-01355-5. Request for other resources should be directed to the corresponding author.

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/2162402X.2024.2340154

Additional information

Funding

This study was supported by the National Natural Science Foundation of China [Grant No. 82072562, 82102808, 82373349], Guangdong Provincial People’s Hospital Young Talent Project [Grant No. KY0120220136] and the Guangdong Provincial Key Lab of Translational Medicine in Lung Cancer [Grant No. 2017B030314120].