Characterization of hypoxia-related molecular clusters and prognostic riskScore for glioma

Background: Gliomas represent a significant burden in the realm of central nervous system (CNS) malignancies, accounting for approximately 30% of all primary brain tumors and a striking 80% of malignant cases. The incidence of gliomas is observed to escalate with advancing age, exhibiting a marginally higher prevalence in the male population. Among these tumors, high-grade gliomas, particularly glioblastoma multiforme (GBM), are characterized by their aggressive nature and dire prognosis. Conventional therapeutic approaches, including surgical intervention, radiotherapy, and chemotherapy, have demonstrated limited efficacy, underscoring an urgent need for the development of targeted therapies and enhanced mechanistic understanding to improve patient outcomes.

Methods: In this study, we aimed to deepen our understanding of the role of hypoxia, a critical factor in cancer progression, within gliomas. Using comprehensive datasets from The Cancer Genome Atlas (TCGA) and Chinese Glioma Genome Atlas (CGGA), we classified gliomas into two distinct subgroups based on hypoxia-related gene expression profiles: C1 and C2. This classification facilitated a comparative analysis of prognostic outcomes and tumor microenvironment characteristics between the two subgroups.

Results: Our findings revealed that patients within the C1 subgroup exhibited significantly poorer prognoses, with an upregulation of genes intricately linked to various tumor progression pathways. Moreover, the immune microenvironment within the C1 subgroup appeared more favorable for tumor survival and growth, coupled with a notable increase in genomic instability compared to the C2 subgroup. A prognostic scoring system developed from key hypoxia-related factors demonstrated substantial predictive value across multiple cohorts.

Conclusion: Ultimately, we identified four core hub genes—SOCS3, CLCF1, PLAUR, and LIF—whose expression was validated under hypoxic conditions via Western blot analysis in glioma cell lines. This study employs bioinformatics to elucidate glioma subtypes, highlighting significant prognostic and functional disparities. The experimental validation of candidate molecules paves the way for future research aimed at unraveling their roles and underlying mechanisms in glioma pathophysiology, potentially guiding novel therapeutic strategies.

Introduction

Glioma is a common cancer that affects the brain and central nervous system, accounting for approximately 30% of all primary brain tumors and 80% of all malignant brain tumors (1–3). Notably, the incidence of glioma increases with age, with slightly higher rates in males than females (4, 5). Unfortunately, the outcomes for these patients are frequently unfavorable. High-grade gliomas, particularly glioblastomas (GBM), are extremely aggressive and are life-threatening. Treatment options for gliomas include surgery, radiation therapy, chemotherapy, or a combination of these approaches (6, 7). However, the limitations of the current treatments underscore the need for continued research and development of new, more targeted treatment approaches. Therefore, improving treatment strategies and understanding the underlying mechanisms of gliomas are vital.

Hypoxia, or low oxygen tension, is an important factor in the growth, development, and progression of cancer (8–12). Several pathways are involved in the response of cancer cells to hypoxia, including the hypoxia-inducible factor (HIF) pathway, which promotes angiogenesis, cell survival, and tumor invasion (9, 13–15). Furthermore, hypoxia affects the immune response, making tumor cells more resistant to immune cell recognition and attack (16, 17). Therefore, targeting the hypoxic pathways in cancer therapy has emerged as an attractive approach for improving cancer treatment outcomes. Importantly, several types of cancers, such as breast cancer, non-small cell lung cancer, and malignant melanoma, have been reported to be affected by hypoxia (18–24). Nevertheless, advancements are crucial for a deeper comprehension of hypoxia’s impact on gliomas. This encompasses elucidating the intricate hypoxic pathways, precisely quantifying hypoxia’s extent and distribution, and innovating effective therapies targeting hypoxia in gliomas.

By employing unsupervised clustering analysis utilizing hypoxia-related gene expression profiles, we successfully categorized gliomas into two distinct subtypes. Each subtype exhibits unique somatic alterations, immune cell composition, metabolic characteristics, and clinical outcomes. Through co-expression network analysis, we identified four crucial hub genes, three of which were upregulated in response to hypoxia treatment in glioma cell lines. These findings have the potential to greatly advance clinical diagnosis and mechanistic research in the field of glioma.

Materials and methods

Patients and samples

