Identification of DELs
2615 genes were differentially expressed between glioma and normal brain tissues in the GEO datasets based on the screening criteria |logFC|>1 and adjusted P < 0.05 to investigate LMGs in gliomas (Fig. 1A-C). Finally, ten DELs—LDHB, PDCD1, SLC16A8, UEVLD, PKM, TLR4, PTH, IL6, ICAM3, PDK1—were identified (Fig. 1D, Supplementary Fig. S1). Importantly, except for the high expression of PDK1 in tumor tissues, the other genes were downregulated (Fig. 1D).
Establishment and validation of lactate metabolism-related risk score model
The risk score of each patient was calculated by weighted regression coefficient and the expression of risk genes to explore their significance. We hypothesized that the risk score determined by the risk score model could predict the prognosis of glioma patients (Table 1). The risk score was calculated using the following formula: Risk score = (-3.03E-05 × LDHB expression) + (7.92E-06 × PKM expression) + (0.014253549 × ICAM3 expression) + (0.000428469 × PDK1 expression). Where LDHB is a beneficial gene, while PKM, ICAM3 and PDK1 are risk genes. According to the median risk score, glioma patients were divided into two groups in the TCGA cohort: low-risk group (n = 342) and high-risk group (n = 352). The Kaplan–Meier survival curve suggested that there is a significant difference in prognosis between the high- and low-risk group (P < 0.001) (Fig. 2A). The prediction performance of the model was evaluated by a time-dependent ROC curve. The area-under-the-curve values of 1-, 3-, and 5-year survival rates were 0.872, 0.866 and 0.803, respectively (Fig. 2B), indicating that the lactate metabolism-related risk score model could accurately predict the prognosis of glioma patients. Moreover, the expression pattern of prognostic genes, the distribution of lactate metabolism-related risk score and the survival status of each patient in the high- and low-risk group were estimated (Fig. 2C-E). Notably, the expected results were obtained after validation with the CGGA database (Supplementary Fig. S2). In addition, we further analyzed the association between risk score and clinical characteristics of patients and found no remarkable difference in risk score among different clinical characteristics (Supplementary Fig. S3).
Table 1
A total of 4 genes significantly associated with overall survival according to multivariate Cox regression analysis
LMGs
|
coefficient
|
HR
|
HR.95L
|
HR.95H
|
p-value
|
LDHB
|
-3.03E-05
|
0.999969721
|
0.999954985
|
0.999984457
|
5.64E-05
|
PKM
|
7.92E-06
|
1.000007919
|
1.000003182
|
1.000012656
|
0.001049881
|
ICAM3
|
0.014253549
|
1.014355615
|
1.009410296
|
1.019325162
|
1.09E-08
|
PDK1
|
0.000428469
|
1.00042856
|
1.000194112
|
1.000663064
|
0.000339587
|
LMGs: Lactate metabolism-related Genes; HR: hazard ratio;
Construction of predictive nomogram
Univariate and multivariate COX regression analyses suggested that lactate metabolism-related risk score was an independent prognostic factor (Fig. 3A-B). To provide a clinical predictive model, we generated a prognostic nomogram to predict the 1-, 3-and 5-year overall survival (OS) rates of glioma patients. Figure 3C shows the survival rate for one of the samples. The correction maps of TCGA and CGGA databases show that the observed consequences were as predicted (Fig. 3D-E), which indicates that the nomogram has a satisfactory prediction performance.
Relationship between expression of LMGs and clinical characteristics
Survival analysis of glioma patients revealed that the mRNA expression levels of ICAM3, PKM, LDHB and PDK1 had a noticeable effect on OS. As depicted in Fig. 4A, in TCGA and CGGA databases, patients with low expression levels of ICAM3, PKM and PDK1 had longer survival times, whereas patients with low expression levels of LDHB were associated with poor OS. IHC staining showed the distribution of LMGs in glioma tissues (Fig. 4B). In light of the low expression level of ICAM3 and PKM in glioma tissues, we then examined whether there are differences between these two genes in different grades of gliomas. The expression levels of ICAM3 and PKM in lower-grade gliomas were substantially lower than those of GBM (Fig. 4B and Supplementary Fig. S4), consistent with the survival analysis results.
Immunocyte infiltration and mutational landscape in the glioma microenvironment
Gliomas promote the infiltration of immune cells (e.g. microglia, peripheral macrophages, myeloid suppressor cells, white blood cells, CD4+T cells and Tregs) into tumors by regulating cell surface molecules and secreting various factors such as cytokines, chemokines and growth factors (Ogden et al., 2008).The immunosuppressive microenvironment established by immune cells infiltrating the tumor area can promote tumor progression and inhibit the immune response, which affects the efficacy of glioma immunotherapy (Audia et al., 2017; Wang et al., 2017a). This is due to the fact that during the interaction of cytokines, chemokines and growth factors secreted by glioma tissues with extracellular matrix components, lactate can serve as signaling molecule involved in intercellular signal transduction, reprogramming the infiltrating immune cells, inducing immune system into an inflammatory or anti-inflammatory state, and contributing to the transformation of tumor microenvironment into an immunosuppressive phenotype, which can promote tumor immune escape(Grabowski et al., 2021; Kes et al., 2020). Therefore, we then explored the infiltration level of immune cells in the glioma microenvironment. According to the degree of immune cell infiltration, glioma samples were divided into low immunity and high immunity cell infiltration. The ESTIMATE algorithm was used to evaluate the glioma microenvironment. We found that there is a drastic difference in the degree of immune infiltration between high and low immune cell infiltration groups (Fig. 5A). The immune, stromal and ESTIMATE scores of the high immunity group were higher than those of the low immunity group (Fig. 5C-E). However, the tumor purity of the high immune cell infiltration group was lower than that of the low immune cell infiltration group (Fig. 5B). The relative infiltration abundance of 22 types of immune cells with different genes in the high- and low-risk groups was calculated by the CIBERSORT algorithm. The results showed that activated CD4+ T memory cells were remarkably infiltrated in the high- risk group, indicating that glioma patients in this group tended to CD4+ T-cell immune response (Fig. 5F). However, successful anti-tumor immune response also requires CD4+T cells to have the ability of self-renewal and effect in the tumor. Following the activation of cytotoxic T cells, the level of extracellular lactate in the tumor microenvironment are increased, leading to accumulation of endogenous lactate in T cells and decreased secretion of pro-inflammatory cytokines, which reduces the anti-tumor immune response of cytotoxic T cells(Mendler et al., 2012). In addition, previous studies have shown that anti-inflammatory cytokines such as TGF-β and IL-10 are the main reasons for the transition from tumor microenvironment to immunosuppressive phenotype(Christofides et al., 2015; Iwami et al., 2011). On this basis, we found that LMGs induce an inflammatory response and are co-expressed with anti-inflammatory cytokines (Supplementary Fig. S5A-B). These results suggest that LMGs are involved in the reprogramming of glioma metabolic pathways, leading to the transformation of glioma microenvironment into immunosuppressive phenotype. Consequently, we explored the mutant genes that caused differences between high- and low-risk groups in the lactate metabolism-related risk model. The 20 genes with the highest mutation frequency in the high- and low-risk groups are shown in Supplementary Fig. S6A-B. By comparing the mutated genes and their mutation frequency between the two groups, we found that TP53 mutations were more prevalent in the high-risk group than in the low-risk group. The mutation rate of the IDH1 gene in the low-risk group was noticeably higher than that in the high-risk group. Moreover, PTEN mutation occurs only in the high-risk group and its mutation frequency was high (21%). We then compared the mutations of DELs, and the results suggested that there were almost no mutations in PKM, ICAM3, PDK1, and LDHB (Supplementary Fig. S6C).
Study on the mechanism of LMGs
To explore the potential upstream and downstream regulatory mechanisms of LMGs in gliomas, we used the STRING database to predict protein-protein interactions (PPI). The PPI network found that PDK1, PKM and LDHB are direct targets of the oncoprotein HIF-1 α (Fig. 6A). Thus, there is a crosstalk among PDK1, PKM, and LDHB. In addition, enrichment pathways were generated by GSVA of the corresponding target genes. LMGs were remarkably correlated with the pathway that promotes tumorigenesis, including NF-κB, PI3K-Akt, mTOR, Wnt, and extracellular matrix- receptor interaction signaling pathways (Fig. 6B).
LMGs may be candidate therapeutic targets in glioma patients
We further investigated the effect of risk score on immune molecules with negative regulation of the anti-tumor immune response to determine whether LMGs have immunosuppressive effects. The genes that negatively regulate anti-tumor immunity and immunosuppressive cytokines were drastically upregulated in the high-risk group (Fig. 6C), indicating that the anti-tumor immune response activity was lower in the high-risk group. We also compared the expression of immunotherapy drug targets in the high- and low-risk groups. As shown in Fig. 6D-F, the expression levels of immunotherapy targets CD274 (PD-L1), CTLA4, HAVCR2, IDO1, LAG3, PDCD1 and PDCD1LG2 were significantly higher in the high-risk group than in the low-risk group, indicating that our risk model may serve as an effective tool to predict the efficacy of therapy in glioma patients. Therefore, targeting LMGs is a viable alternative option for the treatment of glioma patients.