Introduction

The prevalence of cancer is a pressing concern in public health, as both the frequency of new cases and the number of deaths caused by the disease increase each year. The United States recorded a total of 1.9 million new cases of cancer and 609,360 deaths caused by cancer in 20221. Gliomas, characterized by their infiltrative growth pattern resembling “crab claws” within the brain, do not have a distinct boundary separating them from normal brain tissue. This characteristic complicates complete surgical resection, as residual tumor cells can easily lead to recurrence. Furthermore, gliomas have a regenerative and proliferative “leek-like” nature, with the degree of malignancy potentially increasing as the disease progresses, thereby complicating treatment and worsening the prognosis. Another difficulty arises from the intricate nature of the tumor microenvironment (TME) in patients, referred to as interpatient heterogeneity2. Given the unclear molecular pathways involved in the development and spread of tumors, coupled with the limited effectiveness of conventional treatments, it is imperative to clarify the subtypes of glioma to facilitate precision medicine.

The process of aging significantly increases the likelihood of developing a range of diseases, such as diabetes, cardiovascular diseases, neurological disorders, and cancer. The buildup and long-lasting presence of senescent cells in the body are characteristic signs of aging, leading to a steady decline in both regular and disease-related biological functions3. Recent studies have utilized multi-omics techniques to investigate the intricate connection between aging and cancer, given their close association. These studies have combined epigenomic, exomic, and transcriptome data from large-scale datasets to gain a better understanding of this connection4,5,6. These results indicate that aging and cancer share common features that could serve as therapeutic targets, offering important insights into cancer incidence and management. Additionally, some research has revealed differences in molecular patterns across cancer patients of different ages, indicating that age can influence the effectiveness of immune checkpoint therapy and potentially offer new indicators for diagnosis and treatment7,8. Although these studies have focused primarily on the aging process across the entire lifespan, they have failed to address tumor senescence irrespective of patient age.

Apart from elderly individuals, a minority of young people have an uncommon yet lethal condition known as Hutchinson‒Gilford progeria syndrome, which results in fast aging. This finding indicates the existence of regulatory mechanisms that affect the aging process beyond just chronological age9. Furthermore, an abnormal increase in oncogenic signals or deactivation of tumor suppressor genes in cells can trigger a cellular response called oncogene-induced senescence (OIS), which acts as a barrier to tumorigenesis by halting cell proliferation10. Tumor cells that are subjected to exogenous damage, such as radiation therapy and chemotherapy, display dormancy characteristics as a result of therapy-induced senescence (TIS)11. However, recent research has revealed that senescent cells play a key role in promoting tumor phenotypes and are now recognized as novel characteristics of cancer12. These biological features suggest that mapping the senescence landscape in gliomas can assist in accurately categorizing patients and informing treatment choices.

In this study, we collected senescence-related genes and identified two robust clusters of senescence-associated genes based on their gene set characteristics, which we named PMA-high and PMA-low according to their expression profile features. These two subgroups exhibit distinct survival outcomes, clinical traits, tumor microenvironment characteristics, somatic mutation spectra, chromosomal copy number variations, infiltration of hallmark gene sets, and responses to immunotherapy, as well as marked differences in drug treatment responses. The PMA-high subgroup is characterized by poor prognosis, enrichment of oncogenic pathways, immunosuppression, and diverse mutations. Drug prediction suggests that dasatinib may reverse the prognosis associated with the PMA-high signature. Additionally, we extracted a 200-gene signature to quantify our classification. Subsequently, we employed differential analysis, Weighted Gene Co-expression Network Analysis (WGCNA), single-sample Gene Set Enrichment Analysis (ssGSEA), and machine learning to identify biomarkers for the PMA glioma subtype. Ultimately, we determined that CEBPB and LMNA are not only biomarkers of poor prognosis in glioma patients but also key targets in the senescence model of glioma.

Overall, this study proposes the existence of senescence-based subtypes in gliomas, enhancing the understanding of how senescence impacts tumor biology and patient outcomes. These findings highlight the potential of using senescence-related biomarkers for personalized treatment strategies and underscore the importance of further research into senescence-related mechanisms and their role in cancer progression.

Materials and methods

Data acquisition and processing

