Explore the Heterogeneity of Glioblastoma Based on Genes Related to Glutamine Metabolism

Abstract

Objectives:

To establish a glutamine metabolism (GM)-based classification for glioblastoma (GBM) and evaluate its prognostic and immunotherapeutic implications.

Methods:

A total of 237 GBM patients from the Chinese Glioma Genome Atlas (CGGA) database were included as the training set, and 219 patients from the Gene Expression Omnibus (GEO) database served as the validation set. Consensus clustering was performed based on the expression profiles of 13 GM-associated genes to identify robust subgroups. Differences between clusters were analyzed using clinical indices, genomic and transcriptomic biomarkers. Tumor response to immune checkpoint inhibitors (ICIs) was predicted using the tumor immune dysfunction and exclusion (TIDE) algorithm, tumor microenvironment (TME) score, T cell inflammation score, and SubMap algorithm. A GM-based classifier was subsequently developed and validated.

Results:

Consensus clustering of the training set revealed 2 distinct subgroups (cluster 1 and cluster 2) with significant prognostic differences; cluster 2 exhibited poorer overall survival. Immunotherapy response prediction indicated that cluster 2 had a lower likelihood of benefiting from ICIs. The newly developed GM-based classifier demonstrated high accuracy (AUC > 0.9) and maintained strong consistency with the original clustering in terms of subtype classification and immunotherapy prediction across both datasets.

Conclusion:

This study establishes a robust classification system for GBM based on glutamine metabolism-related genes, which effectively stratifies patients into prognostic subgroups and predicts immunotherapy response. The GM-based classifier offers a valuable tool for guiding clinical prognosis and treatment decisions in GBM.

Introduction

GBM is the most common primary malignant brain tumor among adults, ¹ accounting for 27% of all the brain tumors and 80% of malignant brain tumors, ² with a median survival of 14 to 16 months. ³ Despite treatment with surgical resection followed by radiotherapy and adjuvant chemotherapy, the prognosis of GBM remains dismal, with only a 26.5% survival rate at 2 years.^{1 -3} Based on the 2016 report of the World Health Organization (WHO) classification of tumors of the Central Nervous System (CNS), the diagnosis and stratification of the diffuse gliomas were facilitated by recognizing IDH1/2 mutations and 1p/19q codeletion as primary biomarkers. ⁴ However, mutations in IDH1/2 related to carbohydrate metabolism occur in most low-grade gliomas and secondary GBM, and are less frequent in primary GBM.^5,6

Glutamine, a non-essential amino acid that plays a vital role in the metabolism of proliferating cells, including neoplastic cells, ⁶ has been found the most taken up by cancer cells instead of glucose. ⁷ Several studies on cancer cell metabolism have provided evidences that tumor-specific activation of signaling pathways could regulate glutamine uptake, and the metabolism through glutaminolysis could give the cancer cell an alternative energy source.^4,8 In the CNS, GM is particularly relevant because the glutamine-glutamate cycle regulates the production of glutamate-derived neurotransmitters by tightly controlling the bioavailability of the amino acids in a neuron-astrocyte metabolic symbiosis-dependent manner. ⁶ Previous study has demonstrated that inhibiting glioma glutamate release or blockade of glutamate receptors could suppress tumor progression. ⁹ These findings underscore the need for a deeper understanding of glutamine metabolism (GM) in GBM, ¹⁰ which may uncover novel therapeutic targets and facilitate personalized treatment strategies.

Although IDH-based molecular classification has greatly improved the diagnostic and prognostic stratification of GBM, ⁷ the heterogeneity of glutamine metabolism within this tumor type has not been systematically characterized. ¹¹ Specifically, a comprehensive classification of GBM based on GM-related genes, or a robust GM-associated prognostic classifier, remains lacking. In this study, we explored the feasibility of classifying GBM according to GM-related gene expression profiles. We demonstrate that such classification is not only feasible but also significantly associated with patient prognosis and therapeutic response. Moreover, we developed and validated a robust GM-related gene classifier that effectively stratifies GBM patients and holds promise for guiding personalized treatment decisions to improve clinical outcomes.

Methods

The reporting of this study conforms to the REMARK statement ¹² (Supplemental File).

