Shunyao Xia1
Youcheng Xiu- 1Department of Urology, First Affiliated Hospital of Harbin Medical University, Harbin, China
- 2School of Bioinformatics Science and Technology, Harbin Medical University, Harbin, China
Background: Clear cell renal cell carcinoma (ccRCC) is a malignancy with a high incidence rate and poor prognosis worldwide. Copper ionophore–induced death (CID) plays an important role in cancer progression.
Methods: One training and three validation datasets were acquired from TCGA, GEO and ArrayExpress. K-means clustering was conducted to identify the CID subtypes. The ESTIMATE and CIBERSORT algorithms were employed to illustrate the immune microenvironment of ccRCC. LASSO Cox regression was applied to construct the CID feature-based prognostic model. The immunotherapy cohort was acquired from the literature to explore the potential risk scores for predicting immunotherapy responsiveness.
Results: Two CID-related cancer subtypes of ccRCC were identified that performed different immune microenvironment characteristics and prognosis. Based on the identified subtypes, we analyzed the biological heterogeneity and constructed a gene prognostic model. The prognostic model performed well in one training dataset, three validation datasets, and different clinical pathological groups. The prognostic model has a good potential for predicting cancer immune features and immunotherapy responsiveness.
Conclusion: CID plays an important role in the tumor microenvironment progression of ccRCC. The robust gene prognostic model developed can help predict cancer prognosis, immune features, and immunotherapy responsiveness.
Introduction
Renal cell carcinoma (RCC) is the major prevalent urinary system malignancy, with more than 430000 cases diagnosed worldwide in 2020 (Sung et al., 2021). Among RCC subtypes, clear cell renal cell carcinoma (ccRCC) is the most common and lethal form (Jonasch et al., 2014). The current progression of comprehensive therapy strategies for ccRCC, such as tyrosine kinases inhibitors, mTOR inhibitors and immune checkpoint blockades, has significantly improved the prognosis of patients (Motzer et al., 2014, 2015; Choueiri et al., 2017). However, a non-negligible rate of patients remains non-responsive to cancer therapy and with severe side effects (Kennedy and Salama, 2020; Braun et al., 2021). Furthermore, approximately 30% of patients with ccRCC are diagnosed with metastatic cancer, for which effective therapy strategies are limited (Patard et al., 2011). Consequently, ccRCC is a major global public health concern. Therefore, given the high incidence rate and poor prognosis of ccRCC, developing a robust prognostic model to assist patient prognosis evaluation and reveal the underlying heterogeneity mechanism of ccRCC is urgently in demand.
Redox-active copper plays an essential role in maintaining cell homeostasis and takes part in various biological processes, including energy metabolism, biosynthesis, and antioxidant defense (Tsang et al., 2021). Though copper is indispensable for the normal physiological activity of cells, it can be cytotoxic. In 2022, Tsvetkov et al. (2022) revealed an unexpected cell death pattern triggered by copper in a tricarboxylic acid (TCA) cycle metabolism-related mechanism-copper ionophore–induced death (CID). Meanwhile, copper-related drugs, such as copper chelation, have great potential to be developed as a clinic anti-cancer therapy (Yin et al., 2016; Lopez et al., 2019). ccRCC is a cancer type characterized by significant TCA metabolic heterogeneity (Wettersten et al., 2017). As CID may contribute to the heterogeneous ccRCC formation, CID status may be developed as an indicator of the prognosis of patients with ccRCC.
This study conducted an integrated analysis to illustrate the multi-omics features of CID-related genes in ccRCC and identified two CID subtypes of ccRCC. Then, we analyzed the microenvironment heterogeneity across the two CID subtypes. Based on the two identified CID subtypes, we developed a prognostic model to predict patient prognosis, immune characteristics, and immunotherapy responsiveness using a machine-learning method. Our research presented an overview of the regulatory function of CID during ccRCC progression and developed a robust CID-based model to help evaluate the prognosis and immunotherapy suitability during clinical practice.
Methods and materials
Data collection
The genomic data and corresponding clinical information on kidney clear cell carcinoma (KIRC) of The Cancer Genome Atlas (TCGA) were downloaded from the University of California Santa Cruz (UCSC) Xena online tool (https://xenabrowser.net/). LOWESS normalized gene expression profile and quantile normalized, and the log2 transformed gene expression profile of KIRC (GSE29609 and GSE22541) with corresponding clinical information were downloaded from Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo). The log2 quantile normalized expression data of 101 KIRC samples were downloaded from ArrayExpress (https://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-1980/), and clinical information was obtained from Sato et al. (Sato et al., 2013). The expression and clinical data of IMvigor210 trial were accessed with R package “IMvigor210CoreBiologies”. IMvigor210 was a single-arm phase Ⅱ study investigating the anti-PD-L1 antibody agent atezolizumab in patients with metastatic urothelial cancer (mUCC) (Mariathasan et al., 2018). The KIRC samples of TCGA was used as training dataset due to large sample size for statistical accuracy, and complete gene expression and clinical information (eg. stage, grade, survival time). For validation, we chosen the gene expression datasets of KIRC samples with survival time and survival status or disease-free survival time. The samples in IMvigor210 dataset has responsive information for immunotherapy.
Variation and expression correlation
Genes involved in copper ionophore–induced death (CID) were obtained from a study by Tsvetkov et al. (2022), including FDX1, LIAS, LIPT1, DLD, DLAT, PDHA1, PDHB, MTF1, GLS, and CDKN2A. Somatic mutation, copy number variation (CNV) alterations, and differential expression between tumor samples and normal samples of CID genes were demonstrated. The prognostic value of CID genes was analyzed with a univariable Cox proportional hazards regression model (Supplementary Table S1). Co-expression status of CID genes was analyzed by Pearson correlation analysis (Supplementary Table S2), and the correlation network was visualized using Cytoscape software.
Identification of CID subtypes
K-means clustering is an unsupervised learning algorithm that groups data based on each point euclidean distance to a central point called centroid. K-means clustering was performed to identify two CID subtypes based on CID gene expression by R package “pheatmap”. Finally, a total of 197 samples were grouped into “Subtype A” and 329 samples were grouped into “Subtype B” Principle component analysis (PCA) was applied to explore the difference between Subtypes A and B based on CID gene expression. Kaplan-Meier survival analysis and log-rank test were used to analyze the difference in overall survival (OS) among the two subtypes.
Analysis of tumor immune infiltration microenvironment
The ESTIMATE algorithm was applied to evaluate the immune and stromal scores of each KIRC sample in TCGA using R package “estimate.” The proportion of infiltration of 22 immune cells for TCGA KIRC samples was inferred with CIBERSORT algorithm using the web-based analytical tool (https://cibersort.stanford.edu/) (Newman et al., 2015). CIBERSORT estimates the abundances of specific cell types in a mixed cell population using a gene expression-based approach. We focused on mRNA expression of five immune checkpoints, including PD-1, PD-L1, CTLA4, CD47 and BTLA. The immune cytolytic activity (CYT) was calculated as the mean of GZMA and PRF1 expression according to Rooney et al. (Rooney et al., 2015). A one-sided Wilcoxon rank-sum test was used to analyze the differences between subtypes.
Differentially expressed genes between subtypes and functional analysis
A total of 1448 DEGs with |log2FC| > 1 and FDR <0.001 were identified using R package “edgeR” between subtypes (Robinson et al., 2010). Pathway and process enrichment analysis for 344 DEGs with |log2FC| > 2 and FDR <0.001 was performed using the Metascape web-based tool (https://metascape.org/gp/index.html), including many ontology sources such as KEGG Pathway, GO Biological Processes, and Reactome Gene Sets and Canonical Pathways (Supplementary Table S3) (Zhou et al., 2019).
Construction of the prognostic model
A univariable Cox proportional hazards regression model was performed to identify prognostic DEGs. A total of 80 DEGs with p < 0.001 were selected. The least absolute shrinkage and selection operator (LASSO) method was used for significant prognostic DEGs selection in a Cox regression model by fitting a generalised linear model via penalised maximum likelihood. We analysed the lambda value (λ) using the 10 fold cross-validation, between λmin that gives minimum mean cross-validated error or λ1se, that gives a model such that standard error (SE) is within one standard error of the minimum. The process was conducted using R package “glmnet” (Engebretsen and Bohlin, 2019). Finally, a risk score formula was calculated by considering the expression of 17 optimized genes and correlation estimated multivariate Cox regression coefficients using R package “survival” (Supplementary Table S4). The risk score was calculated as follows:
Survival analysis
Patients were classified according to the median of risk score. The log-rank test was used to assess the survival time difference between high-risk and risk score patients using R package “survival.” Additionally, a stratified analysis was performed to determine whether the risk score retained its predictive ability in different subgroups according to gender, age, T stage, N stage, M stage, tumor stage, and tumor grade. Kaplan-Meier plots were used to present the results. Chi-square tests explored the relationships between the risk score and clinical characteristics (Table 1).
TABLE 1. Baseline characteristics of patients in TCGA KIRC cohort
Statistical analysis
A one-sided Wilcoxon rank-sum test was used to test the discrepancy between CID subtypes or high and low-risk groups. Patients were divided into high risk and low risk groups according to the median of risk score. All statistical analyses were performed using R version 4.1.2. p < 0.05 was considered statistically significant.
Results
Multi-omics level alterations of CID genes in KIRC
The analytical process in this study is illustrated in Supplementary Figure S1. We first explored the landscape of variation in CID genes in genome and transcriptome from TCGA KIRC samples. The incidence of somatic mutations of CID genes is shown in Figure 1A. Among them, DLD had the highest mutation frequency (27%), followed by LIAS, MTF1, GLS, PDHA1, and PDHB. The locations of CNV alterations in CID genes on their respective chromosomes and the expression of CID genes are shown in Figure 1B. PDHB showed the highest frequency of CNV deletion, followed by CDKN2A, MTF1, and GLS. GLS showed the highest frequency of CNV amplification, followed by DLD and LIAS. In addition, we explored the expression levels of CID genes between tumor and normal tissues (Figure 1C). In total, 8 (80%) CID genes showed differential expression, CDKN2A showed significant upregulation, and seven CID genes showed significant downregulation in the tumor samples (Figure 1D, p < 0.05, one-sided Wilcoxon rank-sum test).
FIGURE 1. Genetic and transcriptional alterations of CID genes in KIRC. (A) Mutation frequencies of CID genes in KIRC patients of TCGA cohort. (B) Locations of CNV alterations in CID genes on 23 chromosomes and distribution of expression. (C) Frequencies of CNV amplification and deletion of CID genes in TCGA KIRC cohorts. (D) Differential expression of CID genes between tumor and normal samples.
Tumor classification based on CID genes
We explored the prognostic value of CID genes with a univariable Cox proportional hazards regression model. All CID genes were predicted as favorable factors except CDKN2A (Figure 2A). To explore expression correlation among CID genes, we constructed a co-expression network; the thickness of edges means a significance level (Figure 2A). The network indicated a close connection among CID genes.
FIGURE 2. Identification of CID subtypes. (A) Co-expression among CID genes in TCGA KIRC cohorts. The line thickness indicate the strength of the correlation. The color of nodes mean prognostic factors of CID genes. (B) Two heterogeneous subtypes (subtype A and subtype B) were identified according to unsupervised K-means clustering. (C) PCA analysis showing a remarkable difference in expression of CID genes between subtypes. (D) Overall survival analysis between subtype A and subtype B.
To analyze the heterogeneity of KIRC, K-means clustering algorithm was used to identify two CID subtypes based on the expression of CID genes (Figure 2B). The PCA revealed that TCGA KIRC samples had distinctive expression patterns of CID genes between two subtypes (Figure 2C). Next, we explored the difference in prognosis between two subtypes; individuals in Subtype A had significantly worse OS when compared with those in Subtype B (Figure 2D, p = 2.8 E−07, log-rank test).
Characterization of the immune microenvironment between subtypes
The tumor purity distinctions between subtypes, the stromal score, immune score, and ESTIMATE score in Subtype A were significantly higher than those in Subtype B (Figure 3A, p < 0.05, one-sided Wilcoxon rank-sum test). Then, we analyzed the differential expression of five immune checkpoints. PD-L1 expression was significantly higher in subtype A compared with that of Subtype B, and the expression of CTLA4 and PD-1 were significantly lower in Subtype A than those in Subtype B (Figure 3B, p < 0.05). We also evaluated the distinction of immune cells between two subtypes. According to CIBERSORT algorithm, infiltration of “Macrophages M0,” “NK cells activated,” “Plasma cells,” “T cells CD,” “T cells follicular helper,” and “T cells regulatory (Tregs)” were higher in the Subtype A than those in Subtype B (Figure 3C, p < 0.05). Meanwhile, “Dendritic cells resting,” “Eosinophils,” “Macrophages M1,” “Macrophages M2,” “Mast cells resting,” “Monocytes,” and “T cells CD4 memory resting” had significantly lower infiltration in Subtype A compared with Subtype B (Figure 3D, p < 0.05). In addition, we evaluated CYT for KIRC samples, and CYT score in Subtype A was higher than that in Subtype B (Figure 3E, p < 0.05).
FIGURE 3. Distribution of TME between subtype A and subtype B. (A) Distribution of ESTIMATE score in two subtypes. (B) Expression levels of five immune checkpoints between two subtypes. (C,D) Abundance of infiltrating immune cell types between two subtypes. (E) Distribution of immune CYT score between two subtypes.
Identification of DEGs and construction of the prognostic model
To explore the potential biological behavior of CID subtypes, we identified 1448 DEGs between Subtypes A and B (Figure 4A, |log2FC| > 1, FDR <0.001). Pathway and process enrichment analysis for 344 DEGs with |log2FC| > 2 and FDR <0.001 was performed using Metascape tool. DEGs were significantly enriched in “NABA MATRISOME ASSOCIATED,” “acute-phase response,” “Complement and coagulation cascades,” “Transport of small molecules,” and “steroid metabolic process” (Figures 4B,C).
FIGURE 4. Identification of DEGs between subtypes and construction of the prognostic model. (A) Volcano plot showing the differentially upregulated (red points) and downregulated genes (blue points). (B) Pathway and process enrichment analysis has been carried out for 344 DEGs that identified between subtypes. The graphical representation showed top 20 enrichments with p < 0.01. (C) Enrichment terms with a similarity >0.3 are connected by edges. (D) LASSO coefficient profiles of 80 prognostic DEGs. (E) Cross-validation for tuning parameter selection in the LASSO model. (F) Forest plot of the multivariate regression of 17-genes in prognostic model.
A univariable Cox proportional hazards regression model was performed to identify 80 prognostic DEGs with p < 0.001. LASSO method was used for variable selection in a Cox regression model to determine significant prognostic DEGs. One SE above the minimum criteria was chosen, resulting in a model with 17 prognostic genes (Figures 4D,E). Then, based on the expression of the 17 genes, we established a multivariate Cox proportional hazard model (Figure 4F, Supplementary Table S4).
Validation of the prognostic model
According to the formula, the risk score of each patient with KIRC was calculated. Patients were classified into the high- and low-risk score groups using the median as the cutoff value (Figure 5A). The distribution plot of the risk scores revealed that survival time decreased while death rates increased with an increase in risk scores in TCGA cohort (Figure 5B). Figure 5C displays the expression of 17 genes in the prognostic model between high and low-risk groups in TCGA (Figure 5C). Furthermore, patients in the high-risk group had a significantly poorer OS (Figure 5D, p = 1.0 E−15, log-rank test). According to the area under the curve (AUC) of the receiver operating characteristic (ROC) curve, the risk score was able to accurately predict mortality (Figure 5E, AUC = 0.743). Mutations in the tumor suppressor TP53 are associated with various human cancers; consequently, we validated the prognosis power of risk score among TP53 mutation status in TCGA cohort. Patients in the high-risk group had a worse prognosis than patients with TP53 mutation and wild type (Supplementary Figure S2, p = 0.092, p = 2.5E-14).
FIGURE 5. Validation of the prognostic model. (A) Ranked dot shows the distribution of risk score. (B) Scatter plots shows the distribution of patient survival status. (C) Heatmap shows the expression of 17-genes of prognostic model. (D) Assessment of the difference in OS between high risk and low risk samples in TCGA cohort by log-rank test. (E) ROC curves to predict the sensitivity and specificity of 3-years survival according to the risk score. (F) Distribution of DFS between high risk and low risk groups in GSE22541 cohort. (G,H) Kaplan-Meier curves show the independent relevance between overall survival time and risk score in GSE29609 and E-MTAB-1980 cohorts.
Next, we validated the prognosis power of the risk score in independent datasets. In GSE22541 cohort, the disease-free survival (DFS) time of high-risk patients was lower than that of low-risk patients (Figure 5F, p = 0.055, log-rank test). Additionally, survival analysis was carried out in two KIRC cohorts (GSE29606 and E-MTAB-1980), and high-risk scores indicated poor prognosis (Figures 5G,H, p = 0.055, p = 0.047, log-rank test), while the number of surviving patients in the low-risk group were more than those in the high-risk group (Supplementary Figure S3).
In addition, we explored the power of the risk score to predict the outcome of patients within clinical subgroups. The survival analysis revealed that high-risk score patients had a significantly poorer OS compared with that of the low-risk score patients in females, males, age (≥65), age (<65), T1-stage, T2-stage, T3-stage, N0-stage, M0-stage, M1-stage, Stage I, Stage III, Stage IV, Grade 2, Grade 3 or Grade 4 subgroups (Figure 6, p < 0.05, log-rank test).
FIGURE 6. Stratification analysis of the risk score in clinical subgroups. (A)–(P) Survival analysis for high risk score and low risk score patients in sex, age, T-stage, N-stage, M-stage, tumor stage, and tumor grade subgroups.
Correlation of risk scores and immunotherapy
We further investigated the association between risk scores and immune infiltration. The stromal, immune, and ESTIMATE scores of high-risk score samples were higher than those in low-risk score samples (Figure 7A). The Pearson correlation analysis was used to assess the correlation between risk scores and the abundance of immune cells. Infiltration of “T cells regulatory (Tregs),” “T cells CD4 memory activated,” “Plasma cells,” “Macrophages M0,” “Neutrophils,” “T cells CD8,” and “T cells follicular helper” were significantly positively correlated with-risk scores (Figures 7B–H, p < 0.05, Pearson correlation analysis). We also assessed the relationship between the expression of five immune checkpoints and the risk score. PD-1, CTLA4, and BTLA expressions were significantly higher in high-risk samples compared with low-risk samples in TCGA KIRC cohort, and PD-L1 expression was significantly lower in high-risk samples (Figure 7I, p < 0.05).
FIGURE 7. Correlation of risk score and immune cells infiltration. (A) Distribution of ESTIMATE score in high risk and low risk groups. (B)–(H) Positive correlation between risk score and immune cells. (H) Distribution of ESTIMATE score in high risk and low risk groups. (I) Expression of five immune checkpoints in high risk and low risk groups.
To further explore if the risk score can predict patients’ response to immunotherapy, we evaluated the risk score differences among patients in immunotherapy response subgroups. The risk score of patients with progressive disease (PD) or stable disease (SD) was significantly higher compared with patients with partial response (PR) (Figures 8A,B, p = 0.032, p = 0.042). Moreover, the risk score of non-responsive (PD and SD) patients were significantly higher than that of responsive (PR and complete response [CR]) patients (Figure 8C, p < 0.011). The numbers of patients in the high-risk group with PD and SD were more than those in the low-risk group; PR and CR patients were more abundant in the low-risk group (Figure 8D, p = 0.073, hypergeometric test).
FIGURE 8. Prognosis power of risk score in patients with immunotherapy. (A) Distribution of risk score between patients with PD and PR. (B) Distribution of risk score between patients with PR and SD. (C) Distribution of risk score between non-responsive (PD and SD) patients and responsive (PR and CR) patients in IMvigor210 cohort. (D) Distribution of non-responsive and responsive patients between high risk and low risk groups.
Discussion
In this study, we analyzed the role of CID in ccRCC progression, microenvironment alteration, and clinical prognosis. When analyzing the somatic mutation status of CID-related genes in ccRCC, most of them have a high mutation frequency. Meanwhile, the difference in CID-related gene expression is significant in cancer and normal tissues. These results further implied the potential of CID to be developed as a cancer therapy target and prognosis indicator (Ge et al., 2022; Oliveri, 2022).
Accordingly, we then identified ccRCC subtypes based on CID-related genes with a K-means clustering algorithm. Copper and CID plays the essential role in the mudutlation of cancer immune microenvironment. For example, recent research demonstrated that major copper influx transporter copper transporter 1 is correlated with PD-L1 expression across many cancer types (Voli et al., 2020). Meanwhile, copper chelators play the role in the inhibition of STAT3 and EGFR’s phosphorylation and promoted the degradation of PD-L1 (Voli et al., 2020). Further, copper in also correlated with the cancer’s immunogenic cell death in breast cancer (Kaur et al., 2020). Significant immune heterogeneity across the two ccRCC subtypes was observed. In Subtype B, ccRCC has the highest ESTIMATE score, infiltration level of CD8+ T cell and NK cell, and immune cytolytic activity. This result indicated that Subtype B may have higher immune activity. In general, the high immune activity of cancer implies a better prognosis (Chen and Mellman, 2017). However, the immune activity of ccRCC subtype B has a worse prognosis than subtype A. Xu et al. (2019) also indicated that high immune activity relates to poor prognosis. Nakano et al. (2001) also found that a high infiltration level of CD8+ T cells correlates with a poor RCC prognosis. Clonal variation of immune cells of the microenvironment may contribute to this unique characteristic of ccRCC (Borcherding et al., 2021).
Next, we acquired DEGs between Subtypes A and B and conducted the enrichment analysis to reveal the role of CID in ccRCC. According to the results, DEGs concentrated on the immune and metabolic-related processes. Cell toxicity mediated by copper was correlated with glucose metabolism activity (Li et al., 2022). Glucose metabolism alterations of microenvironment components, including cancer cells and immune cells, leading to the formation of different tumor subtypes (Li and Zhang, 2016; Terrén et al., 2019; Zhang et al., 2021). Our research demonstrated that CID may contribute to the immune microenvironment heterogeneity in ccRCC. Consequently, further analysis concentrating on the detailed interaction of CID-mediated metabolism alteration and TME may illustrate the regulatory function of CID in cancer.
Based on the obtained DEGs, we constructed a model to predict the prognosis of patients with ccRCC. Our result demonstrated that the risk score can greatly predict patient prognosis in training and validation datasets. Furthermore, the risk score is effective in different stages (i.e., T, N, and M) and grades. Accordingly, our risk score offers great clinical applicability.
Immunotherapy is widely applied in treating different types of solid cancer (Helmy et al., 2013; Robert et al., 2015; Larkin et al., 2019). Therefore, we further explored the potential of our risk score for predicting the immune features and therapy responsiveness of cancer. High-risk scores predict the lower expression level of PD-L1. Meanwhile, in the anti-PD-L1 cohort, high-risk scores are correlated with a low therapy responsiveness rate. Consequently, anti-PD-L1 therapy may be a suitable choice for low-risk score patients.
It is worth noting the limitations of the research. First, large-scale multi-omics immunotherapy data should be employed to more comprehensively evaluate the potential of the risk score for predicting immunotherapy responsiveness. Due to the lack of high-quality validated data, the enrolled immunotherapy samples are limited. Second, combining the transcriptome analysis of the clinical samples and follow-up data will further test the robustness of the risk score. Third, in vivo and in vitro cell-line and animal models may help explore the potential underlying mechanism of CID in cancer. These shortcomings will be overcome with the rapid progression of big data and our further in-depth research.
In summary, our research revealed the role of CID in ccRCC, identified ccRCC subtypes based on CID features and constructed a robust gene prognostic model to predict patient prognosis. Our research laid a foundation for CID-related analysis and presented a prognostic model which can be potentially applied in the clinical treatment of ccRCC.
Data availability statement
Publicly available datasets were analyzed in this study. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.
Author contributions
YX: manuscript preparation and research conception. SX: manuscript preparation, manuscript revision and data analysis. HJ: manuscript revision. ZQ: bioinformatic analysis and data visualization.
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/fgene.2022.994999/full#supplementary-material
Abbreviations
ccRCC, clear cell renal cell carcinoma; CID, copper ionophore–induced death; CNV, copy number variation; CR, complete response; CYT, cytolytic activity; DEGs, differentially expressed genes; DFS, disease-free survival; TCGA, The Cancer Genome Atlas; GEO, Gene Expression Omnibus; LASSO, Least Absolute Shrinkage and Selection Operator; RCC, renal cell carcinoma; TCA, tricarboxylic acid; KIRC, kidney clear cell carcinoma; mUCC, metastatic urothelial cancer; OS, overall survival; PCA, principle component analysis; SE, standard error; TME, tumor immune microenvironment; Tregs, T cells regulatory; PD, progressive disease; SD, stable disease; PR, partial response; UCSC, the University of California Santa Cruz.
References
Borcherding, N., Vishwakarma, A., Voigt, A. P., Bellizzi, A., Kaplan, J., Nepple, K., et al. (2021). Mapping the immune environment in clear cell renal carcinoma by single-cell genomics. Commun. Biol. 4, 122. doi:10.1038/s42003-020-01625-6
Braun, D. A., Street, K., Burke, K. P., Cookmeyer, D. L., Denize, T., Pedersen, C. B., et al. (2021). Progressive immune dysfunction with advancing disease stage in renal cell carcinoma. Cancer Cell. 39 (5), 632–648.e8. e8. doi:10.1016/j.ccell.2021.02.013
Chen, D. S., and Mellman, I. (2017). Elements of cancer immunity and the cancer-immune set point. Nature 541, 321–330. doi:10.1038/nature21349
Choueiri, T. K., Halabi, S., Sanford, B. L., Hahn, O., Michaelson, M. D., Walsh, M. K., et al. (2017). Cabozantinib versus sunitinib as initial targeted therapy for patients with metastatic renal cell carcinoma of poor or intermediate risk: The alliance A031203 CABOSUN trial. J. Clin. Oncol. 35 (6), 591–597. doi:10.1200/JCO.2016.70.7398
Engebretsen, S., and Bohlin, J. (2019). Statistical predictions with glmnet. Clin. Epigenetics 11, 123. doi:10.1186/s13148-019-0730-1
Ge, E. J., Bush, A. I., Casini, A., Cobine, P. A., Cross, J. R., DeNicola, G. M., et al. (2022). Connecting copper and cancer: From transition metal signalling to metalloplasia. Nat. Rev. Cancer 22, 102–113. doi:10.1038/s41568-021-00417-2
Helmy, K. Y., Patel, S. A., Nahas, G. R., and Rameshwar, P. (2013). Cancer immunotherapy: Accomplishments to date and future promise. Ther. Deliv. 4, 1307–1320. doi:10.4155/tde.13.88
Jonasch, E., Gao, J., and Rathmell, W. K. (2014). Renal cell carcinoma. BMJ 349, g4797. doi:10.1136/bmj.g4797
Kaur, P., Johnson, A., Northcote-Smith, J., Lu, C., and Suntharalingam, K. (2020). Immunogenic cell death of breast cancer stem cells induced by an endoplasmic reticulum-targeting copper(II) complex. Chembiochem 21, 3618–3624. doi:10.1002/cbic.202000553
Kennedy, L. B., and Salama, A. (2020). A review of cancer immunotherapy toxicity. Ca. Cancer J. Clin. 70 (2), 86–104. doi:10.3322/caac.21596
Larkin, J., Chiarion-Sileni, V., Gonzalez, R., Grob, J. J., Rutkowski, P., Lao, C. D., et al. (2019). Five-year survival with combined nivolumab and ipilimumab in advanced melanoma. N. Engl. J. Med. 381, 1535–1546. doi:10.1056/NEJMoa1910836
Li, S. R., Bu, L. L., and Cai, L. (2022). Cuproptosis: Lipoylated TCA cycle proteins-mediated novel cell death pathway. Signal Transduct. Target. Ther. 7, 158. doi:10.1038/s41392-022-01014-x
Li, Z., and Zhang, H. (2016). Reprogramming of glucose, fatty acid and amino acid metabolism for cancer progression. Cell. Mol. Life Sci. 73, 377–392. doi:10.1007/s00018-015-2070-4
Lopez, J., Ramchandani, D., and Vahdat, L. (2019). Copper depletion as a therapeutic strategy in cancer. Metall. Ions Life Sci. 19. doi:10.1515/9783110527872-018
Mariathasan, S., Turley, S. J., Nickles, D., Castiglioni, A., Yuen, K., Wang, Y., et al. (2018). TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature 554, 544–548. doi:10.1038/nature25501
Motzer, R. J., Barrios, C. H., Kim, T. M., Falcon, S., Cosgriff, T., Harker, W. G., et al. (2014). Phase II randomized trial comparing sequential first-line everolimus and second-line sunitinib versus first-line sunitinib and second-line everolimus in patients with metastatic renal cell carcinoma. J. Clin. Oncol. 32 (25), 2765–2772. doi:10.1200/JCO.2013.54.6911
Motzer, R. J., Escudier, B., McDermott, D. F., George, S., Hammers, H. J., Srinivas, S., et al. (2015). Nivolumab versus everolimus in advanced renal-cell carcinoma. N. Engl. J. Med. 373 (19), 1803–1813. doi:10.1056/NEJMoa1510665
Nakano, O., Sato, M., Naito, Y., Suzuki, K., Orikasa, S., Aizawa, M., et al. (2001). Proliferative activity of intratumoral CD8(+) T-lymphocytes as a prognostic factor in human renal cell carcinoma: Clinicopathologic demonstration of antitumor immunity. Cancer Res. 61, 5132–5136.
Newman, A. M., Liu, C. L., Green, M. R., Gentles, A. J., Feng, W., Xu, Y., et al. (2015). Robust enumeration of cell subsets from tissue expression profiles. Nat. Methods 12, 453–457. doi:10.1038/nmeth.3337
Oliveri, V. (2022). Selective targeting of cancer cells by copper ionophores: An overview. Front. Mol. Biosci. 9, 841814. doi:10.3389/fmolb.2022.841814
Patard, J. J., Pignot, G., Escudier, B., Eisen, T., Bex, A., Sternberg, C., et al. (2011). ICUD-EAU international consultation on kidney cancer 2010: Treatment of metastatic disease. Eur. Urol. 60 (4), 684–690. doi:10.1016/j.eururo.2011.06.017
Robert, C., Schachter, J., Long, G. V., Arance, A., Grob, J. J., Mortier, L., et al. (2015). Pembrolizumab versus ipilimumab in advanced melanoma. N. Engl. J. Med. 372, 2521–2532. doi:10.1056/NEJMoa1503093
Robinson, M. D., McCarthy, D. J., and Smyth, G. K. (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinforma. Oxf. Engl. 26, 139–140. doi:10.1093/bioinformatics/btp616
Rooney, M. S., Shukla, S. A., Wu, C. J., Getz, G., and Hacohen, N. (2015). Molecular and genetic properties of tumors associated with local immune cytolytic activity. Cell. 160, 48–61. doi:10.1016/j.cell.2014.12.033
Sato, Y., Yoshizato, T., Shiraishi, Y., Maekawa, S., Okuno, Y., Kamura, T., et al. (2013). Integrated molecular analysis of clear-cell renal cell carcinoma. Nat. Genet. 45, 860–867. doi:10.1038/ng.2699
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. Cancer J. Clin. 71 (3), 209–249. doi:10.3322/caac.21660
Terrén, I., Orrantia, A., Vitallé, J., Zenarruzabeitia, O., and Borrego, F. (2019). NK cell metabolism and tumor microenvironment. Front. Immunol. 10, 2278. doi:10.3389/fimmu.2019.02278
Tsang, T., Davis, C. I., and Brady, D. C. (2021). Copper biology. Curr. Biol. 31 (9), R421–R427. doi:10.1016/j.cub.2021.03.054
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
Wettersten, H. I., Aboud, O. A., Lara, P. N., and Weiss, R. H. (2017). Metabolic reprogramming in clear cell renal cell carcinoma. Nat. Rev. Nephrol. 13 (7), 410–419. doi:10.1038/nrneph.2017.59
Xu, W. H., Xu, Y., Wang, J., Wan, F. N., Wang, H. K., Cao, D. L., et al. (2019). Prognostic value and immune infiltration of novel signatures in clear cell renal cell carcinoma microenvironment. Aging 11, 6999–7020. doi:10.18632/aging.102233
Yin, J. M., Sun, L. B., Zheng, J. S., Wang, X. X., Chen, D. X., and Li, N. (2016). Copper chelation by trientine dihydrochloride inhibits liver RFA-induced inflammatory responses in vivo. Inflamm. Res. 65 (12), 1009–1020. doi:10.1007/s00011-016-0986-2
Zhang, D., Xu, X., and Ye, Q. (2021). Metabolism and immunity in breast cancer. Front. Med. 15, 178–207. doi:10.1007/s11684-020-0793-6
Keywords: ccRCC, machine-learning, bioinformatics, immune microenvironment, precise medicine, CID
Citation: Xia S, Jia H, Qian Z and Xiu Y (2022) Role of copper ionophore–induced death in immune microenvironment and clinical prognosis of ccRCC: An integrated analysis. Front. Genet. 13:994999. doi: 10.3389/fgene.2022.994999
Received: 15 July 2022; Accepted: 29 August 2022;
Published: 03 October 2022.
Edited by:
Lixin Cheng, Jinan University, ChinaReviewed by:
Xinyi Liu, University of Illinois at Chicago, United StatesLiyuan Fan, Shandong University, China
Copyright © 2022 Xia, Jia, Qian and Xiu. 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: Youcheng Xiu, xiuyoucheng0705@sina.com