A prognostic gene signature derived from aging‑related genes predicts survival, immune landscape and therapy response in glioma

Introduction

Gliomas are among the most common and aggressive primary brain tumors, with glioblastoma (GBM) being the most malignant subtype (1). Despite advancements in surgical and therapeutic interventions, GBM prognosis remains poor, with a median survival of ~15 months (2). Low-grade gliomas (LGGs) progress more slowly but often transforme into higher-grade tumors. The heterogeneity and molecular complexity of gliomas pose significant challenges to effective treatment and prognosis (3,4).

Aging is a well-recognized risk factor for numerous cancers, including gliomas. Aging induces cellular and molecular changes that create a favorable environment for tumor development and progression (5). One of the most critical changes is genomic instability, accompanied by significant epigenetic alterations, such as aberrant DNA methylation and histone modifications, which dysregulate gene expression and promote malignancy (6,7).

The tumor microenvironment (TME) also undergoes remodeling with aging, involving changes in the extracellular matrix (ECM), blood vessels and immune cell infiltration, which collectively enhance tumor growth and metastasis (8). In gliomas, these aging-related changes are particularly relevant. Dysregulation of pathways such as the p53 tumor suppressor and PI3K/AKT/mTOR, crucial for cell cycle regulation and survival, is common in gliomas and is exacerbated by aging (9). Moreover, mutations in IDH1 and IDH2, more prevalent in younger patients and associated with improved outcomes, highlight the complex interplay between aging and glioma biology (3,10).

The specific role of aging-related genes in glioma prognosis and their utility in prognostic models remain underexplored. In the present study, an aging-related risk score model that effectively predicts glioma patient prognosis, immunotherapy efficacy, and chemotherapy sensitivity, was developed. However, as the current findings are based on bioinformatics analysis, further experimental validation is required to confirm the model's clinical relevance.

Materials and methods

The present study utilized data from The Cancer Genome Atlas (TCGA) (TCGA-GBM and TCGA-LGG; https://portal.gdc.cancer.gov/) and the Chinese Glioma Genome Atlas (CGGA) (CGGA-mRNAseq_693; http://www.cgga.org.cn/) to analyze gene expression differences between LGGs and GBM. Gene expression profiles and corresponding clinical information were downloaded from these databases. The data acquisition and preprocessing steps are as follows:

Data acquisition and preprocessing

Glioma RNA sequencing and clinical data were obtained from the TCGA website (https://portal.gdc.cancer.gov). Glioma expression data and clinical information were downloaded from the CGGA website (http://www.cgga.org.cn/). Samples with missing or incomplete clinical information were excluded. The gene expression data were normalized and background-corrected via the ‘limma’ (v 3.58.1) package in R.

Differential gene expression analysis

To identify genes differentially expressed between LGG and GBM, TCGA and CGGA expression data were normalized and log2-transformed using the ‘limma’ package in R. |LogFC|>1 was selected and P<0.05 was adjusted to capture biologically meaningful differences while minimizing false positives, consistent with widely accepted thresholds in transcriptome studies.

Senescence-related gene screening

A total of 279 senescence-related genes were downloaded from the Human Aging Genomic Resources (HAGR) database and used Venn diagrams to identify the overlap between these senescence-related genes and the differentially expressed genes (DEGs) from the TCGA and CGGA datasets.

Risk score model construction and validation

From the DEGs, 29 senescence-related genes were selected to construct a risk score model. Univariate Cox regression analysis (P<0.001) identified significant genes, followed by LASSO regression to narrow the selection to 12 genes. Multivariate Cox regression finalized 8 genes for the model. The risk score formula was defined as: risk score=Σ (gene expression x gene coefficient), stratifying patients into high-risk and low-risk groups.

The model's performance was assessed using Kaplan-Meier (K-M) survival and receiver operating characteristic (ROC) curve analyses in TCGA and CGGA datasets. R packages ‘survival’ (v3.2.1) and ‘survminer’ (v3.3.3) were utilized, with P<0.05 considered significant. Area Under Curve (AUC) values >0.9 indicated excellent performance, >0.8 favorable performance, and >0.7 useful discrimination.

Construction and validation of the nomogram

Univariate and multivariate Cox regression analyses were performed to identify potential prognostic factors for overall survival (OS) in patients with glioma. A nomogram was constructed to predict the 1-, 3-, and 5-year OS probabilities. Calibration plots and ROC curves were used to evaluate the nomogram's performance. All analyses were performed via the ‘survival (v3.2.1)’, ‘ggplot2 (v3.3.3)’, and ‘rms (v6.2.0)’ R packages, with P<0.05 considered significant.

Functional enrichment analysis and gene set variation analysis (GSVA)

Using the STRING database (https://string-db.org/) (score ≥0.4, max 50 interactors), 110 genes interacting with the 8 senescence risk score genes were identified. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed using the ‘ggplot2’ (v 3.5.2) and ‘clusterProfiler’ (v 4.10.1) R packages.

GSVA, an unsupervised method, evaluated gene set enrichment based on mRNA expression data. Gene sets from the Molecular Signatures Database (v7.0) (https://www.gsea-msigdb.org/gsea/msigdb) were used to assess functional differences between high- and low-risk groups. Heatmaps were visualized with the ‘ComplexHeatmap’ R package.

Immune microenvironment assessment

The CIBERSORT algorithm estimated differences in immune cell infiltration between high- and low-risk groups, while the ESTIMATE method assessed relationships between the risk score, stromal score, immune score, ESTIMATE score, and tumor purity. CIBERSORT results were analyzed using weighted gene co-expression network analysis (WGCNA). The ‘pickSoftThreshold’ function calculated soft threshold power, and the adjacency matrix was converted into a topological overlap matrix for hierarchical clustering. Co-expressed gene modules were identified using dynamic tree-cutting, and key module genes underwent GO and KEGG enrichment analyses.

Drug sensitivity analysis

Drug sensitivity in high- and low-risk groups was assessed using data from the Genomics of Drug Sensitivity in Cancer (www.cancerrxgene.org) database, combined with gene expression data. Pearson correlation analysis evaluated the relationship between gene expression and drug response.

scRNA-seq data processing and analysis

Raw 10X scRNA-seq data from GSE162631(11) were downloaded from Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/) and processed as follows: i) data were converted into a Seurat object using the ‘Seurat’ R package (12,13); ii) quality control was performed by excluding cells with <200 or >5,000 detected genes and those with >15% mitochondrial content; iii) the ‘FindVariableFeatures’ function identified the top 2,000 highly variable genes; iv) principal component analysis (PCA) and UMAP were applied for dimensionality reduction and clustering (14); v) the ‘FindMarkers’ function identified marker genes for each cluster (fold change >1); vi) clusters were annotated using the ‘SingleR’ package to define cell types (15); and vii) M1 macrophages were identified via markers NOS2 and CD80, while M2 macrophages were marked by CD163, MRC1, IL10, IL6 and ARG1. DEGs between M1 and M2 macrophages were identified using the ‘FindMarkers’ function.

