Mitochondrial metabolic reprogramming-mediated immunogenic cell death reveals immune and prognostic features of clear cell renal cell carcinoma

Background: Mitochondrial metabolic reprogramming (MMR)-mediated immunogenic cell death (ICD) is closely related to the tumor microenvironment (TME). Our purpose was to reveal the TME characteristics of clear cell renal cell carcinoma (ccRCC) by using them.

Methods: Target genes were obtained by intersecting ccRCC differentially expressed genes (DEGs, tumor VS normal) with MMR and ICD-related genes. For the risk model, univariate COX regression and K-M survival analysis were used to identify genes most associated with overall survival (OS). Differences in the TME, function, tumor mutational load (TMB), and microsatellite instability (MSI) between high and low-risk groups were subsequently compared. Using risk scores and clinical variables, a nomogram was constructed. Predictive performance was evaluated by calibration plots and receiver operating characteristics (ROC).

Results: We screened 140 DEGs, including 12 prognostic genes for the construction of risk models. We found that the immune score, immune cell infiltration abundance, and TMB and MSI scores were higher in the high-risk group. Thus, high-risk populations would benefit more from immunotherapy. We also identified the three genes (CENPA, TIMP1, and MYCN) as potential therapeutic targets, of which MYCN is a novel biomarker. Additionally, the nomogram performed well in both TCGA (1-year AUC=0.862) and E-MTAB-1980 cohorts (1-year AUC=0.909).

Conclusions: Our model and nomogram allow accurate prediction of patients’ prognoses and immunotherapy responses.

Background

Renal cell carcinoma (RCC) is one of the most common urological malignancies, especially clear cell renal cell carcinoma (ccRCC), with approximately 431 288 new cases of RCC in 2020, 70% of which were ccRCC (1, 2). For non-metastatic RCC, surgical treatment can achieve curative results, but the high recurrence rate (40%) after surgery is a major factor affecting the long-term survival of patients (3). In addition, although systemic therapy is the standard of care for metastatic RCC, the long-term survival of patients with metastatic RCC has remained suboptimal for many years. With the rise of immunotherapy, combined treatment modalities targeting vascular endothelial growth factor and immune checkpoints offer new hope for improving the prognosis of metastatic RCC and high-risk limited RCC (4, 5). Nevertheless, the response rates and durability of treatments observed in clinical practice have not been satisfactory, mainly due to the lack of guiding and effective biomarkers (6, 7).

Metabolic reprogramming is the phenomenon that occurs when tumor cells undergo a rapid adaptive response to hypoxic and low oxygen conditions to proliferate and migrate in the tumor microenvironment (TME), and it is one of the important hallmarks of cancer (8). Mitochondria are known to play a key role in biosynthesis, energy metabolism, and signaling, and are dynamic organelles that coordinate functions such as metabolism, cell proliferation, and cell survival (9). Tumor cells are no exception, and their malignant transformation and development are dependent on mitochondrial energy metabolism. Interestingly, there is also metabolic heterogeneity in tumor tissues. The adaptive metabolic response of tumor cells not only regulates each other with the TME they are in but also influences the way genes are expressed, leading to great differences in the malignancy and therapeutic response of different individual tumors (10, 11).

For example, in ccRCC, lipid metabolism, and lactate metabolic processes, which are both closely related to the occurrence and progression of ccRCC, as well as the identification of potentially valuable therapeutic targets (FBP1, HADH, and TYMP) based on genes related to lactate metabolism by Sun et al. (11–14). With a deeper understanding of mitochondria in cancer development, cancer therapeutic modalities targeting mitochondrial metabolism are gradually being translated into clinical practice (15, 16). Surprisingly, recent studies have shown that mitochondrial metabolic reprogramming (MMR) can enhance immunogenic cell death (ICD) effects and thus optimize cancer immunotherapy outcomes (17, 18). ICD can activate antitumor immune effects by remodeling the TME. ICD stimulates dendritic cell (DC) maturation under endogenous conditions through the release of tumor-associated antigens and danger-associated molecular patterns, including adenosine triphosphate, high-motility group box 1, and calreticulin. Once DC mature, they deliver antigens to cytotoxic T lymphocytes, which then activate antitumor immunity (19, 20).

Therefore, our study aims to use MMR and ICD-related genes to identify and reveal clear cell renal cell TME and prognostic features to screen individuals suitable for immunotherapy to help clinicians in scientific decision-making.

Materials and methods

Datasets

