A Coordinated Mitochondrial Genome Program Predicts Improved Outcome in Lower-Grade Glioma

1. Introduction

Mitochondria perform numerous functions beyond ATP production, including regulation of cellular signaling, differentiation, and cell death through mechanisms involving reactive oxygen species (ROS), modulation of calcium homeostasis, release of mitochondrial proteins that activate apoptosis, and retrograde communication with the nucleus [1,2,3]. Mitochondria contain a semi-autonomous, self-replicating genome composed of mitochondrial DNA (mtDNA). Human mtDNA is organized as a compact circular genome of 16,569 base pairs (bp) that contains 13 essential protein-coding genes collectively encoding core components of the oxidative phosphorylation (OXPHOS) system, the principal machinery responsible for aerobic ATP production. In genomic order, these mitochondrial protein-coding genes (mtPCGs) are MT-ND1, MT-ND2, MT-CO1, MT-CO2, MT-ATP8, MT-ATP6, MT-CO3, MT-ND3, MT-ND4L, MT-ND4, MT-ND5, MT-ND6, and MT-CYB. Their protein products are incorporated into mitochondrial respiratory chain complexes I, III, IV, and V, where they participate in electron transport, proton gradient generation, and ATP synthesis [4,5,6]. (Figure 1).

Altered mitochondrial function and dysregulated OXPHOS activity have been implicated in multiple aspects of cancer biology. Different tumor types exhibit either increased or reduced mitochondrial content and OXPHOS activity [7,8]. In some cancer cells, activation of the mitochondrial OXPHOS system occurs in parallel with enhanced aerobic glycolysis, representing an adaptive strategy to sustain ATP production and support proliferation under the hypoxic conditions characteristic of tumors. In addition, mtDNA copy number has been associated with tumor stage in several solid malignancies, including hepatocellular carcinoma, pancreatic adenocarcinoma, and melanoma. Abnormal mitochondrial transcription resulting in dysregulated expression of mtDNA-encoded genes has also been implicated in tumor metabolic reprogramming. Furthermore, mtDNA mutations and inhibition of mitochondrial gene transcription may alter intracellular signaling pathways such as mitogen-activated protein kinase (MAPK) and phosphoinositide 3-kinase (PI3K)/AKT/mammalian target of rapamycin (mTOR) signaling [8,9].

Lower-grade gliomas (LGGs) comprise a heterogeneous group of tumors with marked variability in biological and clinical characteristics, as well as patient outcomes, underscoring the need for novel prognostic markers to improve patient stratification and guide clinical trials and therapeutic decisions [10]. Although LGGs generally exhibit a more favorable prognosis than primary glioblastomas (GBMs), these tumors can undergo progressive malignant transformation and still represent a substantial burden in terms of neurological morbidity and mortality [11,12,13]. Mitochondrial DNA represents a promising molecular platform for the development of diagnostic and prognostic biomarkers in cancer due to features such as its short length and relatively simple structure, which facilitate genomic analysis [8]. Therefore, characterizing the expression patterns and potential clinical associations of mtDNA-encoded genes in LGGs may help identify novel prognostic markers and therapeutic opportunities.

Here, we report the expression profiles, associations with patient overall survival (OS), and expression correlation patterns of the 13 mtPCGs in LGG tumors. Our findings demonstrate that increased and coordinated expression of mtPCGs is strongly associated with improved clinical outcomes in patients with LGG.

2. Results

2.1. Expression of mtPCGs Across Different LGG Types

We first investigated mtPCG expression patterns in LGG tumors according to either histological classification or molecular subtype. When The Cancer Genome Atlas (TCGA) and Chinese Gioma Genome Atlas (CGGA) tumors were stratified by histological type, no substantial overall differences in expression profiles were observed among groups, although several individual comparisons reached statistical significance (Supplementary Figure S1, S2). Analyses based on molecular classification revealed a clearer expression pattern, characterized by the highest mtPCG expression levels in LGG harboring IDH mutation with 1p/19q co-deletion (LGG-IDH-mut-codel), the lowest levels in IDH wild-type tumors (LGG-IDH-wt), and intermediate expression in IDH-mutant tumors lacking 1p/19q co-deletion (LGG-IDH-mut-non-codel), in both the TCGA (Figure 2) and CGGA (Figure 3) datasets.