Cell lines and reverse transcription-quantitative (RT-qPCR)

Human astrocytes (HA; cat. no. 1800; ScienCell Research Laboratories), U251 GBM cells (CVCL_0021; Cell Bank of the Chinese Academy of Sciences), and U87 MG GBM of unknown origin (CVCL_0022; American Type Culture Collection) were cultured following the suppliers' protocols. Cell line authentication was not additionally performed in the present study, which represents a limitation.

Total RNA was extracted using TRIzol™ Reagent (Invitrogen; Thermo Fisher Scientific, Inc.). Reverse transcription was performed with the PrimeScript™ RT Reagent Kit (Takara Bio, Inc.) according to the manufacturer's instructions. qPCR was performed using PowerUp™ SYBR™ Green Master Mix (Applied Biosystems; Thermo Fisher Scientific, Inc.) on an ABI 7500 Real-Time PCR System (Applied Biosystems; Thermo Fisher Scientific, Inc.). The thermal cycling program consisted of an initial denaturation at 95˚C for 2 min, followed by 40 cycles of 95˚C for 15 s and 60˚C for 1 min, with a final melt-curve stage to verify amplification specificity. Primer sequences are provided in Table SII. GAPDH was used as the internal control, and relative expression levels were calculated using the 2-ΔΔCq method (16).

Statistical analysis