We obtained the training dataset from the Cancer Genome Atlas (TCGA) database(https://portal.gdc.cancer.gov/). Including TCGA-KIRC TPM and raw counts gene expression data (tumor sample: 541, normal sample: 72) and clinical data of 533 patients. TCGA-KIRC data was processed by log2 (TPM+1) for subsequent analysis. Moreover, we excluded genes with too low expression (more than 25% of samples with counts expression less than 10). From the ArrayExpress database (https://www.ebi.ac.uk/biostudies/arrayexpress), E-MTAB-1980, represented the validation cohort, including 101 patients with ccRCC. The identification of MMR and ICD-related genes was based on the GeneCards database and previously published literature. To begin with, we searched the GeneCards database for “mitochondrial metabolic reprogramming” and “immunogenic cell death” and then screened for protein-coding genes with a correlation score greater than 0.5. We obtained 34 ICD-related genes from the results of Abhishek’s study (21). Eventually, we obtained 3103 genes related to MMR and 1732 genes related to ICD.

Weighted gene co-expression network analysis (WGCNA)

The close association between MMR and ICD-related genes and ccRCC was confirmed using the R package “WGCNA”. Co-expression networks with scale-free topologies were selected using appropriate soft-threshold power. Further, to investigate gene connectivity in this network, a topological overlap matrix (TOM) was converted from the adjacency matrix. A hierarchical clustering tree (dendrogram) of genes was generated based on TOM, gene modules are represented by different branches on the clustering tree, and by different colors.

Clinical and immune features included stage, immune score, and CD8 T-cell infiltration levels. Each sample was assessed using the “estimate” package to determine its immune score, stromal score, estimate score, and tumor purity. The online analysis tool TIMER2.0 (http://timer.cistrome.org/) can download the immune cell infiltration score file for each TCGA-KIRC sample, and we used 22 immune cell infiltration scores based on the CIBERSORT algorithm.

Identification and analysis of target genes (TGs)

Using the “DESeq2” package to screen the differentially expressed genes (DEGs) (tumor VS normal). A differential gene’s screening criteria require the p-value< 0.05, and the |log2fold change|> 1.5. The DEGs, MMR, and ICD-related genes were then intersected to obtain TGs for subsequent analysis.

Somatic mutation data of ccRCC samples were downloaded from the TCGA database, and the “Maftools” package was used to calculate the tumor mutational load (TMB) and visualize the mutation landscape of TGs. In addition, “clusterProfiler” was used for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses to understand which biological functions are associated with the TGs.

Developing a risk model based on prognosis-related genes

Univariate COX regression and K-M survival analysis (screening criteria: p<0.01; GEPIA database, http://gepia.cancer-pku.cn/index.html) were performed using TCGA-KIRC expression and survival data to identify TGs most associated with overall survival (OS), which was then used to perform a least absolute shrinkage and selector operation analysis (LASSO) with the “glmnet” package. Patients’ risk scores were calculated by multiplying their expression value by their coefficient by summing the scores of each OS-related gene. The median risk score was used to distinguish low-risk patients from high-risk patients. In addition, principal component analysis (PCA) was used to assess how well model genes discriminated between samples. Further validating the model’s predictive capability, we analyzed TCGA-KIRC and E-MTAB-1980 cohorts using Kaplan-Meier (K-M) survival analyses. As well, we compared the OS and risk scores of patients with varying clinical characteristics.

The immune landscape across different risk groups

For each sample, the TME score and infiltration level of 22 immune cells were calculated, and then the differences between the two risk groups were observed. Furthermore, immunoinhibitors and immunostimulators are closely related to the TME, we evaluated their expression levels in samples from different risk groups.

Functional and gene set enrichment analysis (GSEA)

DEGs between the two risk groups were identified with the “limma” package. After using the “ClusterProfiler” package, GO and KEGG analyses were performed, and GSEA was to identify differential signaling pathways, biological effects, and gene sets that are enriched in different populations. The “GSVA” package was used to compute normalized enrichment scores for pathway and functional annotations using the gene set variation analysis (GSVA) approach. The “c2.cp.kegg.v7.4.symbols.gmt” annotation file was obtained from the MSigD. To further assess the immune-related functions score for each sample, we used the “GSEABase” and “GSVA” packages and then combined them with the immune gene set annotation file “immune. gmt” for enrichment analysis.

The prediction of immunotherapy response and drug sensitivity

Palmeri’s study indicated that TMB and microsatellite instability (MSI) could predict the effect of immunotherapy (22). Therefore, we assessed the TMB levels and MSI status of the different risk groups. The TMB was calculated using the R package “maftools” and the microsatellite instability score was downloaded using the package “cBioPortalData” in the file “kirc_tcga _pan_can_atlas_2018”. We choose “MSI_SENSOR_SCORE”, if the score is greater than 0.3, it is defined as MSI, and vice versa as MSS. Furthermore, we used the “easier” package (23) to obtain the immune response score file and identify the relationship between risk score and immunotherapy response. Each sample’s sensitivity scores to varying drugs were calculated using the package “oncoPredict” (24), with higher scores indicating greater. The reference file “GDSC2_Expr (RMA Normalized and Log Transformed).rds” was from Cancer Drug Sensitivity Genomics (https://www.cancerrxgene.org/).

Comprehensive analysis of model genes

Our subsequent analysis focused on the three genes with the highest absolute risk coefficients among the twelve model genes. The online analysis tools Gene Set Cancer Analysis (GSCA, http://bioinfo.life.hust.edu.cn/GSCA/#/) and TISIDB (http://cis.hku.hk/TISIDB/) were used to analyze the expression differences, methylation profiles, prognostic significance of these three genes in pan-cancer and immune relevance. In addition, we analyzed their distribution in different cell types using single-cell sequencing data from the online tool Tumor Immune Single-cell Hub 2 (TISCH2, http://tisch.comp-genomics.org/).

Developing and validating nomograms

To screen for independent predictors of OS, univariate and multivariate COX regression analyses were performed for age, gender, stage, and risk score with a screening criterion of p<0.05. Nomograms were then constructed using the “ survival” and “ rms” packages. Furthermore, the model’s predictive performance was validated with calibration plots and receiver operating characteristics (ROC) analyses. The same analysis was performed for the E-MTAB-1980 cohort.

Validation of gene expression and protein expressions

The real-time fluorescence quantitative polymerase chain reaction (RT-qPCR) was used to verify the expression levels of key genes. Human RCC cells are derived from the Procell Life Science&Technology Co., Ltd (Wuhan, China). After RNA was extracted using TRIzol reagent (Invitrogen, Thermo Fisher Scientific, Inc.), reverse transcription was performed using the Takara PrimeScript RT kit (Takara Bio, Inc., Otsu, Japan). Lastly, RT-qPCR was performed using the SYBR premix Ex Taq kit (Takara Bio, Inc., Otsu, Japan) on the Roche LightCycler96 real-time fluorescent quantitative PCR system. Based on the 2^−ΔΔCt method, the FC in mRNA was calculated.

Statistical analysis

Our statistical analyses were performed using R (version 4.2.2) or GraphPad Prism (version 9.0). All statistics were two-sided tests, and p<0.05 was considered statistically significant if not otherwise stated. An analysis of Spearman correlation was performed to determine the correlation between two continuous variables. Comparing two independent samples with Wilcoxon rank sum tests and multiple samples with Kruskal-Wallis tests. Comparing two categorical variables was accomplished with the chi-square test. The significance of differences in K-M survival analysis was assessed using the log-rank test.

Result

Screening of TGs for subsequent analysis

Using the WGCNA algorithm, we evaluated the correlations of MMR and ICD-related genes with ccRCC prognosis and TME. As shown in Figure 1A, we chose 11 as the optimal soft threshold power to construct the co-expression network. Hierarchical clustering analysis was then performed based on weighted correlations, and these genes were classified into different modules (Figure 1B). Figure 1C shows that most of the modules are remarkedly related to the stage, immune infiltration score, and CD8 T cell infiltration abundance, with the red module being the most relevant. Hence, the development of ccRCC appears to be strongly linked to MMR and ICD-related genes. Therefore, we first screened the DEGs between tumor and normal tissues (Figure 1D) and then took the intersection of it and MMR and ICD-related genes to identify 140 target genes for subsequent analysis (Figure 1E). Comparatively to normal tissues, volcano map analysis showed that 108 genes had up-regulation and 32 genes had down-regulation (Figure 1F). In addition, the mutation landscape of TGs showed the highest frequency of mutations in BRCA2 and ERBB4 (Figure 1G). Their close association with immune-related functions and pathways was revealed by enrichment analysis, including regulation of T cell activation, cytokine receptor binding, cytokine activity, and Cytokine-cytokine receptor interaction (Figure 1H). This further corroborates that TGs are closely associated with the TME of ccRCC.

FIGURE 1

Figure 1 Screening process of target genes. (A) Determine the appropriate soft threshold power equal to 11. (B) Generate a hierarchical clustering tree of genes based on TOM, with different colors representing different modules. (C) The modules are closely related to the clinical features and TME of ccRCC, with the red module being the most important. (D) Differential ranking map of DEGs between tumor and normal tissues. (E) Venn diagram of DEGs, MMR, and ICD-related genes taking the intersection. (F) Volcano plot of 140 target genes. (G) Mutation waterfall plot of target genes. H Enrichment bubble map of biological functions and pathways of target genes.

Risk model based on 12 genes

According to the screening criteria (p<0.01), 23 prognostic genes most associated with OS were identified (Figure 2A). On the basis of LASSO regression analysis, a 12-gene risk model was subsequently developed (Figure 2B) [risk score=ABCB1*(-0.114) + CCND1*(-0.002) +CENPA*0.328+CXCL5*0.053+DPEP1*(-0.058) +EPCAM*(-0.028) +FGF1*(-0.021) +KITLG*(-0.094) +MYCN*(-0.232) +OSM*0.052+PLIN2*(-0.047) +TIMP1*0.149]. The results of PCA also showed strong risk differentiation of model genes (Figure 2C). Furthermore, it can be observed from (Figures 2D, E) that the trends of risk scores and survival status of the training set TCGA and validation set E-MTAB-1980 are highly consistent. The heat map of the distribution of model genes with risk scores also remained almost consistent. The results of the K-M survival analysis also showed significant differences in OS between the two risk groups (Figures 2F, G), which indicates the strong predictive performance of our model.

FIGURE 2

Figure 2 Development of risk model. (A) Prognostic genes were screened based on univariate COX regression and K-M survival analysis. (B) LASSO regression was used to construct a 12-gene-based risk model. (C) PCA analysis showed that the model genes had a good discriminatory ability. (D, E) Risk scores of TCGA and E-MTAB-1980 cohorts were correlated with survival status and model genes. (F, G) K-M survival analysis for high- and low-risk groups.

The risk score and clinical variables.

Depending on the patient’s age, gender, and stage, we classify them into different groups. It can be seen from Figure 3A that patients in the high-risk group had a worse prognosis regardless of age (<=65 years, >65 years). Risk scores between the two age groups did not differ significantly, however (Figure 3B). Patients with high-risk outcomes had shorter survival rates in both males and females, and males had a higher percentage of high-risk patients (Figures 3C, D). As well, both early (stage I/II) and late-stage (stage III/IV) patients in the high-risk group had a worse outcome, with a higher proportion of late-stage patients at high risk (Figures 3E, F).

FIGURE 3

Figure 3 Relationship of risk scores to clinical variables. (A, B) K-M survival analysis and comparison of risk scores for high and low-risk groups in patients <=65 years and >65 years. (C, D) K-M survival analysis of high and low-risk groups of male and female patients with comparison of risk scores by gender. (E, F) K-M survival analysis of high- and low-risk groups of early and late-stage patients, with a comparison of risk scores by stage.

Immune landscapes in different risk groups

In the high-risk group, we observed higher estimated scores, immune scores, and tumor purity scores (Figure 4A). Furthermore, higher levels of CD8+ T cells, TFH T cells, and memory-activated CD4+ T cells were found in the high-risk group (Figure 4B). There were also differences in immunoinhibitor and immunostimulator expression levels between the two risk groups. According to Figure 4C, high-risk individuals expressed higher levels of CTLA4, LAG3, PDCD1, TGFB1, and TIGIT. High-risk patients also expressed higher levels of immunostimulators, such as CD27, CD276, CD28, CD70, IL2RA, IL6, and TNFRSF18 (Figure 4D). In addition, we observed the same expression trend in the validation cohort.

FIGURE 4

Figure 4 Immune landscapes in different risk groups. (A) Estimated scores, immune scores, stromal scores, and tumor purity for high- and low-risk groups. (B) 22 levels of immune cell infiltration. (C, D) Expression levels of immunoinhibitor and immunostimulator agents in TCGA and E-MTAB-1980 cohorts. “*”<0.5, "**" <0.01, and "***"<0.001.

Functional and pathway differences between high and low-risk groups

The results of GO and KEGG analysis showed that the differential genes between high and low-risk groups were mainly enriched in immune-related biological functions and pathways, such as immunoglobulin complex and antigen binding (Figure 5A). It appears that pathway enrichment distribution in the two groups differs significantly according to the GSVA results, with the high-risk group mainly enriched in immune-related pathways and the low-risk group enriched in metabolism-related pathways (Figure 5B). As a result of GSEA, the top 5 most significantly enriched pathways are shown. The immune-related pathways that were significantly enriched in the high-risk group included cytokine-cytokine receptor interaction and primary immunodeficiency (Figure 5C). Furthermore, we calculated the immune function score for each sample. High-risk individuals scored higher on all immune-related functions except type II IFN response (Figure 5D). The same trend was maintained in the analysis of the validation cohort. Hence, combining the results of these analyses, it is reasonable to assume that risk scores are closely related to the TME of ccRCC and could potentially guide immunotherapy.

FIGURE 5

Figure 5 Functional and pathway enrichment analysis. (A) GO and KEGG enrichment results of DEGs. (B) GSVA in high and low-risk groups. (C) GSEA of high and low-risk groups. (D) Comparison of differences in immune-related functions between high and low-risk groups in TCGA and E-MTAB-1980 cohorts. “*”<0.5, "**" <0.01, and "***"<0.001.

Immunotherapy response and drug sensitivity

Further exploring immunotherapy response, we compared the difference between the two risk groups according to TMB, MSI, and easier scores, with higher values predicting a greater likelihood of benefiting from immunotherapy. We observed that patients in the high-risk group had higher TMB and easier scores, and that risk scores changed positively with TMB (Figures 6A–C). A higher risk score was also observed in the MSI group (Figure 6B). Consequently, immunotherapy was more likely to succeed in high-risk patients, resulting in an improved outcome. Subsequently, we calculated the sensitivity of each sample to different drugs and then compared the differences in sensitivity of the top 5 drugs in two risk groups by ranking the mean of the sensitivities. Among those at high risk, Carmustine, Nelarabine, MIRA-1, and EPZ5676 appeared to be more sensitive. Perhaps, the combination of these drugs with immune checkpoint inhibitors could lead to surprising therapeutic effects, which provides new ideas for future prospective study designs.

FIGURE 6

Figure 6 Immunotherapy and drug prediction. (A) TMB levels in high and low-risk groups, correlation of risk scores with TMB. (B) Levels of risk scores in MSI and MSS groups, correlation of risk scores with MSI scores. (C) Levels of easier scores in high and low-risk groups, correlation of risk scores with easier scores. (D) Differences in sensitivity to different drugs in high and low-risk groups. “*”<0.5, "**" <0.01, and "***"<0.001

Key model genes affecting risk scores

We in-depth explored the three genes with the highest risk coefficient, including CENPA, TIMP1, and MYCN. Firstly, Figure 7A demonstrates their expression levels in pan-cancer. CENPA and TIMP1 are highly expressed in the tumor tissue of most cancers, while MYCN is usually lowly expressed. We also found that they were closely associated with OS in several cancer types. As with CENPA, TIMP1 was associated with poor prognosis in most cancers. MYCN was mainly associated with good prognosis, especially in ccRCC, where higher MYCN expression was associated with longer OS in patients (Figure 7B). We also evaluated the methylation status of the three genes in tumors and normal tissues. The methylation levels of CENPA and TIMP1 were downregulated in tumor tissues, while MYCN, in contrast to them, had elevated methylation levels (Figure 7C). Moreover, their expression levels were negatively correlated with the methylation levels (Figure 7D), which suggests that methylation could be responsible for their differential expression levels. Using single-cell sequencing data from GSE159115, we further investigated CENPA, TIMP1, and MYCN expression levels in different cell types. As shown in Figure 7E, all cells can be broadly classified into three types, including immune cells, malignant cells, and stromal cells. Compared with normal tissues, CENPA expression in tumor tissues was elevated in all three cell types; TIMP1 was elevated in stromal cells and decreased in malignant cells; MYCN expression levels were not significantly different in different cell types.

FIGURE 7

Figure 7 Comprehensive analysis of key model genes. (A) Pan-cancer analysis of gene expression levels. (B) Pan-cancer survival analysis. (C) Methylation differences between tumor and normal tissues. (D) Relationship between methylation levels and gene expression. (E) Single-cell analysis of key model genes. “****” <0.0001.

Subsequently, we analyzed the relationship of CENPA, TIMP1, and MYCN with the clinical features and TME of ccRCC. To begin with, the expression levels of CENPA and TIMP1 increased with the tumor stage, while MYCN did the opposite (Figure 8A). A similar trend was observed in tumor grade, where the expression levels of CENPA and TIMP1 were positively correlated with grade and that of MYCN was negatively correlated (Figure 8B). Based on these genes, ccRCC can be classified into different immune subtypes (Figure 8C), suggesting a close relationship between them and the tumor’s TME. This is further confirmed by the analysis presented in Figure 8D, where CENPA and TIMP1 were positively correlated with the infiltration level of most immune cells, and MYCN showed the opposite trend (Figure 8D). Therefore, we believe that CENPA, TIMP1, and MYCN are all therapeutic targets with significant potential in ccRCC.

FIGURE 8

Figure 8 Clinical features and immunological correlates. (A) Relationship between gene expression levels and stage. (B) Relationship between gene expression levels and grade. (C) Key model genes classify ccRCC into 6 immune subtypes. (D) The relationship between gene expression levels and the abundance of immune cell infiltration.

Nomogram

An independent predictor of OS in patients with ccRCC was age, stage, and risk score, as identified by univariate and multivariate COX regression analyses (Figures 9A, B). A nomogram was developed to predict survival probability in patients with ccRCC at 1, 3, and 5 years (Figure 9C). As seen in the calibration plots (Figure 9D), the predicted probabilities were almost in a straight line with the actual possibilities for both the TCGA and E-MTAB-1980 cohorts, indicating that the prediction results of the Nomogram were very accurate. Additionally, the results of ROC analysis also show that our nomogram has a high predictive value. The one-year Area Under Curve (AUC) of the nomogram in the TCGA cohort equals 0.862 (one-year AUC for risk score=0.732). The 1-year AUC of the nomogram for the E-MTAB-1980 cohort is 0.909 (0.890 for the risk score). In other words, our nomograms provide accurate prognostic information to patients and clinicians.

FIGURE 9

Figure 9 Construction of Nomogram. (A, B) Univariate and multivariate COX regression analysis. (C) Nomograms predict 1-, 3-, and 5-year survival probabilities. (D) Calibration plots for the TCGA and E-MTAB-1980 cohorts. (E, F) ROC curves of nomograms and risk scores in the TCGA and E-MTAB-1980 cohorts.

Expression validation of key model genes

Wang (25) and Shou (26) have confirmed through detailed in vitro experiments that CENPA and TIMP1 are highly expressed in ccRCC tumor tissues and promote tumor progression through different mechanisms. As a result, we assessed MYCN expression levels in tumors and normal tissues using RT-qPCR. Primer sequence of MYCN: F: CAGGCAAGACAGCAGCA; R: ATGTGCAAAGTGGCAGTGA. The MYCN expression in ACHN and OSRC cell lines was significantly lower than that in HK-2 normal cells (Figure 10B), which was in perfect agreement with the results of TCGA-KIRC analysis (Figure 10A).

FIGURE 10

Figure 10 Expression validation of MYCN. (A) Expression levels of MYCN in TCGA unpaired and paired samples. (B) Expression levels of MYCN in HK-2, ACHN, and OSRC cell lines.

Discussion

As the most common pathological type of RCC, ccRCC is highly immunogenic and heterogeneous. In other words, the TME may be dramatically different for different ccRCC, and this may account for the large differences in response rates to immune checkpoint inhibitors in different patients (27). Furthermore, it is due to such a highly heterogeneous TME that metabolites associated with tumor biological behavior are also increased (28). In numerous studies, metabolic reprogramming has been shown to occur during the development of ccRCC, which controls the tumor’s energetic and biosynthetic metabolisms (29). p M Herst’s study also pointed out that MMR is an important factor in determining tumor fate (30). Furthermore, Bianca’s study showed that MMR can induce the onset of ICD and thus exert tumor-killing effects (31). Therefore, an in-depth investigation of the significance of MMR and ICD in the TME of ccRCC is highly promising. As far as we know, our study is the first to combine these factors to reveal ccRCC prognostic factors and TME characteristics.

Our study integrated MMR and ICD-related genes and identified 140 genes with a significantly different expression between tumor and normal tissues. From these, 12 genes most associated with OS were then screened for use in constructing risk models. PCA and K-M survival analysis showed that risk grouping could distinguish different ccRCC patients well. Subsequently, we analyzed the differences in the TME, TMB, MSI, and functional clustering between the two risk groups to help identify ccRCC patients who have a greater chance of benefitting from immunotherapy. In addition, our comprehensive analysis of key model genes showed that CENPA, TIMP1, and MYCN are all therapeutic targets with significant potential. In the end, we constructed a nomogram with strong predictive performance to predict the survival probability of patients, which makes our study more clinically relevant.

According to our findings, in the high-risk group, CD8 T cells and TFH were more prevalent, as well as higher immune scores. this may seem contradictory, as generally, a high level of immune cell infiltration would inhibit tumor progression and thus improve patient prognosis. Perhaps the heterogeneity between different tumors causes ccRCC to differ from other immunotherapy-responsive solid tumors (32). We also found higher expression of the immune checkpoints PD-1 and CTLA4 in the high-risk group. The study by Pramod et al. points to the abundance of immune cell infiltration in the TME and the expression level of immune checkpoint genes as potential biomarkers for predicting immunotherapeutic response (33), and that abundant CD8 T-cell infiltration and immune checkpoint overexpression are associated with an excellent immunotherapeutic response. High-risk individuals also exhibited a greater abundance of immune-related pathways and functions, and the higher levels of TMB, MSI, and easier scores further validated our findings (22, 23). A risk model associated with the immune landscape of ccRCC was constructed by Zhuo et al. based on lactate metabolism-related genes, and they obtained similar results, inhibitors of immune checkpoints are more likely to work for patients with high-risk conditions (14). However, Wu’s findings were opposed to ours (34). Although they found that the high-risk group also had higher immune checkpoint expression and CD8 T-cell infiltration, they used the Tumor Immune Dysfunction and Exclusion (TIDE) scores to determine the response to immunotherapy in two risk groups. Perhaps, this is the reason for the very different results.

Subsequently, we analyzed the prognostic features and immune correlations of three key genes determining risk scores in ccRCC from multiple perspectives. Pan-cancer analysis revealed that Centromere Protein A (CENPA) and Metallopeptidase Inhibitor 1 (TIMP1) were highly expressed in most cancers and significantly associated with a poorer prognosis, also in ccRCC. CENPA overexpression enhanced ccRCC proliferation and metastasis by activating the Wnt/Kip1 pathway, according to Wang et al. (25). The protein encoded by TIMP1 is not only a natural inhibitor of matrix metalloproteinases but is also closely associated with cell proliferation and apoptosis in several cell types (35). Furthermore, in vitro experiments and bioinformatics analysis conducted by Shou’s study revealed that TIMP1 promotes RCC progression via the epithelial-mesenchymal transition signaling pathway (26). However, they do not elaborate on what roles CENPA and TIMP1 play in the TME of ccRCC. A significant correlation was found between CENPA and TIMP1 expression in tumor tissues and immune cell infiltration, and as expression increased, most immune cell infiltration levels increased. Moreover, ccRCC could be classified into different immune subtypes based on the expression of CENPA and TIMP1. This all suggests that CENPA and TIMP1 are highly promising biomarkers and therapeutic targets. Usually, MYCN is associated with poor prognosis because it is part of the MYC family. In particular, it plays a key role in pediatric neuroblastoma, where it can promote tumor cell growth and appreciation based on metabolic reprogramming (36). Nevertheless, our pan-cancer analysis showed that MYCN appears to be a protective factor for ccRCC and is much less expressed in renal tumor tissues. Moreover, MYCN was closely related to the immune subtype of ccRCC and the infiltration abundance of immune cells. However, no studies have been conducted to analyze in depth the potential mechanisms by which MYCN affects the prognosis of ccRCC, so we do not yet know why MYCN plays a different role in ccRCC than in other cancers. In other words, our study identified a novel and promising biomarker for ccRCC, but more reliable in vitro and in vivo experiments are needed to validate our findings.

In conclusion, our study constructs a 12-gene risk model based on MMR and ICD-related genes, which can accurately distinguish patients at different risks and guide patients for immunotherapy and targeted therapies. Furthermore, a nomogram based on risk scores and clinical variables can also accurately predict a patient’s OS, progress-free survival (PFS), and disease-specific survival (DSS) for 1, 3, and 5 years (Supplementary documents Figure1), and have high clinical application value. However, there are some shortcomings in our study. At first, our study was mainly a retrospective analysis based on public databases, and future high-quality prospective studies are needed to enhance the credibility of the conclusions. Secondly, we did not reveal the specific mechanism of action of key model genes by experimental means.

Conclusion

In this study, we revealed the TME and prognostic features of ccRCC from the perspective of mitochondrial metabolism and programmed cell death and established a new risk model and nomogram to distinguish patients with different molecular and clinical characteristics, which can be used as a reliable tool to predict the prognosis and immunotherapeutic response of ccRCC patients. Furthermore, the identification of the key molecule MYCN will facilitate the development of effective targeted drugs for ccRCC from the perspective of targeting mitochondrial metabolism. To conclude, our study provides both new insights into the regulatory mechanisms of MMR in ccRCC and assists clinicians in making rational decisions.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

Ethics statement

It was approved by Nanchang University’s First Affiliated Hospital’s Research Ethics Committee (Nanchang, China).

Author contributions

LY: Data curation, formal analysis, writing - original draft, writing - review & editing. SL, XL, and WD: Data curation, formal analysis. BF, JX, and WL: Conceptualized research, writing - review & editing. They all agreed to publish. All authors contributed to the article and approved the submitted version.

Funding

This study was supported by the Jiangxi Provincial “Double Thousand Plan” Fund Project (Grant No. jxsq2019201027), the Key Project of Natural Science Foundation of Jiangxi Province (20212ACB206013), and the Youth Project of Natural Science Foundation of Jiangxi Province (20212BAB216037).

Acknowledgments

The authors are grateful to the TCGA and ArrayExpress databases for providing high-quality data.

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/fonc.2023.1146657/full#supplementary-material

Glossary

References

1. Zhu X, Al-Danakh A, Zhang L, Sun X, Jian Y, Wu H, et al. Glycosylation in renal cell carcinoma: mechanisms and clinical implications. Cells (2022) 11(16). doi: 10.3390/cells11162598

CrossRef Full Text | Google Scholar

2. Bukavina L, Bensalah K, Bray F, Carlo M, Challacombe B, Karam JA, et al. Epidemiology of renal cell carcinoma: 2022 update. Eur Urol (2022) 82(5):529–42. doi: 10.1016/j.eururo.2022.08.019

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Larroquette M, Peyraud F, Domblides C, Lefort F, Bernhard J-C, Ravaud A, et al. Adjuvant therapy in renal cell carcinoma: current knowledges and future perspectives. Cancer Treat Rev (2021) 97:102207. doi: 10.1016/j.ctrv.2021.102207

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Ingels A, Campi R, Capitanio U, Amparore D, Bertolo R, Carbonara U, et al. Complementary roles of surgery and systemic treatment in clear cell renal cell carcinoma. Nat Rev Urol (2022) 19(7):391–418. doi: 10.1038/s41585-022-00592-3

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Navani V, Heng DYC. Treatment selection in first-line metastatic renal cell carcinoma-the contemporary treatment paradigm in the age of combination therapy: a review. JAMA Oncol (2022) 8(2):292–9. doi: 10.1001/jamaoncol.2021.4337

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Rosellini M, Marchetti A, Mollica V, Rizzo A, Santoni M, Massari F. Prognostic and predictive biomarkers for immunotherapy in advanced renal cell carcinoma. Nat Rev Urol (2022). doi: 10.1038/s41585-022-00676-0

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Deleuze A, Saout J, Dugay F, Peyronnet B, Mathieu R, Verhoest G, et al. Immunotherapy in renal cell carcinoma: the future is now. Int J Mol Sci (2020) 21(7). doi: 10.3390/ijms21072532

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Yoshida GJ. Metabolic reprogramming: the emerging concept and associated therapeutic strategies. J Exp Clin Cancer Res (2015) 34:111. doi: 10.1186/s13046-015-0221-y

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Ng MYW, Wai T, Simonsen A. Quality control of the mitochondrion. Dev Cell (2021) 56(7):881–905. doi: 10.1016/j.devcel.2021.02.009

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Avolio R, Matassa DS, Criscuolo D, Landriscina M, Esposito F. Modulation of mitochondrial metabolic reprogramming and oxidative stress to overcome chemoresistance in cancer. Biomolecules (2020) 10(1). doi: 10.3390/biom10010135

PubMed Abstract | CrossRef Full Text | Google Scholar

11. di Meo NA, Lasorsa F, Rutigliano M, Loizzo D, Ferro M, Stella A, et al. Renal cell carcinoma as a metabolic disease: an update on main pathways, potential biomarkers, and therapeutic targets. Int J Mol Sci (2022) 23(22). doi: 10.3390/ijms232214360

CrossRef Full Text | Google Scholar

12. Wu Y, Zhang S, Gong X, Tam S, Xiao D, Liu S, et al. The epigenetic regulators and metabolic changes in ferroptosis-associated cancer progression. Mol Cancer (2020) 19(1):39. doi: 10.1186/s12943-020-01157-x

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Heravi G, Yazdanpanah O, Podgorski I, Matherly LH, Liu W. Lipid metabolism reprogramming in renal cell carcinoma. Cancer Metastasis Rev (2022) 41(1):17–31. doi: 10.1007/s10555-021-09996-w

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Sun Z, Tao W, Guo X, Jing C, Zhang M, Wang Z, et al. Construction of a lactate-related prognostic signature for predicting prognosis, tumor microenvironment, and immune response in kidney renal clear cell carcinoma. Front In Immunol (2022) 13:818984. doi: 10.3389/fimmu.2022.818984

CrossRef Full Text | Google Scholar

15. Sainero-Alcolado L, Liaño-Pons J, Ruiz-Pérez MV, Arsenian-Henriksson M. Targeting mitochondrial metabolism for precision medicine in cancer. Cell Death Differ (2022) 29(7):1304–17. doi: 10.1038/s41418-022-01022-y

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Frattaruolo L, Brindisi M, Curcio R, Marra F, Dolce V, Cappello AR. Targeting the mitochondrial metabolic network: a promising strategy in cancer treatment. Int J Mol Sci (2020) 21(17). doi: 10.3390/ijms21176014

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Huang C, Lin B, Chen C, Wang H, Lin X, Liu J, et al. Synergistic reinforcing of immunogenic cell death and transforming tumor-associated macrophages Via a multifunctional cascade bioreactor for optimizing cancer immunotherapy. Adv Mater (2022) 34(51):e2207593. doi: 10.1002/adma.202207593

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Zhou Z, Liu Y, Song W, Jiang X, Deng Z, Xiong W, et al. Metabolic reprogramming mediated pd-L1 depression and hypoxia reversion to reactivate tumor therapy. J Control Release (2022) 352:793–812. doi: 10.1016/j.jconrel.2022.11.004

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Galluzzi L, Vitale I, Warren S, Adjemian S, Agostinis P, Martinez AB, et al. Consensus guidelines for the definition, detection and interpretation of immunogenic cell death. J For Immunotherapy Cancer (2020) 8(1). doi: 10.1136/jitc-2019-000337

CrossRef Full Text | Google Scholar

20. Wang Q, Ju X, Wang J, Fan Y, Ren M, Zhang H. Immunogenic cell death in anticancer chemotherapy and its impact on clinical studies. Cancer Lett (2018) 438:17–23. doi: 10.1016/j.canlet.2018.08.028

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Garg AD, De Ruysscher D, Agostinis P. Immunological metagene signatures derived from immunogenic cancer cell death associate with improved survival of patients with lung, breast or ovarian malignancies: a Large-scale meta-analysis. Oncoimmunology (2016) 5(2):e1069938. doi: 10.1080/2162402X.2015.1069938

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Palmeri M, Mehnert J, Silk AW, Jabbour SK, Ganesan S, Popli P, et al. Real-world application of tumor mutational burden-high (Tmb-high) and microsatellite instability (Msi) confirms their utility as immunotherapy biomarkers. ESMO Open (2022) 7(1):100336. doi: 10.1016/j.esmoop.2021.100336

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Lapuente-Santana Ó, van Genderen M, Hilbers P, Finotello F, Eduati F. Interpretable systems biomarkers predict response to immune-checkpoint inhibitors. Patterns (New York NY) (2021) 2(8):100293. doi: 10.1016/j.patter.2021.100293

CrossRef Full Text | Google Scholar

24. Maeser D, Gruener RF, Huang RS. Oncopredict: an r package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data. Briefings In Bioinf (2021) 22(6). doi: 10.1093/bib/bbab260

CrossRef Full Text | Google Scholar

25. Wang Q, Xu J, Xiong Z, Xu T, Liu J, Liu Y, et al. Cenpa promotes clear cell renal cell carcinoma progression and metastasis Via Wnt/B-catenin signaling pathway. J Transl Med (2021) 19(1):417. doi: 10.1186/s12967-021-03087-8

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Shou Y, Liu Y, Xu J, Liu J, Xu T, Tong J, et al. Timp1 indicates poor prognosis of renal cell carcinoma and accelerates tumorigenesis emt signaling pathway. Front Genet (2022) 13:648134. doi: 10.3389/fgene.2022.648134

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Massari F, Santoni M, Ciccarese C, Santini D, Alfieri S, Martignoni G, et al. Pd-1 blockade therapy in renal cell carcinoma: current studies and future promises. Cancer Treat Rev (2015) 41(2):114–21. doi: 10.1016/j.ctrv.2014.12.013

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Lai Y, Tang F, Huang Y, He C, Chen C, Zhao J, et al. The tumour microenvironment and metabolism in renal cell carcinoma targeted or immune therapy. J Cell Physiol (2021) 236(3):1616–27. doi: 10.1002/jcp.29969

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Wettersten HI, Aboud OA, Lara PN, Weiss RH. Metabolic reprogramming in clear cell renal cell carcinoma. Nat Rev Nephrol (2017) 13(7):410–9. doi: 10.1038/nrneph.2017.59

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Herst PM, Grasso C, Berridge MV. Metabolic reprogramming of mitochondrial respiration in metastatic cancer. Cancer Metastasis Rev (2018) 37(4):643–53. doi: 10.1007/s10555-018-9769-2

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Oresta B, Pozzi C, Braga D, Hurle R, Lazzeri M, Colombo P, et al. Mitochondrial metabolic reprogramming controls the induction of immunogenic cell death and efficacy of chemotherapy in bladder cancer. Sci Transl Med (2021) 13(575). doi: 10.1126/scitranslmed.aba6110

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Braun DA, Bakouny Z, Hirsch L, Flippot R, Van Allen EM, Wu CJ, et al. Beyond conventional immune-checkpoint inhibition - novel immunotherapies for renal cell carcinoma. Nat Rev Clin Oncol (2021) 18(4):199–214. doi: 10.1038/s41571-020-00455-z

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Darvin P, Toor SM, Sasidharan Nair V, Elkord E. Immune checkpoint inhibitors: recent progress and potential biomarkers. Exp Mol Med (2018) 50(12). doi: 10.1038/s12276-018-0191-1

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Wu Y, Zhang X, Wei X, Feng H, Hu B, Deng Z, et al. A mitochondrial dysfunction and oxidative stress pathway-based prognostic signature for clear cell renal cell carcinoma. Oxid Med Cell Longev (2021) 2021:9939331. doi: 10.1155/2021/9939331

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Grünwald B, Schoeps B, Krüger A. Recognizing the molecular multifunctionality and interactome of timp-1. Trends Cell Biol (2019) 29(1). doi: 10.1016/j.tcb.2018.08.006

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Bansal M, Gupta A, Ding H-F. Mycn and metabolic reprogramming in neuroblastoma. Cancers (Basel) (2022) 14(17). doi: 10.3390/cancers14174113

PubMed Abstract | CrossRef Full Text | Google Scholar