Constructing a novel mitochondrial-related gene signature for predicting survival and evaluating the tumor immune microenvironment in clear cell renal cell carcinoma

Abstract

Background:

Mitochondria play an important role in tumors. Cellular energy supply, signaling, metabolism, autophagy, aging, and tumorigenesis are all associated with mitochondria. However, we lack a reliable prognostic model using mitochondrial-related genes in clear cell renal cell carcinoma (ccRCC).

Methods:

A systematic analysis of available TCGA databases and related studies was conducted using the R language and online analysis tools to evaluate the prognostic value of mitochondrial-related genes and the tumor microenvironment in ccRCC.

Results:

We constructed a novel mitochondrial-related gene signature for predicting survival and evaluating the tumor immune microenvironment in ccRCC. The mitochondrial-related gene signature included MICALL2, FKBP10, and ACADSB. According to the risk score of the risk model, ccRCC patients were divided into high- or low-risk groups. The ccRCC high-risk group with a high-risk score is related to poor prognosis and poor efficacy from immune checkpoint inhibitors (ICIs).

Conclusion:

Our mitochondrial-related gene signature, as a risk model, could be a reliable ccRCC prognostic biomarker and could predict the response to immunotherapy. The risk score was correlated with the tumor microenvironment and immune cell infiltration.

1 Introduction

Mitochondria play an important role in tumors (Zong et al., 2016). Cellular energy supply, signaling, metabolism, autophagy, aging, and tumorigenesis are all associated with the mitochondria. The tumor microenvironment affects immune cell evasion or inhibition and drug resistance in malignancies. There are interactions between the tumor microenvironment and mitochondrial function in tumor cells (Yang et al., 2023). Different pathological types of renal tumors have their own characteristics of mitochondrial morphology (Nikolic et al., 2023). Clear cell renal cell carcinoma (ccRCC) exhibits a substantial reduction in mitochondrial mass, paralleled by a reduction in both mtDNA copy number and RNA transcription. Some ccRCC samples display a significant elevation in acyl-carnitines, pointing to a ccRCC-specific disruption of mitochondrially localized beta-oxidation, in line with the reduction in mitochondrial mass (DiNatale et al., 2020). The hypoxic background promoted by tyrosine-kinase inhibitors (TKIs) might be altering the tumor microenvironment in a way that facilitates tumor recognition by immune cells, ultimately leading to improved response in renal cell carcinoma (DiNatale et al., 2020). Mitochondrially targeted anticancer drugs that disrupt the energy-producing systems of cancer are emerging as new potential therapeutics in renal cell carcinoma (Bielcikova et al., 2023). However, we lack a reliable prognostic model in ccRCC that uses mitochondrial-related genes.

2 Materials and methods

2.1 Flow chart and main steps

A flow chart of the work is shown in Supplementary Figure S1. The main steps are to identify differentially expressed genes (DEGs) between TCGA ccRCC normal and tumor samples and then intersect the list of DEGs with mitochondrial-related genes (MRGs) to identify common DEGs. A few genes were further selected by Lasso regression analysis in the common DEGs. Then, a prognostic mitochondrial-related gene signature was constructed and validated based on the selected genes. The TCGA ccRCC samples were divided into low- and high-risk groups based on the prognostic mitochondrial-related gene signature. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were conducted. Finally, the tumor microenvironment analysis and immunotherapy responses of the low- and high-risk groups were predicted.

2.2 Data collection