2.2. High mtPCG Expression Is Associated with Longer OS in Patients with LGG

To verify whether expression of mtDNA-encoded genes is related to patient prognosis, we performed Kaplan-Meier analyses of OS comparing patients with LGG tumors divided into low- or high-expressing groups. We found that high expression of each of the 13 mtPCGs was significantly associated with longer OS in both the TCGA (Supplementary Figure S3) and CGGA (Supplementary Figure S4) datasets. A significant association between gene expression and patient survival was also observed when transcription levels of all 13 mtPCGs were combined into a “mitochondrial gene score” (Figure 4).

To determine whether the mitochondrial transcriptional signature exerted prognostic effects independent of established molecular determinants of LGG outcome, multivariate Cox regression analyses were performed including age, tumor grade, IDH status, 1p/19q codeletion status, and the mitochondrial gene score. Although the mitochondrial score showed a strong univariate association with OS, confidence intervals for the mitochondrial gene score crossed the null hazard ratio, therefore this association did not remain statistically significant after adjustment for established molecular prognostic factors (Supplementary Figure S5). These findings suggest that mitochondrial genome expression is closely linked to favorable molecular glioma subtypes and differentiation states.

2.3. Correlations Among Expression of mtPCGs

We then evaluated correlations among expression of the different mtPCGs. Moderate to strong positive correlations were found among most genes in two different analyses, suggesting a coordinated gene expression program (Table 1; Supplementary Figure S6; Figure 5).

3. Discussion

Large-scale cancer genomics initiatives such as TCGA and CGGA have generated extensive transcriptomic datasets from large cohorts of glioma patients, enabling comprehensive analyses of mitochondrial gene expression patterns in human tumors [9,14,15]. Using this strategy, we identified a consistent transcriptional landscape in which increased expression of the 13 mtPCGs, both individually and as a combined gene set, was strongly associated with improved survival in patients with LGG (our results are summarized schematically in Figure 6). The lack of independent prognostic significance in multivariate models does not necessarily diminish the biological relevance of the mitochondrial transcriptional signature. Instead, these findings suggest that coordinated upregulation of mtPCGs may constitute part of a broader transcriptional and metabolic program associated with favorable glioma molecular subtypes, neuronal differentiation, and reduced tumor aggressiveness.

Early reports of elevated mRNA expression of specific mtDNA-encoded genes in cancer cells [7,16,17,18], were followed by a growing body of work investigating the role of mitochondrial respiratory chain components and mitochondria-associated genes in tumor progression and malignancy [3,7,8]. Recently, an analysis of 515 lung adenocarcinoma samples from the TCGA and Gene Expression Omnibus (GEO) databases demonstrated that patients could be stratified according to mitochondrial-related gene expression profiles, with the subgroup exhibiting enhanced mitochondrial metabolic activity showing a poorer clinical prognosis [19].

Glioma cells show an oxidative metabolic phenotype where mitochondrial respiration is crucial for energy production, and targeting the OXPHOS system has been explored as a potential therapeutic strategy [20]. In GBM cells, depletion of the mitochondria-associated methyltransferase METTL17, which stimulates OXPHOS through mitochondrial RNA methylation, impairs cell proliferation, migration, and invasion, and METTL17 overexpression reversed these effects [21]. In addition, high expression of genes coding glycolytic enzymes, HK2 and PKM2, and low expression of genes involved in mitochondrial oxidative metabolism, SDHB and COX5A, were associated with poorer GBM patient survival [22]. LGG tumors can display dysfunctional mitochondrial functions and increased number of mutations in mitochondrial genome coding genes, possibly related to increased oxidative stress. Variations in mtDNA sequences in these cells may lead to protein malfunction and OXPHOS impairment [23]. In TCGA LGG tumors, Mou et al. [10] carried out an analysis of 200 OXPHOS-related genes and identified two molecular subtypes through consensus clustering, with the subtype termed C2 showing poorer prognosis. In light of the previous literature, our findings are particularly noteworthy, as they suggest a consistent pattern in which enhanced mitochondrial gene expression, specifically mtPCGs, is associated with lower clinical risk across the glioma subtypes investigated.