TCGA (The Cancer Genome Atlas) is a pan-cancer research program initiated jointly by the NCI and NHGRI in the United States, which includes genomic, transcriptomic, and clinical data from 33 types of cancer and over 11,000 patients. The Chinese Glioma Genome Atlas (CGGA) is a specialized database led by Beijing Tiantan Hospital, focusing on multi-omics research of primary and recurrent brain gliomas. Its core integrates genomic, transcriptomic, epigenetic data, and complete clinical follow-up information, with a sample size exceeding 2,000 cases (as of 2023). Data repositories, including TCGA (LGG and GBM) and CGGA (mRNA_325 and mRNA_301) databases, were searched for available glioma genomics, transcriptomics, and clinical information. Our research scope includes patients from WHO stage II to IV. The aforementioned datasets were downloaded from various websites (https://xenabrowser.net/datapages/ and http://www.cgga.org.cn/).

Identification of hypoxia cluster

The WINTER_HYPOXIA_METAGENE gene set obtained from the Molecular Signatures Database (MSigDB) (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp) was selected for k-means clustering. Normalized expression data were subjected to k-means clustering using the Consensus ClusterPlus R package. The optimal number of clusters was determined by selecting the k-value at which the magnitude of the cophenetic correlation coefficient decreased. Additionally, principal component analysis (PCA) was employed to gain a better understanding of variations between clusters.

Differentially expressed gene analysis and functional enrichment

Differential expression analyses were performed using the “DESeq2” package. In this analysis, genes with an absolute log2 fold change greater than 1 and a false discovery rate (FDR) <0.05 were considered as differentially expressed genes between two clusters. Gene ontology and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was conducted using the “ClusterProfiler” R package, with a cutoff value of FDR < 0.05 (25). To investigate the difference in biological process terms between two clusters, GSEA and GSVA were performed using the R package, considering the gene sets of “h.all.v7.5. symbols” and “c2.cp.kegg.v7.5. symbols” downloaded from MSigDB.

Estimation of immune infiltration and tumor purity

The CIBERSORT algorithm was used to estimate the relative fractions of the 22 immune cell types in each sample (26). ESTIMATE was used to evaluate the immune cells and stromal contents of each sample (27).

Somatic mutation identification

Somatic mutation data of all patients categorized under the “Masked Somatic Mutation” category in TCGA were processed using VARSCAN software (https://portal.gdc.cancer.gov/). Mutation analysis and visualization were performed using the “maftools” package. To detect differences in copy number alterations between the subtypes, GISTIC2.0 analysis was conducted. Here, loci with a GISTIC value greater than 1 or less than -1 were defined as amplification or deletion, respectively.

Weighted co-expression network

A normalized expression matrix was used to construct a weighted co-expression network (WGCNA) using R package (28). A co-expression network was created using the blockwiseModules function with default parameters. The modules were evaluated for significance based on the correlation between module eigengenes and information using Pearson’s test. Hub genes were then chosen based on the modular connectivity of each gene and its relationship with the phenotypic traits in the hub module.

Cell lines and culture

LN229 and U118 cells were procured from the American Type Culture Collection (Rockville, MD) and maintained in Dulbecco’s modified Eagle’s medium supplemented with high glucose, sodium pyruvate, 10% fetal bovine serum, and 1% penicillin–streptomycin. To simulate hypoxic conditions, cells were cultured in a hypoxia chamber with 94% N₂, 5% CO₂, and 1% O₂ at 37°C.

Western blotting

Protein extraction and western blot analysis were performed as previously described (29). The antibodies used were anti-HIF-1 (1:5000, 20960-1-AP), anti-SOCS3 (1:2000, 14025-1-AP), anti-PLAUR (1:1000, 10286-1-AP), and anti-LIF (1:500, 26757-1-AP), all of which were obtained from Proteintech Group (Wuhan Sanying, China). The antibodies were validated using a commercial vendor. The bands on each membrane were detected using an ECL kit from Beyotime Biotechnology (Beijing, China), and ImageJ software (National Institutes of Health, Bethesda, MD, USA) was used for quantification.

Bioinformatic and statistical analyses

To detect survival differences between clusters, Kaplan–Meier analysis with log-rank tests was performed. Unpaired Student’s t-tests were used for normally distributed variables, and Wilcoxon rank-sum tests were conducted for non-normally distributed variables when comparing the two groups. The Benjamini–Hochberg method was used to adjust the P-value. For the univariate analysis, we selected factors known to affect outcomes and patient characteristics. Cox regression analysis was then conducted with age, sex, and grade as factors. Cox proportional hazard models were used to examine the impact of various risk factors on event outcomes, and the reliability of the model was assessed using the Schoenfeld residual. To adjust for explanatory confounding variables that were prognostic in the univariate analysis, Cox multivariate analysis was conducted. A hypoxia-related risk score was established by including normalized gene expression values weighted by their Least Absolute Shrinkage and Selection Operator (LASSO) Cox coefficients (30). All statistical analyses were performed using R software, and statistical significance was considered at a P-value < 0.05.

Results

Consensus clustering reveals two hypoxia clusters in glioma

To identify the heterogeneity of hypoxia within gliomas, 240 hypoxia-related genes were obtained from MSigDB for clustering analysis. We conducted univariate Cox regression and residual analyses to identify genes with prognostic significance, resulting in a total of 86 candidate genes for subsequent clustering analysis in both TCGA and CGGA cohorts (Figure 1A, Supplementary Table). Using an unsupervised consensus clustering analysis, we identified two clusters, C1 and C2 (Figure 1B). Survival analysis revealed that patients with the C1 subtype were likely to have worse outcomes than those with the C2 subtype (P <0.0001, log-rank test, Figure 1C). Furthermore, we performed PCA to confirm the assignment of clusters and robust differences in expression patterns between the two clusters (Figure 1D).

Figure 1

Figure 1. Hypoxia gene profile of glioma yielded two clusters in the TCGA cohort. (A) Venn diagram showing 86 candidate gene. (B) Heatmap displaying consensus clustering with the robust classification (k = 2). (C) Survival analysis of C1 and C2 based on OS. The P-value is calculated by the log-rank test between clusters. (D) Principal component analysis (PCA) of two clusters using whole transcriptome data.

Identification of differentially expressed genes and functional analysis

We performed differential analysis using C1 as the experimental group. Based on the cutoff criteria of |log2 (fold change) | >1.0 and FDR <0.01 using R package “DESeq2”, we identified a total of 10,482 differentially expressed genes (7,622 upregulated and 2,860 downregulated), as shown in the volcano plot (Figure 2A). To explore the functional status, we calculated KEGG signaling scores using the GSVA method and ran a GSEA analysis for the HALLMARKER pathways. Our results revealed that tumor malignancy-related signaling pathways, such as angiogenesis, epithelial–mesenchymal transition, and the cell cycle, were enriched in C1. In contrast, gap and tight junctions showed lower enrichment scores (Figure 2B). Furthermore, C1 exhibited higher enrichment of the Janus Kinase/Signal Transducer and Activator of Transcription (JAK-STAT), Tumor Necrosis Factor alpha/Nuclear Factor kappa-light-chain-enhancer of activated B cells (NF-κB), and P53 pathways (Figures 2C-E). In summary, our findings suggest that C1 glioma cells exhibit a higher degree of hypoxia and a more significant malignant phenotype.

Figure 2

Figure 2. Functional insights of distinct hypoxia clusters in glioma. (A) Volcano plots showing genes that are differentially expressed in C1 or C2 samples. Red dots indicate genes upregulated in C1; blue dots indicate genes upregulated in C2. (B) Comparison of KEGG functions scores(calculated by GSVA analyze) between C1 and C2 (the color of the squares indicates the high and low average scores of all samples within each cluster). (C-E) GSEA analyses of HALLMARK gene set between C1 and C2, including JAK-STAT pathway, NF-κb pathway and P53 signaling.

Immune characteristics and hypoxia clusters

Next, we utilized several previously reported immune-related tools to decode the immune infiltration of the immune subtypes. First, we used the CIBERSORT method to compare the composition of infiltrating immune cells between the groups. C2 had a higher percentage of monocytes and plasma cells compared to C1, whereas C1 had a higher percentage of macrophages, including M0,M1 and M2(Figure 3A). On the other hand, the ratio of CD4+ T cells and CD8+ T cells in C1 is slightly higher than that in C2, although the total amount is relatively low. Subsequently, we determined stromal and immune scores using the ESTIMATE algorithm. The C1 cluster showed higher immune and stromal scores and lower tumor purity than did cluster C2 (Figures 3B-D). Furthermore, we use the EaSIeR package, a tool based on system tumor microenvironment characteristics, to quantify immune cell composition and predict different features of immune response through intracellular and intercellular communication,We discovered that the C1 cluster might have a higher immunotherapy response score compared to cluster C2 (Figure 3E) (31). In addition, most classical checkpoint genes, such as PDL1, CTLA4, CD86, HAVCR2, LGALS9, and CD48, were highly expressed in C1 tumors (Figure 3F), indicating an increased level of immunosuppression in these tumors.

Figure 3

Figure 3. Investigation of the immunologic intertumoral heterogeneity. (A) Relative proportion of 22 infiltrating immune cells estimated by CIBERSORT algorism between C1 and C2. (B-D). Comparison of immune, stromal, and tumor purity scores from ESTIMATE for C1 and C2. (E) Immune therapy response scores(calculated by EaSIeR algorism) between C1 and C2. (F) Heatmap displaying the average expression levels of inhibitive checkpoint genes. *P < 0.05; **P < 0.01;***P < 0.001.

Genomic alterations of hypoxia clusters

It has been reported that tumor genomic mutations correlate with oxygen supply levels (32, 33). We analyzed the differences in gene mutations between the two clusters and specifically looked at the top 20 genes with different mutation frequencies in gliomas. C2 showed significantly higher frequencies of IDH1 and ATRX (93% and 40%, respectively) than did C1 (25% and 23%, respectively). The mutation rate of TP53 in C1 was similar to that in C2 (43%–48%) (Figures 4A, B). In addition, there was significant heterogeneity in the CNV profiles between the two clusters. C1 had more amplified and deleted variant samples as well as a higher tumor mutational burden (TMB) than did C2 (Figures 4C-E). These results indicate that gene mutations may be associated with the hypoxic cluster phenotype.

Figure 4

Figure 4. Comparison of genomic alterations between C1 and C2. (A, B) The top 20 mutated genes in glioma between C1 and C2. The colors of rectangles in the body of the heatmap indicate different types of somatic mutations. (C-E) The number of mutations and copy number aberrations including tumor mutation burdens, amplifications and deletions *P < 0.05; **P < 0.01;***P < 0.001.

Construction of the prognostic model

To refine prognostic biomarkers from the initial 86 hypoxia-related genes, we employed a multi-step filtering approach. First, genes were subjected to univariate Cox proportional hazards regression (P < 0.05), followed by multivariate Cox regression and Schoenfeld’s residuals test (P > 0.05 for proportionality assumption), yielding 60 candidate genes. Subsequently, a Least Absolute Shrinkage and Selection Operator (LASSO) Cox proportional hazards regression with 10-fold cross-validation was applied to this candidate set to mitigate overfitting and identify core prognostic features. The optimal penalization coefficient lambda (lambda.min= –2.7894) minimized the cross-validation error, resulting in a 9-gene prognostic signature (ANXA5, DDIT3, FABP5, LOX, PLAUR, SLC16A1, SLC20A1, TBPL1, and TFRC). Their prognostic significance was statistically meaningful and consistent with the PH assumption, and they were also free from multicollinearity (Figures 5A–C). A prognostic index was calculated for all cancer samples using the following formula: risk score = ANXA5×0.15648040+DDIT3×0.01273471+FABP5×0.21023267+LOX×0.0112542+PLAUR×0.22763644+SLC16A1×0.02450256+SLC20A1×0.10289532+TBPL1×-0.01552882+TFRC×0.1555086.

Figure 5

Figure 5. Construction of the hypoxia-associated prognostic model. (A, B) Partial likelihood deviance for the LASSO regression and Lasso regression analysis. (C) Forest plot for nine robust candidate genes. (D). Nomogram obtained using multivariate Cox regression analysis for predicting the proportion of patients with OS. (E) Plots depicting the calibration of the model in terms of the agreement between the predicted and observed OS. The model performance is shown in the plot relative to the 45-degree line, which represents a perfect prediction. (F) Concordance indices (C-index) of the nomogram-based signature in TCGA, CGGA_325, and CGGA301 datasets. (G) ROC curve AUC plotted for different durations of OS for nomogram-based signatures in TCGA, CGGA_325, and CGGA301 datasets.

In the multivariate Cox proportional hazards analysis (Table 1), the riskScore was identified as an independent prognostic biomarker (HR = 1.92, 95% CI = 1.69–2.19, p <0.001). To improve prognostic prediction in patients with glioma, we developed a nomogram by integrating three independent predictors – risk score, age, and grade – into a multivariate Cox regression model. We evaluated and validated this nomogram using data from TCGA, CGGA_325, and CGGA_301 databases. The nomogram generated a score that predicted the 1-, 3-, and 5-year overall survival (OS) rates for individual patients (Figure 5D). The performance of the nomogram in predicting patient OS was evaluated using a calibration plot, which demonstrated that it accurately predicted patient survival according to an ideal model (Figure 5E). The concordance indices for predicting OS with the nomogram model were 0.857, 0.775, and 0.724 in TCGA, CGGA_325, and CGGA_301 cohorts, respectively (Figure 5F), indicating good predictive accuracy. Furthermore, for the 1- to 5-year OS prediction, the time-dependent receiver operating characteristic curve demonstrated excellent predictive performance (area under the curve >0.7) across all three cohorts (Figure 5G).

Table 1

Table 1. Univariate and multivariate Cox regression analyses of riskScore with WHO grade, age, and sex.

Landscape of riskScore associations: clinical subgroups, functional activities, and metabolic signatures

Next, we analyzed the relationship between risk score and clinical parameters. As shown in Figure 6A, the risk score of the majority of the C1 cluster was higher than that of the C2 cluster. Additionally, there was a negative correlation between the risk score and the Proneural score of glioblastoma (r = -0.4, p <0.001), whereas a positive correlation was observed with the mesenchymal score (r = 0.33, p <0.001) Figures 6B, C. Moreover, the risk score was higher in patients with GBM, World Health Organization grade IV, wild-type IDH-1, and non-methylated MGMT promoter Figures 6D–G.

Figure 6

Figure 6. Relationship between clinical features of gliomas and nomogram-based signature. (A) Distribution of riskScore and hypoxia-related clusters. (B, C) The correlation between riskScore and proneural and mesenchymal scores. (D-G) Boxplot indicating riskScore in different histological type, WHO grades, IDH-1 with 1p19q status, and MGMT promoter methylation from TCGA dataset. (H) Correlations of the riskScore with different metabolism marker gene expression, the color of the scatter points represents the significance of the Pearson correlation analysis, with red indicating a positive correlation, blue indicating a negative correlation, and gray indicating no correlation. *P < 0.05; **P < 0.01;***P < 0.001.

Hypoxia can also alter the expression of genes involved in energy metabolism, leading to changes in protein and enzyme levels. To identify differences in energy metabolism, we examined seven central metabolic pathways: glycolysis, tricarboxylic acid cycle, oxidative phosphorylation, lipid, and glutaminolysis. As shown in Figure 6H, most markers of glycolysis and oxidative phosphorylation were positively correlated with the risk score (r >0.2, P <0.001). However, further studies are required to investigate the global effects of hypoxia on metabolism.

Identification of hub genes through WGCNA analysis

Weighted gene co-expression network analysis (WGCNA) was performed to identify hypoxia-regulated hub genes using transcriptomic data from the TCGA samples (n=661 samples). After quality control and normalization, 19,583 protein-coding genes were retained for network construction. A soft thresholding power of β=10 was selected to achieve scale-free topology (scale-free fit index >0.8);. Hierarchical clustering with dynamic tree cutting (minModuleSize = 50) identified 22 distinct co-expression modules. Among these modules, the blue one was found to have the highest correlation with hypoxia (r = 0.92, p <0.0001) (Figures 7A, B). Biological process analysis revealed that genes in the blue module were enriched in wound healing, angiogenesis regulation, and histone modification. Furthermore, KEGG pathway analysis demonstrated that the genes were associated with the JAK-STAT and NF-κB signaling pathways (Figure 7C). Finally, after identifying hub genes with a high degree of connectivity in the blue module, we found that CLCF1, PLAUR, SOCS3, and LIF were the top four hub genes (Figure 7D). HIF-1 is the most important transcription factor that regulates gene expression under hypoxic conditions. To investigate its role in the glioma cell lines LN229 and U118, we exposed the cells to prolonged hypoxic conditions for 48 h and analyzed the expression of several hub genes using western blotting. Our results, shown in Figure 8, demonstrate that the expression of HIF-1, PLAUR, SOCS3 and LIF was significantly upregulated under prolonged hypoxic conditions. The above evidence proves that hub genes can be modulated under low levels of oxygen by HIF-1.

Figure 7

Figure 7. Identification of four hub genes of Hypoxia cluster. (A, B) WGCNA analysis identified the blue module as the highest correlated module with hypoxia (r = 0.92, p <0.001). (C) GO and KEGG analysis of blue module genes. (D) Hub genes of blue modules filtered by degree of connectivity. *P < 0.05.

Figure 8

Figure 8. Western blotting validation of HIF-1 and hub genes in glioma cell lines under hypoxic conditions. (A, B) Western blotting showing that the protein expressions of HIF-1, PLAUR, SOCS3, and LIF in U118MGand LN229 celllines are upregulated under hypoxic conditions. (C-J) Using Image J to perform semi-quantification of the PLAUR, SOCS3, and LIF protein bands in Figure 8, calculating the relative expression of each protein in the cells based on the quantitative values. All data were analyzed using the t-test and were shown as the mean ± SD (two independent experiments). *P < 0.05; **P < 0.01;***P < 0.001.

Discussion

Tumor hypoxia, characterized by an insufficient supply of oxygen to solid tumors, plays a pivotal role in the development and progression of cancer. This is due to an imbalance between the oxygen demand and supply caused by the rapid growth of tumors. Hypoxia triggers diverse biological processes, including metabolic alterations, angiogenesis, and metastasis, thereby contributing to treatment resistance and unfavorable patient outcomes (34–36). Notably, extensive research has highlighted the profound effect of hypoxia on cancer cells and its potential as an attractive therapeutic target. Moreover, several hypoxia gene signatures have been developed to predict the prognosis of various cancer types, including head and neck, breast, prostate, and bladder cancers. Nonetheless, prior studies have encountered limitations, such as oversimplification of hypoxia gene networks and insufficient integration with conventional prognostic systems. Nevertheless, the significance of tumor hypoxia in prognosis and its potential as a therapeutic target continues to be actively investigated, with substantial promise for enhancing cancer treatment outcomes.

This study comprehensively explored the heterogeneity and molecular characteristics of hypoxic clusters in gliomas. Consensus clustering analysis revealed two distinct clusters, C1 and C2, with C1 being associated with a more aggressive phenotype and poorer prognosis. Differential expression analysis identified genes linked to tumor malignancy, angiogenesis, and epithelial–mesenchymal transition that were upregulated in C1. Immune infiltration analysis showed differences in immune cell composition, indicating higher immunosuppression in C1. Genomic analysis revealed distinct mutation patterns, with C2 exhibiting a higher frequency of IDH1 and ATRX mutations. The prognostic model constructed using LASSO regression and nomogram integration yielded accurate survival predictions. Energy metabolism analysis revealed correlations between hypoxia and glycolysis/oxidative phosphorylation markers. WGCNA identified the blue module as being highly correlated with hypoxia and enriched in wound healing, angiogenesis regulation, and histone modification processes. Hub genes within the blue module (CLCF1, PLAUR, SOCS3, and LIF) were found to be modulated by HIF-1 under hypoxic conditions. Overall, this assay enhances our understanding of hypoxic cluster heterogeneity and the molecular characteristics of gliomas. These findings may have implications for targeted therapies and for improving patient outcomes.

Our research focused on investigating the involvement of hub genes in tumor hypoxia, particularly in gliomas. Although these hub genes have been extensively studied in various cancers, their roles in tumor hypoxia remain largely unexplored. For instance, CLCF1 expression has been identified as an independent prognostic factor and a potential target for immunotherapy in hepatocellular carcinoma (37). Additionally, PLAUR promoted anoikis resistance and metastasis in cholangiocarcinoma (34). LIF, a pleiotropic cytokine with diverse roles in different systems, has been implicated in hematopoietic differentiation and has recently emerged as a biomarker and therapeutic target in pancreatic ductal adenocarcinoma (38–40). However, it is important to note that the specific biological functions of these hub genes in the context of tumor hypoxia still require further investigation. Further studies are needed to elucidate the precise roles and mechanisms of tumor hypoxia, potentially opening new avenues for therapeutic interventions and biomarker development.

Conclusion

In summary, our study introduces a novel hypoxia classification for gliomas, encompassing two well-defined clusters exhibiting distinct prognoses, somatic variations, immune infiltration, and metabolic phenotypes. Our model identified three hub genes that were significantly upregulated under hypoxic conditions. Additional in vitro and in vivo experiments are required to gain a comprehensive understanding of the underlying molecular mechanisms driving these findings. These studies will provide valuable insights and shed light on the intricate pathways and processes involved in the observed phenomena.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.

Author contributions

FX: Funding acquisition, Supervision, Investigation, Conceptualization, Software, Writing – original draft, Resources, Validation, Project administration, Data curation, Formal analysis, Visualization, Methodology. XW: Funding acquisition, Validation, Writing – review & editing, Project administration, Formal analysis, Supervision, Conceptualization, Resources, Methodology, Software, Data curation, Investigation, Visualization. CX: Validation, Writing – review & editing, Funding acquisition, Formal analysis, Supervision, Data curation, Software, Project administration, Investigation, Conceptualization, Resources, Writing – original draft, Visualization, Methodology.

Funding

The author(s) declare that no financial support was received for the research, and/or publication of this article.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2025.1605949/full#supplementary-material

References

1. Alberto FD, Francesca DL, and v.W. Danielle FD. Anna, Arterial spin labeling MR imaging for differentiation between high- and low-grade glioma-a meta-analysis. Neuro Oncol. (2018) 20. doi: 10.1093/neuonc/noy095

PubMed Abstract | Crossref Full Text | Google Scholar

2. Bill H D, Heng L, Rui Y, Landon J H, Alexis L Z, Fangping Z, et al. Sensitive and rapid detection of TERT promoter and IDH mutations in diffuse gliomas. Neuro Oncol. (2018) 21. doi: 10.1093/neuonc/noy167

PubMed Abstract | Crossref Full Text | Google Scholar

3. Lei Z, Anshun H, Bohan C, Jinfang B, Jun C, Dianhao G, et al. A HOTAIR regulatory element modulates glioma cell sensitivity to temozolomide through long-range regulation of multiple target genes. Genome Res. (2020) 30. doi: 10.1101/gr.251058.119

PubMed Abstract | Crossref Full Text | Google Scholar

4. Antje W, Tobias K, E. Andrew EH, W. Frank B, Tracy T, Michael P, et al. Glioblastoma in elderly patients: solid conclusions built on shifting sand? Neuro Oncol. (2017) 20. doi: 10.1093/neuonc/nox133

PubMed Abstract | Crossref Full Text | Google Scholar

5. Quinn T O, Warren C, William H, Joshua B R, Justin D L, Michael E B, et al. Sex-specific gene and pathway modeling of inherited glioma risk. Neuro Oncol. (2018) 21. doi: 10.1093/neuonc/noy135

PubMed Abstract | Crossref Full Text | Google Scholar

6. Yushen Q, Satoshi M, Haju K, Erqi L P, Kiran A K, Alexander L C, et al. Cost-effectiveness of radiation and chemotherapy for high-risk low-grade glioma. Neuro Oncol. (2017) 19. doi: 10.1093/neuonc/nox121

PubMed Abstract | Crossref Full Text | Google Scholar

7. Yuxiang W, Aaron T W, Sevin T, Wei H W, Carlie S, David S K, et al. Targeting therapeutic vulnerabilities with PARP inhibition and radiation in IDH-mutant gliomas and cholangiocarcinomas. Sci Adv. (2020) 6. doi: 10.1126/sciadv.aaz3221

PubMed Abstract | Crossref Full Text | Google Scholar

8. Thakur A, Qiu G, Xu C, Han X, Yang T, NG SP, et al. Label-free sensing of exosomal MCT1 and CD147 for tracking metabolic reprogramming and Malignant progression in glioma. Sci Adv. (2020) 6. doi: 10.1126/sciadv.aaz6119

PubMed Abstract | Crossref Full Text | Google Scholar

9. Bernadette A C, Inge VDW, Catriona J, and S. Robert AJ. Post-transcriptional regulation of homeostatic, stressed, and Malignant stem cells. Cell Stem Cell. (2020) 26. doi: 10.1016/j.stem.2020.01.005

PubMed Abstract | Crossref Full Text | Google Scholar

10. Edward L L and Amato J G. The ever-expanding role of HIF in tumour and stromal biology. Nat Cell Biol. (2016) 18. doi: 10.1038/ncb3330

PubMed Abstract | Crossref Full Text | Google Scholar

11. Qing Y, Mohd Parvez K, Christophe M, Edward L L, Zachary T, Laura M, et al. Suppressing mitochondrial respiration is critical for hypoxia tolerance in the fetal growth plate. Dev Cell. (2019) 49. doi: 10.1016/j.devcel.2019.04.029

PubMed Abstract | Crossref Full Text | Google Scholar

12. Rafal P, Liam W, B. Lauren CJ, S. Maria A, W. Kee H, R. Angela M, et al. Noninvasive imaging of cycling hypoxia in head and neck cancer using intrinsic susceptibility MRI. Clin Cancer Res. (2017) 23. doi: 10.1158/1078-0432.Ccr-16-1209

PubMed Abstract | Crossref Full Text | Google Scholar

13. Fei H, Chien-Feng L, Zhen C, Xian Z, Guoxiang J, Wei-Na Z, et al. The critical role of AMPK in driving Akt activation under stress, tumorigenesis and drug resistance. Nat Commun. (2018) 9. doi: 10.1038/s41467-018-07188-9

PubMed Abstract | Crossref Full Text | Google Scholar

14. Keiko K, Wataru S, Takatoshi Y, Kohei M, Katsuhiko I, Kazuki T, et al. An anionic phthalocyanine decreases NRAS expression by breaking down its RNA G-quadruplex. Nat Commun. (2018) 9. doi: 10.1038/s41467-018-04771-y

PubMed Abstract | Crossref Full Text | Google Scholar

15. Timothy D OB, Peilin J, Neil E C, Maria Teresa L, and Zhongming Z. Weak sharing of genetic association signals in three lung cancer subtypes: evidence at the SNP, gene, regulation, and pathway levels. Genome Med. (2018) 10. doi: 10.1186/s13073-018-0522-9

PubMed Abstract | Crossref Full Text | Google Scholar

16. Erinn B R and Amato J G. Hypoxic control of metastasis. Science. (2016) 352. doi: 10.1126/science.aaf4405

PubMed Abstract | Crossref Full Text | Google Scholar

17. Min-Zu W, Wei-Chung C, Su-Feng C, Shin N, Carolyn OC, Chia-Lin L, et al. miR-25/93 mediates hypoxia-induced immunosuppression by repressing cGAS. Nat Cell Biol. (2017) 19. doi: 10.1038/ncb3615

PubMed Abstract | Crossref Full Text | Google Scholar

18. Bo H, Kira H W, Zhi-Jie L, Judith A E, Ji L, Edward R B, et al. Selective tubulin-binding drugs induce pericyte phenotype switching and anti-cancer immunity. EMBO Mol Med. (2025). doi: 10.1038/s44321-025-00222-6

PubMed Abstract | Crossref Full Text | Google Scholar

19. Danyang C, Xiaohui C, Mingqiang Y, Qiunuo L, Shaojuan W, Siyue K, et al. H3K36me2 methyltransferase NSD2/WHSC1 promotes triple-negative breast cancer metastasis via activation of ULK1-dependent autophagy. Autophagy. (2025). doi: 10.1080/15548627.2025.2479995

PubMed Abstract | Crossref Full Text | Google Scholar

20. Jialin Z, Xiaoyu X, Ziwei Z, Kanghui G, Yi X, and Hualei D. Magea13 attenuates myocardial injury in acute myocardial infarction by inhibiting the cAMP-PKA signaling pathway. Apoptosis. (2025) 30. doi: 10.1007/s10495-025-02078-0

PubMed Abstract | Crossref Full Text | Google Scholar

21. Muhammad Zubair S, Ruyi H, Yingying H, Xin G, Yang L, Miao G, et al. Targeting TRAP1-dependent metabolic reprogramming to overcome doxorubicin resistance in quiescent breast cancer. Drug Resist Update. (2025) 81. doi: 10.1016/j.drup.2025.101226

PubMed Abstract | Crossref Full Text | Google Scholar

22. Na C, Zhanfeng L, Heyuan L, Aimin J, Liqiang Z, Siqi Y, et al. Enhancing PD-1 blockade in NSCLC: Reprogramming tumor immune microenvironment with albumin-bound statins targeting lipid rafts and mitochondrial respiration. Bioact Mater. (2025) 49. doi: 10.1016/j.bioactmat.2025.03.003

PubMed Abstract | Crossref Full Text | Google Scholar

23. Seung Hee S, Ji Hae L, Eun Kyung C, Seung Bae R, and Kyungsil Y. C/EBPβ Regulates HIF-1α-driven invasion of non-small-cell lung cancer cells. Biomolecules. (2025) 15. doi: 10.3390/biom15010036

PubMed Abstract | Crossref Full Text | Google Scholar

24. Tongxuan S, Ziqi J, Jiayi L, Heng C, Hengyi X, Lin C, et al. Unraveling the triad of hypoxia, cancer cell stemness, and drug resistance. J Hematol Oncol. (2025) 18. doi: 10.1186/s13045-025-01684-4

PubMed Abstract | Crossref Full Text | Google Scholar

25. Guangchuang Y. Thirteen years of clusterProfiler. Innovation (Camb). (2024) 5. doi: 10.1016/j.xinn.2024.100722

PubMed Abstract | Crossref Full Text | Google Scholar

26. Aaron M N, Chih Long L, Michael R G, Andrew J G, Weiguo F, Yue X, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. (2015) 12. doi: 10.1038/nmeth.3337

PubMed Abstract | Crossref Full Text | Google Scholar

27. Etienne B, Nicolas A G, Laetitia L, Bénédicte B, Nabila E, Florent P, et al. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol. (2016) 17. doi: 10.1186/s13059-016-1070-5

PubMed Abstract | Crossref Full Text | Google Scholar

28. Langfelder P and Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. (2008) 9:559. doi: 10.1186/1471-2105-9-559

PubMed Abstract | Crossref Full Text | Google Scholar

29. Chengran X, Jinhai H, Yi Y, Lun L, and Guangyu L. Increased expression of homeobox 5 predicts poor prognosis: A potential prognostic biomarker for glioma. Int J Gen Med. (2022) 15. doi: 10.2147/ijgm.S350454

PubMed Abstract | Crossref Full Text | Google Scholar

30. Kenneth TJ, Balasubramanian N, and Trevor H. Elastic net regularization paths for all generalized linear models. J Stat Softw. (2023) 106. doi: 10.18637/jss.v106.i01

PubMed Abstract | Crossref Full Text | Google Scholar

31. Óscar L-S, v.G. Maisa H, Peter AJ, Francesca F, and Federica E. Interpreta ble systems biomarkers predict response to immune-checkpoint inhibitors. Patterns (N Y). (2021) 2. doi: 10.1016/j.patter.2021.100293

PubMed Abstract | Crossref Full Text | Google Scholar

32. Bhuvic P, Sangami P, Collin W E, Vijay N, Shree S P, Tatenda M, et al. NF2 loss-of-function and hypoxia drive radiation resistance in grade 2 meningiomas. J Natl Cancer Inst. (2025) 117. doi: 10.1093/jnci/djaf022

PubMed Abstract | Crossref Full Text | Google Scholar

33. Frank S L. Under (Genetic selection) pressure: human tumors and human populations in hypoxia. Cancer Discov. (2025) 15. doi: 10.1158/2159-8290.Cd-25-0175

PubMed Abstract | Crossref Full Text | Google Scholar

34. Khalia R P, T. Joanne TM, S. Lauren N, Victoria A, Liam G S, Emma L S, et al. Reconstituted high-density lipoproteins rescue diabetes-impaired endothelial cell metabolic reprograming and angiogenic responses to hypoxia. Arterioscler Thromb Vasc Biol. (2025). doi: 10.1161/atvbaha.124.320110

PubMed Abstract | Crossref Full Text | Google Scholar

35. Ming B, Pengzhi X, Rong C, Na L, Sunrun C, Qiqiang G, et al. ROS-ATM-CHK2 axis stabilizes HIF-1α and promotes tumor angiogenesis in hypoxic microenvironment. Oncogene. (2025). doi: 10.1038/s41388-025-03336-w

PubMed Abstract | Crossref Full Text | Google Scholar

36. Xiaoqing C, Xing Y, Yu W, Pengtao S, Ruichao H, Chen Y, et al. Hypoxia-mimicking microenvironment scaffold for enhanced tendon regeneration. ACS Appl Mater Interfaces. (2025) 17. doi: 10.1021/acsami.4c18082

PubMed Abstract | Crossref Full Text | Google Scholar

37. Mengjia S, Junyi H, Qiu-Zhong P, Jieying Y, Jingjing Z, Yao-Jun Z, et al. Cancer-associated fibroblast-mediated cellular crosstalk supports hepatocellular carcinoma progression. Hepatology. (2021) 73. doi: 10.1002/hep.31792

PubMed Abstract | Crossref Full Text | Google Scholar

38. Gu-Wei J, Zheng-Gang X, Shuo-Chen L, Shu-Ya C, Chen-Yu J, Ming L, et al. Radiogenomics of intrahepatic cholangiocarcinoma predicts immunochemotherapy response and identifies therapeutic target. Clin Mol Hepatol. (2025). doi: 10.3350/cmh.2024.0895

PubMed Abstract | Crossref Full Text | Google Scholar

39. Jaime C-G, Nadia H, Magda C, María Ángeles S, Aitziber U, Antonio L-P, et al. Prognostic impact of tumor location and gene expression profile in sporadic desmoid tumor. Eur J Cancer. (2024) 209. doi: 10.1016/j.ejca.2024.114270

PubMed Abstract | Crossref Full Text | Google Scholar

40. Jianming W, Chun-Yuan C, Xue Y, Fan Z, Juan L, Zhaohui F, et al. Leukemia inhibitory factor, a double-edged sword with therapeutic implications in human diseases. Mol Ther. (2022) 31. doi: 10.1016/j.ymthe.2022.12.016

PubMed Abstract | Crossref Full Text | Google Scholar