All the statistical analyses were conducted via R software (version 4.0.5; https://www.r-project.org/). P<0.05 was considered to indicate a statistically significant difference.

Results

Screening of DEGs and model construction

In the TCGA dataset, 12,178 DEGs were identified (7,604 upregulated, 4,574 downregulated) (|LogFC|>1, adjusted P<0.05) (Fig. 1A), while 1,890 DEGs were identified in the CGGA dataset (930 upregulated, 960 downregulated) (Fig. 1B). A total of 279 senescence-related genes were obtained from the HAGR. By intersecting DEGs from TCGA and CGGA with these genes, 29 overlapping genes were identified (Fig. 1C).

Figure 1 Screening of DEGs and model construction. (A) Volcano plot of significantly upregulated (red) and downregulated (blue) genes in the TCGA dataset (|LogFC|>1, adjusted P<0.05). (B) Volcano plot of DEGs in the CGGA dataset (|LogFC|>1, adjusted P<0.05). (C) Venn diagram showing 29 shared senescence-related genes from TCGA and CGGA datasets. (D) Heatmap of hierarchical clustering of 29 genes, separating LGG and GBM samples. (E) Partial likelihood deviance plot from LASSO regression identifying 12 genes at the optimal lambda. (F) Coefficient profiles of 12 genes across lambda values in LASSO analysis. (G) Forest plot from multivariate Cox regression of 8 genes significantly associated with prognosis, displaying HRs and P-values. DEGs, differentially expressed genes; TCGA, The Cancer Genome Atlas; CGGA, Chinese Glioma Genome Atlas; LGG, low-grade glioma; GBM, glioblastoma; HR, hazard ratio; CI, confidence interval.

To evaluate the prognostic potential of these genes, univariate Cox regression was used (Table SI). Hierarchical clustering separated samples into GBM-dominated (23 upregulated genes) and LGG-dominated (6 downregulated genes) groups (Fig. 1D). Using LASSO regression, 12 genes were selected by minimizing partial likelihood deviance (Fig. 1E), with the coefficient paths shown in Fig. 1F. Multivariate Cox regression further identified 8 of these genes as independently prognostic (Fig. 1G). Based on these 8 genes, a risk score model was constructed.

Evaluation of the prognostic predictive ability of the risk score model

The prognostic predictive ability of the 8-gene risk score model was further evaluated in the TCGA dataset. K-M analysis demonstrated that high-risk patients had significantly lower survival rates compared with low-risk patients in the training set (HR=6.84, 95% CI=4.44-10.53, P<0.001; median survival: 16.6 vs. 95.5 months; Fig. 2A), test set (HR=6.86, 95% CI=3.46-13.60, P<0.001; 11.9 vs. 50.5 months, Fig. 2B), and the entire dataset (HR=6.676, 95% CI=4.70-9.73, P<0.001; 17.7 vs. 95.5 months, Fig. 2C). ROC curves confirmed strong model performance, with AUC values of 0.877, 0.891 and 0.831 for 1-, 3- and 5-year predictions in the training set (Fig. 2D), 0.847, 0.932 and 0.886 in the test set (Fig. 2E), and 0.870, 0.904 and 0.853 in the entire dataset (Fig. 2F). These results indicated high sensitivity and specificity in predicting patient survival at various time points. Risk score distribution and survival analysis revealed clear stratification between high- and low-risk groups. High-risk patients, primarily distributed on the right side, showed shorter survival times and higher mortality in the training (Fig. 2G) and test sets (Fig. 2H). Similar trends in the entire dataset validated the model's generalizability and effectiveness (Fig. 2I).

Figure 2 Evaluation of the prognostic predictive ability of the 8-gene risk score model. (A-C) Kaplan-Meier survival curves showing significantly lower survival in the high-risk group across the (A) training set, (B) test set and (C) entire dataset (P<0.001). (D-F) Receiver operating characteristic curves for 1-, 3- and 5-year survival predictions in the (D) training set, (E) test set and (F) entire dataset, demonstrating high predictive accuracy. (G-I) Risk score distribution and survival status plots illustrating shorter survival and higher mortality in high-risk patients in the (G) training set, (H) test set and (I) entire dataset.

Validation of the risk score model in the CGGA dataset

The prognostic performance of the senescence-related risk score model was validated in the CGGA dataset. K-M analysis showed significantly lower survival rates in the high-risk group compared with the low-risk group (HR=3.72, 95% CI=3.00-4.60, P<0.001; 15.6 vs. 94.4 months) (Fig. 3A). The median survival time was notably shorter in the high-risk group, confirming the model's effectiveness. ROC analysis demonstrated AUC values of 0.750, 0.788 and 0.785 for 1-, 3- and 5-year prognoses, respectively, indicating favorable sensitivity and specificity (Fig. 3B). The risk score distribution revealed a clear survival time distinction, with early-stage deaths more frequent in the high-risk group, while the low-risk group showed longer survival (Fig. 3C). A heatmap highlighted differential expression of the 8 key genes, with ribonucleotide reductase regulatory subunit (RRM2) and WEE1 upregulated in the high-risk group and glutamate ionotropic receptor kainate type subunit 2 (GRIK2) and growth arrest and DNA damage inducible gamma (GADD45G) downregulated (Fig. 3C).

Figure 3 Validation of the 8-gene risk score model in the CGGA dataset. (A) Kaplan-Meier survival curves showing lower survival in the high-risk group (Hazard ratio=3.72; P<0.001). (B) Receiver operating characteristic curves for 1-, 3-, and 5-year survival predictions with area under the curve values of 0.750, 0.788 and 0.785, indicating favorable performance. (C) Risk score distribution and survival status plot illustrating shorter survival and higher mortality in the high-risk group. CGGA, Chinese Glioma Genome Atlas.

Prognostic ability of the risk score model in different clinical and molecular subgroups

The prognostic predictive ability of the senescence-related risk score model was assessed across clinical and molecular subgroups (Fig. 4A-R).

Figure 4 Prognostic ability of the risk score model in clinical and molecular subgroups. (A-C) Kaplan-Meier curves for age groups under 40 years (A), 40-60 years (B), and over 60 years (C), showing worse survival in high-risk patients. (D and E) Survival stratification by sex, with high-risk groups having worse outcomes for both males and females. (F-H) Stratification by tumor grade (G2-G4), where high-risk patients consistently had shorter survival. (I-L) Stratification by histological type, showing significant differences for astrocytoma (I), oligoastrocytoma (K) and glioblastoma (L), but not for oligodendroglioma (J). (M and N) (M and N) Kaplan-Meier curves for IDH wild-type (M) and IDH-mutant (N) patients, revealing worse survival in the high-risk group. (O and P) Stratification by 1p/19q status, with significant differences in patients with (O) and without (P) codeletion. (Q and R) Stratification by the MGMT methylation status, with worse survival in the high-risk group for both methylated (Q) and unmethylated (R) patients.

Age stratification. The model demonstrated strong discriminatory power across all age groups. Low-risk patients had significantly improved survival compared with high-risk patients, with HR values of 4.04 (95% CI=1.74-9.40, P<0.001) for all ages (Fig. 4A), 10.19 (95% CI=5.65-20.92, P<0.001) for ages 40-60 (Fig. 4B), and 2.41 (95% CI=1.23-4.70, P=0.010) for patients over 60 (Fig. 4C). Regarding sex stratification, both male and female low-risk patients exhibited significantly improved survival than high-risk patients (P<0.001, Fig. 4D and E).

Tumor grade stratification. High-risk patients had consistently shorter survival across tumor grades G2 to G4 (Fig. 4F-H).

Histological type stratification. Low-risk patients with astrocytoma, oligoastrocytoma and GBM had improved survival compared with high-risk patients (Fig. 4I, K and L). For oligodendroglioma, the model showed limited discriminatory ability, with a non-significant difference (HR=2.55, 95% CI=0.91-7.17, P=0.076; Fig. 4J).

Molecular characteristic stratification. IDH status: High-risk patients had worse survival in both IDH wild-type (HR=6.77, 95% CI=3.35-13.67, P<0.001, Fig. 4M) and IDH-mutant groups (HR=2.53, 95% CI=1.40-4.56, P=0.002, Fig. 4N).

1p/19q status. Low-risk patients with 1p/19q codeletion showed improved survival (HR=4.21, 95% CI=1.27-13.89, P=0.018, Fig. 4O), while high-risk patients without codeletion had significantly poorer outcomes (HR=3.62; 95% CI=1.52-7.25; P<0.001; Fig. 4P).

MGMT methylation status. Low-risk patients with MGMT methylation had improved survival (HR=5.77, 95% CI=3.36-9.91, P<0.001; Fig. 4Q), with significant differences also observed in unmethylated patients (HR=4.32, 95% CI=2.05-9.10, P<0.001; Fig. 4R).

Construction and evaluation of the nomogram

Univariate and multivariate Cox regression analyses were performed, incorporating the senescence-related risk score and clinicopathological factors (Fig. 5A). Multivariate analysis identified tumor grade, age, IDH status and risk score as independent prognostic factors. G3 (HR=2.153, P=0.004) and G4 (HR=4.418, P<0.001) grades, along with ages 40-60 (HR=2.943, P<0.001) and over 60 (HR=5.623, P<0.001), were associated with poorer survival. IDH wild-type (HR=2.017, P=0.036) and high-risk scores (HR=1.824, P=0.048) also indicated increased mortality. A prognostic nomogram was developed to predict 1-, 3-, and 5-year survival probabilities (Fig. 5B). G4 grade and age over 60 contributed the highest points, with total scores mapped to survival probabilities. The nomogram was validated using the C-index (0.853, P<2x10-16), calibration plots and ROC curves. Calibration plots closely matched the diagonal line (Fig. 5C). ROC curves showed AUCs of 0.882, 0.934 and 0.902 for 1-, 3- and 5-year predictions, confirming excellent accuracy (Fig. 5D).

Figure 5 Construction and evaluation of the nomogram. (A) Multivariate Cox regression analysis identifying tumor grade, age, IDH status, and risk score as independent prognostic factors. (B) Nomogram predicting 1-, 3- and 5-year survival probabilities based on clinicopathological factors and risk score. (C) Calibration plot showing strong alignment between predicted and actual survival outcomes. (D) Time-dependent receiver operating characteristic curves with area under the curve values of 0.882, 0.934 and 0.902 for 1-, 3- and 5-year survival predictions, demonstrating high model accuracy.

Biological function differences between patients with high-risk and low-risk glioma

Using the STRING database, 110 genes interacting with the 8 senescence-related risk score genes were identified (threshold: score ≥0.4, interactors ≤50). GO analysis revealed enrichment in processes like mitotic cell cycle transition, nuclear division, spindle assembly and DNA binding, while KEGG analysis highlighted pathways such as the cell cycle, senescence, p53 signaling and FoxO signaling (Fig. 6A and B). GSVA showed that high-risk groups had upregulated pathways, including KRAS signaling, MYC targets, G2M checkpoint, and glycolysis, whereas low-risk groups exhibited upregulation in pathways like oxidative phosphorylation, DNA repair and apical junctions, reflecting differences in cell cycle, metabolism and immune responses (Fig. 6C).

Figure 6 Biological Function Differences between patients with high-risk and low-risk glioma. (A) Gene Ontology enrichment analysis showing involvement in the cell cycle, nuclear division, and protein kinase activity. (B) Kyoto Encyclopedia of Genes and Genomes pathway analysis highlighting enrichment in the cell cycle, senescence, p53 and FoxO signaling pathways. (C) Gene set variation analysis showing the upregulation of KRAS signaling and the inflammatory response in the high-risk group and oxidative phosphorylation and DNA repair in the low-risk group.

Evaluation of the immune microenvironment in patients with high-risk and low-risk glioma

The senescence-related risk score was analyzed with the immune microenvironment using CIBERSORT. In the high-risk group, M2 macrophages (31.11%) and monocytes (18.79%) were most abundant, while monocytes (25.19%) and M2 macrophages (24.72%) dominated in the low-risk group (Fig. 7A). High-risk tumors showed greater infiltration of immune cells, including naïve B cells, CD8+ T cells, memory CD4+ T cells, Tregs, gamma delta T cells, and macrophages (M0, M1 and M2). Low-risk tumors had higher levels of memory B cells, naïve CD4+ T cells, activated natural killer (NK) cells, monocytes, mast cells, eosinophils and neutrophils (Fig. 7B). Immune, stromal and ESTIMATE scores were higher in the high-risk group (P<0.001), while tumor purity was lower (P<0.001) (Fig. 7C-F). Greater infiltration of M2 macrophages and neutrophils in high-risk patients correlated with worse survival (P=0.003), while higher eosinophil levels indicated improved survival (P=0.025). In low-risk patients, M2 macrophage infiltration was linked to poorer survival (P=0.001) (Fig. 7G-J).

Figure 7 Immune microenvironment in patients with high-risk and low-risk glioma. (A) Proportions of immune cells, with M2 macrophages and monocytes being the most abundant in both groups. (B) High-risk group showed increased CD8+ T cells, Tregs and macrophage infiltration, while low-risk group had higher memory B cells, monocytes and neutrophils. (C-F) High-risk group exhibited higher ImmuneScore and StromalScore, but lower tumor purity. (G-J) Survival analysis showed worse outcomes with greater M2 macrophage and neutrophil infiltration, and improved outcomes with higher eosinophil infiltration (*P<0.05, **P<0.01, ***P<0.001).

WGCNA analysis of CIBERSORT results and functional enrichment analysis

WGCNA was performed on CIBERSORT results to identify gene modules associated with immune cell infiltration. At a soft threshold of 12, the R² value stabilized near 1, indicating scale-free topology with sparse yet moderate connectivity (Fig. 8A). The gene dendrogram identified 20 co-expression modules (Fig. 8B), and correlation analysis revealed functional roles of specific modules in the immune microenvironment (Fig. 8C). The MEblue module was positively correlated with immunosuppressive cells (M2 macrophages, Tregs) and negatively correlated with antitumor cells (activated NK cells), suggesting its role in promoting an immunosuppressive environment.

Figure 8 Weighted gene co-expression network analysis and functional enrichment. (A) Soft threshold at 12 ensures scale-free network topology. (B) Gene dendrogram identifying 20 co-expression modules. (C) Correlation heatmap showing that the MEblue module is positively associated with M2 macrophages and Tregs and negatively associated with activated NK cells. (D) The risk score was positively correlated with the MEblue module (R=0.718, ***P<0.001). (E) Gene Ontology analysis of the MEblue module highlights enrichment in extracellular matrix organization and cell adhesion. (F) Kyoto Encyclopedia of Genes and Genomes analysis revealed enrichment of the PI3K-Akt, focal adhesion and MAPK signaling pathways.

A strong positive correlation was observed between the risk score and the MEblue module (R=0.718, P<0.001, Fig. 8D). Functional enrichment analysis of MEblue genes revealed GO terms related to ECM organization, cell adhesion and actin binding (Fig. 8E). KEGG analysis highlighted enrichment in PI3K-Akt signaling, proteoglycans in cancer, focal adhesion and MAPK signaling pathways (Fig. 8F).

Association of the risk score with immune checkpoint molecules, immunotherapy and chemotherapy sensitivity

Immune checkpoint gene expression and its association with immune evasion, suppression and the immune microenvironment were analyzed in patients with glioma. Immune checkpoint genes, including PDCD1 (PD-1), CD274 (PD-L1), PDCD1LG2 (PD-L2), CTLA4, HAVCR2 (TIM-3), LAG3, CEACAM1 and IDO1, were significantly upregulated in high-risk patients compared with low-risk patients (P<0.001, Fig. 9A). The risk score demonstrated a strong positive correlation with these molecules (P<0.05, Fig. 9B).

Figure 9 Associations of the risk score with immune checkpoints, immunotherapy and chemotherapy sensitivity. (A) Increased expression of immune checkpoint molecules (PD-1, PD-L1 and CTLA4) in the high-risk group. (B) Positive correlation between the risk score and immune checkpoint expression. (C-F) High-risk patients have higher TIDE, immune exclusion and microsatellite instability scores, indicating stronger immune evasion. (G-K) High-risk patients have reduced sensitivity to chemotherapy, with higher IC50 values for multiple drugs. (L-N) Higher risk scores are linked to IDH wild-type, 1p/19q non-codeletion and MGMT unmethylated status (***P<0.001).

Using the TIDE algorithm, high-risk patients exhibited higher immune evasion and exclusion scores (P<0.001, Fig. 9C and D), while dysfunction scores did not differ significantly (P>0.05, Fig. 9E). Higher microsatellite instability scores were observed in the high-risk group (P<0.001, Fig. 9F).

Chemotherapy sensitivity analysis revealed reduced sensitivity in high-risk patients to temozolomide, carmustine, vincristine, cisplatin and palbociclib, as indicated by higher IC50 values (P<0.001, Fig. 9G-K). High-risk patients with wild-type IDH, 1p/19q non-codeletion, and unmethylated MGMT had significantly higher risk scores (P<0.001, Fig. 9L-N), explaining their lower drug sensitivity or resistance.

Single-cell differential expression analysis of GADD45G in M1 and M2 macrophages

Single-cell analysis was performed on four GBM samples from GSE162631. After filtering cells with mitochondrial gene expression >15%, 55,757 high-quality cells were retained (Fig. S1A). The ‘FindVariableFeatures’ function identified the top 2,000 most variable genes (Fig. S1B). Dimensionality reduction using PCA and UMAP resulted in 17 cell clusters (Fig. 10A), which were annotated with ‘SingleR’ to identify subpopulations, including monocytes, macrophages, endothelial cells, neutrophils and NK cells (Fig. 10B). Expression analysis of the eight model genes showed retinol-binding protein 1(RBP1) was highly expressed in endothelial cells, while GADD45G was expressed in both macrophages and endothelial cells (Fig. 10C). UMAP further confirmed GADD45G's predominant expression in macrophages (Fig. 10D). M1 and M2 macrophages were manually annotated (Fig. 10E), and differential expression analysis revealed significantly higher GADD45G expression in M2 macrophages compared with M1 macrophages (Fig. 10F).

Figure 10 Single-cell differential expression analysis of GADD45G in M1 and M2 macrophages. (A) UMAP plot showing 17 cell clusters identified from glioblastoma single-cell data. (B) Cell type annotation via ‘SingleR’, identifying subpopulations such as monocytes, macrophages, endothelial cells, neutrophils and natural killer cells. (C) Expression distribution of eight model genes across subpopulations, with GADD45G highly expressed in macrophages and endothelial cells. (D) UMAP plot showing predominant GADD45G expression in macrophages. (E) Manual annotation of M1 and M2 macrophages. (F) GADD45G expression is significantly greater in M2 macrophages than in M1 macrophages. GADD45G, growth arrest and DNA damage inducible gamma (***P<0.001).

Validation of 8-gene expression by RT-qPCR in normal and GBM cells

RT-qPCR was performed to validate the expression of eight genes [netrin-4 (NTN4), RBP1, Twist Family BHLH Transcription Factor 1 (TWIST1), GADD45G, NUAK2, GRIK2, WEE1 and RRM2] in normal astrocytes (HA) and GBM cell lines (U251 and U87).

Regarding NTN4, significantly higher expression was observed in HA compared with U251 and U87 (Fig. 11A; P<0.05). RBP1 was significantly elevated in U251 compared with HA (P<0.05; Fig. 11B). TWIST1 was strongly upregulated in U251 and U87 compared with HA (P<0.05; Fig. 11C). GADD45G demonstrated higher expression in HA compared with U251, but differences between HA and U87 were not significant (Fig. 11D). NUAK2 was increased in HA compared with U251 (P<0.05), while no significant difference was observed between HA and U87 (Fig. 11E). GRIK2 exhibited significantly higher expression in U251 and U87 compared with HA (P<0.05; (Fig. 11F). WEE1 revealed increased expression in both U251 and U87 compared with HA (P<0.05; (Fig. 11G). Regarding RRM2, no significant differences were observed across HA, U251 and U87 groups (Fig. 11H).

Figure 11 RT-qPCR validation of 8-gene expression in normal astrocytes (HA) and glioblastoma cell lines (U251 and U87). (A-H) RT-qPCR analysis of eight genes: (A) NTN4, (B) RBP1, (C) TWIST1, (D) GADD45G, (E) NUAK2, (F) GRIK2, (G) WEE1 and (H) RRM2. *P<0.05. RT-qPCR, reverse transcription-quantitative PCR; NTN4, netrin-4; RBP1, retinol-binding protein 1; TWIST1, Twist Family BHLH Transcription Factor 1; GADD45G, growth arrest and DNA damage inducible gamma; RRM2, ribonucleotide reductase regulatory subunit; GRIK2, glutamate ionotropic receptor kainate type subunit 2; ns, not significant.

Discussion

Aging-related genes play critical roles in tumors by regulating the cell cycle, DNA repair, metabolism and immune evasion (6,17). In the present study, 29 aging-related genes critical to GBM progression were identified and an eight-gene risk score model was developed with strong prognostic and discriminative performance in the TCGA and CGGA datasets, demonstrating its broad applicability and reliability. Unlike prior prognostic signatures based on senescence or immune-related genes, our model integrates aging-related genes curated from HAGR with multi-cohort transcriptomic validation, scRNA-seq profiling and RT-qPCR assays and extends to immune microenvironment and drug resistance analyses. These genes participate in diverse biological processes relevant to glioma. NTN4 is associated with poor prognosis, likely due to its role in ECM regulation and tumor invasion (18); RBP1 may affect tumor differentiation and survival via the retinoic acid signaling pathway (19); TWIST1 plays a pivotal role in cancer progression, primarily by driving epithelial-mesenchymal transition, which in turn enhances the invasive and metastatic potential of tumor cells (20,21); Accumulating evidence suggests that GADD45G proteins participate in the regulation of tumor progression, including that of gliomas, particularly under conditions of oncogenic stress (22-24); NUAK2 is involved in neural tube development, and its aberrant expression is strongly associated with tumor progression and poor patient prognosis, including in GBM (25-27); GRIK2 has been reported to participate in the development of gastric (28,29); however, its potential involvement in glioma has not yet been clearly elucidated; WEE1, a key cell cycle regulator, has shown therapeutic potential in GBM clinical trials (30); and RRM2 has been reported to facilitate the proliferation, migration and invasive behavior of human GBM cells, while simultaneously suppressing apoptosis (31). Collectively, these genes converge on pathways of cell cycle regulation, DNA damage response, invasion and metabolic adaptation, consistent with their reported functions in glioma biology. Their integration into a unified signature highlights their combined prognostic value.

The risk score model showed strong predictive power across clinical and molecular subgroups, demonstrating its utility in both overall and subgroup prognosis. High-risk patients with wild-type IDH and 1p/19q non-codeletion had significantly lower survival rates, consistent with these molecular features being key glioma prognostic factors. IDH mutations are linked to improved prognosis and treatment response, while 1p/19q codeletion enhances sensitivity to chemotherapy and radiotherapy (10). To improve clinical applicability, the risk score model was integrated with clinicopathological features to construct a nomogram, which exhibited high predictive accuracy and consistency. This validates the model's utility in providing valuable prognostic insights across diverse molecular subtypes and clinical contexts.

GO and KEGG analyses highlighted key biological functions and pathways related to the risk score model, including the cell cycle, DNA repair and metabolism. WEE1, a critical regulator of cell cycle checkpoints, prevents damaged cells from entering mitosis (32). KEGG analysis linked the model genes to pathways such as the cell cycle, cellular senescence, p53 signaling and FoxO signaling, which are central to tumorigenesis and progression. For instance, p53 suppresses glioma by regulating the cell cycle and inducing apoptosis (33,34). Cellular senescence, triggered by DNA damage, telomere shortening, or oncogenic stress, plays a dual role in glioma by suppressing tumor formation while potentially promoting progression via proinflammatory factor secretion (35).

GSVA showed significant upregulation of cancer-related pathways, including KRAS signaling, MYC targets and E2F targets, in the high-risk group. KRAS drives sustained proliferation and survival via the MAPK and PI3K/AKT pathways, contributing to poor prognosis and therapy resistance in GBM (36,37). MYC overexpression promotes malignant transformation and high proliferative capacity, with MYC-regulated genes playing critical roles in cancer progression (38,39). In the low-risk group, pathways related to oxidative phosphorylation, DNA repair, tight junctions and interferon responses were upregulated. Enhanced oxidative phosphorylation, indicating a reliance on mitochondrial energy production, is linked to reduced tumor proliferation and improved prognosis (40). Upregulated DNA repair-related genes promote genomic stability, reduce carcinogenesis risk, and improve chemotherapy response (41,42).

The TME plays a vital role in tumor growth, metastasis and treatment response, consisting of immune cells, fibroblasts, vascular cells and stromal components (43). In the high-risk group, elevated infiltration of M2 macrophages, an immunosuppressive subtype linked to poor prognosis, was observed. M2 macrophages promote glioma growth and immune evasion by secreting anti-inflammatory cytokines and stimulating angiogenesis (44). Regulatory T cells (Tregs), also abundant in this group, suppress antitumor responses by inhibiting effector T-cell and NK cell activity via cytokines like IL-10 and TGF-β (45). High-risk patients exhibited significantly higher immune, stromal and ESTIMATE scores (P<0.001), indicating a greater abundance of nontumor components such as fibroblasts and stromal cells, which promote tumor expansion and invasion by secreting growth factors and matrix-degrading enzymes (46). These components, along with immunosuppressive cells, create barriers that reduce the effectiveness of chemotherapy and radiotherapy, potentially explaining the poorer therapeutic responses in patients with high-risk GBM.

WGCNA analysis of CIBERSORT results identified the MEblue module as positively correlated with immunosuppressive cells (M2 macrophages and Tregs) and negatively correlated with antitumor immune cells (activated NK cells), suggesting an association in immune evasion and tumor progression. GO and KEGG analyses showed that MEblue genes regulate ECM organization and cell adhesion, which are known to support the immunosuppressive tumor environment. ECM remodeling, a key step in cancer progression, involves proteins and enzymes like collagen, matrix metalloproteinases, and integrins (47). In gliomas, ECM remodeling promotes tumor invasion and metastasis through cell-matrix interactions and focal adhesions. Integrins bind ECM components, such as fibronectin, and activate downstream pathways, including PI3K-Akt and MAPK, which regulate cell proliferation, survival and migration, driving GBM progression (48,49).

In the high-risk group, key immune checkpoint genes, including PDCD1 (PD-1), CD274 (PD-L1) and PDCD1LG2 (PD-L2), were significantly upregulated, consistent with an immunosuppressive state. PD-1 binding to PD-L1 or PD-L2 reduces T-cell proliferation and cytotoxicity, allowing tumor cells to evade immune surveillance (50). High CTLA-4 expression further supports immune evasion by inhibiting T-cell activation through competitive binding to B7 molecules (51). The high-risk group also exhibited significantly higher TIDE scores (P<0.001), consistent with stronger immune evasion potential and poorer responses to immune checkpoint inhibitors. Elevated immune exclusion scores in this group were associated with reduced T-cell infiltration and cytotoxic activity through physical and chemical barriers, aligning with findings of increased immune and stromal components in the high-risk group.

The glioma immune microenvironment is not static but evolves with tumor stage, progression and recurrence. GBM generally exhibits higher levels of M2 macrophages and exhausted T cells than LGG, and recurrent tumors often show stromal fibrosis and stronger immunosuppressive signaling (52-54). In the present study, the computational analysis could not capture these stage-specific dynamics, which should be addressed in future studies using longitudinal sampling, spatial transcriptomics and multiplex immunohistochemistry.

Drug sensitivity analysis showed that high-risk patients were less responsive to chemotherapeutic agents, including temozolomide, carmustine, vincristine, cisplatin and palbociclib, compared with low-risk patients. In the present study, GSVA and KEGG enrichment analyses indicated that the high-risk group was enriched in PI3K/AKT/mTOR, E2F target, and KRAS signaling pathways, which are known to promote tumor survival, inhibit apoptosis, and drive therapeutic resistance in glioma (37,55,56). Key genes in our signature, including RRM2 and WEE1, regulate DNA synthesis and G2/M checkpoint control, respectively, mechanisms that facilitate DNA repair and contribute to resistance to temozolomide and radiotherapy (57,58). Furthermore, immune deconvolution revealed that the high-risk subgroup exhibited increased infiltration of M2 macrophages and Tregs, which secrete immunosuppressive cytokines such as IL-10 and TGF-β, thereby limiting antitumor immunity and drug efficacy (59). These changes are consistent with the immunosuppressive remodeling and stromal fibrosis observed in advanced gliomas (60). By contrast, the low-risk group was enriched for IDH mutation, 1p/19q codeletion, and MGMT promoter methylation, well-established molecular features associated with improved chemosensitivity and prognosis (61-63). Taken together, these findings suggest that the diminished drug sensitivity in the high-risk group results from the combined effects of cell-cycle dysregulation (RRM2 and WEE1), PI3K/AKT-mediated survival signaling, and an immunosuppressive TME, consistent with previous mechanistic studies in glioma biology.

In summary, the present study highlights the critical role of aging-related genes in glioma and introduces a gene expression-based risk score model for predicting patient survival. These findings offer new insights and potential targets for personalized treatment and prognosis. However, this study has several limitations. Functional validation experiments (for example, knockdown/overexpression assays and in vivo models) were not performed, and drug sensitivity and immune analyses relied on computational inference without experimental confirmation. The scRNA-seq dataset was relatively small, and the RT-qPCR validation was conducted only in established glioma cell lines (U87, U251 and HA). Although these models are widely used and provide preliminary support for the expression patterns of the signature genes, they do not fully recapitulate the molecular heterogeneity and TME interactions observed in patient tissues. Long-term cultured cell lines may acquire genetic drift and altered signaling pathways, potentially affecting gene expression and treatment response. Therefore, future studies should include validation using patient-derived glioma tissues, organoid cultures, or xenograft models to better reflect the in vivo tumor context and confirm clinical relevance. In addition, grade-specific analyses for LGG and GBM were not conducted. These limitations warrant further validation in future studies.

Supplementary Material

Data processing and identification of variable genes. (A) Cells with >15% mitochondrial gene expression were filtered, leaving 55,757 high-quality cells. (B) The top 2,000 most variable genes identified via the ‘FindVariableFeatures’ function.

Each gene was evaluated for its association with overall survival in 698 patients with glioma.

Sequences of primers used for reverse transcription-quantitative PCR.

Acknowledgements

Not applicable.

Funding

Funding: No funding was received.

Availability of data and materials

The data generated in the present study may be found at the following URL: (portal.gdc.cancer.gov) (https://www.cgga.org.cn/download.jsp).

Authors' contributions

QZ conceptualized the study, and wrote, reviewed and edited the manuscript. ZL developed methodology and wrote the original draft. LZ and PL performed data visualization. XLi and XLiang conducted formal analysis. TW, JY and QY curated data. QZ and ZL confirm the authenticity of all the raw data. All authors read and approved the final version of the manuscript.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References