RNA-seq data and clinical information about the ccRCC TCGA training cohort were downloaded from the GDC data portal (https://portal.gdc.cancer.gov/), and 521 ccRCC tumor samples with enough clinical information were used in the survival-related analysis and the tumor immune microenvironment-related analysis. The RNA array assay data and clinical information of the 101 ccRCC tumor samples in the E-MTAB-1980 validation cohort were downloaded from the ArrayExpress database (https://www.ebi.ac.uk/biostudies/arrayexpress/) of EMBL-EBI (European Bioinformatics Institute of the European Molecular Biology Laboratory). The list of 2030 MRGs was identified in a published study (Chang et al., 2023), and the data were collected from the MitoCarta 3.0 database and the Gene Set Enrichment Analyses (GSEA). Protein level expression data, including figures and immunohistochemical (IHC) images, were obtained from CPTAC through UALCAN (https://ualcan.path.uab.edu/) and HPA (https://www.proteinatlas.org/).

2.3 Identification of differentially expressed genes

We used the “DESeq2” package of R to identify DEGs between normal and tumor samples, or between high- and low-risk tumor groups from the TCGA training cohort. Only genes that were expressed >0 in all samples were retained. The criteria for defining DEGs were |Log₂ fold change|>1 and adjusted P < 0.05. For visualizing the DEGs, we used the “ggplot2” package of R to create the volcano plot, which used |Log₂ fold change|>1 and adjusted P < 0.05 as boundaries. The list of 2,955 DEGs was also obtained from GEPIA2 (http://gepia2.cancer-pku.cn/) (Tang et al., 2017) by using the ANOVA method with |Log₂ fold change|>1 and q-value<0.01 for further screening of DEGs. A Venn plot was created using the EVenn online tool (Chen et al., 2021) to display the DEGs and mitochondrial-related genes (MRGs) common to both DEG groups.

2.4 Construction and validation of prognostic mitochondrial-related gene signature

The common DEGs were further screened by Lasso regression analysis using the “glmnet” and “survival” packages of R and multivariable Cox regression analysis (the coxph function in the “survival” package of R). The coefficient index was achieved by using the cph function in the “rms” package of R. The specific multivariable Cox regression cph statement is as follows: “cph (Surv (OS.time, OS)∼gene 1+gene 2+…+gene n, data name).” Risk score = ∑gene n × βn. The “gene” represents the log₂ (TPM+1) value or other normalized data value of genes selected into the mitochondrial-related risk score signature. The “β” represents the coefficient index value. The “n” represents the serial number. The samples were divided into low- and high-risk groups based on the risk score (median cut-off value). To analyze the survival condition for the prognosis signature, a Kaplan–Meier (K-M) survival curve was constructed by using the “survival” and “survminer” packages of R. The predictive performance was presented by the ROC curve. The prognostic mitochondrial-related risk score signature was further validated in the E-MTAB-1980 cohort, and the data of this cohort were normalized through the robust multi-array average (RMA) algorithm.

2.5 GO and KEGG pathway analyses

DEGs between the high- and low-risk groups were input into the DAVID website (https://david.ncifcrf.gov/) for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses. Results were visualized by using the “ggplot2” and “GOplot” packages of R.

2.6 Tumor microenvironment analysis and immunotherapy response prediction

The stromal score, immune score, and ESTIMATE score were calculated by using the “estimate” package of R. Tumor purity = cos (0.6049872018 + 0.0001467884×ESTIMATE score) (Yoshihara et al., 2013). A box plot was created by using the “ggpubr” package of R. A correlation scatter plot was achieved by using the “ggpubr” and “ggExtra” packages of R.

The abundance of 22 tumor-infiltrating immune cells (TIICs) was calculated by using the CIBERSORT algorithm (Newman et al., 2015) in R (https://cibersortx.stanford.edu/), and the results were visualized as a box plot by using the “ggpubr” package of R. The Wilcoxon test was used for assessing between-group differences.

The list of 79 immune checkpoint genes (ICGs) was obtained from a published study (Hu et al., 2021). The statistically significant results of 39 genes were visualized as a box plot by using the “ggpubr” package of R. A t-test was used to determine statistically significant correlations between groups.

The immunotherapy response prediction was analyzed by the tumor immune dysfunction and exclusion (TIDE) score (http://tide.dfci.harvard.edu/login/). A violin plot for TIDE was achieved by using the “ggpubr” package of R. A bar chart to visualize the immunotherapy response prediction results was created by using the “ggplot2,” “ggthemes,” and “ggprism” packages of R.

2.7 Statistics analyses

R (version 4.3.2) software was used for all statistical analyses. Some test methods are annotated in the figures and tables. The t-test was used to analyze the expression and distribution of TIDE score, exclusion score, dysfunction score, and microsatellite instability (MSI) score in different groups. Pearson’s chi-squared test was used to evaluate the immunotherapy response difference in TIDE prediction. The correlation was evaluated using the Spearman method. P < 0.05 was defined as statistically significant. A count value was used in the “DESeq2” package analysis. The TPM value was used in the CIBERSORT algorithm. Data were normalized according to the requirements in the TIDE algorithm. In other analyses, the log₂ (TPM+1) value was used.

3 Results

3.1 Identification of DEGs related to mitochondria

By using the “DESeq2” package of R in the TCGA training cohort, 3,111 DEGs (986 down; 2,125 up) were obtained from 14,449 genes expressed in the normal and tumor groups. The DEGs are visualized by a volcano plot (Figure 1A). A total of 1,658 common DEGs between the 3,111 DEGs from the “DESeq2” package analysis and the 2,955 DEGs from the ANOVA method analysis in GEPIA2 are shown in the Venn plot. We found 164 common DEGs between the 1,658 DEGs and the 2030 MRGs (Figure 1B). In the 164 common DEGs, we ultimately selected 41 genes with |Log₂ fold change|>2 for Lasso regression analysis and multivariable Cox regression analysis. Through Lasso regression analysis (Figure 1C), four genes were filtered out: MICALL2, FKBP10, ACADSB, and ALDH6A1. Through multivariable Cox regression analysis (coxph function in the “survival” packages of R, Figure 2), three mitochondrial-related DEGs, including MICALL2, FKBP10, and ACADSB, were finally selected to build a ccRCC prognostic model. The RNA level expression (Figure 1D) and protein level expression (Figures 1E,F) of three DEGs between normal samples and tumor samples are shown in Figure 1. In renal cell carcinoma (unknown pathological subtype) IHC images from HPA (Figure 1F), eight with MICALL2 medium staining and four with MICALL2 low staining; six with FKBP10 medium staining and six without FKBP10 staining; nine without ACADSB staining, one with ACADSB low staining, one with ACADSB medium staining, and one with ACADSB high staining. The information about the three genes is shown in Table 1.

FIGURE 1

FIGURE 2

The effects of age, ethnicity, gender, race, stage, ALDH6A1, FKBP10, MICALL2, and ACADSB were all involved in the forest plot (Figure 2), and we obtained a list of more variable coefficient index values from this analysis. Age, ethnicity, stage, and MICALL2 showed statistically significant differences. The P-values of FKBP10 and ACADSB were less than 0.1. Meanwhile, we also obtained a list of multivariable coefficient index values by using the cph function in the “rms” package of R, and only the effects of MICALL2, FKBP10, and ACADSB were involved in this analysis.

3.2 Construction and validation of prognostic mitochondrial-related gene signature

The more variable coefficient index value and the multivariable coefficient index value were validated separately through a Kaplan–Meier (K-M) survival curve and the ROC curve. Risk score (more variable coefficient index value, Supplementary Figure S2) = 0.379045×MICALL2 + 0.127374×FKBP10 − 0.249544×ACADSB. Risk score (multivariable coefficient index value, Figures 3A,B) = 0.3089×MICALL2 + 0.1995×FKBP10 − 0.4429×ACADSB. The samples were divided into low- and high-risk groups based on the risk score (median cut-off value). In Supplementary Figure S2, we show the intersection between the K–M survival curve and the lower AUC values of the ROC curve. Finally, a multivariable coefficient index value was selected. Then, the risk score for each ccRCC patient in both the training and the validation cohorts was computed based on the following formula:

FIGURE 3

The stability of the prognostic risk model was further validated in the E-MTAB-1980 cohort, and similar results are shown in the K–M survival curve and the ROC curve (Figures 3C,D).

Clinical characteristics of the low- and high-risk groups in the TCGA 521 ccRCC cohort are shown in Table 2, which was created using the “gtsummary” package of R with univariate Cox regression analysis. The high-risk group had a higher proportion of men, T3+T4 stage, M1 stage, stage III + IV, and dead status.

Multivariate Cox regressions of various prognostic parameters in the TCGA 521 ccRCC cohort are shown in Table 3. Through multivariable Cox regression analysis (coxph function in “survival” packages of R) and the “gtsummary” package of R, the effects of age, gender, race, ethnicity, detailed stage subgroups, and risk group were examined. Stage IV and high-risk groups are related to an unfavorable prognosis. Non-Hispanic or non-Latino ethnicity and high age might be related to unfavorable prognosis; however, the left 95% confidence interval (CI) of age is near 1, and the P-value of race is not statistically significant.

3.3 GO and KEGG pathway analyses

A total of 484 DEGs between the high- and low-risk groups were selected in the TCGA 521 ccRCC cohort and were input into the DAVID website. GO and KEGG pathway analyses results are shown in Figure 4. ECM–receptor interaction and focal adhesion pathway, which are involved in tumor invasion and tumor microenvironment, were enriched in KEGG analysis. JAK-STAT, PI3K-Akt, NF-kappa B, and HIF-1 signaling pathways were also enriched in KEGG analysis. Some extracellular-related functions, some immune-related functions, and cytokine-mediated signaling pathways were enriched in the GO analysis.

FIGURE 4

3.4 Tumor microenvironment analysis and immunotherapy response prediction

In the high-risk group, we found higher stromal scores, immune scores, and ESTIMATE scores and lower tumor purity. The risk score and stromal score, immune score, and ESTIMATE score were positively correlated (Figures 5A,F). This means the infiltration of immune cells and stromal cells was higher in the high-risk group.

FIGURE 5

In 22 TIICs, the contents of regulatory T cells (Treg), M0 macrophages, and CD8⁺ T cells were remarkably higher in the high-risk group. The numbers of M1 macrophages, M2 macrophages, and monocytes were higher in the low-risk group (Figure 5B). M2 macrophages might be correlated with pro-tumor activity, whereas M1 macrophages might be correlated with antitumor activity. CD8⁺ T cells are the major antitumor cells. Treg cells can inhibit the anti-tumor immune effect (Dou and Fang, 2021; Peña-Romero and Orenes-Piñero, 2022). In summary, the difference between high- and low-risk groups in the tumor immune microenvironment was complex.

A total of 39 of 79 ICGs showed a statistically significant difference (Figure 5C). PDL1 and PDL2 were lower in the high-risk group, which indicates poor immunotherapy response.

In TIDE immunotherapy response prediction, the TIDE score and dysfunction score were higher in the high-risk group. The MSI (microsatellite instability) score was lower in the high-risk group (Figure 5E), and the response percentage was lower in the high-risk group (Figure 5D), which indicates poor immunotherapy response.

4 Discussion

MICALL2 is considered a typical cell mobility-promoting factor. Overexpression of MICALL2 could promote the proliferation, migration, and invasion of ccRCC cell lines (Zeng et al., 2023). FKBP10 knockdown could cause cell cycle arrest at the G0/G1 phase and reduce cell proliferation, invasion, and migration in 786-O and A-704 cell lines (Ge et al., 2017). FKBP10 might participate in the process of type I collagen synthesis in ccRCC via regulating crosslinking of pro-collagen chains (Zhang et al., 2021). FKBP10 might promote ccRCC progression and regulate sensitivity to HIF2α blockade by facilitating LDHA phosphorylation (Liu et al., 2024). ACADSB is involved in the metabolism of fatty acids and branched-chain amino acids and plays an important role in glioma, colorectal cancer, and hepatocellular carcinoma (Liu et al., 2022).

The mitochondrial-related gene signature based on MICALL2, FKBP10, and ACADSB could provide prognostic predictive value in ccRCC. The ccRCC patients in the high-risk group had an unfavorable prognosis. DEGs between the high- and low-risk groups were involved in pathways about extracellular-related function, immune-related function, or cytokine-mediated signaling pathways. Those pathways were associated with tumor invasion, tumor microenvironment, tumor immunotherapy response, or tumor metabolism. ESTIMATE analysis showed that the infiltration of immune cells and stromal cells was higher in the high-risk group. More Treg cells, fewer M1 macrophages, lower PDL1 expression levels, and lower PDL2 expression levels in the high-risk group might be associated with the poorer immunotherapy response. In addition, more elaborate analyses, such as single-cell data analyses, might be more informative. New methods of single-cell data analysis are showing potential (Yuan et al., 2024a; Yuan et al., 2024b; Yuan et al., 2025).

5 Conclusion

Our results suggested that our mitochondrial-related gene signature based on MICALL2, FKBP10, and ACADSB, as a risk model, could be a reliable ccRCC prognostic biomarker and could predict the response to immunotherapy. The risk score was correlated with the tumor microenvironment and immune cell infiltration.

References

• 1

  BielcikovaZ.WernerL.StursaJ.CernyV.KrizovaL.SpacekJ.et al (2023). Mitochondrially targeted tamoxifen as anticancer therapy: case series of patients with renal cell carcinoma treated in a phase I/Ib clinical trial. Ther. Adv. Med. Oncol.15, 17588359231197957. 10.1177/17588359231197957

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 2

  ChangJ.WuH.WuJ.LiuM.ZhangW.HuY.et al (2023). Constructing a novel mitochondrial-related gene signature for evaluating the tumor immune microenvironment and predicting survival in stomach adenocarcinoma. J. Transl. Med.21 (1), 191. 10.1186/s12967-023-04033-6

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 3

  ChenT.ZhangH.LiuY.LiuY. X.HuangL. (2021). EVenn: easy to create repeatable and editable Venn diagrams and Venn networks online. J. Genet. Genomics48 (9), 863–866. 10.1016/j.jgg.2021.07.007

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 4

  DiNataleR. G.SanchezA.HakimiA. A.ReznikE. (2020). Metabolomics informs common patterns of molecular dysfunction across histologies of renal cell carcinoma. Urol. Oncol.38 (10), 755–762. 10.1016/j.urolonc.2019.04.028

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 5

  DouA.FangJ. (2021). Heterogeneous myeloid cells in tumors. Cancers (Basel)13 (15), 3772. 10.3390/cancers13153772

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 6

  GeY.XuA.ZhangM.XiongH.FangL.ZhangX.et al (2017). FK506 binding protein 10 is overexpressed and promotes renal cell carcinoma. Urol. Int.98 (2), 169–176. 10.1159/000448338

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 7

  HuF. F.LiuC. J.LiuL. L.ZhangQ.GuoA. Y. (2021). Expression profile of immune checkpoint genes and their roles in predicting immunotherapy response. Brief. Bioinform22 (3), bbaa176. 10.1093/bib/bbaa176

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 8

  LiuX.ZhangW.WangH.ZhuL.XuK. (2022). Decreased expression of ACADSB predicts poor prognosis in clear cell renal cell carcinoma. Front. Oncol.11, 762629. 10.3389/fonc.2021.762629

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 9

  LiuR.ZouZ.ChenL.FengY.YeJ.DengY.et al (2024). FKBP10 promotes clear cell renal cell carcinoma progression and regulates sensitivity to the HIF2α blockade by facilitating LDHA phosphorylation. Cell Death Dis.15 (1), 64. 10.1038/s41419-024-06450-x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 10

  NewmanA. M.LiuC. L.GreenM. R.GentlesA. J.FengW.XuY.et al (2015). Robust enumeration of cell subsets from tissue expression profiles. Nat. Methods12 (5), 453–457. 10.1038/nmeth.3337

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 11

  NikolicG.ZivoticM.CirovicS.DespotovicS.DundjerovicD.Radojevic SkodricS. (2023). The utility of mitochondrial detection methods applied as an additional tool for the differentiation of renal cell tumors. Diagn. (Basel)13 (14), 2319. 10.3390/diagnostics13142319

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 12

  Peña-RomeroA. C.Orenes-PiñeroE. (2022). Dual effect of immune cells within tumour microenvironment: pro- and anti-tumour effects and their triggers. Cancers (Basel)14 (7), 1681. 10.3390/cancers14071681

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 13

  TangZ.LiC.KangB.GaoG.LiC.ZhangZ. (2017). GEPIA: a web server for cancer and normal gene expression profiling and interactive analyses. Nucleic Acids Res.45 (W1), W98–W102. 10.1093/nar/gkx247

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 14

  YangY.AnY.RenM.WangH.BaiJ.DuW.et al (2023). The mechanisms of action of mitochondrial targeting agents in cancer: inhibiting oxidative phosphorylation and inducing apoptosis. Front. Pharmacol.14, 1243613. 10.3389/fphar.2023.1243613

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 15

  YoshiharaK.ShahmoradgoliM.MartínezE.VegesnaR.KimH.Torres-GarciaW.et al (2013). Inferring tumour purity and stromal and immune cell admixture from expression data. Nat. Commun.4, 2612. 10.1038/ncomms3612

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 16

  YuanL.SunS.JiangY.ZhangQ.YeL.ZhengC. H.et al (2024a). scRGCL: a cell type annotation method for single-cell RNA-seq data using residual graph convolutional neural network with contrastive learning. Brief. Bioinform26 (1), bbae662. 10.1093/bib/bbae662

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 17

  YuanL.ZhaoL.JiangY.ShenZ.ZhangQ.ZhangM.et al (2024b). scMGATGRN: a multiview graph attention network-based method for inferring gene regulatory networks from single-cell transcriptomic data. Brief. Bioinform25 (6), bbae526. 10.1093/bib/bbae526

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 18

  YuanL.XuZ.MengB.YeL. (2025). scAMZI: attention-based deep autoencoder with zero-inflated layer for clustering scRNA-seq data. BMC Genomics26 (1), 350. 10.1186/s12864-025-11511-2

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 19

  ZengX.WangH.YangJ.HuJ. (2023). Micall2 is responsible for the malignancy of clear cell renal cell carcinoma. Yonago Acta Med.66 (1), 171–179. 10.33160/yam.2023.02.021

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 20

  ZhangY.YinY.LiuS.YangL.SunC.AnR. (2021). FK506 binding protein 10: a key actor of collagen crosslinking in clear cell renal cell carcinoma. Aging (Albany NY)13 (15), 19475–19485. 10.18632/aging.203359

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 21

  ZongW. X.RabinowitzJ. D.WhiteE. (2016). Mitochondria and cancer. Mol. Cell.61 (5), 667–676. 10.1016/j.molcel.2016.02.011

  • Pubmed Abstract
  • CrossRef
  • Google Scholar