Skip to main content
Download PDF

Development and validation of a novel mitophagy-related gene prognostic signature for glioblastoma multiforme

BMC Cancer volume 22, Article number: 644 (2022) Cite this article

Abstract

Background

Glioblastoma multiforme (GBM) is one of the most malignant tumors in brain with high morbidity and mortality. Mitophagy plays a significant role in carcinogenesis, metastasis, and invasion. In our study, we aim to construct a mitophagy-related risk model to predict prognosis in GBM.

Methods

RNA-seq data combined with clinical information were downloaded from TCGA. The 4-gene risk model and nomograph was then constructed and validated in external cohort. Evaluation of immune infiltration, functional enrichment and tumor microenvironment (TME) were then performed.

Result

A mitophagy-related risk model was established and patients in TCGA and CGGA were classified into low-risk and high-risk groups. In both cohorts, patients in low-risk group had improved survival, while high-risk group had poor prognosis. Also, the risk model was identified as an independent factor for predicting overall survival via Cox regression. Furthermore, a prognostic nomogram including mitophagy signatures was established with excellent predictive performance. In addition, the risk model was closely associated with regulation of immune infiltration as well as TME.

Conclusion

In conclusion, our study constructed a mitophagy-related risk model, which can be utilized for the clinical prognostic prediction in GBM.

Peer Review reports

Introduction

Gliomas are the most common brain tumors in central nervous systems and represent 75% malignant brain cancer in adults. Glioblastoma (GBM), also known as glioblastoma multiforme, is the most lethal type in gliomas (World Health Organization [WHO] grade IV), with a 5-year relative survival of 3 ~ 7% [1, 2]. Despite of the progress in multimodality therapies including surgical resection, radiotherapy, systemic therapies (chemotherapy, targeted therapy) and supportive care, the median survival time of GBM is less than 2 years [3]. These conditions suggest that the conventional stage system for predicting prognosis such as WHO grade is insufficient to cover the clinical diversity of GBM. Thus, a novel prognostic model for GBM is needed to be established.

The role of autophagy in malignant metastasis and response to treatment has been elucidated in many types of tumors, which is characterized by degrading dysfunctional and cellular materials or damaged organelles through lysosomal system to maintain cellular homeostasis [4, 5]. Mitophagy, as a specific type of autophagy, is a fundamental process of removing damaged or excessive mitochondria via autophagolysosomes. This process is often triggered by oxidative stress or increasing need of bioenergy, which is necessary for cancer formation and invasion [6]. Some studies have discovered the role of mitophagy in gliomas [7, 8]. Huang et al. found that inhibition of mitophagy partially reverted cannabidiol-induced glioma cell death, suggesting the positive role of mitophagy on anti-tumor [9]. Two studies demonstrated that induction of mitophagy by FOXO3a protect the gliomas from temozolomide-induced cytotoxicity, indicating double-sword effect of mitophagy on glioma [10, 11]. Thus, whether mitophagy is correlated to the prognosis of GBM, and the possible involvement of mitophagy-related genes remain to be explored.

Given the existing findings, we constructed a mitophagy-related prognostic model based on mRNA expression and clinical data of GBM patients from TCGA, and we also explored the correlations between mitophagy and the tumor immune microenvironment.

Materials and methods

Datasets

