The gene expression data and clinical characteristics of patients with glioma in the training set were obtained from the Chinese Glioma Genome Altas (CGGA)-RNA-seq dataset (693) (http://www.cgga.org.cn/). The Cancer Genome Atlas (TCGA) dataset (https://portal.gdc.cancer.gov) and the CGGA RNA-seq dataset (325) were used as two independent datasets. The included patients and detailed metabolism-related genes are listed in Supplementary Table 1. Genetic mutation profiles were obtained from the cBioPortal dataset (http://www.cbioportal.org/). Notably, 26 LMGs were downloaded from the Molecular Signature Database version 7.0 (MSigDB) (http://www.broad.mit.edu/gsea/msigdb/). The representative immunohistochemical (IHC) staining images of LMGs were obtained from the Human Protein Atlas (HPA) database (https://www.proteinatlas.org/).

To determine the expression of LMGs between tumor tissues and paracancerous tissues of glioma, six samples from Changhai Hospital were collected for quantitative real-time polymerase chain reaction (qRT-PCR). The primer of the six LMGs is shown in Supplementary Table 2. Ethical approval was granted by the Ethics Committee of Shanghai Changhai Hospital (Ethics approval number: CHEC2020-164). All participants were fully informed of the research and provided the informed consent forms.

The STRING functional protein association network (https://string-db.org) (13) was used to construct the interaction network. The input genes consisted of LMGs to a high confidence (0.4) of active interaction. Moreover, functional enrichment analysis was carried out on the LMGs.

Univariate cox regression analysis on the LMGs in the training cohort was performed to identify the association between the expression levels of the genes and overall survival (OS) time using the survival package. Significant genes with a p-value of < 0.05 identified by univariate cox regression were further screened by least absolute shrinkage and selection operator (LASSO) cox regression. Then, six optimal LMGs were used to construct a prognostic risk score model using the following formula: risk score = ∑Coefi·Expi (Coef is the regression coefficient calculated by the LASSO). According to the calculated prognostic risk score, all patients were divided into a high-risk and a low-risk group. Kaplan-Meier survival curve and time-dependent receiver operating characteristic (ROC) curve were used to evaluate the prognostic performance of the constructed risk model.

The abundance of immune cells in tumor tissue was calculated using the Cibersort algorithm (http://CIBERSORT.stanford.edu/), which transformed the normalized gene expression matrix into the composition of infiltrating immune cells. The LM22 signature matrix defined 22 types of immune cell components and was used as a reference expression signature with 1,000 permutations.

Hallmark gene sets and curated gene sets were downloaded from the MsigDB dataset. Gene set variation analysis (GSVA) was utilized to evaluate the relative enrichment of the 50 hallmark pathways across samples using a non-parametric approach (14). The correlation between the LMG score and the 50 cancer Hallmark Pathway score was calculated by Spearman analysis and was visualized in a heatmap. Gene set enrichment analysis (GSEA) for the curated gene sets was performed using the R package ClusterProfiler (15) based on the differentially expressed genes. Significant gene sets were defined by false discovery rate (FDR) < 0.05 and normalized enrichment score (NES) > 2. The results were visualized by R package enrichplot.

The clinical features and LMG signature were used to develop the nomogram for glioma prognostic prediction using the “rms” package (16). Calibration curves and ROC curves for 1, 2, 3, and 5 years were used to assess the accuracy and discrimination ability of the nomogram.

The single nucleotide variation (SNV) and copy number variation (CNV) data of patients with glioma were obtained from the TCGA dataset. In brief, the TCGAbiolinks R package was used to download the SNV and CNV data from the TCGA database. The maftools R package was utilized to analyze SNV data and visualize the mutation waterfall plots. CNV data from the two different risk groups were analyzed using the GISTIC2.0 module on GenePattern (https://cloud.genepattern.org/gp/pages/index.jsf) with default parameters.

Genomics of Drug Sensitivity in Cancer (GDSC; https://www.cancerrxgene.org/) (17) is the largest public pharmacogenomics database, which is used for exploring molecular cancer therapy and mutation. The R package pRRophetic (18) was used to screen chemotherapeutic agents in the high- and low-risk groups.

First, RNA from glioma and paracancerous tissues was extracted using Trizol (Invitrogen, USA) extraction reagents. Isolated RNA was converted into cDNA via a cDNA reverse transcription kit (TaKaRa, Cat: RR047A, Japan). Following this, the relative expression of cDNA was visualized using the SYBR qPCR Master Mix (Takara, Cat: RR820A, Japan) on a Real-Time PCR system (Applied Biosystems, Foster City, USA). β-Actin was used as a reference gene, and relative gene expression was quantified using the formula: 2–^ΔΔCt.

Statistical analyses were mainly performed using R version 3.6.2. Kaplan–Meier and log-rank analyses were performed to evaluate the survival differences between different groups of patients. Student's t-test and one-way ANOVA analysis were used to estimate the differences between two groups and more than two groups. The correlation analysis was calculated using the “Spearman” method. Two-sided p < 0.05 was regarded as statistically significant.

To identify the significance of lactate metabolism in glioma, we collected 26 LMGs from the MSigDB dataset and analyzed their expression in glioma with different clinical characteristics. It was found that the expression levels of LDHB, LDHD, SLC16A7, SLC25A12, PER2, and TP53 in 1p/19q Codel glioma were higher than those in non-Codel glioma. The expression levels of EMB, LDHA, SLC16A1, SLC16A3, SLC16A8, PARK7, and PFKFB2 were lower in 1p/19q Codel glioma than those in non-Codel glioma (Figure 1A). Furthermore, the expressions of HAGH, LDHB, LDHD, SLC16A7, SLC25A12, and TP53 were upregulated, and the expressions of EMB, LDHA, SLC16A3, and SLC16A8 were downregulated in isocitrate dehydrogenase (IDH) wild-type glioma compared with mutant glioma (Figure 1B). To investigate the protein interactions between LMGs and their biological functions, we constructed the interactions between genes related to lactate metabolism and the enriched their gene ontology (GO) functions by these interactions from the STRING database. The results showed that there were interactions between LMGs and that LDHA, LDHB, and SLC16A3 had more interactions relative to other genes (Figure 1C). These genes were enriched in the lactate transmembrane transport, lactate metabolic process, and biosynthetic process (Figure 1C). Correlation analysis showed a strong positive correlation between some of the LMGs. For example, TP53, HIF1A, MRS2, DRAS2, and SLC16A1 were positively correlated, and LDHB, PARK7, HAGH, and PNKD were also positively correlated, suggesting that they may play a synergistic role in lactate metabolism (Figure 1D). We also examined the mutational panorama of LMGs and observed that the frequencies of mutation change with the LMGs were quite low in the TCGA cohort, except for TP53, which may imply that their transcriptional levels were relatively stable in performing biological functions (Figure 1E).

To explore the prognostic significance of LMGs in glioma, we performed univariate cox proportional hazards regression analysis and evaluated the relationship between the gene expression value of LMGs and the OS status of patients with glioma. We found that 11 genes were significantly associated with OS time (p < 0.05, Figure 2A). Among the 11 genes, 6 genes (HIF1A, LDHA, SLC16A1, MRS2, PFKB2, and TP53) were identified as high-risk factors (hazard ratio >1) and 5 genes (HAGH, LDHB, SLE16A1, SLC25A12, and PER2) were identified as protective factors (hazard ratio < 1). These genes were subjected to Kaplan-Meier analysis, and the representative genes are shown in Figure 2B. Moreover, the IHC profiles generated from the HPA datasets also provided information about the location of LMGs and their expression status (Supplementary Figures S1A–L). We used qPCR to further validate the LMG expression. The results showed that glioma tissues have higher expression of LDHA, HLF1A, SLC16A1, and MRS2 and lower expression of LDHB and SLC25A12 compared with paracancerous tissues (Figure 2C). These results suggested that LMGs were correlated with outcomes in patients with glioma and might play a role in glioma progression.

Considering the prognostic value of some LMGs in glioma, we intended to construct a risk score model to evaluate the prognosis status of glioma samples more accurately. The 11 prognosis-related genes were further performed using the LASSO regression analysis to screen the most valuable predictive genes (Figures 3A, B). Five genes were selected with a p < 0.05 and the risk signature included MRS2, LDHA, SLC16A1, LDHB, and SLC25A12 (Figure 3C). Based on the risk score, we divided the patients with glioma into high-risk and low-risk groups (Figure 3D). It was found that the LMG signature could effectively differentiate cancer prognosis. Patients with glioma with high-risk scores had a shorter survival time and poorer prognosis as compared with those with low-risk scores (Figures 3E, F). In the high-risk group, the expression of MRS2, LDHA, and SLC6A1 was higher, and the expression of LDHB and SLC25A12 was lower than that in the low-risk group (Figure 3G). The risk score exhibited a high prognostic validity, with an area under curve (AUC) value of 0.729 (Figure 3H). To further validate the robustness of the LMG model, we utilized two independent sources of glioma to test the robustness of our model. The result of the Kaplan-Meier survival analysis showed that the LMG model was effective in both datasets (Figures 3I, J). The above results suggest that our LMG-based prognostic model could effectively predict glioma prognosis.

We further compared the risk score between patients with different clinical characteristics and found that the LMG signature was associated with age, tumor grade, IDH mutation, radiotherapy, and the O-methylguanine-DNA-methyltransferase (MGMT) status (Figure 4A). However, the LMG signature was not associated with gender (Figure 4A). To further identify the role of the LMG signature for clinical subgroups, patients in the training cohort were divided into subgroups according to the clinical features, including the age, sex, IDH status, 1p/19q status, and WHO classification of the CGGA-693 cohort to compare the survival curves between the subgroup based on the median risk score. The results showed that the LMG signature also had prognostic value in each subgroup, except for patients with grade II in the training cohort (Figures 4B–H).

To test whether the LMG signature was an independent predictor for glioma, multivariate analysis was performed using a cox proportional hazards model with the predictive risk score and clinical factors. The results showed that IDH mutation, 1p19q codel, and the LMG-based risk score were independent prognostic predictors of OS (Figure 5A). To evaluate the clinical application of the LMG signature, we constructed a nomogram model by combining the clinical features (age, sex, and neoadjuvant) with the risk score of patients with glioma (Figure 5B). The nomogram showed a favorable predictive ability for OS rates, with a high AUC value of 0.894 (Figure 5C). The calibration curve showed that the predictions of our nomogram model were in good agreement with the 1-, 2-, 3-, and 5-year actual observations (Figure 5D). Thus, our nomogram was more accurate and clinically valuable than any single independent recurrence factor.

To evaluate the relationship between these TME indicators and the LMG signature, we used the CIBERSORT algorithm (19) to calculate the proportion of the tumor-infiltrating immune cells in the TME in patients with glioma. The proportion of M0 macrophages, monocytes, CD4⁺ T memory resting cells, and Treg cells was high and that of memory B cells, NK cells, and CD8⁺ T cells was low in patients of the high-risk group (Figure 6A), suggesting a lack of the cell killing activity within tumors of patients in the high-risk group, most probably due to the reduced number of NK cells and less CD8⁺T infiltration. Meanwhile, we compared the expression level of immune-checkpoint genes between the low- and high-risk LMG signature groups and found that CD80, CD86, PDL1, PDL2, HAVCR2, and TIGIT were highly expressed in the high-risk groups (Figure 6B). The upregulation of PDL1, PDL2, HAVCR2, and TIGIT could be in an immunosuppressive status, thus facilitating tumor progression.

To evaluate the relationship between the LMG signature and cancer hallmark pathways, we collected the hallmark gene set from the MsigDB dataset and calculated the pathway activity in each sample using GSVA (14). The GSVA scores of these samples were then subjected to Spearman correlation analysis with the risk coefficients of the samples. It was found that the risk score was highly correlated with the “IL6 JAK STAT3 pathway,” “epithelial-mesenchymal response,” and “TNF-α signaling via NFKB” pathway (Figure 6C). We further used ESTIMATE (20) to calculate the immune score of each patient and then analyzed the correlation between the risk score and the immune score. It was found that there was a significant positive correlation between the risk score and the immune score, and the r-value calculated by Pearson was 0.35 (Figure 6D). GSEA was further performed to identify the involvement of pathways regulating tumorigenesis in the high-risk group. We found that multiple classic tumor-related pathways were enriched in the high-risk group, such as “hypoxia,” “MYC targets,” and “cancer poor survival” pathways (Figure 6E). In contrast, metabolism-related pathways were enriched in the low-risk group, including “Oxidative phosphorylation,” “selenoamino acid metabolism,” and “pyruvate metabolism” pathways (Figure 6F). These results indicate that the LMG signature was associated with the immune response and tumorigenesis pathway. Lactic acid metabolism may affect immune cell function to affect glioma progression.

To investigate the underlying genetic alteration of the risk level defined by LMG signature in glioma, we obtained TCGA-LGG/GBM somatic mutation profiles and analyzed the mutation landscape for patients with high- and low-risk (Figures 7A, B). We explored the top 20 mutated genes in two groups, respectively. The gene with the most mutation frequency was IDH1 (80%) in the low-risk group and that in the high-risk group was TP53 (51%). We also calculated the significant differentially mutated genes between the two risk groups. As shown in Figures 7C, D, IHD1, CIC, FUBP1, and NOTCH1 were found with a much higher mutation rate in the low-risk group than in the high-risk group. In contrast, PTEN, EGFR, TTN, and RB1 showed a much lower mutation rate in the low-risk group than in the high-risk group. Additionally, the CNV alteration landscapes of the high- and low-risk groups were further compared. The two groups had distinct CNV alteration profiles (Figures 7E, F). The above results suggested that genetic alteration may affect the expression of LMG, which further affect tumor progression.

To explore potential therapeutic agents for patients with glioma, GDSC analysis was used to predict the chemotherapy response between the two groups (Figures 8A–L). Several chemotherapeutic agents were significantly different between high- and low-risk groups with regard to their sensitivity. EGFR inhibitor Lapatinib and Gefitinib had higher IC50 in the low-risk group than that in the high-risk group (Figures 8A, B), suggesting EGFR inhibitors could be potential pharmacological therapies for patients with high-risk glioma. Other chemotherapeutic agents, such as Vinblastine, Nilotinib, Vorinostat, Axitinib, and Sorafenib also had distinct sensitivity between the two risk groups. The IC50 values of these agents were significantly higher in patients with high-risk glioma than the patients with low-risk glioma (Figures 8C–L), suggesting that patients with high-risk glioma may have tolerance for these drugs and ultimately be able to resist.