We have recently begun exploring the hypothesis that increased expression of genes associated with excitatory synaptic signaling and synaptic plasticity may contribute to a more favorable clinical behavior by conferring a more mature, neuronal-like phenotype to LGG tumors. Higher expression of DLG2, DLG3, and DLG4, which encode synaptic scaffolding proteins of the membrane-associated guanylate kinase (MAGUK) family, PSD-93 (SAP-93), SAP-102, and PSD-95, respectively, was strongly correlated with the expression of a synaptic gene signature and was associated with significantly longer survival in both the TCGA and CGGA LGG cohorts [24]. The mitochondrial OXPHOS system is increasingly recognized as an integral component of the synaptic and neuronal plasticity machinery [25,26,27,28,29]. Synaptic transmission and activity-dependent plasticity are highly dependent on mitochondrial function, and pharmacological or metabolic inhibition of mitochondrial activity has been shown to impair synaptic signaling and long-term potentiation (LTP) [30,31,32,33]. Notably, neuronal activity associated with hippocampal LTP induction promotes increased expression of mitochondrial genes [34]. Within this framework, the association between elevated expression of the mitochondrial protein-coding genome and improved clinical outcome may reflect a broader biological program in which coordinated upregulation of genes involved in excitatory synaptic transmission, neuronal differentiation, and synaptic plasticity contributes to a less aggressive phenotype in LGGs.

An important limitation of the present study is that transcriptomic associations with patient survival, as well as correlations among different gene expression signatures, should not be interpreted as evidence of direct biological causation. Rather, these findings are best viewed as a framework for the generation of mechanistic hypotheses and the identification of potential prognostic transcriptional state. Functional studies will be necessary to determine whether mtPCGs actively contribute to the modulation of LGG aggressiveness. Experimental approaches involving gene silencing or genetic manipulation in LGG models may help clarify whether mtPCGs exert direct tumor-suppressive effects.

Another relevant consideration is that mRNA abundance does not necessarily correlate with protein expression or functional protein activity. Although several studies across different biological systems have reported substantial concordance between transcriptomic and proteomic measurements, others have demonstrated only modest correlations between mRNA and protein levels [35,36,37]. Consequently, validation at the protein and functional levels will be important to further establish the biological significance of the transcriptional patterns identified in this study.

4. Materials and Methods

4.1. Datasets and Gene Expression Analyses

RNA-sequencing data together with corresponding clinical annotations for patients with lower-grade glioma (LGG) were obtained from two independent cohorts comprising World Health Organization (WHO) grade 2 and grade 3 gliomas: The TCGA https://www.cancer.gov/tcga) and Chinese Glioma Genome Atlas (CGGA, http://www.cgga.org.cn) [14,15]. For the TCGA cohort, gene-level expected count data were retrieved from the TOIL repository through the UCSCXenaTools package (version 1.7.0; rOpenSci project, Auckland, New Zealand). Clinical information was obtained from the TCGA Pan-Cancer Atlas through the cBioPortalData package. For the CGGA cohort, RNA-sequencing and clinical datasets from the CGGA-325 and CGGA-693 projects were downloaded and integrated by retaining genes common to both datasets. To minimize potential technical variability between datasets, batch origin (CGGA-325 versus CGGA-693) was included as a covariate in the DESeq2 design formula during differential expression analyses. Samples lacking survival information or incomplete clinical annotations were excluded from downstream analyses.

TCGA expected counts were converted back to integer-compatible count values to enable count-based modeling with DESeq2. For CGGA, raw count matrices from both datasets were merged after harmonization of shared genes. When duplicate gene symbols were identified, the transcript displaying the highest average expression across samples was retained. Variance-stabilizing transformation was then applied to normalized counts using DESeq2, followed by gene-wise z-score standardization within each cohort.