RNA sequencing (RNA-seq) data of 169 GBM patients with 5 normal controls were extracted from TCGA database on 1 August 2021 (https://portal.gdc.cancer.gov/repository). The validation data were extracted from CGGA database mRNAseq_693 on 1 August 2021 with a total of 134 GBM patients [12, 13] (http://www.cgga.org.cn/).

Identification of mitophagy-related differential expressed genes (DEGs)

Fifty-one mitophagy-related genes were extracted from Gene Set Enrichment Analysis (http://www.gsea-msigdb.org/gsea/index.jsp) combined with GeneCards (https://www.genecards.org/). Then mRNA expression data were compared utilizing “limma” package between 169 GBM patients and 5 controls [14]. DEGs were screened out with the following cut off: P value < 0.05. PPI network of DEGs was constructed using Search Tool for the Retrieval of Interacting Genes (STRING), version 11.0 (https://string-db.org/).

Establishment and validation of risk model for prognosis

To further screen out DEGs with high predicting value, we performed Cox regression analysis to evaluate the relation between DEGs and survival status in TCGA cohort. 5 mitophagy-related genes were firstly identified for further analysis with P < 0.1. LASSO regression analysis was performed via R package “glmnet” to identify prognosis-related DEGs and developed risk model [15]. Finally, 4 genes were selected as the optimal gene to construct risk score. The calculation of risk score was as follows: Risk Score = ∑ \({\sum}_i^4 Xi\times Yi\) (X: coefficients, Y: gene expression level). According to median value of the risk score, patients in TCGA were well classified into low-risk and high-risk groups and the survival status of 2 groups were further compared using Kaplan–Meier analysis. PCA of 4-genes risk model was constructed by R package “Rtsne”. The time-dependent ROC curve was applied to evaluate the predicting efficacy of risk score utilizing R “survivalROC” package [16]. The images of immunohistochemistry (IHC) staining of prognosis-related genes were extracted from The Human Protein Atlas database [17]. To analyze independent prognosis value, univariate and multivariable Cox regression models were employed in TCGA cohort and CGGA cohort. To validate the risk model, CGGA cohort (mRNAseq_693) was applied. 4-genes risk signature was then calculated by the same method in TCGA cohort and patients in CGGA were divided into low-risk and high-risk groups and the difference was further analyzed according to method mentioned above. The nomograph was constructed using R package “rms”, “foreign”, “survival” and evaluated by ROC curve. The Calibration plots were applied to evaluate the discriminative ability of the nomogram.

Functional enrichment based on GO and KEGG

Gene Ontology analysis was performed to enrich the biological processes, molecular functions and cellular components based on DEGs between low-risk and high-risk groups. And the Kyoto Encyclopedia of Genes and Genomes pathway analysis were conducted to determine the pathway corelated to the Mitophagy-related risk signature. GO and KEGG analysis both used R package “clusterProfiler” [18].

Estimation of immune infiltration and immune related pathway

The scores of infiltrating immune cells and immune-related pathways were evaluated using R package “gvsa” in TCGA cohort and CGGA cohort [19]. Furthermore, the relationships between immune cells and mitophagy-related genes were analyzed through TIMER database [20] (https://cistrome.shinyapps.io/timer/).

Tumor microenvironment analysis

The 22 types of immune cells composition in subgroups were analyzed using CIBERSORT algorithm (https://cibersort.stanford.edu/). The correlations between 22 immune cells composition were evaluated via “corrplot” R package.

Statistical analysis

R software (version 4.1.0), SPSS (version 23.0) and R studio (version 1.1.463) were used to perform statistical analysis. Single-factor analysis was utilized to compare the different expression of genes. Person chi-square test was used to evaluate the categorical variables. Kaplan-Meier method combined with log-rank test was used to analyze the overall survival (OS) of GBM patients. Univariate and multivariate Cox regression were used to assess the independent prognostic value of the risk model. The analysis of immune cell infiltration and immune pathways were conducted via Wilcoxon test.

Results

Identifying prognostic-related DEGs that were associated with mitophagy

A total of 5 normal and 169 GBM samples with mRNA expression profile and clinical data were extracted from TCGA. The 51 mitophagy-related genetic expressions were compared between normal group and tumor group via “limma” package. Twenty-seven mitophagy-related DEGs were identified with P < 0.05, among which 15 genes (GABARAPL1, PRKN, OPTN, PINK1, MAP1LC3A, ULK1, DNM1L, AMBRA1, TOMM20, ATG13, MFN2, PTEN, TOMM70, RNF41, VPS13C) were downregulated while 12 genes (TOMM40, GABARAP, CSNK2B, TOMM22, TOMM5, ROCK1, TOMM7, RIPK2, MTERF3, UBA52, PHB2, RPS27A) were upregulated. The heatmap of 51 DEGs was illustrated in Fig. 1A. To further explore the interaction between DEGs, we performed protein–protein interaction (PPI) on mitophagy-related gene (Fig. 1B). We identified 28 hub genes with the minimum required interaction score of 0.900 (highest confidence), among which 10 genes (ATG 13, CSNK2B, DNM1L, GABARAP, MAP1LC3A, MFN2, MTERF3, OPTN, PINK1, PTEN) were DEGs (Fig. 1B). The correlation network of 33 mitophagy-related genes expressions with significant difference is further shown in Fig. 1C (caption: red = positive correlations; blue = negative correlations).

Fig. 1
figure 1

Expression pattern of 51 mitophagy-related genes between GBM and normal patients in TCGA cohort. A Heatmap indicating different expression genes (DEGs) between tumor (red) and normal sample (blue). * P < 0.05, **P < 0.01, *** P < 0.001. B PPI network for the interaction between mitophagy-related genes. (Interaction score = 0.9). C The correlation network of the 33 DEGs. (Red line: positive relation, blue line: negative relation)

Full size image

Consensus clustering analysis of GBM based on different expression pattern

To further classify GBM subtypes based on their distinct expression pattern of mitophagy-related DEGs, we constructed consensus clustering analysis in GBM patient in the TCGA cohort. By using ConsensusClusterPlus package based on 27 mitophagy-related DEGs, we identified 2 different regulation patterns (k = 2) including 139 cases in cluster 1 and 30 in cluster 2 with the highest intragroup correlations (Fig. 2A). Then we integrated the cluster with the mRNA expression level and clinical informations including age (≤60 or > 60 years), sex (male or female) and survival status (dead or alive), which are illustrated in heatmap (Fig. 2B).

Fig. 2
figure 2

Classification of tumor samples based on mitophagy-related DEGs. A 159 GBM patients were separated into two clusters based on the consensus clustering matrix (k = 2). B Heatmap of two clusters with clinical information

Full size image

Construction of risk signature based on the TCGA cohort

After deleting duplicate and missing value, a total of 159 GBM samples were match with clinical information. To screen out candidates for constructing risk signature, we performed univariate Cox regression analysis on DEGs and identified 5 genes (MAP1LC3A, TOMM20, TOMM22, PHB2, UBA52) with the criteria of P < 0.1 (Fig. 3A). Among them, 1 gene (MAP1LC3A) was related to increased risk with HRs > 1, while the other 4 genes (TOMM20, TOMM22, PHB2, UBA52) were associated with protective effect with HRs < 1 (Fig. 3A). Next, least absolute shrinkage and selection operator (LASSO) Cox regression analysis was used to identified 4 survival-related genes (MAP1LC3A, TOMM20, PHB2, UBA52) to construct risk model for prognosis according to the optimum λ value (Fig. 3B, C). The risk score formula was as follows: risk score = (0.0272 * MAP1LC3A exp.) + (− 0.00546 * PHB2 exp.) + (− 0.008623 * TOMM20 exp.) + (− 0.002654 * UBA52 exp.). The downregulation of MAP1L3A, TOMM20 and upregulation of PHB2, UBA52 in GBM tissue were further validated by The Human Protein Atlas database (Fig. 3D).

Fig. 3
figure 3

Construction of risk model in TCGA cohort based on mitophagy-related genes. A Univariate cox regression analysis of OS for 5 mitophagy-relation genes with P < 0.1. B least absolute shrinkage and selection operator (LASSO) regression of 5 OS-related genes. Black: MAP1LC3A, Blue: TOMM 22, Baby blue: PHB2, Red: TOMM 20, Green: UBA52. C Cross-validation in LASSO regression. D Immunohistochemistry staining of normal and tumor sample. Positive area was measured by image J

Full size image

The 159 samples were classified into high-risk group (79) and low-risk group (80) (Fig. 4A). Principal component analysis (PCA) showed that GBM patients were well divided into 2 clusters (Fig. 4B). As shown in Fig. 4C, patients in high-risk group exhibited more death (red dot) and lower survival time compared to low-risk group. Furthermore, the difference of OS between high-risk and low-risk groups was prominent (P = 0.0033) (Fig. 4D). Time-dependent ROC curves were utilized to evaluate the predictive performance of the model and the area under the curve (AUC) was 0.635 for 1-year, 0.75 for 2-year, and 0.724 for 3-year survival (Fig. 4E).

Fig. 4
figure 4

Predicting value of risk model in TCGA cohort. A Distribution of patients according to risk model. B Principal component analysis (PCA) plot of GBM patient based on risk model. C The survival status and time of patient with GBM. (Low-risk: on the left side of dotted line. High-risk on the right side of dotted line). D Kaplan-meier (KM) survival analysis of patients in low-risk and high-risk groups. E ROC curve for predicting efficiency of the risk model. F Univariate Cox regression analysis of risk score in TCGA cohort. G Multivariate Cox regression analysis of risk score in TCGA cohort. H Heatmap of patients in different risk group with clinical information. (Red: positive expression. Blue: negative expression)

Full size image

To evaluate whether risk score was an independent risk factor, we utilized univariate and multivariable Cox regression analyses in TCGA cohort. In univariate Cox regression, the risk model had a significant relationship with OS (HR = 1.7566, 95%CI (1.1988–2.5739), P < 0.005) (Fig. 4F). The multivariate analysis also indicated the risk model was an important factor correlated to OS (HR = 1.6734, 95%CI (1.1368–2. 4633), P < 0.01) (Fig. 4G). Furthermore, the heatmap of clinical features showed that the survival status was better in the low-risk group (Fig. 4H).

Validation of the prognostic model in the Chinese Gliomas Genome Atlas (CGGA) databases

After screening out, 134 GBM patient’s samples with clinical informations were extracted from CGGA mRNAseq_693. According to the median risk score built in TCGA cohort, patients were classified into high-risk (66) and low-risk (67) groups (Fig. 5A). PCA analysis indicated that patients were separated well into two groups (Fig. 5B). The dot line plot showed that patients with low-risk exhibited better survival time and survival rate (Fig. 5C). Furthermore, Kaplan–Meier analysis also showed significant difference in survival status between two groups (P = 0.0175) (Fig. 5D). Time-dependent ROC curves indicated that the 1-year, 2-year, 3-year AUC were 0.603, 0.709, 0.653 respectively, suggesting satisfactory predicting efficacy of the risk score model (Fig. 5E).

Fig. 5
figure 5

Validation of risk model in CGGA cohort. A Distribution of patients according to risk model in CGGA cohort. B Principal component analysis (PCA) plot of GBM patient based on risk model. C The survival status and time of patient with GBM. (Low-risk: on the left side of dotted line. High-risk on the right side of dotted line). D Kaplan-meier (KM) survival analysis of patients in low-risk and high-risk groups. E ROC curve for predicting efficiency of the risk model. F Univariate Cox regression analysis of risk score in CGGA cohort. G Multivariate Cox regression analysis of risk score in CGGA cohort

Full size image

To evaluate whether risk score was an independent risk factor, we utilized univariate and multivariable Cox regression analyses in CCGA cohort. In univariate Cox regression, the risk model had a significant relationship with OS (HR = 1.5757, 95%CI (1.0786–2.3019), P = 0.0187) (Fig. 5F). The multivariate analysis also indicated the risk model was an important factor correlated to OS (HR = 1.5702, 95%CI (1.0599–2. 3262), P = 0.0245) (Fig. 5G).

Establishment and validation of nomograph

To construct a clinical-based method for predicting the prognosis of GBM patients, we established a nomograph based on 4 prognostic parameters including age, sex, radiation therapy, riskscore (Fig. 6A). The calibration plots indicated excellent agreement between the predicted and actual observation in both training and validation cohort (Fig. 6B, C). The AUC of nomograph in predicting 1-year, 2-year, 3-year survival were 0.693, 0.726, 0.731 respectively in TCGA cohort and were 0.637, 0.704, 0.681 respectively in CGGA cohort, indicating a favorable predictive ability of nomograph (Fig. 6D, E).

Fig. 6
figure 6

Establishment of nomograph for predicting 1-year, 2-year and 3-year survival probability. A Construction of prognostic nomograph to predict survival of GBM in TCGA cohort. B Calibration curves of the prognostic nomogram for predicting 1-year, 2-year and 3-year survival probability in the TCGA cohort. C Calibration curves of the prognostic nomogram for predicting 1-year, 2-year and 3-year survival probability in the CGGA cohort. The prognostic value of the nomogram evaluated by ROC curve in training (D) and validation cohort (E)

Full size image

The GO and KEGG enrichment analysis based on risk model

To further clarify the pathway and the biological function of DEGs according to risk model, we performed Gene ontology (GO) enrichment analysis and Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathway analysis on DEGs between low-risk and high-risk groups in TCGA cohort. First, we applied limma package in R language to identify 187 DEGs between two groups with the criteria FDR < 0.05 and |log2FC | ≥ 1. In 187 DEGs, 27 genes were downregulated in high-risk group, while the other 160 genes were upregulated. GO and KEGG enrichment analysis were applied on DEGs. According to GO analysis, DEGs based on risk model were mainly enriched in neurotransmitter transport, presynapse, ion channel activity (P < 0.005) (Fig. 7A, B). KEGG analysis showed that DEGs were closely associated with neuroactive ligand−receptor interaction and calcium signaling pathway (P < 0.01) (Fig. 7C, D).

Fig. 7
figure 7

Functional enrichment of DEGs based on low-risk and high-risk groups. A, B Gene Ontology (GO) enrichment analysis of DEGs. C, D Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of DEGs

Full size image

Evaluation of the immune activity between subgroups

As immune micro-environmental abnormality plays an important part in oncogenesis, invasion and metastasis, we further evaluated the enrichment of 16 types of immune cells and 13 immune-related pathways in TCGA and CGGA in different risk score via single-sample gene set enrichment analysis (ssGSEA) package. As shown in TCGA boxplot, patients in high-risk group showed lower level of immune cells (CD8+ T cells, natural killer (NK) cells, follicular helper T cell (Tfh) cells, Th2 cell) than those in low-risk group (P < 0.05) (Fig. 8A), while the difference in immune pathway is not prominent (Fig. 8B).

Fig. 8
figure 8

Differences of immune cells infiltration and immune pathway based on risk score. A, B ssGSEA scores of immune cells (A.16 immune cells) and immune pathway (B.13 immune pathways) in TCGA cohort. C, D single-sample gene set enrichment analysis (ssGSEA) scores of immune cells (A.16 immune cells) and immune pathway (B.13 immune pathways) in CGGA cohort. * P < 0.05, ** P < 0.01, *** P < 0.001

Full size image

In CGGA cohort, high-risk group had lower infiltration of immune cells (CD8+ T cells, NK cells, activated dendritic cells (aDCs), Mast cells, Neutrophils), which was almost the same as TCGA (Fig. 8C). Additionally, APC co stimulation, cytolytic, inflammation promoting, Type I IFN response were enriched in low-risk group, indicated the low immune activity in high-risk group (Fig. 8D). To further analyze the relationship between mitophagy-related DEGs (MAP1LC3A, TOMM20, PHB2, UBA52) and immune cells, the TIMER 2.0 was applied (Fig. 9). As shown in Fig. 9, TOMM20 was positively associated with CD8+ T cells (P < 0.0001) and negatively related to B cells and CD4+ cells (P = 0.043, 0.0052). UBA52 had positively relationship with macrophage, neutrophil, dendritic cell (DC) (P = 0.00076, 0.0045, 0.017). PHB2 had negatively relationship with CD4+ cells (P = 0.016). These findings suggested that the risk model genes were closely associated with immune infiltration.

Fig. 9
figure 9

TIMER analysis of correlations between prognostic genes and immune cell infiltration. A, B, C, D Relationship between MAP1LC3A, TOMM20, PHB2, UBA52 and immune cell infiltration

Full size image

Evaluation of the tumor microenvironment between subgroups

Due to the close relationship between mitophagy-related risk model and immune activity, CIBERSORT was utilized to analysis the tumor microenvironment (TME) between low-risk and high-risk groups. The overview of 22 immune cell compositions in GBM samples was shown in Fig. 10A. The high-risk group had higher proportion of resting NK cells as well as plasm cells, while low-risk group had increased proportion of activated NK cells, indicating that the activation of NK cells regulated anti-tumor effects in low-risk group (Fig. 10B, C).

Fig. 10
figure 10

CIBERSORT analysis of immune cell composition in TCGA. A Overview of immune cell composition in GBM sample. B immune cell composition in high-risk and low-risk groups

Full size image

Discussion

The heterogeneity of GBM and the lack of effective stage make it important to develop stable prognostic model. Our current study firstly analyzed the relationship between 51 mitophagy-related genes and the prognosis of GBM in TCGA cohort and constructed a mitophagy-related prognostic predicting model containing 4 genes (MAP1LC3A, TOMM20, PHB2, UBA52), which was validated well in external databases. Also, we then constructed a mitophagy-based nomograph with potent predicting value. Furthermore, GO and KEGG pathway enrichment analysis uncovered that the mitophagy-related genes are associated with synaptic activity, ion channel, neuroactive ligand-receptor interaction and calcium signaling. In addition, the mitophagy genes were closely related to immune infiltration and TME.

Mitophagy, an evolutionarily conserved programmed cell death essential for cellular homeostasis, is an autophagic response targeting damaged, dysfunction mitochondria. Despite numerous studies have pointed out that the autophagy plays a vital part in tumorigenesis, metastasis and drug resistance, limited studies reported about the role of mitophagy in cancer especially gliomas [6, 21]. In our research, we systematically combined mitophagy-related genes with survival data and constructed a 4-gene prognostic predicting model to predict OS in GBM patients, suggesting a strong relationship between GBM and mitophagy.

Four important mitophagy-related genes (MAP1LC3A, TOMM20, PHB2, UBA52) were identified in our study. Microtubule Associated Protein 1 Light Chain 3 Alpha (MAP1LC3A) encodes 2 different isoforms: MAP1A and MAP1B. MAP1A and MAP1B are two microtubule-related proteins mediating the physical interactions between microtubules and components of the cytoskeleton. Compelling studies reported that the expression of MAP1LC3A was suppressed in many tumor cells including GBM, indicating that it might be involve in the tumorigenesis of various cancers [22, 23]. Wang reported that upregulation of MAP1LC3A in GBM could predict poor prognosis [24]. Consistent with the studies above, we observed that the MAP1LC3A was downregulated in GBM samples compared to normal samples and it was enriched in high-risk group. Further Cox regression analysis showed that high expression of MAP1LC3A was related to poor survival, suggesting that dysregulation of mitophagy promoted the tumor progression. Translocase Of Outer Mitochondrial Membrane 20 (TOMM20) is a pre-protein receptor on the translocation complex of the mitochondrial outer membrane. According to previous studies, upregulation of TOMM20 could be observed in many cancers such as hepatocarcinoma [25,26,27,28]. However, the expression of TOMM20 was lower in GBM compared to normal samples in our study and it was highly expressed in low-risk group. One possible explanation is that TOMM20 upregulated the CD8+ T cell, providing tumor-suppressive effect in our study based on TIMER analysis (Fig. 9). Prohibitin 2 (PHB2) is highly conserved protein mainly in mitochondria, nucleus and plasma membrane, which is required for Parkin-induced mitophagy in mammalian cells and cancer cell proliferation and adhesion [29]. Overexpression of PHB2 had been reported in hepatocarcinoma, breast cancer as well as lung cancer [30,31,32]. PHB2 was upregulated in GBM patients and enriched mainly in low-risk group in our study. The current discrepancy still urged further exploration on the two-edge sword effect of mitophagy. Ubiquitin A-52 Residue Ribosomal Protein Fusion Product 1 (UBA52) is a ubiquitin coding gene encoding a ubiquitin fusion protein which is comprised of ribosomal protein L40 at C-terminus and ubiquitin at the N-terminus [33]. Ubiquitin is closely related to cell cycle regulation and lysosomal degradation. Also, UBA52 was observed in tumor tissues [34, 35]. It is reported that over-expression of UBA52 induced cellular apoptosis in tumor tissue [36] and UBA52 was mainly enriched in low-risk group in our study.

Until now, the role of mitophagy in tumor has not been fully understood. On the one hand, inhibition of mitophagy suppressed the growth of glioblastoma cells [7]. On the other hands, lethal mitophagy can inhibited the proliferation of glioma [9]. According to data from TCGA cohort, about 56% mitophagy-related DEGs were downregulated with 44% upregulated in tumor sample, indicating the two-edges sword effect of mitophagy, which is consistent with the findings above. Furthermore, the TOMM20 can facilitate ROS-induced pyroptosis and PHB2 is associated with apoptosis, indicating the interaction and coexistence of different programmed cell death as tumor grows [37].

Through LASSO and Cox regression analysis, we then constructed a mitophagy-related risk model and nomograph to predict prognosis of GBM. Based on risk model, patients were divided into high-risk and low-risk group. As patients with high-risk developed poor prognosis, more aggressive methods and closer follow-up time interval are required, indicating the risk model offer precise individualized treatment in clinical practice. Nomograph is characterized by intuitive visual presentation in guiding clinical practice and we established a mitophagy-related nomograph for the first time with excellent predictive performance, which is superior than conventional WHO stage.

To further analyze the functional enrichment between low-risk and high-risk groups, GO and KEGG pathway analysis were conducted. Go analysis suggested that the DEGs mainly involved in synaptic activity and ion channel, which is consistent with previous findings that synaptic activity drives the progression of gliomas and ion channel is closely related to the proliferation, metastasis, invasion of GBM [38, 39]. KEGG pathway analysis indicated that the DEGs between 2 subgroups mainly enriched in neuroactive ligand-receptor interaction as well as calcium signaling. Previous studies demonstrated that neuroactive ligand-receptor interaction plays an important part in development of GBM [40, 41]. A comprehensive analysis conducted by Pal J found that GBM patients with defective neuroactive ligand-receptor interaction had poor prognosis [40]. Calcium signaling was closely related to the tumorigenesis, progression of GBM [42]. These findings suggested that the 4-gene risk model might regulate tumor progression through the pathways or processes above.

As previous studies suggested that immune infiltration played a crucial role in the prognosis of patients with glioma, we further analyzed and demonstrated that the immune filtration and immune pathway were statistically different between 2 subgroups [43]. It is reported by recent studies that GBM patient with numerous CD8+ T cells tend to have better survival, while activated NK cells predict better prognosis in GBM [44]. In our study, the level of CD8+ T cells, NK+ cells, Tfh cells and Th2 T cells were all downregulated in high-risk group, suggesting a suppressive immune infiltration in patients with poor prognosis. In addition, the prognostic gene TOMM20 is positively corelated to CD8+ T cells in TIMER analysis (Fig. 9), suggesting a close link between mitophagy and immune infiltration. The enigmatic and sophisticated association linking immunity with mitophagy are being gradually uncovered with the progress of the experimental research. One study from Paul K revealed that mitophagy could trigger the proliferation of CD8+ cell to improve prognosis in cancer [45]. Another study from Alejandro López-Soto suggested that the function of NK cell can be modulated by mitophagy [46]. Furthermore, the immune pathway analysis mainly focused on APC costimulation, inflammation-promoting, type I IFN response pathway. Finally, CIBERSORT analysis revealed that the primary immune cells enriched in the low-risk group were activated NK cells, suggesting the crucial role of activation of NK cells in mediating better prognosis. Taken all together, our studies indicated that the mitophagy-related genes modulated immune microenvironment and affected the prognosis.

There are some limitations in our study. First, the clinical data from TCGA are incomplete such as therapy and therapeutic effect, which might provide clues on biomarker of treatment. Second, validation of risk signature is best carried out in vivo or well-established study, as external databases often have bias on race or area.

Conclusion

In summary, we utilized 4 mitophagy-related genes to construct a risk model that accurately predicts the prognosis of GBM patients. Our findings suggested the crucial role of mitophagy in GBM, which might be related to tumor immune microenvironment modulation. Further studies are needed to verify these results in vitro and in vivo.

Availability of data and materials

The datasets used and/or analysed during the current study available from the corresponding author on reasonable request.

Abbreviations

GBM:

Glioblastoma multiforme

DEGs:

Differentially expressed genes

LASSO:

Least absolute shrinkage and selection operator

TME:

Tumor microenvironment

PPI:

Protein protein interaction

PCA:

Principal component analysis

AUC:

Area under the curve

KEGG:

Kyoto Encyclopaedia of Genes and Genomes

ssGSEA:

Aingle-sample gene set enrichment analysis

GO:

Gene ontology

STRING:

Search Tool for the Retrieval of Interacting Genes

NK:

Natural killer

Tfh:

Follicular helper T cell

MAP1LC3A:

Microtubule Associated Protein 1 Light Chain 3 Alpha

TOMM20:

Translocase Of Outer Mitochondrial Membrane 20

PHB2:

Prohibitin 2

UBA52:

Ubiquitin A-52 Residue Ribosomal Protein Fusion Product 1

References

  1. Wesseling P, Capper D. WHO 2016 Classification of gliomas. Neuropathol Appl Neurobiol. 2018;44(2):139–50.

    Article  CAS  PubMed  Google Scholar 

  2. Ostrom QT, Cioffi G, Gittleman H, Patil N, Waite K, Kruchko C, et al. CBTRUS Statistical Report: Primary Brain and Other Central Nervous System Tumors Diagnosed in the United States in 2012-2016. Neuro-Oncology. 2019;21(Suppl 5):v1–v100.

    Article  PubMed  PubMed Central  Google Scholar 

  3. Wang Y, Jiang T. Understanding high grade glioma: molecular mechanism, therapy and comprehensive management. Cancer Lett. 2013;331(2):139–46.

    Article  CAS  PubMed  Google Scholar 

  4. Yu L, Chen Y, Tooze SA. Autophagy pathway: Cellular and molecular mechanisms. Autophagy. 2018;14(2):207–15.

    Article  CAS  PubMed  Google Scholar 

  5. Onorati AV, Dyczynski M, Ojha R, Amaravadi RK. Targeting autophagy in cancer. Cancer. 2018;124(16):3307–18.

    Article  PubMed  Google Scholar 

  6. Bernardini JP, Lazarou M, Dewson G. Parkin and mitophagy in cancer. Oncogene. 2017;36(10):1315–27.

    Article  CAS  PubMed  Google Scholar 

  7. Maiti P, Scott J, Sengupta D, Al-Gharaibeh A, Dunbar GL. Curcumin and solid lipid curcumin particles induce autophagy, but inhibit mitophagy and the PI3K-Akt/mTOR pathway in cultured glioblastoma cells. Int J Mol Sci. 2019;20(2):399.

    Article  PubMed Central  CAS  Google Scholar 

  8. Liu C, Peng Z, Li P, Fu H, Feng J, Zhang Y, et al. lncRNA RMST suppressed gbm cell mitophagy through enhancing FUS SUMOylation. Mol Ther Nucleic Acids. 2020;19:1198–208.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Huang T, Xu T, Wang Y, Zhou Y, Yu D, Wang Z, et al. Cannabidiol inhibits human glioma by induction of lethal mitophagy through activating TRPV4. Autophagy. 2021:1–15.

  10. He C, Lu S, Wang XZ, Wang CC, Wang L, Liang SP, et al. FOXO3a protects glioma cells against temozolomide-induced DNA double strand breaks via promotion of BNIP3-mediated mitophagy. Acta Pharmacol Sin. 2021;42(8):1324–37.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Zhang J, Liu L, Xue Y, Ma Y, Liu X, Li Z, et al. Endothelial monocyte-activating polypeptide-II induces BNIP3-mediated mitophagy to enhance temozolomide cytotoxicity of glioma stem cells via down-regulating MiR-24-3p. Front Mol Neurosci. 2018;11:92.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  12. Wang Y, Qian T, You G, Peng X, Chen C, You Y, et al. Localizing seizure-susceptible brain regions associated with low-grade gliomas using voxel-based lesion-symptom mapping. Neuro-Oncology. 2015;17(2):282–8.

    Article  PubMed  Google Scholar 

  13. Liu X, Li Y, Qian Z, Sun Z, Xu K, Wang K, et al. A radiomic signature as a non-invasive predictor of progression-free survival in patients with lower-grade gliomas. Neuroimage Clin. 2018;20:1070–7.

    Article  PubMed  PubMed Central  Google Scholar 

  14. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):e47.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  15. Friedman J, Hastie T, Tibshirani R. Regularization paths for generalized linear models via coordinate descent. J Stat Softw. 2010;33(1):1–22.

    Article  PubMed  PubMed Central  Google Scholar 

  16. Blanche P, Dartigues JF, Jacqmin-Gadda H. Estimating and comparing time-dependent areas under receiver operating characteristic curves for censored event times with competing risks. Stat Med. 2013;32(30):5381–97.

    Article  PubMed  Google Scholar 

  17. Uhlen M, Fagerberg L, Hallstrom BM, Lindskog C, Oksvold P, Mardinoglu A, et al. Proteomics. Tissue-based map of the human proteome. Science. 2015;347(6220):1260419.

    Article  PubMed  CAS  Google Scholar 

  18. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Hanzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013;14:7.

    Article  PubMed  PubMed Central  Google Scholar 

  20. Li T, Fu J, Zeng Z, Cohen D, Li J, Chen Q, et al. TIMER2.0 for analysis of tumor-infiltrating immune cells. Nucleic Acids Res. 2020;48(W1):W509–14.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Ziegler PK, Bollrath J, Pallangyo CK, Matsutani T, Canli O, De Oliveira T, et al. Mitophagy in intestinal epithelial cells triggers adaptive immunity during tumorigenesis. Cell. 2018;174(1):88–101.e116.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Bai H, Inoue J, Kawano T, Inazawa J. A transcriptional variant of the LC3A gene is involved in autophagy and frequently inactivated in human cancers. Oncogene. 2012;31(40):4397–408.

    Article  CAS  PubMed  Google Scholar 

  23. Giatromanolaki A, Sivridis E, Mitrakas A, Kalamida D, Zois CE, Haider S, et al. Autophagy and lysosomal related protein expression patterns in human glioblastoma. Cancer Biol Ther. 2014;15(11):1468–78.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Wang Z, Gao L, Guo X, Feng C, Lian W, Deng K, et al. Development and validation of a nomogram with an autophagy-related gene signature for predicting survival in patients with glioblastoma. Aging (Albany NY). 2019;11(24):12246–69.

    Article  CAS  Google Scholar 

  25. Johnson JM, Lai SY, Cotzia P, Cognetti D, Luginbuhl A, Pribitkin EA, et al. Mitochondrial metabolism as a treatment target in anaplastic thyroid cancer. Semin Oncol. 2015;42(6):915–22.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Mikkilineni L, Whitaker-Menezes D, Domingo-Vidal M, Sprandio J, Avena P, Cotzia P, et al. Hodgkin lymphoma: A complex metabolic ecosystem with glycolytic reprogramming of the tumor microenvironment. Semin Oncol. 2017;44(3):218–25.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Grillon E, Farion R, Reuveni M, Glidle A, Remy C, Coles JA. Spatial profiles of markers of glycolysis, mitochondria, and proton pumps in a rat glioma suggest coordinated programming for proliferation. BMC Res Notes. 2015;8:207.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  28. Park SH, Lee AR, Choi K, Joung S, Yoon JB, Kim S. TOMM20 as a potential therapeutic target of colorectal cancer. BMB Rep. 2019;52(12):712–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Sievers C, Billig G, Gottschalk K, Rudel T. Prohibitins are required for cancer cell proliferation and adhesion. PLoS One. 2010;5(9):e12735.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  30. Zhang H, Yin C, Liu X, Bai X, Wang L, Xu H, et al. Prohibitin 2/PHB2 in Parkin-Mediated mitophagy: a potential therapeutic target for non-small cell lung carcinoma. Med Sci Monit. 2020;26:e923227.

    CAS  PubMed  PubMed Central  Google Scholar 

  31. Cheng J, Gao F, Chen X, Wu J, Xing C, Lv Z, et al. Prohibitin-2 promotes hepatocellular carcinoma malignancy progression in hypoxia based on a label-free quantitative proteomics strategy. Mol Carcinog. 2014;53(10):820–32.

    Article  CAS  PubMed  Google Scholar 

  32. Yang J, Li B, He QY. Significance of prohibitin domain family in tumorigenesis and its implication in cancer diagnosis and treatment. Cell Death Dis. 2018;9(6):580.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  33. Kobayashi M, Oshima S, Maeyashiki C, Nibe Y, Otsubo K, Matsuzawa Y, et al. The ubiquitin hybrid gene UBA52 regulates ubiquitination of ribosome and sustains embryonic development. Sci Rep. 2016;6:36780.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Kanayama H, Tanaka K, Aki M, Kagawa S, Miyaji H, Satoh M, et al. Changes in expressions of proteasome and ubiquitin genes in human renal cancer cells. Cancer Res. 1991;51(24):6677–85.

    CAS  PubMed  Google Scholar 

  35. Barnard GF, Mori M, Staniunas RJ, Begum NA, Bao S, Puder M, et al. Ubiquitin fusion proteins are overexpressed in colon cancer but not in gastric cancer. Biochim Biophys Acta. 1995;1272(3):147–53.

    Article  PubMed  Google Scholar 

  36. Han XJ, Lee MJ, Yu GR, Lee ZW, Bae JY, Bae YC, et al. Altered dynamics of ubiquitin hybrid proteins during tumor cell apoptosis. Cell Death Dis. 2012;3:e255.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  37. Zhou B, Zhang JY, Liu XS, Chen HZ, Ai YL, Cheng K, et al. Tom20 senses iron-activated ROS signaling to promote melanoma cell pyroptosis. Cell Res. 2018;28(12):1171–85.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Venkataramani V, Tanev DI, Strahle C, Studier-Fischer A, Fankhauser L, Kessler T, et al. Glutamatergic synaptic input to glioma cells drives brain tumour progression. Nature. 2019;573(7775):532–8.

    Article  CAS  PubMed  Google Scholar 

  39. Griffin M, Khan R, Basu S, Smith S. Ion channels as therapeutic targets in high grade gliomas. Cancers (Basel). 2020;12(10):3068.

    Article  CAS  Google Scholar 

  40. Pal J, Patil V, Kumar A, Kaur K, Sarkar C, Somasundaram K. Loss-of-function mutations in Calcitonin Receptor (CALCR) identify highly aggressive glioblastoma with poor outcome. Clin Cancer Res. 2018;24(6):1448–58.

    Article  CAS  PubMed  Google Scholar 

  41. Wei B, Wang L, Du C, Hu G, Wang L, Jin Y, et al. Identification of differentially expressed genes regulated by transcription factors in glioblastomas by bioinformatics analysis. Mol Med Rep. 2015;11(4):2548–54.

    Article  CAS  PubMed  Google Scholar 

  42. Leclerc C, Haeich J, Aulestia FJ, Kilhoffer MC, Miller AL, Neant I, et al. Calcium signaling orchestrates glioblastoma development: Facts and conjunctures. Biochim Biophys Acta. 2016;1863(6 Pt B):1447–59.

    Article  CAS  PubMed  Google Scholar 

  43. Gieryng A, Pszczolkowska D, Walentynowicz KA, Rajan WD, Kaminska B. Immune microenvironment of gliomas. Lab Investig. 2017;97(5):498–518.

    Article  CAS  PubMed  Google Scholar 

  44. Alexiou GA, Vartholomatos G, Karamoutsios A, Batistatou A, Kyritsis AP, Voulgaris S. Circulating progenitor cells: a comparison of patients with glioblastoma or meningioma. Acta Neurol Belg. 2013;113(1):7–11.

    Article  PubMed  Google Scholar 

  45. Lobinger D, Gempt J, Sievert W, Barz M, Schmitt S, Nguyen HT, et al. Potential role of Hsp70 and activated NK cells for prediction of prognosis in glioblastoma patients. Front Mol Biosci. 2021;8:669366.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Lopez-Soto A, Bravo-San Pedro JM, Kroemer G, Galluzzi L, Gonzalez S. Involvement of autophagy in NK cell development and function. Autophagy. 2017;13(3):633–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

We thank TCGA for providing dataset. And all methods were carried out in accordance with relevant guidelines and regulations from TCGA.

Funding

This work is supported by the National Natural Science Foundation of China (82100238, 82171698, 82170561, 81300279, 81741067), the Natural Science Foundation for Distinguished Young Scholars of Guangdong Province (2021B1515020003), the Climbing Program of Introduced Talents and High-level Hospital Construction Project of Guangdong Provincial People’s Hospital (DFJH201923, DFJH201803, KJ012019099, KJ012021143, KY012021183).

Author information

Author notes

  1. Jinghua Wang, Xinqi Qiu, Jiayu Huang and Zewei Zhuo contributed equally to this work.

Authors and Affiliations

  1. Department of Hematology, Guangdong Provincial People’s Hospital, Guangdong Academy of Medical Sciences, Guangzhou, Guangdong, China

    Jinghua Wang

  2. Zhuguang Community Healthcare Center, Guangzhou, 510080, China

    Xinqi Qiu

  3. Department of Otorhinolaryngology-Head and Neck Surgery, Huizhou Municipal Central People’s Hospital, Huizhou, 516001, People’s Republic of China

    Jiayu Huang

  4. School of Bioscience and Bioengineering, South China University of Technology, Guangzhou, 510006, China

    Zewei Zhuo & Huihuan Wu

  5. Department of Gastroenterology, Guangdong Provincial People’s Hospital, Guangdong Academy of Medical Sciences, Guangzhou, Guangdong, 510080, People’s Republic of China

    Zewei Zhuo, Hao Chen, Ruijie Zeng, Huihuan Wu, Kehang Guo, Qi Yang & Yujun Luo

  6. Department of General Medicine, The First Affiliated Hospital of Guangzhou Medical University, Guangzhou, 510120, China

    Huiling Ye

  7. Department of Pathology and Laboratory Medicine, University of Cincinnati College of Medicine, Cincinnati, OH, 45267, USA

    Wei Huang

Authors
  1. Jinghua Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Xinqi Qiu
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Jiayu Huang
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Zewei Zhuo
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Hao Chen
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Ruijie Zeng
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Huihuan Wu
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Kehang Guo
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Qi Yang
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Huiling Ye
    View author publications

    You can also search for this author in PubMed Google Scholar

  11. Wei Huang
    View author publications

    You can also search for this author in PubMed Google Scholar

  12. Yujun Luo
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

WJH, CH, ZWZ, QXQ, JYH acquired the data, performed the analysis and wrote the manuscript. ZWZ, RJZ, JYH, YQ, WHH, GKH participated in data analysis. YJL, WH, YHL were involved in study design, supervision and acquiring funding. All authors contributed to editorial changes in the manuscript. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Huiling Ye, Wei Huang or Yujun Luo.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare no conflict of interest.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Wang, J., Qiu, X., Huang, J. et al. Development and validation of a novel mitophagy-related gene prognostic signature for glioblastoma multiforme. BMC Cancer 22, 644 (2022). https://doi.org/10.1186/s12885-022-09707-w

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12885-022-09707-w

Keywords

  • Mitophagy
  • Risk model
  • Glioblastoma multiforme
  • Prognosis
  • Risk signature

BMC Cancer

ISSN: 1471-2407

Contact us