Yi-Pan Zhu1†
Michal A. Rahat
Qiang-Zhe Zhang- 1State Key Laboratory of Medicinal Chemical Biology and College of Pharmacy, Tianjin Key Laboratory of Molecular Drug Research, Nankai University, Tianjin, China
- 2Tianjin Key Laboratory of Extracorporeal Life Support for Critical Diseases, Institute of Hepatobiliary Disease, Tianjin Third Central Hospital, Nankai University, Tianjin, China
- 3Department of Neurosurgery, Tianjin First Central Hospital, School of Medicine, Nankai University, Tianjin, China
- 4Immunotherapy Laboratory, Carmel Medical Center, and the Ruth and Bruce Rappaport Faculty of Medicine, Technion-Israel Institute of Technology, Haifa, Israel
- 5Department of Neurosurgery, Tianjin Huanhu Hospital, School of Medicine, Nankai University, Tianjin, China
Background: T cell exhaustion (TEX) heterogeneity leads to unfavorable immunotherapeutic responses in patients with cancer. Classification of TEX molecular phenotypes is pivotal to overcoming TEX and improving immunotherapies in the clinical setting. Cuproptosis is a novel form of programmed cell death associated with tumor progression. However, the relation between cuproptosis-related genes (CuRGs) and the different TEX phenotypes has not been investigated in lung adenocarcinoma (LUAD).
Methods: Unsupervised hierarchical clustering and principal component analysis (PCA) algorithm were performed to determine CuRGs-related molecular subtypes and scores for patients with LUAD. The tumor immune microenvironment (TIME) landscape in these molecular subtypes and scores was estimated using ESTIMATE and ssGSEA algorithms. Furthermore, TEX characteristics and phenotypes were evaluated in distinct molecular subtypes and scores through GSVA and Spearman correlation analysis. Finally, TIDE scores, immunophenoscore, pRRophetic, GSE78220, and IMvigor210 datasets were employed to appraise the distinguishing capacity of CuRGscore in immunotherapy and pharmacotherapy effectiveness.
Results: We identified three CuRGclusters, three geneClusters, and CuRGscore based on 1012 LUAD transcriptional profiles from five datasets. Compared with other molecular subtypes, CuRGcluster B, geneCluster C, and low-CuRGscore group with good prognosis presented fewer TEX characteristics, including immunosuppressive cells infiltration and TEX-associated gene signatures, signal pathways, checkpoint genes, transcription and inflammatory factors. These molecular subtypes were also responsive in distinguishing TEX phenotype in the terminal, GZMK+, and OXPHOS- TEX subtypes, but not the TCF7+ TEX subtype. Notably, copper importer and exporter, SLC31A1 and ATP7B, were remarkably associated with four TEX phenotypes and nine checkpoint genes such as PDCD1, CTLA4, HAVCR2, TIGIT, LAG3, IDO1, SIGLEC7, CD274, PDCD1LG2, indicating that cuproptosis was involved in the development of TEX and immunosuppressive environment in patients with LUAD. Moreover, CuRGscore was significantly related to the TIDE score, immunophenoscore, and terminal TEX score (Spearman R = 0.62, p < 0.001) to effectively predict immunotherapy and drug sensitivity in both training and external validation cohorts.
Conclusion: Our study demonstrated the extensive effect of cuproptosis on TEX. CuRGs-related molecular subtypes and scores could illuminate the heterogeneity of TEX phenotype as reliable tools in predicting prognosis and directing more effective immunotherapeutic and chemotherapeutic strategies for patients with LUAD.
1 Introduction
Lung cancer is the leading cause of cancer-related deaths worldwide (Sung et al., 2021), with a 5-year survival rate of only 22% (Siegel et al., 2022). Approximately 85% of lung cancer cases are non-small cell lung cancer (NSCLC), of which lung adenocarcinoma (LUAD) are the major subtypes (Molina et al., 2008; Herbst et al., 2018). The carcinogenesis of LUAD is associated with smoking, alcohol consumption, metabolic disorders, and other factors (Alberg et al., 2013), which cause molecular and pathological alterations, genetic mutations, and immune escape in tumor (Zhang et al., 2021). The death rate of LUAD is high due to lack of effective diagnostics and treatment.
Immune checkpoint blockade (ICB) therapy has revolutionized the treatment of advanced NSCLC. It has improved the overall survival of NSCLC patients, but only a few patients have benefited in the clinical setting (Jiang et al., 2021). The objective response rate to PD-1 blockade (Pembrolizumab) is 19.4% for all the patients with a median overall survival of 12.0 months (Garon et al., 2015). In another monotherapy of anti-PD-1, median overall survival of Nivolumab was 9.9 months in 129 patients with an objective response rate of 17% (Gettinger et al., 2015). Functional phenotypes of T cells are considered as key deciders for immunotherapy efficiency.
Tumor-infiltrating CD8+ T cells serve as primary executors that can recognize and kill tumor cells (Durgeau et al., 2018). They encounter immunosuppressive signals in the tumor microenvironment (TME) respond by becoming exhausted, and multiple subsets of T cell exhaustion (TEX) were identified (Crespo et al., 2013; Cai et al., 2020). Four TEX subsets have been defined based on Ly108 and CD69 expression: TEX progenitor 1 (Ly108+ CD69+) and 2 (Ly108+ CD69−), TEX intermediate (Ly108- CD69−), and TEX terminal (Ly108- CD69+) (Beltra et al., 2020). Progenitor and intermediate TEX selectively respond to anti-PD-1 therapy, while terminal TEX cells do not (Blackburn et al., 2008; Miller et al., 2019; Beltra et al., 2020). The heterogeneity of TEX could obstruct the favorable therapeutic response in patients with cancer (Dolina et al., 2021). Therefore, a deeper understanding and microdissection of TEX heterogeneity that would point to possible new biomarkers are crucial to alleviate TEX and improve tumor immunotherapy efficacy and clinical outcomes.
Cuproptosis is recently identified as a new copper-dependent form of regulated cell death, which is different from other known regulated cell death (Tang et al., 2022; Tsvetkov et al., 2022; Wang et al., 2022). Accumulation of intracellular copper induces mitochondrial stress through the aggregation of lipoylated tricarboxylic acid cycle enzymes and the loss of Fe–S cluster proteins, and ultimately causes cuproptosis (Chen et al., 2022; Tsvetkov et al., 2022). Copper (Cu) ionophores such as elesclomol raising intracellular Cu concentration induced cell death and executed antineoplastic activity (Tsvetkov et al., 2022; O'Day et al., 2013). Similarly, artificially overexpressing copper importer SLC31A1 resulted in Cu over-accumulation, thereby triggered cuproptosis in A549 lung cancer cells (Tsvetkov et al., 2022). Whereas, depleting cupric ion with copper chelators reduced PD-L1 expression by both inhibiting PD-L1 transcription and promoting PD-L1 ubiquitination degradation in xenograft mouse (Voli et al., 2020). The evidence suggested that cuproptosis was tightly associated with tumor progression. However, the relationship between cuproptosis-related genes (CuRGs) dissecting molecular subtypes and TEX phenotypes in patients with LUAD was largely unknown.
In the study, we chose thirteen previously identified cuproptosis-related genes (CuRGs) and four gene expression omnibus (GEO) datasets with the same platform GPL570 as candidate genes and data to hold the accuracy of model building. The principal component analysis (PCA) algorithm considers the contribution of each gene in a large gene set rather than several genes to define the CuRGscore for each patient; hence, we employed it together with unsupervised hierarchical clustering to construct three CuRGclusters, three geneClusters, and CuRGscore based on 1012 patients with LUAD in the training cohort. Two independent GEO datasets and two pre-treatment immunotherapy response datasets were used to validate the accuracy of the results. Finally, our results concluded that these molecular subtypes and CuRGscores could distinguish the heterogeneity of TEX phenotype to predict prognosis and immuno- and chemotherapeutic responses reliably in patients with LUAD.
2 Materials and methods
2.1 LUAD data sets and preprocessing
The flow diagram of the study is illustrated in Figure 1. The transcriptome data and clinical information of patients with LUAD were obtained from The Cancer Genome Atlas (TCGA) repositories (https://portal.gdc.cancer.gov). Six microarray datasets, namely, GSE30219, GSE31210, GSE37745, GSE50081, GSE72094, and GSE68465, were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). Information about these datasets is described in detail in Supplementary Table S1. Four GEO datasets (GSE30219, GSE31210, GSE37745, and GSE50081 from the same platform GPL570) and the TCGA-LUAD (Illumina HiSeq) dataset were chosen as training cohort, while two independent GEO datasets (GSE72094 from the GPL96 platform and GSE68465 from the GPL15048 platform) were used as external validation cohorts. Samples without documented survival information were excluded. The fragments per kilobase million (FPKM) values in the TCGA-LUAD dataset were converted to transcripts per kilobase million (TPM) values for our analysis (Conesa et al., 2016; Gustafsson et al., 2021). For GEO datasets, the expression of genes with multiple probes was averaged. Then the “Combat” algorithm of the SVA package was used to correct the batch effect of samples and merge the training data into one dataset (Zhang et al., 2020). Based on the above procedures, 1012 patients with LUAD in the training dataset and 840 in the testing cohort were collected for subsequent analysis.
FIGURE 1. The flow diagram of the study.
2.2 Unsupervised hierarchical clustering for CuRGclusters and geneClusters
Firstly, 13 reported genes were chosen as CuRGs in the study, including 3 genes of copper transporter protein (CuPT; ATP7A, ATP7B, and SLC31A1), 4 components in lipoic acid (LA) pathway (FDX1, LIAS, LIPT1, and DLD), 3 pyruvate dehydrogenase (PDH) complex (DLAT, PDHA1, and PDHB) and other three proteins (MTF1, GLS, and CDKN2A) (Tsvetkov et al., 2022). Based on the expression level of the 13 CuRGs, 1012 patients with LUAD were categorized into distinct molecular subtypes using the “ConsensusClusterPlus” R package. To ensure cluster stability, 1,000 iterations were performed (Wilkerson and Hayes, 2010). The clusters were identified according to the optimal k value of the consensus matrix in the training cohort and were named CuRGclusters. Subsequently, gene expression profiles were analyzed between these CuRGclusters using an empirical Bayesian approach in the “limma” R package (Ritchie et al., 2015). The differentially expressed genes (DEGs) from different CuRGclusters with adjusted p-values <0.001 were considered statistically significant. The univariate Cox regression analysis was used to process the DEGs to extract the expression profile of significant prognostic genes. Based on these prognostic DEGs, 1012 patients with LUAD were dissected again by unsupervised hierarchical clustering. The clusters were identified according to the optimal k value and were termed geneClusters.
2.3 Construction of the CuRGscore
CuRGs-related gene signature was developed using the PCA algorithm based on the prognostic-related DEGs, named CuRGscore. Principal components 1 and 2 were extracted to serve as the signature score using a method similar to GGI (Sotiriou et al., 2006; Zhang et al., 2020). The CuRGscore for each patient was calculated using: CuRGscore = Σ(PC1i + PC2i), where is the expression of these significant prognostic genes. The patients were divided into high- and low-CuRGscore groups according to the optimal cutoff value derived from the Survminer algorithm. Kaplan-Meier survival analysis was executed to evaluate the overall survival (OS) of patients with LUAD (log-rank test, p < 0.001).
2.4 Evaluation of TEX molecular characterization in different molecular subtypes
Six gene signatures of CD8 TEX were chosen from previous studies to investigate the molecular characterization of TEX in different molecular subtypes, including NSCLC (Guo et al., 2018), colorectal cancer (Zhang et al., 2018), hepatocellular carcinoma (Zheng et al., 2017), melanoma (Sade-Feldman et al., 2018), metastatic melanoma (Tirosh et al., 2016) and chronic infection (Bengsch et al., 2018). Four gene signatures of the TEX phenotype were collected from pan-cancer single-cell sequencing data (Zheng et al., 2021). Eight signaling pathways and hallmark genesets associated with TEX were extracted from the Molecular Signatures Database (MSigDB, V7.5.1) (Liberzon et al., 2015). These TEX signatures were further used to conduct gene set variation analysis (GSVA) to evaluate the TEX stage in CuRGclusters, geneClusters, and CuRGscore. Subsequently, checkpoint genes, transcription, and inflammatory factors associated with TEX were further analyzed in different molecular subtypes. Nine checkpoint genes such as PDCD1, CTLA4, HAVCR2, TIGIT, LAG3, IDO1, SIGLEC7, CD274, PDCD1LG2 were chosen as TEX-associated checkpoint genes (Yang et al., 2022). Furthermore, eleven transcription factors, including TOX, TBX21, EOMES, BATF, NFATC1-4, NR4A1-3, as well as seven inflammatory factors such as IL1RN, IL4, IL10, CCL18, TGFB1, IL1B, CCL2 were employed to evaluate the TEX progress (Zheng et al., 2021; Zhang et al., 2022). Finally, the difference in these genes was compared between distinct subtypes using the “limma” R package.
2.5 Correlation analysis for CuRGscore and TEX characterization
To uncover the correlation between CuRGscore and TEX progress, Spearman correlation analysis was conducted between CuRGscore and four TEX phenotypes. A p-value less than 0.05 was considered statistically significant. Moreover, CuRGs were extracted and normalized by the limma R package. The correlation of CuRGs with TEX phenotypes and checkpoint genes was estimated with the Spearman analysis.
2.6 Stratification analyses between CuRGscore and TEX phenotype or immune checkpoint genes
Patients were divided into subgroups according to the four TEX phenotypes or the nine immune checkpoint genes. Subsequently, the Kaplan-Meier survival curve of patients in each subgroup was analyzed to evaluate the stability of CuRGscore prognostic efficacy (log-rank test, p < 0.01).
2.7 Evaluation of TIME in different molecular subtypes
The cell infiltration of LUAD samples with the immune, stromal, and ESTIMATE scores, as well as tumor purity, were evaluated using the ESTIMATE R package (Yoshihara et al., 2013). A higher stromal score or immune score indicates a higher relative content of stromal or immune cells in the TME, while an ESTIMATE score indicates an aggregated TME stromal or immune score. Tumor purity represents the amount of tumor cells in TME. Subsequently, the ssGSEA algorithm was used to assess the abundance of 23 infiltrated immune cells in TIME for each LUAD sample (Barbie et al., 2009).
2.8 Assessment of immunotherapeutic response of CuRGscore
To test the predictive ability of CuRGscore in immunotherapy, the immunophenoscore was downloaded from the cancer immunome atlas (TCIA, https://tcia.at) (Charoentong et al., 2017). Tumor immune dysfunction and exclusion (TIDE) score, another robust predictor for immunotherapy, was estimated by uploading the data normalized to the average to the TIDE web (http://tide.dfci.harvard.edu/) (Jiang et al., 2018; Fu et al., 2020). High immunophenoscore represents better responses to immune checkpoint blockade (ICB) therapy, while a high TIDE score implies more immune evasion. Then, the Spearman correlation was calculated between CuRGscore and immunophenoscore or TIDE score. The distribution of TIDE score or immunophenoscore was also analyzed in the high- and low-CuRGscore groups.
Moreover, the GSE78220 dataset with pre-treatment melanomas undergoing anti-PD-1 therapy and the IMvigor210 dataset with pre-treatment urothelial cancer samples underlying anti-PD-L1 treatment were employed to evaluate the predictive performance of CuRGscore in immunotherapeutic response. GSE78220 (Illumina HiSeq 2000) was downloaded from the GEO database and normalized using the limma package. IMvigor210 dataset (Illumina HiSeq 2500) was obtained from http://research-pub.gene.com/IMvigor210CoreBiologies/. The raw data were normalized by DEseq2 and converted to TPM value to calculate the CuRGscore. The distribution of distinct anti-PD-L1 clinical responses in the CuRGscore subgroups was analyzed with chi-square or Wilcoxon test.
2.9 Assessment of chemotherapeutic response of CuRGscore
To predict the therapeutic performance in the high- and low-CuRGscore patients with LUAD, the R package “pRRophetic” was employed to determine the half maximal inhibitory concentration (IC50) in different CuRGscore groups by ridge regression. The cutoff criteria for the spearman correlation and drug sensitivity were p-value <0.001, and for the |Spearman Cor| was >0.1.
2.10 Statistical analysis
Statistical analyses of data were performed in the R software (version 4.1.3). For comparisons between two groups, unpaired Student’s t-tests were used to analyze variables with a normal distribution, while Wilcoxon tests were used for non-normal distribution. One-way ANOVA and Kruskal–Wallis tests were employed to compare the differences between three or more groups (Hazra and Gogtay, 2016). Spearman correlation analysis then estimated the correlation coefficients. Log-rank tests were used to determine the significance of Kaplan-Meier analysis-generated OS curves. Unless noted otherwise, statistical significance in two-tailed tests was considered when p < 0.05.
3 Results
3.1 TEX characterization of CuRGs-related molecular subtypes in LUAD
Three distinct molecular subtypes were identified according to the optimal consensus matrix k = 3, which included 377 patients in cluster A, 325 in cluster B, and 310 in cluster C (Supplementary Figure S1A and Supplementary Table S2). The survival curve revealed that patients in CuRGcluster B had a prominent survival advantage, while patients in CuRGcluster C exerted an unfavorable prognosis (log-rank test, p < 0.001; Figure 2A). The ESTIMATE analysis showed that CuRGcluster C and A had higher immune scores than CuRGcluster B (Figure 2B). Further, we investigated the proportions of 23 types of infiltrating immune cells using ssGSEA. Consistent with ESTIMATE results, most innate and adaptive immune cells were high in CuRGcluster C, while these cells were infiltrated lowly in CuRGcluster B (Figure 2C). According to the above findings, we speculated that TEX was prevalent in high levels of immune cell infiltration CuRGcluster C with poor prognosis. Hence, TEX-associated signatures and signaling pathways were analyzed to verify this hypothesis. Six TEX signatures derived from published articles were found to be highly expressed in CuRGcluster C compared with CuRGcluster B (Figure 2D). The signaling pathways associated with TEX in CuRGcluster C, including IFNα, IFNγ, TNFα, IL2, IL6, and glycolysis, except TGFβ, were also significantly upregulated compared to those in CuRGcluster B (Figure 2E). The canonical immune checkpoint genes of TEX, such as PDCD1, CD274, PDCD1LG2, CTLA4, HAVCR2, TIGIT, LAG3, IDO1, and SIGLEC7 were highly expressed in CuRGcluster C (Figure 2F). We further investigated the inflammatory and transcription factors involved in TEX development in the three CuRGclusters. While TBX21, BATF, and NFATC2 were significantly upregulated, NFATC3, NFATC4, NR4A1, and NR4A2 were significantly downregulated in CuRGcluster C (Figure 2G). Meanwhile, inflammatory factors, including IL1RN, IL10, CCL18, TGFB1, IL1B and CCL2, were also highly expressed in CuRGcluster C (Figure 2H). The collective data demonstrated that patients in CuRGcluster C with high TEX levels result in poor OS. In contrast, CuRGcluster B patients presented a low TEX level with a better prognosis.
FIGURE 2. TEX and TIME characteristics of each CuRGcluster. (A) Survival analyses for the three CuRGclusters based on 1012 patients from the training cohort (GSE30219, GSE31210, GSE37745, GSE50081, TCGA-LUAD). (B) Immune score analyses using the ESTIMATE algorithm in three CuRGclusters. (C) The abundance of each tumor microenvironment infiltrating cell in CuRGcluster A-C. (D–H) GSVA enrichment analysis between CuRGcluster A-C with TEX-related signatures (D), signal pathways (E), immune checkpoint genes (F), transcription factors (G), or inflammatory factors (H). The line in the box represents the median value, and the asterisk represents the p-value. *p < 0.05; **p < 0.01; ***p < 0.001.
3.2 Identification and TEX characterization of CuRGs-related gene subtypes in LUAD
Based on the optimal consensus matrix k = 3 from the unsupervised clustering of 150 prognostic-related DEGs identified by univariate Cox regression, the whole training cohort was classified into three gene clusters, defined as CuRGs-related gene clusters A-C, respectively (Supplementary Figure S1B and Supplementary Table S3). The relationship between CuRGs-related gene cluster A-C and the clinical characteristics were presented in Figure 3A. The expression of CuRGs was notably different in the three gene clusters, similar to the results with CuRGcluster A-C (Figure 3B). These results demonstrated the existence of three distinct CuRGs-classified patterns in LUAD. The Kaplan-Meier survival analysis suggested that 375 patients of gene cluster C presented the best prognosis, while 254 patients of gene cluster B appeared to have the worst OS (log-rank test, p < 0.001; Figure 3C). Consistent with the ESTIMATE results of CuRGclusters, geneCluster C showed low immune scores compared to geneCluster B and A (Figure 3D). Moreover, the ssGSEA results of infiltrating immune cells reconfirmed the immune score in gene cluster C, which had a low level in most innate and adaptive immune cells (Figure 3E). Together, we found that geneCluster B shows a similar result to that of CuRGcluster C, with the worst prognosis in highly immune-infiltrated subtypes in patients with LUAD. Further, geneCluster B was significantly enriched in infiltrating immune cells associated with TEX, including regulatory T cells (Treg), myeloid-derived suppressor cells (MDSC), and tumor-associated macrophages (TAM) (Figure 3E). All TEX-related signatures, signaling pathways (except TGFβ), and checkpoint genes were highly expressed in CuRGcluster B (Figures 3F–H). TBX21, EOMES, and BATF were considerably upregulated, whereas TOX, NFAT, and NR4A were remarkably downregulated in CuRGcluster B (Figure 3I). Besides IL4, other inflammatory factors were highly expressed in CuRGcluster B compared with CuRGcluster C (Figure 3J). Collectively, these results indicated that CuRGs-related geneCluster could potentially direct T cell exhaustion class to predict the prognosis of patients with LUAD.
FIGURE 3. TEX and TIME characterization of the three geneClusters in LUAD. (A) Unsupervised clustering for 150 prognostic-related DEGs identified by univariate Cox regression classifying patients of training cohort into geneClusters A-C, respectively. The geneClusters, CuRGclusters, projects, tumor stage, survival status, relapse, smoking, gender, and age were used as clinical characteristics. (B) Expression of 13 cuproptosis genes in the three geneClusters. (C) Survival analyses for the three geneClusters based on 1,012 patients. (D) Immune score analyses using the ESTIMATE algorithm in the three geneClusters. (E) The abundance of each TME infiltrating cell in the three geneClusters. (F–J) GSVA enrichment analysis between geneClusters A-C with TEX-related signatures (F), signal pathways (G), checkpoint genes (H), transcription factors (I), or inflammatory factors (J). The line in the box represents the median value, and the asterisk represents the p-value. *p < 0.05; **p < 0.01; ***p < 0.001.
3.3 Construction and TEX characterization of the CuRGscore in LUAD
To accurately quantify the grade of LUAD in individual patient, we performed a PCA algorithm to compute the CuRGscore of each tumor. Based on the cutoff value (1.034) calculated by the survminer algorithm, 1012 patients in the training cohort were divided into low- or high-CuRGscore groups (low = 562 and high = 450; Supplementary Table S4). The alluvial diagram illustrated the attribute changes of patients in CuRGclusters, geneClusters, CuRGscore, and survival status (Figure 4A). CuRGcluster B had the lowest median score, whereas the median score of subtype C was the highest (Figure 4B). On the contrary, geneCluster B also had the highest median score, while the median score of geneCluster C was the lowest (Figure 4C). Kaplan-Meier survival curves indicated that low CuRGscore signifies a better prognosis in patients with LUAD (log-rank test, p < 0.001; Figure 4D). This conclusion was confirmed in the external validation cohorts, including GSE72094, GSE68465, and GSE72094 + GSE68465 datasets (log-rank test, p < 0.001; Figures 4E, F, and Supplementary Figure S2A). We observed high immune, stromal, and ESTIMATE scores and low tumor purity in the high-CuRGscore group compared with the low-CuRGscore group (Figure 4G and Supplementary Figure S2B). Further, the high-CuRGscore group exhibited abundant TEX-related infiltrating immune cells in both training and external validation cohorts (Figure 4H and Supplementary Figure S2C). Moreover, all gene signatures, signaling pathways, and checkpoint genes of TEX were significantly increased in the high-CuRGscore group (Figures 4I–K and Supplementary Figures S2D–F). Notably, TBX21, EOMES, and NFATC2 were significantly upregulated, while TOX, NFATC3, and NR4A were downregulated in the high-CuRGscore group (Figure 4L and Supplementary Figure S2G). Except for IL4, other inflammatory factors were highly expressed in the high-CuRGscore group compared with the low-CuRGscore group (Figure 4M and Supplementary Figure S2H). These findings indicated that CuRGscore could evaluate TEX phenotype and clinical outcomes of patients with LUAD in both training and external validation cohorts.
FIGURE 4. TEX and TIME Characterisation of the CuRGscore in training and external validation cohorts. (A) Alluvial diagram showing the relation of patients in CuRGcluster, geneCluster, CuRGscore, and survival status. (B, C) Differences of CuRGscore among three CuRGclusters (B) and geneClusters (C) in the training cohort (GSE30219, GSE31210, GSE37745, GSE50081, TCGA-LUAD; Kruskal–Wallis test). (D–F) Survival analyses for the high- and low-CuRGscore groups in the training cohort (D) and two independent external cohorts, GSE72094 (E) and GSE68465 (F). (G) Immune, stromal, ESTIMATE, and tumor purity scores in high- and low-CuRGscore groups. (H) The abundance of each TME infiltrating cells in high- and low-CuRGscore groups. (I–M) GSVA enrichment analysis in high- and low-CuRGscore groups with TEX-related signatures (I), signaling pathways (J), checkpoint genes (K), transcription factors (L), and inflammatory factors (M). The line in the box represents the median value, and the asterisk represents the p-value. *p < 0.05; **p < 0.01; ***p < 0.001.
3.4 CuRGscore was positively associated with TEX in LUAD
To further illustrate the correlation between TEX phenotypes and molecular subtypes, TEX scores were calculated for pan-cancer gene signatures of four TEX phenotypes (TCF7+, GZMK+, and OXPHOS- and terminal TEX) (Zheng et al., 2021) for 1012 patients with LUAD. Kruskal–Wallis test showed that terminal, GZMK+, and OXPHOS- TEX were significantly higher in CuRGcluster C and geneCluster B, while TCF7+ TEX was the lowest (Figures 5A, B). Similarly, terminal, GZMK+, and OXPHOS- TEX were higher in the high-CuRGscore group compared to the low-CuRGscore group (Figure 5C). Spearman correlation analyses revealed that CuRGscore was positively associated with the terminal, GZMK+, and OXPHOS- TEX scores, but not with the TCF7+ TEX score. The correlation coefficient between the terminal TEX score and CuRGscore was 0.62 (p < 0.001; Figures 5D–G). These results demonstrated that CuRGscore was tightly related to the terminal, GZMK+, and OXPHOS-TEX cell phenotypes. To uncover the hub genes in the terminal, GZMK+, and OXPHOS- TEX cell phenotypes caused by CuRGs, we further performed Spearman correlation analyses between the CuRGs and four TEX phenotypes or nine checkpoint genes (Figure 5H and Supplementary Table S5). We observed a significant positive correlation between SLC31A1, a copper importer, and the nine checkpoint genes, PDCD1, CTLA4, HAVCR2, TIGIT, LAG3, IDO1, SIGLEC7, CD274, PDCD1LG2. Terminal, GZMK+, and OXPHOS- TEX scores were also positively associated with SLC31A1. HAVCR2 and terminal TEX scores had the strongest positive correlation with SLC31A1 among the nine checkpoint genes and four TEX scores, respectively (Figures 5I, J). Another copper exporter, ATP7B, showed a significant negative association with these TEX markers. HAVCR2 and GZMK+ TEX scores had the highest negative correlation with ATP7B among the nine checkpoint genes and four TEX scores, respectively (Figures 5K, L). Similar results also were found in the external validation cohort (Supplementary Figure S3 and Supplementary Table S6). Hence, the data indicated that cuproptosis, especially SLC31A1 and ATP7B as copper ion importers and exporters, play an important role in the development of TEX in patients with LUAD.
FIGURE 5. Correlation analyses between CuRGscore and TEX in the training cohort. (A, B) Differences in the four TEX phenotype scores among CuRGclusters A-C (A) or geneClusters A-C (B). (C) Differences in the four TEX phenotype scores between high- and low-CuRGscores. (D–G) Spearman correlation analysis between CuRGscores and GZMK+ TEX (D), OXPHOS- TEX (E), Terminal TEX (F), or TCF7+ TEX (G) scores. (H) Spearman correlation analyses between CuRGs and TEX markers, including TEX scores of the four subtypes and nine checkpoint genes. (I) Spearman correlation analysis between HAVCR2 and SLC31A1. (J) Spearman correlation analysis between SLC31A1 and terminal TEX score. (K) Spearman correlation analysis between HAVCR2 and ATP7B. (L) Spearman correlation analysis between ATP7B and GZMK+ TEX score.
3.5 Impact of CuRGscores, TEX phenotypes, and immune checkpoint genes on clinical outcome in LUAD
We performed a Kaplan-Meier survival analysis to examine the effect of TEX phenotypes on the OS of patients with LUAD. We found that the two intermediate phenotypes, GZMK+ and OXPHOS- TEX subtypes, do not exhibit significantly different OS between the high- and low-score groups (log-rank test, p > 0.05; Figures 6A, B). The early phenotype, TCF7+ TEX subtype, showed a favorable prognosis in the high-score group (log-rank test, p = 0.001; Figure 6C), whereas the terminal TEX phenotype demonstrated poor prognosis in the high-score group (log-rank test, p < 0.001; Figure 6D). We further determined whether the CuRGscores could influence the prognosis of patients with LUAD at the same TEX phenotype or similar expression level of immune checkpoint genes. The low-CuRGscore implied a significantly longer OS than the high-CuRGscore, regardless of the TEX phenotypes (log-rank test, p < 0.001; Figures 6E, F, H). Notably, low-CuRGscore indicated a prolonged OS of patients with a high-TCF7+ TEX score than those with a low-TCF7+ TEX score (Figure 6G). Low-CuRGscores were also predictive of a good prognosis for the TEX-related immune checkpoint genes compared with high-CuRGscores, regardless of the expression of TIGIT (log-rank test, p < 0.001; Figure 6I). Similar patterns of clinical outcomes were also found among the four groups divided by CuRGscore and CTLA4, IDO1, PDCD1, LAG3, HAVCR2, SIGLEC7, CD274, or PDCD1LG2 in the training cohort (log-rank test, p < 0.001; Figures 6J–P and Supplementary Figure S4J). We observed similar results in the external validation cohort (GSE72094), that low-CuRGscores indicated better prognosis in patients with LUAD, regardless of the expression level of immune checkpoint genes (log-rank test, p < 0.01; Supplementary Figures S4A–I). These findings further confirmed that CuRGscore was a reliable biological signature for predicting the OS of patients with LUAD.
FIGURE 6. Survival analyses for subgroups of patients stratified by CuRGscores and immune checkpoint genes in the training cohort. (A–D) Survival analysis to examine the OS of patients with LUAD in four TEX subtypes, including GZMK+ (A), OXPHOS- (B), TCF7+ (C), and terminal TEX phenotypes (D). (E–H) Prognosis analysis for patients with LUAD patients stratified into four groups by CuRGscore and four TEX subtypes, including GZMK+ (E), OXPHOS- (F), TCF7+ (G), and terminal TEX phenotypes (H). (I–P) Kaplan-Meier curves analyzed the overall survival of patients with LUAD, stratified into four groups by CuRGscores and checkpoint genes, such as TIGIT (I), CTLA4 (J), IDO1 (K), LAG3 (L), PDCD1 (M), HAVCR2 (N), SIGLEC7 (O) and CD274 (P).
3.6 CuRGscore is a reliable predictor of immunotherapeutic responses in patients with LUAD
Since only 20% of patients with cancer benefit from ICB therapy, an accurate alternate predictor is vital (Mamdani et al., 2022). Hence, we first calculated the TIDE score to evaluate the predictive ability of CuRGscore for immunotherapy in LUAD. We found that the patients in the high-CuRGscore group had a high TIDE score (Figure 7A). Also, these patients had higher microsatellite instability (MSI) and cancer-associated fibroblast (CAF) scores (Figures 7B, C), which is suppressive for the function of immune cells. The external testing cohort exhibited similar results (Figures 7D–F). Moreover, Spearman correlation analysis showed a positive correlation between the CuRGscores and the TIDE scores in both training and external testing cohort (Figures 7G, H). These findings indicated that patients in the low-CuRGscore group respond better to ICB therapy than those in the high-CuRGscore group. We further found a negative correlation between the CuRGscores and immunophenoscores in the TCGA-LUAD dataset (Figure 7I). Group-wise comparison also demonstrated a higher immunophenoscore in the low-CuRGscore group (Figure 7J), strengthening the hypothesis that patients of the low-CuRGscore group respond better to ICB immunotherapy. To test the predictive potential of the CuRGscores in the clinical setting, we determined the relationship between the anti-PD-1/PD-L1 immunotherapeutic response and CuRGscore in the GSE78220 dataset with pre-treatment melanomas undergoing anti-PD-1 therapy (Hugo et al., 2016). We found that patients with complete response (CR) and partial response (PR) were mainly fell into low score area as well as patients with progressive disease (PD) were mostly in high score area (Figure 7K). The Chi-square test indicated that 71% of patients in the low-rated CuRGScore group benefited from anti-PD-1 therapy (CR or PR), compared with the high-rated group (36%, Figure 7L). Finally, we employed pre-treatment urothelial cancer samples from a large phase II trial (IMvigor210) with anti-PD-L1 to evaluate the prognosis of patients. Kaplan-Meier survival curve revealed that patients in the low-CuRGscore group exhibited a better OS from the immune therapy than those in the high-CuRGscore group (Figure 7M). Patients with CR/PR showed a significantly lower CuRGscore than those with stable Disease (SD)/PD (Figure 7N). The percentage of patients with CR/PR in the low-CuRGScore group was about 2.5 times higher than that in the high-CuRGscore group (Figure 7O). Collectively, these results indicated that patients in the low-CuRGscore group have a better chance of benefiting from ICB immunotherapy.
FIGURE 7. The role of CURGscore in ICB immunotherapy. (A–C) The distribution of TIDE (A), MSI (B), and CAF (C) scores between patients in the high- and low-CuRGscore groups in the training cohort. (D–F) The distribution of TIDE (D), MSI (E), and CAF (F) scores between patients in the high- and low-CuRGscore groups in the testing cohort. (G–I) Spearman correlation analysis for TIDE scores or immunophenoscore with CuRGscore in the training (G), external testing (H), and TGCA-LUAD (I) cohorts. (J) Difference of immunophenoscore between patients in the high- and low-CuRGscore groups in the TGCA-LUAD cohort. (K) Distribution of patients with a distinct anti-PD-1 clinical response based on CuRGscores in the GSE78220 dataset. (L) Distribution of distinct anti-PD-1 clinical response in CuRGscore subgroups in the GSE78220 dataset (chi-square test). (M) Kaplan-Meier plot of overall survival for patients of CuRGscore subgroups in the IMvigor210 cohort (log-rank test). (N) The difference in the CuRGscores in anti-PD-L1 clinical response subgroups in the IMvigor210 cohort (Wilcoxon test). (O) The distribution of distinct anti-PD-L1 clinical response in high- and low-CuRGscore subgroups in the IMvigor210 cohort (chi-square test). PR, Partial Response; PD, Progressive Disease; SD, Stable Disease; and CR, Complete Response.
3.7 CuRGscore directs chemotherapeutics for patients with LUAD
To evaluate the predictive potential of CuRGscore on drug sensitivity, pRRophetic analysis was performed to estimate the IC50 of chemotherapeutics with a ridge regression model for high- and low-CuRGscore patients. We found 129 drugs coexisted in the training and external testing cohorts, significantly associated with CuRGscores, including 2 resistant and 127 sensitive drugs (Figure 8A and Supplementary Table S7). Notably, two receptor tyrosine kinase inhibitors, MP470 and SB52334 (ALK5 inhibitor), exhibited drug resistance in patients of the high-CuRGscore group (Figures 8B, C). While patients in the high-CuRGscore group were more sensitive to ALK F1174L and ALK R1275Q inhibitor, CH5424802 (Figure 8D), indicating that these patients had more ALK mutation. An ATP-competitive and selective CDK1 inhibitor, RO-3306, had the highest drug sensitivity and Spearman correlation for patients of the high-CuRGscore group among all 127 sensitivity drugs in both the training and testing cohorts (Figure 8E). Vinblastine, paclitaxel, and docetaxel, disturbing the microtubule cytoskeleton in mitosis, showed better chemotherapeutic potential for patients in the high-CuRGscore group (Figures 8F–H). Protein stability and degradation regulator, lenalidomide, as well as three DNA replication inhibitors, vinblastine, paclitaxel, and Cisplatin, also revealed better chemotherapeutic sensitivity in the high-CuRGscore group (Figures 8I–L). Additionally, the metabolism regulator, AICAR, was also more sensitive in the high-CuRGscore group than in the low-CuRGscore group (Figure 8M). Together, our results suggested that CuRGscore could be used for chemotherapeutic evaluation in patients with LUAD.
FIGURE 8. Chemotherapeutic responses in high- and low-CuRGscore patients with LUAD. (A) Common drug sensitivity and resistance for high- and low-risk patients with LUAD using pRRophetic analysis in the training and external validation cohorts. (B, C) Drug resistance for receptor tyrosine kinases inhibitor, MP470 (B), and ALK5 inhibitor, SB52334 (C) in high-risk patients with LUAD. (D, E) Drug sensitivity for ALK F1174L and ALK R1275Q inhibitor, CH5424802 (D), and CDK1 inhibitor, RO-3306 (E), in high-risk patients with LUAD. (F–H) Chemotherapeutic sensitivity for mitosis inhibitors, such as vinblastine (F), paclitaxel (G) and docetaxel (H) in high-risk patients with LUAD. (I) Chemotherapeutic sensitivity for protein stability and degradation (lenalidomide) in high-risk patients with LUAD. (J–L) Chemotherapeutic sensitivity for DNA replication, including vinblastine (J), paclitaxel (K), and Cisplatin (L) in high-risk patients with LUAD. (M) Drug sensitivity for metabolism regulator, AICAR, in high-risk patients with LUAD.
4 Discussion
Although ICB immunotherapy is a milestone for cancer treatment, only a minority of patients with cancer could benefit from it. For example, PD-L1 expresses in at least 50% of tumor cells in patients with advanced NSCLC, only 45.2% of these patients benefited from the anti-PD-1 antibody, pembrolizumab (Garon et al., 2015). The heterogeneity of TEX was considered an essential determinant of the immunotherapeutic response of cancer (Liu et al., 2020; Inflammatory et al., 2021). Hence, developing a virtual microdissection analytical approach to resolve the heterogeneity of TEX is eminent. Herein, based on 1012 LUAD transcriptional data from five datasets, we employed unsupervised hierarchical clustering and PCA algorithm to construct three CuRGclusters, three geneClusters, and CuRGscore. These molecular subtypes could distinguish the phenotypes of TEX and malignant grades of TME to reliably predict prognosis as well as immunotherapeutic and chemotherapeutic responses in patients with LUAD.
Multiple signatures of TEX have been identified from cancers, including NSCLC (Guo et al., 2018). In our study, identified CuRGs-related clusters and CuRGsores could significantly distinguish these TEX signatures. Furthermore, four gene signatures of TEX phenotype, terminal, GZMK+, OXPHOS- and TCF7+ TEX subtypes, were identified from pan-cancer single-cell sequencing data (Zheng et al., 2021). As progenitors of TEX, the TCF7+ phenotype had stem cell-like capacity, and they differentiated into intermediate phenotype with a loss of TCF1. Studies have shown that the terminal phenotype cannot be functionally reinvigorated through ICB therapy (Beltra et al., 2020). Compared with other molecular subtypes, CuRGcluster C, geneCluster B, and high-CuRGscore group with poor prognoses had a higher fraction in terminal TEX, GZMK+ TEX, and OXPHOS- TEX. Especially, the TEX terminal phenotype was remarkably related to the CuRGscore (R = 0.62 and p < 0.001). These demonstrated that CuRGscore was a good signature to direct TEX phenotype and predict prognosis in LUAD.
Dysfunction in IFN, TNF, IL2, IL6, and glycolysis signaling pathways, as well as high expression of immunosuppressive receptors such as PDCD1, CTLA4, HAVCR2, TIGIT, LAG3, and SIGLEC7, are classic characteristics of TEX (Blank et al., 2019; Saleh et al., 2020). Our study found that CuRGcluster C, geneCluster B, and high-CuRGscore groups represented TEX characteristics with poor prognosis. TEX is also regulated by a network of transcription factors and multiple inflammatory factors (Collins and Henderson, 2014; Fang et al., 2022). For example, high expression levels of EOMES could accelerate the exhaustion of anti-tumor CD8+ T cells (Li et al., 2018). We found that EOMES was also highly expressed in geneCluster B and high-CuRGscore groups with a short OS. Moreover, inflammatory factors, L1RN, IL10, CCL18, TGFB1, IL1B, and CCL2, were enhanced in CuRGcluster C, geneCluster B, and high-CuRGscore groups compared with the other groups. Interestingly, copper importer SLC31A1 was significantly and positively associated with TEX phenotypes and checkpoint genes, while copper exporters, ATP7A and ATP7B, were negatively correlated with TEX phenotypes and checkpoint genes. Recently, two pieces of research also demonstrated that SLC31A1-dependent copper absorption increases the malignancy grade of pancreatic and breast cancers (Yu et al., 2019; Li et al., 2022). We also observed that three anti-cuproptosis genes (CDKN2A, GLS and MTF1) were positively associated with TEX characteristics. Additionally, seven pro-cuproptosis genes (DLAT, DLD, FDX1, LIAS, LIPT1, PDHA1, and PDHB) were also positively associated with TEX characteristics. These results were more significant in the external validation dataset, which reconfirmed the conclusion that cuproptosis was involved in the development of TEX of LUAD.
Besides TEX, other immunosuppressive characteristics are also essential to determine the degree of malignancy of TME, which influences immuno- and chemotherapeutic responses. Studies have shown that Treg, MDSC, and TAM extensively infiltrate advanced NSCLC (Liang et al., 2022). In this study, CuRGcluster C, geneCluster B, and high-CuRGscore subtypes were significantly enriched in infiltrating Treg, MDSC, and TAM. Moreover, CuRGscore was significantly related to the TIDE score and immunophenoscore, as well as can predict immunotherapeutic responses and drug sensitivity. Our study showed that two receptor tyrosine kinase inhibitors, MP470 (c-Kit, PDGFα, Flt3 inhibitor) and SB52334 (ALK5 inhibitor), showed drug resistance in patients with high CuRGscore. Although the inhibitors have not been used in clinic, they are promising antitumor drugs (Mahadevan et al., 2007; Qi et al., 2009). Vinblastine (Velban), paclitaxel (Taxol), docetaxel (Taxotere) and cisplatin (Platinol AQ) as clinical antitumor agents, are more sensitive in high CuRGscore group, which further indicates that CuRGscore can guide the use of drugs. Collectively, these results indicate that cuproptosis is involved in the development of immunosuppressive and malignant TME in patients with LUAD.
5 Conclusion
Our study demonstrated the extensive regulatory mechanisms of cuproptosis on TEX. Our results indicate that CuRGs-related molecular subtypes and scores could distinguish the phenotypes of TEX to reliably predict prognosis and direct more effective immunotherapeutic and chemotherapeutic strategies for patients with LUAD.
Data availability statement
Publicly available datasets were analyzed in this study. This data can be found here: https://portal.gdc.cancer.gov, https://www.ncbi.nlm.nih.gov/geo/, http://research-pub.gene.com/IMvigor210Core Biologies/.
Author contributions
Y-PZ and H-TD collected, analyzed and generated data. Y-PZ and Q-ZZ drafted the manuscript. H-TD, Q-ZZ, XW, MR, and SS reviewed and edited the manuscript. Q-ZZ initiated the idea and designed the experiments. All authors have read and approved the submitted version of this manuscript.
Funding
This work was supported by the National Natural Science Foundation of China (No. 82073233), the Natural Science Foundation of Tianjin (NO. 20JCQNJC00960, 21JCQNJC00130, and 20JCQNJC01880), and the Science and Technology Talents Development program of Tianjin Health Commission (NO. KJ20129).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphar.2023.1146468/full#supplementary-material
References
Alberg, A. J., Brock, M. V., Ford, J. G., Samet, J. M., and Spivack, S. D. (2013). Epidemiology of lung cancer: Diagnosis and management of lung cancer, 3rd ed: American college of chest physicians evidence-based clinical practice guidelines. Chest 143, e1S–e29S. doi:10.1378/chest.12-2345
Barbie, D. A., Tamayo, P., Boehm, J. S., Kim, S. Y., Moody, S. E., Dunn, I. F., et al. (2009). Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. Nature 462, 108–112. doi:10.1038/nature08460
Beltra, J. C., Manne, S., Abdel-Hakeem, M. S., Kurachi, M., Giles, J. R., Chen, Z., et al. (2020). Developmental relationships of four exhausted CD8(+) T cell subsets reveals underlying transcriptional and epigenetic landscape control mechanisms. Immunity 52, 825–841.e8. doi:10.1016/j.immuni.2020.04.014
Bengsch, B., Ohtani, T., Khan, O., Setty, M., Manne, S., O'Brien, S., et al. (2018). Epigenomic-guided mass cytometry profiling reveals disease-specific features of exhausted CD8 T cells. Immunity 48, 1029–1045.e5. doi:10.1016/j.immuni.2018.04.026
Blackburn, S. D., Shin, H., Freeman, G. J., and Wherry, E. J. (2008). Selective expansion of a subset of exhausted CD8 T cells by alphaPD-L1 blockade. Proc. Natl. Acad. Sci. U. S. A. 105, 15016–15021. doi:10.1073/pnas.0801497105
Blank, C. U., Haining, W. N., Held, W., Hogan, P. G., Kallies, A., Lugli, E., et al. (2019). Defining 'T cell exhaustion. Nat. Rev. Immunol. 19, 665–674. doi:10.1038/s41577-019-0221-9
Cai, M. C., Zhao, X., Cao, M., Ma, P., Chen, M., Wu, J., et al. (2020). T-cell exhaustion interrelates with immune cytolytic activity to shape the inflamed tumor microenvironment. J. Pathol. 251, 147–159. doi:10.1002/path.5435
Charoentong, P., Finotello, F., Angelova, M., Mayer, C., Efremova, M., Rieder, D., et al. (2017). Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade. Cell Rep. 18, 248–262. doi:10.1016/j.celrep.2016.12.019
Chen, L., Min, J., and Wang, F. (2022). Copper homeostasis and cuproptosis in health and disease. Signal Transduct. Target Ther. 7, 378. doi:10.1038/s41392-022-01229-y
Collins, M. H., and Henderson, A. J. (2014). Transcriptional regulation and T cell exhaustion. Curr. Opin. HIV AIDS 9, 459–463. doi:10.1097/COH.0000000000000091
Conesa, A., Madrigal, P., Tarazona, S., Gomez-Cabrero, D., Cervera, A., McPherson, A., et al. (2016). A survey of best practices for RNA-seq data analysis. Genome Biol. 17, 13. doi:10.1186/s13059-016-0881-8
Crespo, J., Sun, H., Welling, T. H., Tian, Z., and Zou, W. (2013). T cell anergy, exhaustion, senescence, and stemness in the tumor microenvironment. Curr. Opin. Immunol. 25, 214–221. doi:10.1016/j.coi.2012.12.003
Dolina, J. S., Van Braeckel-Budimir, N., Thomas, G. D., and Salek-Ardakani, S. (2021). CD8(+) T cell exhaustion in cancer. Front. Immunol. 12, 715234. doi:10.3389/fimmu.2021.715234
Durgeau, A., Virk, Y., Corgnac, S., and Mami-Chouaib, F. (2018). Recent advances in targeting CD8 T-cell immunity for more effective cancer immunotherapy. Front. Immunol. 9, 14. doi:10.3389/fimmu.2018.00014
Fang, L., Liu, K., Liu, C., Wang, X., Ma, W., Xu, W., et al. (2022). Tumor accomplice: T cell exhaustion induced by chronic inflammation. Front. Immunol. 13, 979116. doi:10.3389/fimmu.2022.979116
Fu, J., Li, K., Zhang, W., Wan, C., Zhang, J., Jiang, P., et al. (2020). Large-scale public data reuse to model immunotherapy response and resistance. Genome Med. 12, 21. doi:10.1186/s13073-020-0721-z
Garon, E. B., Rizvi, N. A., Hui, R., Leighl, N., Balmanoukian, A. S., Eder, J. P., et al. (2015). Pembrolizumab for the treatment of non-small cell lung cancer. N. Engl. J. Med. 372, 2018–2028. doi:10.1056/NEJMoa1501824
Gettinger, S. N., Horn, L., Gandhi, L., Spigel, D. R., Antonia, S. J., Rizvi, N. A., et al. (2015). Overall survival and long-term safety of Nivolumab (Anti-Programmed death 1 antibody, BMS-936558, ONO-4538) in patients with previously treated advanced non-small cell lung cancer. J. Clin. Oncol. 33, 2004–2012. doi:10.1200/JCO.2014.58.3708
Guo, X., Zhang, Y., Zheng, L., Zheng, C., Song, J., Zhang, Q., et al. (2018). Global characterization of T cells in non-small cell lung cancer by single-cell sequencing. Nat. Med. 24, 978–985. doi:10.1038/s41591-018-0045-3
Gustafsson, J., Robinson, J., Nielsen, J., and Pachter, L. (2021). Butterfly: Addressing the pooled amplification paradox with unique molecular identifiers in single-cell RNA-seq. Genome Biol. 22, 174. doi:10.1186/s13059-021-02386-z
Hazra, A., and Gogtay, N. (2016). Biostatistics series module 3: Comparing groups: Numerical variables. Indian J. dermatology 61, 251–260. doi:10.4103/0019-5154.182416
Herbst, R. S., Morgensztern, D., and Boshoff, C. (2018). The biology and management of non-small cell lung cancer. Nature 553, 446–454. doi:10.1038/nature25183
Hugo, W., Zaretsky, J. M., Sun, L., Song, C., Moreno, B. H., Hu-Lieskovan, S., et al. (2016). Genomic and transcriptomic features of response to anti-PD-1 therapy in metastatic melanoma. Cell 165, 35–44. doi:10.1016/j.cell.2016.02.065
Inflammatory, C., Wandmacher, A. M., Mehdorn, A. S., and Sebens, S. (2021). The heterogeneity of the tumor microenvironment as essential determinant of development, progression and therapy response of pancreatic cancer. Cancers (Basel) 13, 4932. doi:10.3390/cancers13194932
Jiang, P., Gu, S., Pan, D., Fu, J., Sahu, A., Hu, X., et al. (2018). Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat. Med. 24, 1550–1558. doi:10.1038/s41591-018-0136-1
Jiang, Z., Zhou, Y., and Huang, J. (2021). A combination of biomarkers predict response to immune checkpoint blockade therapy in non-small cell lung cancer. Front. Immunol. 12, 813331. doi:10.3389/fimmu.2021.813331
Li, J., He, Y., Hao, J., Ni, L., and Dong, C. (2018). High levels of eomes promote exhaustion of anti-tumor CD8(+) T cells. Front. Immunol. 9, 2981. doi:10.3389/fimmu.2018.02981
Li, L., Li, L., and Sun, Q. (2022). High expression of cuproptosis-related SLC31A1 gene in relation to unfavorable outcome and deregulated immune cell infiltration in breast cancer: An analysis based on public databases. BMC Bioinforma. 23, 350. doi:10.1186/s12859-022-04894-6
Liang, J., Bi, G., Shan, G., Jin, X., Bian, Y., and Wang, Q. 2022, Tumor-associated regulatory T cells in non-small cell lung cancer: Current advances and future perspectives. J. Immunol. Res. 2022, 4355386. doi:10.1155/2022/4355386
Liberzon, A., Birger, C., Thorvaldsdottir, H., Ghandi, M., Mesirov, J. P., and Tamayo, P. (2015). The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst. 1, 417–425. doi:10.1016/j.cels.2015.12.004
Liu, F., Liu, W., Sanin, D. E., Jia, G., Tian, M., Wang, H., et al. (2020). Heterogeneity of exhausted T cells in the tumor microenvironment is linked to patient survival following resection in hepatocellular carcinoma. Oncoimmunology 9, 1746573. doi:10.1080/2162402X.2020.1746573
Mahadevan, D., Cooke, L., Riley, C., Swart, R., Simons, B., Della Croce, K., et al. (2007). A novel tyrosine kinase switch is a mechanism of imatinib resistance in gastrointestinal stromal tumors. Oncogene 26, 3909–3919. doi:10.1038/sj.onc.1210173
Mamdani, H., Matosevic, S., Khalid, A. B., Durm, G., and Jalal, S. I. (2022). Immunotherapy in lung cancer: Current landscape and future directions. Front. Immunol. 13, 823618. doi:10.3389/fimmu.2022.823618
Miller, B. C., Sen, D. R., Al Abosy, R., Bi, K., Virkud, Y. V., LaFleur, M. W., et al. (2019). Subsets of exhausted CD8(+) T cells differentially mediate tumor control and respond to checkpoint blockade. Nat. Immunol. 20, 326–336. doi:10.1038/s41590-019-0312-6
Molina, J. R., Yang, P., Cassivi, S. D., Schild, S. E., and Adjei, A. A. (2008). Non-small cell lung cancer: Epidemiology, risk factors, treatment, and survivorship. Mayo Clin. Proc. 83, 584–594. doi:10.4065/83.5.584
O'Day, S. J., Eggermont, A. M., Chiarion-Sileni, V., Kefford, R., Grob, J. J., Mortier, L., et al. (2013). Final results of phase III SYMMETRY study: Randomized, double-blind trial of elesclomol plus paclitaxel versus paclitaxel alone as treatment for chemotherapy-naive patients with advanced melanoma. J. Clin. Oncol. 31, 1211–1218. doi:10.1200/JCO.2012.44.5585
Qi, W., Cooke, L. S., Stejskal, A., Riley, C., Croce, K. D., Saldanha, J. W., et al. (2009). MP470, a novel receptor tyrosine kinase inhibitor, in combination with Erlotinib inhibits the HER family/PI3K/Akt pathway and tumor growth in prostate cancer. BMC cancer 9, 142. doi:10.1186/1471-2407-9-142
Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic acids Res. 43, e47. doi:10.1093/nar/gkv007
Sade-Feldman, M., Yizhak, K., Bjorgaard, S. L., Ray, J. P., de Boer, C. G., Jenkins, R. W., et al. (2018). Defining T cell states associated with response to checkpoint immunotherapy in melanoma. Cell 175, 998–1013.e20. doi:10.1016/j.cell.2018.10.038
Saleh, R., Taha, R. Z., Toor, S. M., Sasidharan Nair, V., Murshed, K., Khawar, M., et al. (2020). Expression of immune checkpoints and T cell exhaustion markers in early and advanced stages of colorectal cancer. Cancer Immunol. Immunother. 69, 1989–1999. doi:10.1007/s00262-020-02593-w
Siegel, R. L., Miller, K. D., Fuchs, H. E., and Jemal, A. (2022). Cancer statistics, 2022. CA a cancer J. Clin. 72, 7–33. doi:10.3322/caac.21708
Sotiriou, C., Wirapati, P., Loi, S., Harris, A., Fox, S., Smeds, J., et al. (2006). Gene expression profiling in breast cancer: Understanding the molecular basis of histologic grade to improve prognosis. J. Natl. Cancer Inst. 98, 262–272. doi:10.1093/jnci/djj052
Sung, H., Ferlay, J., Siegel, R. L., Laversanne, M., Soerjomataram, I., Jemal, A., et al. (2021). Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA a cancer J. Clin. 71, 209–249. doi:10.3322/caac.21660
Tang, D., Chen, X., and Kroemer, G. (2022). Cuproptosis: A copper-triggered modality of mitochondrial cell death. Cell Res. 32, 417–418. doi:10.1038/s41422-022-00653-7
Tirosh, I., Izar, B., Prakadan, S. M., Wadsworth, M. H., Treacy, D., Trombetta, J. J., et al. (2016). Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq. Sci. (New York, N.Y.) 352, 189–196. doi:10.1126/science.aad0501
Tsvetkov, P., Coy, S., Petrova, B., Dreishpoon, M., Verma, A., Abdusamad, M., et al. (2022). Copper induces cell death by targeting lipoylated TCA cycle proteins. Sci. (New York, N.Y.) 375, 1254–1261. doi:10.1126/science.abf0529
Voli, F., Valli, E., Lerra, L., Kimpton, K., Saletta, F., Giorgi, F. M., et al. (2020). Intratumoral copper modulates PD-L1 expression and influences tumor immune evasion. Cancer Res. 80, 4129–4144. doi:10.1158/0008-5472.CAN-20-0471
Wang, Y., Zhang, L., and Zhou, F. (2022). Cuproptosis: A new form of programmed cell death. Cell. Mol. Immunol. 19, 867–868. doi:10.1038/s41423-022-00866-1
Wilkerson, M. D., and Hayes, D. N. (2010). ConsensusClusterPlus: A class discovery tool with confidence assessments and item tracking. Bioinforma. Oxf. Engl. 26, 1572–1573. doi:10.1093/bioinformatics/btq170
Yang, M., Lin, C., Wang, Y., Chen, K., Zhang, H., and Li, W. (2022). Identification of a cytokine-dominated immunosuppressive class in squamous cell lung carcinoma with implications for immunotherapy resistance. Genome Med. 14, 72. doi:10.1186/s13073-022-01079-x
Yoshihara, K., Shahmoradgoli, M., Martínez, E., Vegesna, R., Kim, H., Torres-Garcia, W., et al. (2013). Inferring tumour purity and stromal and immune cell admixture from expression data. Nat. Commun. 4, 2612. doi:10.1038/ncomms3612
Yu, Z., Zhou, R., Zhao, Y., Pan, Y., Liang, H., Zhang, J. S., et al. (2019). Blockage of SLC31A1-dependent copper absorption increases pancreatic cancer cell autophagy to resist cell death. Cell Prolif. 52, e12568. doi:10.1111/cpr.12568
Zhang, B., Wu, Q., Li, B., Wang, D., Wang, L., and Zhou, Y. L. (2020). m(6 A regulator-mediated methylation modification patterns and tumor microenvironment infiltration characterization in gastric cancer. Mol. Cancer 19, 53. doi:10.1186/s12943-020-01170-0
Zhang, D., Jiang, Q., Ge, X., Shi, Y., Ye, T., Mi, Y., et al. (2021). RHOV promotes lung adenocarcinoma cell growth and metastasis through JNK/c-Jun pathway. Int. J. Biol. Sci. 17, 2622–2632. doi:10.7150/ijbs.59939
Zhang, L., Yu, X., Zheng, L., Zhang, Y., Li, Y., Fang, Q., et al. (2018). Lineage tracking reveals dynamic relationships of T cells in colorectal cancer. Nature 564, 268–272. doi:10.1038/s41586-018-0694-x
Zhang, Z., Chen, L., Chen, H., Zhao, J., Li, K., Sun, J., et al. (2022). Pan-cancer landscape of T-cell exhaustion heterogeneity within the tumor microenvironment revealed a progressive roadmap of hierarchical dysfunction associated with prognosis and therapeutic efficacy. EBioMedicine 83, 104207. doi:10.1016/j.ebiom.2022.104207
Zheng, C., Zheng, L., Yoo, J. K., Guo, H., Zhang, Y., Guo, X., et al. (2017). Landscape of infiltrating T cells in liver cancer revealed by single-cell sequencing. Cell 169, 1342–1356.e16. doi:10.1016/j.cell.2017.05.035
Keywords: cuproptosis, T cell exhaustion, tumor immune microenvironment, lung adenocarcinoma, immunotherapy
Citation: Zhu Y-P, Deng H-T, Wang X, Rahat MA, Sun S and Zhang Q-Z (2023) Cuproptosis-related molecular subtypes direct T cell exhaustion phenotypes and therapeutic strategies for patients with lung adenocarcinoma. Front. Pharmacol. 14:1146468. doi: 10.3389/fphar.2023.1146468
Received: 17 January 2023; Accepted: 29 March 2023;
Published: 11 April 2023.
Edited by:
Zhijie Xu, Xiangya Hospital, Central South University, ChinaReviewed by:
Allen K. L. Cheung, Hong Kong Baptist University, Hong Kong SAR, ChinaChangbin Ke, Hubei University of Medicine, China
Copyright © 2023 Zhu, Deng, Wang, Rahat, Sun and Zhang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Qiang-Zhe Zhang, zhangqiangzhe@nankai.edu.cn
†These authors have contributed equally to this work