Expression patterns of mtPCGs were initially examined according to the previous histopathological classification of gliomas [14,38], including astrocytoma (TCGA, n = 188; CGGA, n = 239), oligodendroglioma (TCGA, n = 179, CGGA, n = 145), and oligoastrocytoma (TCGA, n = 123, CGGA, n = 23). Analyses were subsequently performed using the updated molecular classification system [39] in both TCGA and CGGA cohorts, consisting of LGG harboring IDH mutation with 1p/19q co-deletion (LGG-IDH-mut-codel; TCGA, n = 164; CGGA, n = 113), LGG with IDH mutation without 1p/19q co-deletion (LGG-IDH-mut-non-codel; TCGA, n = 235; CGGA, n = 142), and IDH wild-type LGG (LGG-IDH-wt; TCGA, n = 91; CGGA, n = 88). Groupwise comparisons were performed using the Wilcoxon rank-sum test, and resulting p-values were corrected for multiple comparisons using the Benjamini–Hochberg procedure.

4.2. Survival Analysis

Associations between (OS) and expression levels of selected genes were investigated using Kaplan–Meier survival analyses. Patients were stratified into high- and low-expression groups according to optimal expression thresholds identified with the surv_cutpoint function from the survminer package (Supplementary Table S1). Differences between survival curves were assessed using the log-rank test, and p-values were adjusted for multiple testing across genes and cohorts using the Benjamini–Hochberg correction.

All statistical analyses were conducted in the R environment (version 4.5.1). Variance-stabilizing transformation and batch correction procedures were performed using DESeq2 (version 1.48.1). Survival analyses were carried out with the survival package (version 3.8.3), and Kaplan–Meier plots were generated using survminer. Data visualization was performed with ggplot2 (version 4.0.2). Data retrieval and preprocessing were conducted using UCSCXenaTools (version 1.7.0) and cBioPortalData (version 2.20.0).

4.2.1. Multivariate Survival Analysis

Multivariate Cox proportional hazards regression analyses were performed to evaluate whether the mitochondrial gene score was independently associated with overall survival after adjustment for established clinical and molecular prognostic variables. Covariates included age, tumor grade, IDH mutation status, and 1p/19q co-deletion status. Hazard ratios (HRs) and 95% confidence intervals (CIs) were calculated using the survival package in R.

4.3. Gene Expression Correlation Analyses

Co-expression relationships among mtPCGs were examined using both the Evergene platform [40] and analyses conducted within the R statistical environment. Correlation coefficients (r) ≥ 0.5 combined with p-values < 0.05 were considered indicative of moderate-to-strong correlations. For correlation analyses performed in R, gene expression values were standardized across samples through z-score normalization. Pairwise associations between genes were calculated using Spearman’s rank correlation coefficient in order to evaluate monotonic relationships without assuming normal data distribution. A symmetric correlation matrix was subsequently generated to represent gene co-expression patterns. Hierarchical clustering was performed using Ward’s minimum variance method, and heatmaps were produced with the pheatmap package in R. Correlation strength was represented by a color gradient ranging from blue (negative correlations) to red (positive correlations). Missing values were handled through pairwise deletion, and no data imputation procedures were applied.

5. Conclusions

Our findings indicate that increased expression of mtPCGs is consistently associated with improved survival in LGG. This association was observed both for individual genes and for the combined mitochondrial transcriptional signature across two independent large-scale transcriptomic cohorts, and may be, associated with glioma molecular subtypes with better prognosis. These results contrast with observations reported in other malignancies, in which enhanced mitochondrial metabolic activity has been linked to tumor aggressiveness and poor prognosis. These findings highlight the prognostic potential of mitochondrial transcriptional states in glioma biology and suggest that mtDNA-encoded genes may serve as informative biomarkers of tumor differentiation and clinical behavior. Further mechanistic studies will be necessary to determine whether enhanced mitochondrial gene expression directly contributes to reduced glioma aggressiveness or instead reflects underlying lineage-specific and metabolic states associated with improved patient outcome.