Data Collection

We retrospectively collected the gene expression profile and corresponding clinical information (gender, age, and survival information) of 237 patients from CGGA_693 of the CGGA database as the training set. Similarly, as validation sets, the gene expression profiles and clinical information of 219 GBM patients from GSE16011 and GSE74187 of the GEO database were obtained. In addition, 13 GM-related genes were selected based on the classic glutamine metabolism pathway from Molecular Signature Database v7.1 (MSigDB; http://www.broad.mit.edu/gsea/msigdb/). The detailed metabolism-related genes are listed in Table S1.

Consensus Clustering

The unsupervised consensus clustering was taken separately in the 3 databases using “ConsensusClusterPlus” R package ¹³ based on the expression patterns of 13 GM-related genes. The Euclidean distance was applied to calculate the similarity distance between samples, and K-means methods were utilized for clustering. Through resampling analysis, 80% of the samples were sampled 100 times. ³ The cumulative distribution function (CDF) determined the optimal number of clusters and was validated in the GSE16011 and GSE74187 datasets. Principal component analysis (PCA) was also carried out to understand better the variations between clusters using the R package “princomp.”

Differentially Expressed Gene Analysis (DEGs) and Functional Enrichment Analysis

DEGs between the primary GBM subtypes were identified through fitting the gene expression profile of GSE16011 database to a linear model using the limma package in R, with an absolute Log 2 (fold change) ⩾0.585^14,15 and P ⩽ .05. Subsequently, univariate Cox regression were performed separately on the up and down-regulated DEGs to select the optimal prognostic gene set. Gene set enrichment analysis (GESA) was conducted for the gene expression profile of GBM samples from CGGA/GSE16011/GSE74187 datasets using the Kyoto Encyclopedia of Genes and Genomes (KEGG) and NetworkAnalyst 3.0. Subsequently, ssGSEA (“GSVA” function in R) was performed to calculate the enrichment scores of each classical oncogenic pathway in each case using transcriptomics.^{16 -18}

Immune Infiltration Estimation, Immune Checkpoint Analysis, and Immunotherapy Assessment

Based on the clustering and expression profile of each gene, CIBERSORT was applied to estimate immune cell infiltration using the default signature gene matrix LM22, which contains 547 genes distinguishing 22 human immune cell subtypes(http://cibersortx.stanford.edu/).^{19 -21} The differences in expression values for each immune checkpoint gene between the 2 clusters in GSE16011 were calculated using the Wilcoxon test and visualized. According to the CGP2014 database, the IC50 for drug analysis of each subtype was estimated through the pRRophetic R package. TIDE (http://tide.dfci.harvard.edu/), a computational method based on the transcriptome to model the induction of T cell dysfunction and the prevention of T cell infiltration in tumor immune evasion, calculate T cell dysfunction and exclusion scores and merged as the TIDE score to predict tumor response to ICIs. ²² In addition, the TME score and T cell inflammation score were obtained based on the T cell inflammatory marker genes and their corresponding coefficient. The submap algorithm (https://cloud.genepattern.org/gp/pages/index.jsf) were also utilized to evaluate the response due to the immunotherapy treatment.

Mutation Spectrum Analysis

The gene mutation data was downloaded from the CGGA database. Mutation situation of each gene in different GBM samples was tested using Maftools package. ²³ Genes are sorted by mutational frequency, differentiating the mutational spectrum between the subtypes of GBM. The top 30 mutated genes with the highest mutation frequencies were selected to display.

Survival Prognostic Analysis

Kaplan-Meier (KM) survival curve is defined as the probability of surviving in a given length of time while considering time in many small intervals. ²⁴ The KM univariate survival analysis was undergone for the differences in overall survival (OS) between the 2 clusters in the training and validation databases using the survival package of the R software. A log-rank test P < .05 was used to define differences in survival time.

The Classifier Construction

The top 30 upregulated DEGs with Log2 (fold change) ⩾0.585 and significant univariate cox regression in the GSE16011 database were selected as seed genes to construct the least absolute shrinkage and selection operator (LASSO) Cox regression model through the R package “glmnet” (binomial). ²⁵ The accurate value of the classifier was tested through the receiver operating characteristic (ROC) curves with R package “pROC.” Risk scores were evaluated by multiplying gene expression using the regression coefficient acquired upon LASSO regression.^1,20 Patients were categorized into high and low-risk score groups (predicted cluster 2 and cluster 1) based on the median risk score.

Statistical Analysis

All statistical analyses were conducted using R software (version 4.0.3, http://www.r-project.org) Student’s t-test was applied to identify the DEGs and IC50 drug analysis between the 2 clusters. The Wilcoxon test were used to access differences in the immune cell infiltration, immune checkpoint analysis, and immunotherapy response across subgroups. The chi-square test was applied to compare clinical indices, such as age and gender, between the 2 clusters. P < .05 was regarded as statistically significant.

Results

Molecular Cluster Identification and Validation

All GBM samples in the CGGA cohort were categorized into k (k = 2, 3, 4, 5) different subtypes based on the expression profiles of 13 GM-related genes. Based on the CDF curves of the consensus score, k = 2 was selected as the optimal division (Figures 1A-C and S1). PCA revealed that the 2 clusters could be separated from each other (Figure 1D). These results indicated significant differences in metabolic phenotypes between the 2 clusters. To validate the stability of the molecular subtypes, we further performed clustering analysis using the GSE16011 and GSE74187 datasets. The clustering results obtained from GSE16011 and GSE74187 were largely consistent with those from the CGGA cohort, as summarized in Figures S2 and S3 respectively. Thus, we identified 2 stable clusters of GBM depending on the expression of 13 GM-related genes.

Figure 1.

GM-related genes could distinguish GBM patients in CGGA database: (A) consensus clustering matrix heatmap plots of 237 samples from CGGA dataset for k = 2, (B) consensus clustering CDF for k = 2 to k = 5, (C) relative change in area under CDF curve for k = 2 to k = 5, and (D) PCA analysis of the GM-related genes expression when k = 2.

Characteristic Differences of the Clustering

The GBM samples from CGGA/GSE16011/GSE74187 databases were divided into 2 groups (cluster 1 and cluster 2). The cluster 2 showed higher-level GM compared to cluster 1, with a poorer prognosis. The differences between the 2 clusters in each database were analyzed based on clinical indices (including gender, age, and OS), functional gene enrichment, immune cell infiltration, and gene mutation frequency.

Clinical Characteristics

KM survival analysis revealed that the OS of the 2 clusters in the GSE16011 and GSE74187 databases showed a significant difference (P < .05; Figure 2A and B), indicating that the GM-related genes could participate in the malignant progression of GBM and had vital implications on patient prognosis. However, no survival discrepancy was found in the CGGA cohort, which may be related to its smaller sample size or shorter follow-up duration. Importantly, the prognostic trend was consistent across all 3 datasets, with cluster 2 showing a worse outcome and a hazard ratio > 1 in each validation set (Figure S4A). Besides, males had a higher likelihood of developing GBM than females in both clusters, but the difference was not significant, as was the case for age (Figure S4B-E).

Figure 2.

Validation of inter-cluster differences after clustering: (A and B) Kaplan-Meier analysis of patients between the 2 clusters in GSE16011 and GSE74187 datasets, (C) functional enrichments between the 2 clusters of CGGA by GSEA analysis, (D) heatmap depicting the normalized enrichment scores of 10 oncogenic pathways between clusters in CGGA, (E) quantification of oncogenic pathways enriched in 2 clusters, (F) boxplots of the CCP and EMT score for 2 clusters in CGGA dataset, (G) a comparison of the abundance of tumor-infiltrating immune cells between the 2 clusters was shown, (H) heatmap showing the association of immune cell infiltration and clinical indicators based on the immune cells. Cluster 1(I) and cluster 2 (J) from the CGGA dataset had some similar mutated genes. The mutation frequency of TP53 were the highest in the 2 clusters.

Different Biological Processes at the Transcription Levels

GSEA was peformed on the gene expression data of GBM patients to explore the differences in biological processes between the 2 clusters in the 3 cohorts. As shown in the CGGA database (Figure 2C), GSEA indicated that graft-versus-host disease, malaria, and toll-like receptor signaling were significantly enriched in cluster 2 patients. In contrast, cluster 1 cases showed enrichment of ribosome, steroid hormone biosynthesis, porphyrin, and chlorophyll metabolism. The significant pathways were selected based on the nominal P < .05 and false discovery rate (FDR) < 0.25. In the other 2 databases, significant differences in biological processes between the 2 GBM subgroups were also observed (Figures S5A and S6A).

We then referred to established signatures to evaluate the enrichment scores of 10 classical oncogenic pathways ¹⁸ (Figure 2D). Oncogenic pathways such as NFR2, MYC, phosphatidylinositol 3-kinase (PI3K), and NOTCH-related pathways had higher enrichment scores in cluster1 but lower scores in cluster 2 (Figure 2E), consistent with a previous study. ¹⁸ However, abnormal TGF-β, cell cycle, TP53, and WNT signaling was enriched in cluster 2. Additional validation was also performed using the meteoric cohort (Figures S5B-C and S6B-C).

The 2 factors, cell cycle progression (CCP) and epithelial-mesenchymal transition (EMT), determining tumor proliferation and metastasis, were also compared between the 2 clusters. Both the factors were significantly enriched in cluster 2 (Figures 2F, S5D, and S6D). Moreover, the box plots and heatmaps of the immune cell fractions in glioma tissues showed significant differences between the 2 subgroups (Figures 2G and H, S5E-F, and S6E-F).

Differences in the Gene Mutation Frequency

The top 30 mutated genes with the highest mutation frequency in the 2 clusters of the CGGA cohort were separately screened using the Maftools package, depicting that the mutation frequency of TP53 gene was the highest in both clusters and was higher in cluster 2 than that in cluster 1 (49% vs 45%; Figure 2I and J). The TP53 tumor suppressor gene is frequently mutated in human cancers. A mutant TP53 RNA expression signature shows significant correlation with reduced survival in 11 cancer types. ²⁶ Moreover, the IDH1 gene also had a high mutation frequency in the cluster 2 subgroup. Mutation of IDH1 gene at R132 can convert the α-ketoglutarate (α-KG) into α-hydroxyglutaric acid (2-HG) and the accumulation of 2-HG will further lead to the tumorigenesis. Along with TP53 and IDH1 mutations, ATRX gene also showed a higher mutation frequency in cluster 2. The 3 mutated genes could promote the malignant development of GBM together.

Screening of the DEGs

Total 17 438 GM-related DEGs were exploited between the 2 clusters with an absolute Log2 (fold change) ⩾0.585 and P ⩽ .05, including 455 upregulated and 281 downregulated genes (Figure 3A). As shown in Figure 3B, the heatmap depicted the top 20 up- and downregulated genes, and the top 5 genes were further screened based on the order of the absolute Log2 (fold change; Table 1). The results indicated 4 out of the 5 genes were upregulated in cluster 2, promoting the GM and indicating a poorer prognosis. Specifically, F13A1 is associated with tumor-associated macrophages and poor prognosis in GBM ²⁷ ; CXCL14 promotes GBM cell migration and proliferation ²⁸ ; TOX3 is involved in transcription regulation and correlates with GBM survival ²⁹ ; DKK1 acts as a Wnt inhibitor but promotes glioma malignancy^30,31; and NNMT is linked to mesenchymal GBM subtype and poor outcome.^32,33 In addition, we identified the top 30 GM-related genes significantly associated with the survival of GBM by conducting a univariate Cox regression analysis with P < .05 (Table S2).

Figure 3.

DEGs screening and the guidance of clustering for clinical prognosis and immunotherapy in GSE16011 dataset: (A) the volcano plots of DEGs, including 455 upregulated and 281 downregulated genes, (B) heatmap showing the top 20 up- and downregulated genes in the 2 clusters, (C) cluster 1 showed higher OS than cluster 2 in GBM patients by gender (P > .05), (D) the expression of 26 immune-checkpoints in 2 clusters, (E) the bar plot of the estimated IC50 of the top 50 ICIs screened based on the T value. The T cell inflammation score (F), TIDE algorithm (G), TME score (H), and submap analysis (I) all showing the poorer response to immunotherapy in cluster 2 subgroup.

Table 1.

The Top 5 GM-Related DEGs Associated with GBM Survival.

mRNA	Log2 fold change	P-value	Q-value
F13A1	1.807	2.31E−09	1.47E−07
CXCL14	1.707	3.40E−07	7.92E−06
TOX3	−1.680	6.70E−09	3.35E−07
DKK1	1.680	9.05E−07	1.74E−05
NNMT	1.639	5.58E−08	1.85E−06

The Guiding Significance for the Clinical Prognosis and Treatment

Although the OS analysis revealed that the GBM patients (both males and females) in cluster 1 had a better prognosis than those in cluster 2, the difference was not significant (P > .05; Figure 3C). Immunotherapy, such as ICIs, still exhibits a low complete response rate (~10%) and requires further investigation. Notably, cluster 2 had significantly higher gene expressions than cluster1 in several immune checkpoints, including BTLA, CD274, CD28, CEACAM1, HAVCR2, TNFRSF18, TNFRSF9, and TNFSF14 (Figure 3D). We then selected the top 50 ICIs based on the above immune checkpoints. As shown in Figure 3E, there was a significant difference between the 2 clusters with the absolute t > 5 (P < .05). However, cluster 1 exhibited a greater number of significantly effective ICI options than cluster 2. Moreover, the T cell inflammation score, TIDE algorithm, and TME scores were used to were used to predict the response to immunotherapy among different subtypes (Figure 3F-H). According to the results, the cluster 1 might have fewer clinical side effects and better response to immunotherapy than cluster 2. The SubMap algorithm also depicted the similar results (Figure 3I). Cluster 1 displayed a higher response rate to ICIs, indicating that cluster 1 may benefit from immunotherapy. Therefore, these data suggest that immunometabolism subtypes can help identify potential patient groups that would benefit from immune checkpoint inhibitors.

Concordance Between the Classifier and Original Clustering Results

The top 30 upregulated DEGs associated with the survival of GBM were identified as candidate genes for model training. We performed LASSO Cox regression with cross-validation (Figure 4A and B). Then, the GBM samples from the GSE16011 cohort were then randomly divided into training and testing subsets at a 1:1 ratio to validate the model’s accuracy. The results indicated that the AUC values for the training and testing samples were 0.903, 0.942, and 0.901 (Figure 4C), suggesting that the model was successfully established with high accuracy. In addition, there was a high agreement between the classification results predicted by the classifier and the original clustering assignments (Figure 4D).

Figure 4.

(A) Robust GM-based classifier was defined to evaluate immunotherapy response, (A) tenfold cross validation for tuning parameter selection in the LASSO model, (B) LASSO coefficient profiles of the top 30 survival-associated GM-related genes from the upregulated DEGs in GSE16011 dataset, (C) accuracy-dependent ROC curves validation at the training and testing samples of GSE16011 database, (D) a decision tree showing the high consistency of the 2 subtyping methods, (E) Kaplan-Meier analysis of patients between the 2 predicted clusters, (F) the expression of 26 immune-checkpoints in 2 predicted clusters, (G) the bar plot of the estimated IC50 of the top 50 ICIs screened based on the T value. The T cell inflammation score (H) and TIDE algorithm (I) showing the poorer response to immunotherapy in cluster 2 subgroup, keeping a consistent with the results of the original clustering.

Furthermore, the clinical utility of the classifier was evaluated in terms of prognosis predictions and immunotherapy efficacy. In Figure 4E, the OS of predicted cluster 2 was lower than that of predicted cluster 1, but no significant difference was observed. However, the results for some immunotherapy responses remained largely consistent with the original clustering. The predicted cluster 2 subgroup also had significantly higher gene expressions than predicted cluster1 in some immune checkpoints, including CD274, CD28, CEACAM1, HAVCR2, TNFRSF18, and TNFSF14 (Figure 4F). Based on the ICIs in Figure 4G, predicted cluster1 was more sensitive to immunotherapy than predicted cluster 2, as validated by the T cell inflammation score, TIDE algorithm, and other metrics (Figure 4H and I).

Discussion

GM Subtyping Discovery and Validation

Aberrant reprograming of metabolism is a hallmark in various malignant tumors. Previous studies had described the Warburg effect in GBM and other tumor types.^3,34 According to the Warburg hypothesis, cancers are partly caused by impaired mitochondrial function and OxPhos, with cancer cells producing energy through glucose fermentation (such as aerobic glycolysis) while exhibiting limited OxPhos capacity.^3,34 However, Reinfeld et al ⁷ observed that myeloid cells had the greatest capacity to take up intratumoral glucose, followed by T cells and cancer cells. In contrast, the cancer cells showed the highest uptake of glutamine, indicating that targeting GM could hamper the growth of cancer cells. Currently, the somatotype of GBM mainly depends on IDH1 mutation status, which plays a vital role in glucose metabolism by converting isocitrate to α-KG. ³⁵ In our study, we attempted to subtype the GBM using GM-related genes and to validate this classification by defining a robust classifier for clinical prognosis and treatment.

We reviewed the literature on glutamine transport, decomposition, and synthesis for the first time and obtained 13 GM-related genes from the MSigDB. ³ Based on the 13 genes, the somatotype of GBM samples from 3 databases was analyzed. The training set comprised gene expression profiles and corresponding clinical information from 237 GBM patients in the CGGA database, while another 219 patients from the GEO (GSE16011 and GSE74187) database were used for validation. First, the 237 samples were clustered into 2 distinct subclasses (k = 2). For k > 2, clustering stability did not improve significantly, a finding validated in the other 2 databases. Thus, k = 2 was selected as the optimal number of clusters, indicating that GBM can be subtyped based on the GM-related genes.

Differences in Pathways, Immunity, and Clinical Characteristics Between the Two GM Subtypes

Next, we explored differences between the 2 clusters regarding the clinical indices, genomic and transcriptomic biomarkers. OS differed significantly between the 2 clusters in all 3 databases except for CGGA, suggesting that GM-related genes are involed in the malignant GBM progression and had prognostic implications. No significant differences were observed in terms of age and gender. We found that males had a higher incidence of GBM than females in both cluster 1 and cluster 2, consistent with previous studies.^36,37 GSEA suggested significant differences in the gene expression of patients between the 2 subgroups, with pathways related to poorer prognosis (eg, graft-vs-host disease, malaria) enriched in cluster 2. Recent studies have shown that oncogenic pathway activation could upregulate certain metabolic pathways. ³⁸ Indeed, we found that GBM samples in cluster 1 exhibited several glycolysis-related oncogenic pathways, including NFR2, MYC, PI3K, and NOTCH pathways. ¹⁸ Bollong et al ³⁹ found that inhibition of the glycolytic enzyme PGK1 leads to accumulation of the reactive metabolite methylglyoxal, which selectively modifies KEAP1 to form a methylimidazole crosslink between proximal cysteine and arginine residues (MICA). This posttranslational modification results in the KEAP1 dimerization, NRF2 accumulation and activation of the NRF2 transcriptional program, demonstrating direct crosstalk between glycolysis and the KEAP1–NRF2 axis. MYC is a potent oncogene that drives multiple cancer-related phenotypes, including altered metabolism to support rapid cell growth and proliferation. It increases glycolytic flux in GBM cells, creating a vulnerability to glycolytic inhibition. ⁴⁰ PI3K/Akt regulates fructose 2,6-bisphosphatase (PFKFB2) expression and enhances glycolysis, closely associating with a variety of enzymatic biological effects and glucose metabolism. ⁴¹ Genes encoding proteins involved in glucose uptake, glycolysis, lactate-to-pyruvate conversion, and repression of the tricarboxylic acid cycle are direct transcriptional targets of Notch signaling. ⁴²

In contrast, abnormal TGF-β, cell cycle, TP53, and WNT signaling were enriched in cluster 2. TGF-β is a pleiotropic cytokine associated with poor prognosis in multiple tumor types, including GBM ²⁵ ; it promotes immunosuppression, angiogenesis, metastasis, EMT, fibroblast activation, and desmoplasia. ³⁵ Moreover, the p53, the RB, and the receptor tyrosine kinase (RTK) signaling pathways are frequently implicated in GBM. ⁴³ In particular, the p53 pathway is affected by genomic alterations (CDKN2 A deletion, MDM2/MDM4 amplification, TP53 mutation/deletion) in 87% of samples. ⁴⁴ We found that the mutation frequency of TP53 was relatively higher in cluster 2. TP53, an independent predictor of OS regardless of histology, ⁴⁵ was associated with a poor prognosis. Activated WNT/β-catenin pathway also drives EMT and predict glioma recurrence and poor prognosis. ⁴⁶ Both proliferation- and metastasis-related factors were significantly enriched in cluster 2, correlating with a poorer prognosis. ³⁸

We then compared immune cell infiltration between clusters. Cluster 2 showed higher infiltration of immunosuppressive cells, such as tumor-associated macrophages (TAMs) and Tregs, than cluster 1. Contrary to the current polarization dogma, glioma-infiltrating TAMs exhibit a continuum phenotype between M1- and M2-like states, resembling M0 macrophages. ⁴⁷ M0-like macrophages, a feature of GBM malignancy, revealed distinct assembly in glioma with high level of EGF containing fibulin-like extracellular matrix protein 2 (EFEMP2), which closely correlats with OS. ⁴⁷ Additionally, a cytokine storm induced by M2 macrophages upon neoadjuvant treatment promotes tumor growth and progression. ⁴⁸ Conversely, effect immune cells (NK cells, CD4⁺ T cells) were decreased in cluster 2, transforming the tumor microenvironment into a microenvironment conducive to tumor survival, growth, and uncontrollable proliferation, eventually causing immune escape of tumor cells. ⁴⁹ Consistently, EVA1B correlates with immune cell infiltration and prognosis in glioma, ⁵⁰ and the immunosuppressive microenvironment of cluster 2 mirrors the EVA1B-associated signature.

Clinical Application Value

We further explored the clinical prognosis and treatment value of clustering. First, we found that GBM patients in cluster 2 had shorter OS than those in cluster 1 in both males and females, providing clues for its potential clinical application. Next, we examined 26 immune checkpoints in GBM from the GSE16011 database. The expression of BTLA, CD274, CD28, CEACAM1, HAVCR2, TNFRSF18, TNFRSF9, and TNFSF14—especially HAVCR2—was significantly higher in cluster 2 than in cluster 1. HAVCR2 (TIM-3) can inhibit antitumor immunity by mediating T cell depletion. ⁵¹ Interestingly, these results might suggest better efficacy and greater sensitivity to immunotherapy in cluster 2. However, more immune checkpoint inhibitors (ICIs) were selected based on the above checkpoints in cluster 1 than in cluster 2. As Wang et al ¹ found, heterogeneity of inhibitors targeting integrin (ITG) subunits may explain why ITG-targeted inhibitors show only limited efficacy in a small group of lung cancer patients. The heterogeneity of expression and distinctly different biological roles of immune checkpoint subunits between the 2 clusters may be 1 reason for the differential treatment responses. Moreover, the T cell inflammation score, TIDE algorithm, and TME scores all indicated a poor response to immunotherapy in cluster 2, suggesting a propensity for drug resistance and the need for developing new immunosuppressants.

Based on these findings, this GM-based subtyping system has several potential clinical applications. First, prognostic stratification: patients in cluster 2 consistently showed shorter OS, identifying a high-risk subgroup that may benefit from more aggressive treatment. Second, guidance for immunotherapy: given that cluster 2 shows elevated immune checkpoints but poor predicted response, these patients may be more suitable for combination strategies (eg, checkpoint blockade plus glutamine metabolism inhibitors), whereas cluster 1 patients with a glycolytic-dominant phenotype might benefit from glycolysis-targeting agents. Third, combination metabolic therapy: the glutaminase inhibitor CB-839 has entered clinical trials ¹⁶ ; cluster 2 patients, with higher glutamine dependency, might derive greater benefit from CB-839-based combinations. Therefore, clustering based on GM-related genes may become a promising classification for the diagnosis and treatment of GBM patients in future clinical practice.

However, the clinical application of this subtyping system must consider ethnic and population heterogeneity. To enhance geographical representativeness, we incorporated the latest molecular typing data from the Chinese population. ⁵² Qu et al established a large-scale Chinese primary glioma cohort (798 cases, including 355 GBM) and demonstrated that Chinese glioma patients have significantly longer survival than CGGA/TCGA cohorts, highlighting the need for population-specific reference data. Another study on East Asian GBM (367 cases) identified 4 molecular subtypes, including a metabolic subtype (active lipid metabolism) and an immunoregulatory subtype. Importantly, the metabolic subtype is closely related to our GM-based classification, as glutamine metabolism is a core component of tumor metabolic reprograming. Furthermore, East Asian GBM patients show a lower frequency of EGFR mutations and lack a typical EGFR-driven subtype, which may influence the expression patterns of GM-related genes. Therefore, the distribution of GM-based subtypes may differ across ethnic populations, and the applicability of our classifier in East Asian patients requires direct validation using Chinese cohorts. Future studies should integrate multi-omic data from both Western and Chinese populations to establish a more generalizable GM subtyping system.

Moreover, extending GM-based subtyping to other glioma grades and age groups, such as adolescents and young adults (AYA), may reveal additional heterogeneity. A recent population-based study highlighted the importance of molecular classification in redefining treatment paradigms for younger glioma patients. ⁵³ The AYA population (ages 15-39) exhibits distinct molecular features and clinical outcomes compared to older adults. Because glutamine metabolism is critical for rapidly proliferating cells and metabolic reprograming may differ across age groups, applying our GM-based classifier to lower-grade gliomas and AYA cohorts could identify subtype-specific vulnerabilities. Future studies should validate our findings in these populations.

Classifier Construction and Performance

We defined a classifier for GBM subtyping by selecting the top 30 upregulated DEGs associated with the survival of GBM from GSE16011 database. Satisfactorily, our model achieved high accuracy (AUC > 0.9) and maintained a high degree of consistency with the original clustering results in GBM subtyping and immunotherapy prediction, outperforming other classification models. ⁴⁴ This is the first clinical model to classify GBM from the perspective of GM, and it has guiding significance for clinical immunotherapy. However, the prognosis value of the model was not significant due to insufficient samples in the database, and further exploration with larger datasets will be needed in the future.

Limitations

First, tumor heterogeneity is an essential feature of gliomas, ⁵⁴ with different microenvironment characteristics across tumor sites. ²⁰ All data and corresponding information were collected from public databases, making it impossible to detect immune status in the same or different tumor regions. Therefore, this gene signature should be validated in well-designed, multicenter, prospective studies. Second, the clinical information in the 3 databases only included age, gender, and OS. Moreover, one of the validation datasets lacked some clinical data, so the prognostic value should be further validated using more external datasets.

Conclusion

Glutamine is an essential nutrient for tumor growth. In clinical practice, the classification of GBM mainly depends on IDH1 status, which is associated with glucose metabolism. In this study, we successfully demonstrated that GBM (garde IV) could be newly subtyped into 2 distinct clusters based on GM-related genes tusing consensus clustering. This classification had significant implications for the prognosis and clinical treatment of GBM and indicates that GM is closely related to the malignant progression of GBM. Moreover, we defined a robust GM-related classifier for GBM subtyping, which may further guide clinical prognosis and treatment. In summary, we have revealed the metabilic heterogenicity of GBM and provided a novel subtyping method to better understand the importance of GM in tumors.

Supplemental Material

sj-docx-1-cix-10.1177_11769351261449032 – Supplemental material for Explore the Heterogeneity of Glioblastoma Based on Genes Related to Glutamine Metabolism

Supplemental material, sj-docx-1-cix-10.1177_11769351261449032 for Explore the Heterogeneity of Glioblastoma Based on Genes Related to Glutamine Metabolism by Ling Wang, Tong Chen, Jiajia Chen, Yiwen Yang, Linglin Sun, Chao Gu and Chenliang Cao in Cancer Informatics

Supplemental Material

sj-docx-2-cix-10.1177_11769351261449032 – Supplemental material for Explore the Heterogeneity of Glioblastoma Based on Genes Related to Glutamine Metabolism

Supplemental material, sj-docx-2-cix-10.1177_11769351261449032 for Explore the Heterogeneity of Glioblastoma Based on Genes Related to Glutamine Metabolism by Ling Wang, Tong Chen, Jiajia Chen, Yiwen Yang, Linglin Sun, Chao Gu and Chenliang Cao in Cancer Informatics