Transcriptome expression data and clinical data were obtained from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/) and the Chinese Glioma Genome Atlas (CGGA) database (http://www.cgga.org.cn)13,14,15,16,17,18. The external dataset is sourced from Gliovis, and the drug treatment data is derived from the GEO database (GSE159609)19. Additionally, genomic mutation data for gliomas were acquired from the GDC Data Portal (https://www.example.com).

Clustering analysis of glioma patients

A total of 368 aging-associated genes were identified (Table S1). Unsupervised clustering analysis was performed via the R package CancerSubtypes to select the most suitable number of clusters, and the k-means algorithm was used to perform clustering for a range of 2–6 clusters with a cutoff of 0.00120. The delta area and consensus cumulative distribution function (CDF) determined that the most suitable clustering solution was k = 2. The associations of the two identified clusters with clinical phenotypes were examined, highlighting significant distinctions between the subtypes. The R package limma was used to conduct differential gene expression analysis, identifying genes that were differentially expressed between the two clusters. The Gene Ontology (GO) database was used to evaluate biological process (BP), molecular function (MF), and cellular component (CC) categories. The Kyoto Encyclopedia of Genes and Genomes (KEGG) database21,22,23 was utilized to analyze pathways associated with the identified genes, and specific gene sets were analyzed to provide a more detailed understanding of the functions and interactions of the identified SRGs.

WGCNA

WGCNA was conducted to find modules linked to the two identified clusters24. The computation of the ideal value for the soft-thresholding power was led by the scale-free topology criterion, with the minimum module size set to 50 genes. Modules were identified via the dynamic tree cut method, with the MEDissThres parameter set at 0.25.

Somatic mutation analysis

Differences in somatic mutations between clusters were investigated25. The R package maftools were used to visualize the mutation landscape, estimate mutational allele tumor heterogeneity (MATH), determine the tumor mutation burden (TMB), and evaluate the occurrence of coexistence and mutual exclusivity events. The chi-square test was used to evaluate whether the distribution of specific gene mutations differed significantly between the two clusters, with a P value < 0.05 considered significant.

TME characterization

The R software package ESTIMATE26 was used to predict immune and stromal scores for individual tumor samples. The ESTIMATE score was computed to provide an overall estimate of nontumor cell infiltration within the TME. Features from the TISIDB portal were used to elucidate tumor-immune cell interactions. ssGSEA was performed to quantify the relative abundance of various immune cell types within the TME27. The R package IOBR was used to assess features indicative of immune exhaustion, suppression, and exclusion within the TME28.

Machine learning model construction and gene intersection analysis

Machine learning techniques, including RF, LASSO regression, and SVM, were employed to analyze the data. Differential analysis was initially performed on the two subtypes, followed by an intersection with module-related genes from the WGCNA results to identify differentially expressed genes (DEGs) associated with traits. Survival analysis was conducted via the R packages survival and survminer, with coxPfilter = 0.05 as the filtering criterion, to identify prognosis-related genes. The intersection of these genes yielded core prognostic genes.

Immunohistochemistry

The tissue sections embedded in paraffin were dewaxed and rehydrated, followed by microwave antigen retrieval. The tissues were then blocked with 0.3% hydrogen peroxide solution for 30 min. Primary antibodies against CEBPB (Proteintech, Cat No. 23431–1-AP, 1:200) and LMNA (Servicebio, GB152400-100, 1:500) were added dropwise to the tissues and incubated overnight at 4 °C. Subsequently, the tissue samples were incubated with biotinylated secondary antibodies at room temperature for 10 min. After incubation, the samples were stained with diaminobenzidine (DAB) and hematoxylin. Fluorescent images were captured using a microscope.

Statistical analysis

Statistical analyses were performed via R software version 4.3.0 (https://www.R-project.org/), and graphical plots were generated via GraphPad 8.0. The Wilcoxon test was used to compare two paired groups, whereas analysis of variance (ANOVA) was used for comparisons among three or more groups. The log-rank test was performed via the R package Survminer to compare the survival distributions of the two groups. P values less than 0.05 were considered statistically significant.

Results

SRG expression defines distinct glioma subgroups with prognostic implications

We identified a total of 368 senescence-related genes from the literature, forming our final aging gene cluster. Using unsupervised clustering methods on the merged datasets (CCGA_325, and TCGA_693), glioma patients were stratified into two distinct clusters, the C1 and C2 clusters (Fig. 1A–D). Patients in the C2 cluster had a better prognosis than those in the C1 cluster (Fig. 1E). We designated C1 as PMA-high and C2 as PMA-low. The clinical characteristics of the samples indicated that the PMA-high subgroup had a greater proportion of grade IV gliomas, with 38 cases of WHO II (9%), 104 cases of WHO III (24%), and 299 cases of WHO IV (68%), in contrast to the PMA-low subgroup, which had 234 cases of WHO II (44%), 218 cases of WHO III (42%), and 75 cases of WHO IV (14%) (Fig. 1F). Additionally, an increase in the WHO grade was correlated with poorer patient prognosis (Fig. 1G). We further assessed the distribution of PMA in relation to different IDH statuses. In the PMA-low subgroup, 398 cases (81%) were IDH_Mutant, and 92 cases (19%) were IDH_Wildtype, whereas the PMA-high subgroup consisted of 102 cases (24%) IDH_Mutant and 329 cases (76%) IDH_Wildtype (Fig. 1H). Patients with IDH_Wild-type status had a significantly worse prognosis than those with IDH_Mutant status (Fig. 1I), with all the observed differences reaching statistical significance (P < 0.01). Subsequently, in the study of high-grade gliomas, we conducted an in-depth analysis based on the classification system proposed by Verhaak et al. in 201029, namely Classical, Mesenchymal, Proneural, and Neural subtypes, as well as the two DNA methylation (DNAm) subgroups identified by Drexler et al. in 202430. We performed a one-to-one correspondence study between the subgroups from these two studies and our PMA_High and PMA_low subtypes and found that the PMA-high subtype often exhibits consistency with the subgroup associated with poorer prognosis (Figure S1A).

Fig. 1
figure 1

SRG expression defines distinct glioma subgroups with prognostic implications. (A)–(C) Identification of glioma prognosis subgroups using SRGs; (D) Heatmap of prognosis subgroups; (E) Subgroup C1 (PMA_H) exhibits poor prognosis. (F, G) Proportion and prognosis of PMA in WHO grades. (H, I) Proportion and prognosis of PMA in IDH status.

Full size image

Immune characteristics of glioma subclusters

The tumor microenvironment (TME) encompasses the interplay between malignant tumor cells, immune cells, and stromal cells, which collectively promote or inhibit tumor progression. We employed the ESTIMATE algorithm to assess the overall TME status across two aging clusters. Notably, the immune scores, stromal scores, and ESTIMATE scores were all highest in the PAM-high cluster, indicating a more abundant composition of immune and stromal cells within this subtype (Fig. 2A). Analysis of immune function and immune cell profiles revealed enhanced malignant proliferation in the PMA-high cluster. Specifically, certain pathways and immune signatures, particularly human leukocyte antigen (HLA), C–C chemokine receptor (CCR), and T cell co-inhibition, were significantly more active or prevalent in the PMA-high cluster. However, TH1 cells, which play a crucial role in anti-tumor immunity by directly killing tumor cells, promoting adaptive immune responses, and modulating the TME, exhibited subdued activity in the PAM-high cluster, suggesting a higher degree of malignancy and cachexia in PAM-high tumors. To better evaluate the immune status between subgroups, we calculated Cytolytic Scores for the PMA subtypes (Fig. 2C). Cytolytic Score is a biomarker that reflects the cancer immune status by assessing the cytolytic activity of immune cells in the TME, calculated as the geometric mean of mRNA expression levels of granzyme A (GZMA) and perforin 1 (PRF1). Our study found that higher CYT scores were associated with an immunosuppressive TME and poorer survival prognosis, consistent with previous research. Additionally, we conducted a CIBERSORT analysis of immune cell infiltration in the subgroups, which revealed a higher enrichment of Macrophages M2 in the PMA-HIGH cluster (Fig. 2D). TIDE immunotherapy prediction indicated that PMA-low patients were more likely to benefit from immunotherapy (Fig. 2E PMA_L: PMA_H = 56.8%:31%). Furthermore, the TIDE total score, exclusion score (Exclusion), dysfunction score (Dysfunction), CD274 (also known as programmed death-ligand 1, PD-L1), and Cancer-Associated Fibroblasts (CAF) were all higher in PMA_H (Fig. 2E–J), suggesting weakened T cell activity and cytolytic capacity, severe dysfunction, and the exclusion of tumor cells from the tumor tissue, which renders them ineffective in attacking and killing tumor cells, leading to immunotherapy insensitivity. Subsequently, we conducted drug prediction on the two subpopulations, focusing primarily on the efficacy of senescence-targeting drugs and broad-spectrum anticancer drugs. The results revealed that the PMA_H is more likely to benefit from treatment with the senescence-targeting drug dasatinib and the broad-spectrum anticancer drug fluorouracil (Figure S1B-1C). To validate our findings, we retrieved and obtained glioma sample data related to dasatinib treatment studies (GSE159609) from the GEO database. After performing differential analysis to identify differentially expressed genes (Figure S1D-1E), despite the limitation imposed by the sample size, we observed a reversal in pathways that were originally upregulated in the PMA_H, such as HALLMARK_MYC_SIGNALING, PI3K_AKT_MTOR_SIGNALING, and KEGG_MAPK_SIGNALING_PATHWAY, following drug treatment (Figure S1F-1G). Additionally, the infiltration level of M2 macrophages decreased after drug treatment (Figure S1H).

Fig. 2
figure 2

Immune molecular functional landscape of glioma subtypes. (A) Differences in immune microenvironment-related scores among subtypes; (B) score distribution of immune cells and immune functions among subtypes; (C) differences in cytolytic scores among subtypes; (D) CIBERSORT immune cell infiltration; (E) distribution of TIDE immunotherapy response between two subtypes; (F) TIDE Scores of Two Subtypes; (G) EXCLUSION scores of two subtypes; (H) dysfunction scores of two subtypes; (I) proportion of CD247 between two subtypes; (J) proportion of CAFs between two subtypes.

Full size image

Pathway activation and functional differences across glioma subtypes

Through differential analysis of the subtypes, we identified a total of 1579 differential genes (Fig. 3A, Table s2). Further analysis utilizing the Kyoto Encyclopedia of Genes and Genomes (KEGG) database, Gene Set Enrichment Analysis (GSEA), and specific gene set analysis revealed that the PAM-high subtype specifically activated the MAPK signaling pathway, cell cycle, FoxO signaling pathway, and Wnt signaling pathway (Fig. 3B). These signaling pathways are generally closely related to the initiation and progression of tumors as well as the aging process.

Fig. 3
figure 3

Pathway Activation and Functional Differences between Glioma Subtypes. (A) Volcano plot of differentially expressed genes between subgroups; (B) KEGG pathway analysis of differentially expressed genes identifying prognosis-related pathways in subgroups; (C)–(E) An overview of the upregulation of 50 hallmark gene sets; (F) The aging-related pathway is upregulated in PMA-high; (G) An overview of the downregulation of 50 hallmark gene sets; (H, I) GO analysis differences between subgroups identified by differentially expressed genes.

Full size image

In addition, through hallmark gene set GSEA analysis, we observed that the PAM-high subtype upregulated gene sets related to epithelial-mesenchymal transition (HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION), interferon-gamma response (HALLMARK_INTERFERON_GAMMA_RESPONSE), and tumor necrosis factor (TNF) signaling via the NF-κB pathway (HALLMARK_TNFA_SIGNALING_VIA_NFKB) (Fig. 3C–E). Concurrently, the activation of the Notch signaling pathway, JAK/AKT signaling pathway, and NF-κB signaling pathway, which are closely associated with aging features, were also prominent characteristics of the PAM-high subtype (Fig. 3F).

In contrast, in the PMA-high subtype, we observed downregulation of HALLMARK_KRAS_SIGNALING_DN (downregulation of KRAS signaling pathway) and HALLMARK_PANCREAS_BETA_CELLS (gene set related to pancreatic beta cells) (Fig. 3G).

Further, GO (Gene Ontology) analysis suggested that these differential genes are closely related to biological processes such as DNA-binding transcription factor binding, GTPase regulatory activity, proteasome-mediated ubiquitin-dependent protein catabolic process, and regulation of nervous system development, implying that the PAM-high subtype may possess proliferation-related characteristics (Fig. 3H, I).

Multi-omics characteristics of glioma subtypes: genomic variations, mutational burden, and chromosomal copy number variations

We aimed to construct a scoring system to accurately quantify the characteristics of aging genes in glioma. By downloading and merging all samples from TCGA_LGG and TCGA_GBM (TCGA RNAseq) through CGGA, and excluding those samples with missing survival time and status, we retained data from 603 glioma patients as our validation dataset. We identified the 100 genes with the most significant differences between the two subtypes and quantified their expression in the validation set (Fig. 4A). Analysis of the TCGA database indicated that the PMA-high subtype was associated with poor prognosis (Fig. 4B). We subsequently stratified the data via the CMOIC, NTP, and PMA clustering methods. The consistency between different clustering methods was quantified via the kappa value and P value from consensus matrix analysis (Fig. 4C–E). High consistency was observed among CMOIC, NTP, and PMA clustering (Kappa values of 0.834, 0.892, and 0.924, respectively, with P values < 0.001), indicating considerable accuracy and reliability in these classifications and predictions.

Fig. 4
figure 4

Genomic variations and mutation burden in glioma subtypes. (A) Heatmap of PMA subgroups validated in the TCGA dataset using 200 genes. (B) Prognosis validation. (C)–(E): Clustering validation and display of CMOIC, NTP, and PMA clusters. (F) Display of tumor mutation burden and clustering prognosis; (G) High tumor mutation burden in PMA_H. (H, I) Presentation of mutation burden genes in different clusters.

Full size image

We analyzed the TMB to account for genomic variations. The overall mutation spectrum revealed significant differences in TMB (Fig. 4G). The PMA-high subgroup had a higher TMB. Among the two subtypes, a greater mutation burden was correlated with a worse prognosis, with H-TMB + PMA-high patients having the poorest prognosis and L-TMB + PMA-low patients having the best prognosis (Fig. 4F). The top 25 highly mutated genes are presented for each cluster. Both subtypes frequently exhibited mutations in IDH1, TP53, and ATRX, but IDH1 mutations were present in only 18% of the PMA-high patients compared with 89% of the PMA-low patients. The PMA-low subtype presented a greater frequency of mutations in TP53 (48%) and ATRX (40%) (Fig. 4H, I). Additionally, various low-frequency mutations, such as those in EGFR (1%) and PTEN (1%), were observed, which may affect different signaling pathways and the biological characteristics of glioma. In the PMA-high subgroup, the mutation spectrum was more complex and diverse. IDH1 mutations were present in 18% of the patients, with reduced frequencies of TP53 (34%) and ATRX (14%) mutations. New high-frequency mutations, such as CIC (21%) and NOTCH1 (7%), have emerged, reflecting different genetic backgrounds and disease processes. The presence of multiple low-frequency mutations in PMA-high patients suggests a more complex genetic background and greater heterogeneity, potentially promoting glioma progression and deterioration. Furthermore, gene copy number variation (CNV) analysis revealed significant differences in CNVs between the two subtypes. These differences were primarily observed on chromosomes 7, 9, and 10. Specifically, in the PMA_H subtype, we observed deletions in the regions 10q26.3, 10q23.31, and 9p21.3, as well as gains in the 7p11.2 region. In contrast, chromosomal variations were generally not prominent in the PMA_L subtype, with the most notable deletion occurring in the 19q13.43 region and the most significant gain in the 7p32.3 region (Figure S1I-1 J).

In addition, we obtained the public cohort studies from Gravendel and Rembrandt as additional validation sets (Figure S2A-2H). The studies demonstrated that in both validation sets, PMA-H exhibited poor prognosis, high scores across various tumor microenvironment indicators, and immune cell and immune function scores that were consistent with those in the validation sets. Furthermore, the hallmark gene set scores were also found to be consistent with our test set.

Co-expression modules associated with glioma subtypes

The soft threshold power was determined as β = 11 (Fig. 5A) by applying the scale-free topology criterion. We next established co-expression modules with a cutting height ≥ 0.25. Hierarchical clustering analysis of modules within the same branch revealed comparable gene expression patterns (Fig. 5B–D). As a result, these gene modules with comparable characteristics were organized into eleven separate co-expression modules, each represented by a different color, such as gray, brown, blue, or yellow (Fig. 5E). We subsequently employed visualization techniques to represent the gene clusters and analyzed them to determine the relationships between different modules. The light green module was significantly correlated with the PMA-high subtype (Fig. 5F, G). This finding emphasizes the importance of further investigating the light green module, as shown in Table S3.

Fig. 5
figure 5

Identification of Special Gene Groups by WGCNA. (A) Determination of the optimal scale-free index (β) and mean connectivity across various soft-thresholding powers; (B) gene dendrogram and module colors; (C) hierarchical clustering analysis; (D) gene dendrogram based on clustering results; (E) heatmap of correlations between module genes and two isoforms; (F, G) correlation between genes and PMA_H, and gene clustering diagram.

Full size image

Core prognostic genes identified in glioma subtypes

A total of 898 module-related genes and 1,579 DEGs were identified. By intersecting these genes with the aging gene set, we identified 12 DEGs strongly associated with the PMA-high subtype (Fig. 6A). Univariate analysis of these DEGs revealed that each was an adverse prognostic factor in glioma (Fig. 6B). The clinical trait heatmap revealed the highest expression levels in the clinical group with poor prognosis (Fig. 6C), and high expression levels were particularly evident in the PMA-high group (Fig. 6D–A1). High expression was significantly associated with poor prognosis. We subsequently analyzed the 12 genes via LASSO regression, RF, and SVM methods. LASSO regression identified eight genes (Fig. 7A–C), SVM identified five genes, and RF identified seven genes. Considering the intersection of these results, we identified two core genes (Fig. 7D). Expression and survival analyses of CEBPB and LMNA, which use both TCGA and CCGA data, consistently indicated that these genes are not only biomarkers for poor prognosis in glioma patients but also key targets in the glioma aging model (Fig. 7E–H).

Fig. 6
figure 6

Identification of core prognostic genes. (A) SRG genes associated with the differential module in PMA_H. (B) Univariate analysis of SRG genes associated with the differential module in PMA_H. (C) Heatmap display of SRG genes associated with the differential module in PMA_H. (D)–(O) Expression of genes in glioma PMA. (P)–(A1): Prognosis of genes in glioma.

Full size image
Fig. 7
figure 7

Machine learning-based screening of core genes for cluster prognosis. (A) Lasso regression analysis; (B) support vector machine analysis; (C) random forest analysis; (D) intersection of genes identified by three machine learning methods; (E, F) Expression of genes in glioma PMA (TCGA); (G, H) Prognosis of genes in glioma (TCGA).

Full size image

Expression characteristics of CEBPB and LMNA in glioma and their verification

To further delve into the expression characteristics and potential roles of CEBPB and LMNA in gliomas, we systematically analyzed the RNA and protein expression levels of these two genes in brain tumors using the TCGA database and the Human Protein Atlas database. The results indicated that the expression of CEBPB and LMNA is lower in normal brain tissue compared to gliomas, and they exhibit high expression levels in 80 brain tumor cell lines, particularly in common glioma cell lines such as U87, LN229, and U251MG, where their expression is especially prominent (Fig. 8A–D). In gliomas of different grades, higher-grade gliomas exhibit higher expression levels of CEBPB and LMNA, which is confirmed by both CGGA and TCGA data (Fig. 8E, H). High expression of CEBPB and LMNA is more likely to benefit from the senescence-eliminating drug dasatinib (Fig. 8I, J). After drug treatment, the expression of both genes decreased compared to the untreated state. Although limited by sample size, the P value for LMNA did not reach statistical significance (P < 0.05), whereas CEBPB exhibited a significant response (Fig. 8K, L). Furthermore, we conducted a detailed analysis of gliomas of different grades and found that the expression levels of CEBPB and LMNA increased with the grade of gliomas, reaching the highest levels in glioblastoma. To validate this finding, we performed immunohistochemistry experiments, which further confirmed that the expression levels of CEBPB and LMNA were indeed higher in WHO grade III gliomas compared to WHO grade II gliomas (Fig. 8M–R). These findings provide important clues for our in-depth understanding of the roles of CEBPB and LMNA in the development and progression of gliomas.

Fig. 8
figure 8

Expression characteristics of CEBPB and LMNA in glioma and their verification. (A) The expression of CEBPB in glioma is higher than that in normal tissue. (B) CEBPB expression in brain tumor cell lines; (C) the expression of LMNA in glioma is higher than that in normal tissue. (D) LMNA expression in brain tumor cell lines; (E) CEBPB expression in different grades of glioma (CGGA); (F) CEBPB expression in different grades of glioma (TCGA); (G) LMNA expression in different grades of glioma (CGGA); (H) LMNA expression in different grades of glioma (TCGA); (I) Drug prediction of dasatinib in different CEBPB groups; (J) Drug prediction of dasatinib in different LMNA groups; (K) Decreased CEBPB expression after dasatinib treatment; (L) Decreased LMNA expression after dasatinib treatment; (M)–(O) LMNA expression in normal and different grades of glioma; (P)–(R) CEBPB expression in normal and different grades of glioma.

Full size image

In gliomas of different grades, higher-grade gliomas exhibit higher expression levels of CEBPB and LMNA, which is confirmed by both CGGA and TCGA data (Fig. 8E–H). High expression of CEBPB and LMNA is more likely to benefit from the senescence-eliminating drug dasatinib (Fig. 8I, J). After drug treatment, the expression of both genes decreased compared to the untreated state. Although limited by sample size, the P value for LMNA did not reach statistical significance (P < 0.05), whereas CEBPB exhibited a significant response (Fig. 8K, L). Furthermore, we conducted a detailed analysis of gliomas of different grades and found that the expression levels of CEBPB and LMNA increased with the grade of gliomas, reaching the highest levels in glioblastoma. To validate this finding, we performed immunohistochemistry experiments, which further confirmed that the expression levels of CEBPB and LMNA were indeed higher in WHO grade III gliomas compared to WHO grade II gliomas (Fig. 8M–R). These findings provide important clues for our in-depth understanding of the roles of CEBPB and LMNA in the development and progression of gliomas.

Discussion

Glioma, the most common primary malignant tumor originating within the cranial cavity, is characterized by poor prognosis, high disability rates, and frequent recurrence, imposing a significant burden on patients, families, and society. Glioblastoma is the most malignant type of glioma. Even with the most extensive and safe removal of tumors through standard surgical procedures combined with concurrent temozolomide (TMZ) chemoradiotherapy, patients still have a median survival time of no more than 15 months. Unlike tumors in other organs, gliomas exhibit more complex characteristics, with malignant glioma cells demonstrating robust invasive capabilities early on and infiltrating surrounding normal tissues. Consequently, even if the tumor is completely resected, recurrence at the surgical margins is difficult to avoid31,32,33.

Previous research has revealed the dual role of cellular senescence in cancer. In the early stages, it inhibits tumorigenesis and limits malignant progression. Nevertheless, the persistence of the tumor can lead to the promotion of tumor growth, stimulation of invasion and migration, and inhibition of immune surveillance3,12,34,35. Despite notable advancements in the examination of senescence markers such as senescence-associated β-galactosidase, p16, and p21, identifying specific and definitive markers for tumor senescence remains a major challenge.

The hallmarks of aging encompass genomic instability, epigenetic alterations, chronic inflammation, and gut microbiota dysbiosis, among others. The “driving factors” of cancer include the generation of sustained proliferative signals, evasion of growth suppression, resistance to cell death, acquisition of replicative immortality, induction of angiogenesis, and activation of invasion and metastasis12,36,37. During the process of human aging, cellular DNA accumulates various damages, such as point mutations, insertions, and deletions, which lead to genomic instability and subsequently trigger cellular senescence or death38. Notably, genomic instability is also a prominent feature of cancer cells, which continuously accumulate genetic mutations and exhibit chromosomal abnormalities, thereby gaining the ability to proliferate indefinitely.

Previous studies have revealed that senescence-associated secretory phenotypes play a crucial role in intercellular communication and exert a profound impact on the dynamic evolution of the tumor microenvironment (TME). Specifically, hypoxic adaptation, metabolic dysregulation, and phenotypic shifts in immune cells regulated by cellular senescence synergistically promote the development of an immunosuppressive microenvironment and chronic inflammation, thereby driving tumor progression.

Given that senescence is a characteristic feature of cancer and that the molecular classification of cancer aids in predicting specific targets for treatment, thereby driving the advancement of precision medicine39, investigating the role of senescence in glioma classification to guide clinical applications holds significant research importance.

In this study, we delved into the distinct roles played by 368 SRGs in gliomas. Drawing upon the foundational work by Verhaak et al. in 201029, which classified gliomas into classic, mesenchymal, proneural, and neural subtypes, and the recent study by Drexler et al. in 202430 that identified two DNA methylation (DNAm) subgroups within glioblastomas, our investigation revealed the existence of two distinct subtypes: PMA-high and PMA-low, based on glioma prognosis. The PMA-high subtype was found to be associated with poor prognosis and demonstrated varying distributions of glioma patients across WHO grades and IDH statuses in our sample analysis. Specifically, among 374 WHO grade IV patients, 299 (79.95%) were classified as PMA-high, whereas only 75 (20.05%) were PMA-low. Similarly, in the cohort of 421 IDH_Wild-type patients, 329 (78.14%) were identified as PMA-high, with 92 (21.85%) falling into the PMA-low category. Further analysis revealed significant variations in expression profiles, immunity, and hallmark gene set scores between the two subtypes, highlighting the heterogeneous nature of PMA-high patients. Moreover, the DEGs activated proliferation and immune function-related characteristics as well as various tumor-related pathways, including the MAPK signaling pathway, cell cycle regulation, the FoxO signaling pathway, and the Wnt signaling pathway. Specific gene set analysis further demonstrated that PMA-high upregulated signaling pathways such as HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION, HALLMARK_INTERFERON_GAMMA_RESPONSE, and HALLMARK_TNFA_SIGNALING_VIA_NFKB. Shedding light on the underlying molecular mechanisms contributing to the distinct prognostic outcomes. These pathways are closely linked to tumors and have been reported to be activated in various types of cancer40,41,42,43,44. Additionally, studies have suggested that the MAPK signaling pathway, NF-KB signaling pathway, PI3K/AKT signaling pathway, JAK/STAT signaling pathway, and NOTCH signaling pathway are associated with cellular senescence. Therefore, we investigated the activation status of these pathways in the two subgroups. Our research indicated that these pathways were activated in the PMA-HIGH subgroup. Moreover, the PAM-high subgroup harbors more chromosomal copy number variations, which aligns with an EMT-based subclassification. In cases with worse prognosis, there is often a higher tumor mutational burden and a greater occurrence of copy number variations.

To enhance the generalizability and reliability of our findings, we went beyond initial discoveries by validating our conclusions across multiple independent datasets. This approach confirmed the stability and reproducibility of our results, providing stronger support for our research. Mutation analysis further substantiated our subtype classification. Our data revealed a greater TMB in the PMA-high subgroup. Importantly, while mutations in IDH1, TP53, and ATRX were common in both subgroups, the IDH1 mutation rate was significantly different: 18% in the PMA-high subgroup compared with 89% in the PMA-low subgroup. Furthermore, TP53 and ATRX exhibited distinct mutation patterns in the PMA-high subgroup, contributing to a more complex and diverse mutation spectrum. Despite the presence of IDH1 mutations (18%), the mutation frequencies of TP53 (34%) and ATRX (14%) were significantly lower. New high-frequency mutations, such as CIC (21%) and NOTCH1 (7%), were also identified. This mutation diversity may indicate different genetic backgrounds and disease processes, with multiple low-frequency mutations in PMA-high patients suggesting a more intricate genetic landscape and greater heterogeneity. The coexistence of these mutations could drive glioma progression and deterioration. Moreover, the PAM-high subgroup harbors more chromosomal copy number variations, which aligns with an EMT-based subclassification45. In cases with worse prognosis, there is often a higher tumor mutational burden and a greater occurrence of copy number variations.

We utilized various algorithms, including machine learning, to identify characteristic genes differentiating the two subtypes. Our studies consistently revealed that the genes encoding CEBPB and LMNA are not only biomarkers of poor prognosis in glioma patients but also critical targets in glioma senescence models. It has been reported in the literature that the expression of C/EBPbeta and STAT3 in human gliomas is associated with mesenchymal differentiation and portends poor clinical outcomes46, which is consistent with our findings. Specifically, CEBPB is highly expressed in tumors, particularly in the PMA-HIGH subgroup, where the expression of both CEBPB and LMNA is even more pronounced. Furthermore, CCL20 secreted by CRC cells can recruit Tregs through the FOXO1/CEBPB/NF-κB signaling pathway47, thereby promoting chemotherapy resistance. This suggests that the FOXO1/CEBPB/NF-κB/CCL20 axis may represent a promising therapeutic target for CRC treatment. In breast cancer, tumor glycolysis coordinates a molecular network involving the AMPK-ULK1, autophagy48, and CEBPB pathways to influence myeloid-derived suppressor cells (MDSCs) and maintain tumor immune suppression. In RCC cells, CEBPB enhances the expression of IL6, which in turn promotes the phosphorylation of STAT3 and the expression of its downstream target genes49. This mechanism ultimately facilitates tumor progression. Additionally, the NAT10/CEBPB/vimentin axis plays a crucial role in the malignant progression of salivary gland adenoid cystic carcinoma (SACC), and these findings could be leveraged to improve the treatment of SACC50. Collectively, these observations underscore the pivotal role of CEBPB in multiple types of tumors. Molecular characterization has revealed a novel LMNA-NTRK1 rearrangement on chromosome 1 with oncogenic potential. A patient harboring this rearrangement achieved a partial response to treatment with the pan-TRK inhibitor entrectinib, with the target liver lesions decreasing in size from 6.8 and 8.2 cm in longest diameter to 4.7 and 4.3 cm, respectively51.

It is noteworthy that we conducted TIDE immunotherapy prediction and drug sensitivity analysis on different patient groups. The results showed that the immunotherapy response in the PMA_H was significantly inferior to that in the PMA_L. Specifically, patients in the PMA_L were more likely to benefit from immunotherapy, the immunotherapy response rates for PMA_L and PMA_H were 56.8% and 31%, respectively). In addition, the TIDE total score, exclusion score, dysfunction score, CD274 (also known as programmed death ligand 1, PD-L1) level, and cancer-associated fibroblasts (CAF) level were all higher in the PMA_H. This indicates that in the PMA_H, T-cell activity and lytic capacity were reduced, dysfunction was severe, and tumor cells were excluded from the tumor tissue, leading to immune cells being unable to effectively attack and kill tumor cells. Therefore, immunotherapy was not sensitive for patients in the PMA_H.

Drug sensitivity analysis also revealed that although patients in the PMA_H had a poorer prognosis, they seemed to benefit from treatment with the anti-aging drug dasatinib and the broad-spectrum antitumor drug fluorouracil. Multiple studies have shown that CEBPB and LMNA are highly expressed in malignant gliomas and are closely related to poor prognosis. Patients with high expression of CEBPB and LMNA are more likely to benefit from the senescence-eliminating drug dasatinib, and the expression levels of these two genes decreased after drug treatment. Although limited by sample size, the change in LMNA expression did not reach statistical significance (P  < 0.05), while the response of CEBPB was very significant.

Furthermore, we collected clinical samples for immunohistochemical experiments. The experimental results showed that the expression levels of CEBPB and LMNA increased with the grade of glioma in normal tissue and WHO grade II and III gliomas. This result is consistent with the findings from database analysis.

In this study, we conducted an in-depth analysis of the specific roles of 368 SRGs in gliomas. Using unsupervised clustering methods, we identified two distinct senescence-related clusters that exhibited significant differences in clinical traits, immune landscapes, mutation profiles, and TME infiltration. To accurately quantify the senescence gene sets and pinpoint molecular markers for the two subtypes, we conducted differential analysis, WGCNA, and survival analysis and used various machine learning techniques. These multifaceted analyses allowed us to identify key genes that not only explain the prognosis of different subpopulations but also enable precise and consistent classification of glioma samples, laying a strong foundation for future research.

While our multi-cohort integrative analysis and IHC validation strongly associate CEBPB/LMNA expression with distinct molecular subtypes and clinical outcomes, we acknowledge that functional validation of their mechanistic roles remains to be explored. Specifically, in vitro and in vivo model experiments are needed to dissect whether CEBPB and LMNA act as drivers or bystanders in shaping glioma phenotypes. In addition, the interactions between these biomarkers and canonical oncogenic pathways also require systematic investigation. These issues are currently the focus of our ongoing work. Importantly, the correlation patterns and preliminary protein-level evidence presented in this paper provide a testable framework for such hypothesis-driven research. We anticipate that elucidating these mechanisms will refine subtype-specific therapeutic strategies.

In summary, our molecular classification method based on SRGs provides a more accurate prediction of survival outcomes for glioma patients. This study offers insights into the impact of cellular senescence on gliomas, revealing two subtypes with significantly different prognostic characteristics and highlighting the inherent heterogeneity of gliomas. Moreover, our data provide a resource for understanding the complex relationship between senescence and cancer, potentially opening new avenues for research and therapeutic strategies.