The clinical and gene expression data of patients with LGGGBM were obtained from the TCGA database (https://cancergenome.nih.gov/) and the CGGA database (http://www.cgga.org.cn/). Normal brain RNA sequencing data were obtained from GTEx and utilized as normal control data. In this research, data from 609 glioma samples and 1152 normal brain samples were analyzed. Among them, the data used in the gene expression data were normalized. Next, DEGs were chosen using the R software’s “limma” package with the absolute value of the log2-transformed fold change (FC) > 1 and the adjusted p-value (adj. P) < 0.05 as the cutoff. In addition, the ComBat method was used to reduce batch effects with the R package “sva”.

The profiles of lncRNAs were acquired from the RNAseq dataset. The log2 transformation was used to normalize the total RNA expression data. The list of cellular senescence-related genes was obtained from the Human Aging Genomic Resources (HAGR) (https://genomics.senescence.info). The Pearson correlation was used to calculate the relationship between lncRNAs and cellular senescence-related genes. Cellular senescence-related lncRNAs (SRlncRNAs) were screened by a correlation coefficient R ² > 0.4 and a p < 0.001. Finally, coexpression networks were visualized using Cytoscape software 3.8.0.

To construct a validated prognostic model, 609 glioma patients were randomly divided into training and testing cohorts. Ultimately, 306 patients were enrolled in the training cohort, and 303 patients were enrolled in the testing cohort. The key characteristics of each cohort are shown in Table 1. The SRlncRNA signature was derived based on the training cohort, and its potential to predict patient survival was validated utilizing the testing cohort and the whole cohort. We also confirmed the prognostic signature in the TCGA-LGG cohort, the TCGA-GBM cohort and the CGGA cohort (Supplementary Tables S1, S2).

The prognostic significance of cellular senescence-related lncRNAs was initially determined using univariate Cox regression. Least absolute shrinkage and selection operator (LASSO) regression was used to integrate the cellular senescence-related lncRNAs with p < 0.05 in univariate analysis. The LASSO results were then included in a multivariate Cox model to generate a risk score. A risk score was calculated using a linear combination of cellular senescence-related lncRNA expression levels multiplied by a regression coefficient (β): risk score = ∑ i = 1 n β i × (expression of lncRNA_i). Based on the median risk score, the patients were categorized into high-risk and low-risk groups. The log-rank test was used to compare the survival differences between the two groups.

To construct an independent prognostic model, Cox regression was used. Patient survival was predicted utilizing a nomogram. The correctness of the model was assessed by utilizing the concordance index (C-index), calibration curves, and receiver operating characteristic (ROC) curves. We used multivariate Cox regression with demographic data to determine whether the risk score was an independent predictor of patient outcomes.

The functional enrichment of the gene expression data was interpreted by employing Gene Set Enrichment Analysis (GSEA, http://www.broadinstitute.org/gsea/index.jsp). The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways associated with cellular senescence are displayed, as well as the functional enrichment of cell senescence-related lncRNAs with predictive significance.

To assess the riskscore and the level of correlation of the immune Microenvironment, immune-infiltrating scores of 16 immune cells, 13 immune function and immune checkpoint were calculated by singlesample gene set enrichment analysis (ssGSEA) using the “gsva” and “GSEABase” packages in R software. Additionally examined was the relationship between riskscore and immunological subtypes (C1–C6). At the same time, we calculated the immune infiltration between the high- and low-risk groups using the CIBERSORT, TIMER and ESTIMATE databases. Pearson correlation was used to determine the relationship between risk scores and immune infiltration.

All procedures for this research were approved by the Ethics Committee of Harbin Medical University (Harbin, China). Glioma samples and adjacent normal brain tissue were collected from 24 patients who were resected at Harbin Medical University Cancer Hospital. Patients had not been treated with radiotherapy or chemotherapy before surgery. All tissue samples were pathologically confirmed and immediately frozen in liquid nitrogen for storage until the RNA was extracted. Each patient signed a written informed consent form before the tissue samples were used for research purposes. The clinical characteristics of the 24 glioma patients are shown in Supplementary Table S4.

The human GBM cell lines (T98G and U251) and the normal human astrocyte cell line (NHA) were obtained from our laboratory (Li et al., 2019a). All cell lines were cultured in Dulbecco’s modified Eagle’s medium (DMEM; SIGMA) supplemented with 10% fetal bovine serum (FBS; PAN SERATECH) at 37°C in a humidified chamber containing 5% CO₂. All glioma patients from our hospital were informed and ethical approval for the current research was obtained through the Ethics Committee of Harbin Medical University (KY2021-42).

The RNAsimple Total RNA Kit (DP419, Tiangen Biotech Co., Ltd. Beijing, China) was used to extract total RNA from tissue and cells. The lnRcute lncRNA First-Strand cDNA Kit (KR202, Tiangen Biotech Co., Ltd.) was used to construct complementary DNA. The lnRcute lncRNA qPCR Kit (FP402, Tiangen Biotech Co., Ltd.) was used to perform RT-qPCRs according to the manufacturer’s instructions. To evaluate relative gene expression, the 2^−ΔΔCt method was utilized. GAPDH was used as an endogenous control. The primer sequences were as follows: SNAI3-AS1, 5-TGA​GGT​GCT​CCT​CCG​AGA​AT-3 (Forward) and 5-GTA​TAG​CTC​CCT​GGC​AGA​GTT​CA-3 (Reverse); GAPDH, 5-CAG​GAG​GCA​TTG​CTG​ATG​AT-3 (Forward) and 5GAA​GGC​TGG​GGC​TCA​TTT-3 (Reverse).

All tissues were fixed overnight in formalin solution, dehydrated in ethanol, embedded in paraffin and sectioned at 5 mm. To remove the paraffin, specimens were treated with xylene and ethanol. Slides were blocked with 5% normal goat serum and incubated overnight at 4°C with anti-YPEL3 (Proteintech). After washing with PBS, slides were incubated with goat anti-rabbit horseradish peroxidase for 30 min at room temperature. Immunohistochemical (IHC) reactions were detected using the DAB kit. Slides were examined under a phase-contrast light microscope (Nikon).

The survival curves were constructed using the Kaplan-Meier method and evaluated using the log-rank test. The predictive impact of the cellular SRlncRNAS signature and clinicopathological data was estimated by combining Cox regression and LASSO regression. The statistical analyses were carried out using the R programming language (version 4.1.3). The differences between groups were calculated utilizing a two-tailed Student’s t-test. The statistical tests were conducted bilaterally, with a significance level of p < 0.05.

The flowchart of this study is shown in Supplementary Figure S1. A total of 279 cellular senescence-related genes were obtained from HAGR, among which 273 genes were expressed in glioma, and a total of 483 cellular SRlncRNAs were obtained by constructing a cellular senescence-related mRNA and lncRNA coexpression network. Applying LASSO regression, 31 lncRNAs associated with cell senescence were identified (Figures 1A, B; Supplementary Table S3). In total, 13 lncRNAs with a low risk (hazard ratio (HR) < 1) and 18 lncRNAs with a high risk (hazard ratio (HR) > 1) were found. Furthermore, multivariate Cox analysis revealed 6 lncRNAs with prognostic significance and their coefficients among the aforementioned 31 cellular SRlncRNAs (Table 2). The names of the 6 lncRNAs were SNAI3-AS1, CRNDE, AGAP2-AS1, HOXD-AS2, PAXIP1-AS2, WAC-AS1. These 6 lncRNAs were exploited to construct an SRlncRNA signature. The risk score was calculated using the following formula: (0.54993 × CRNDE) + (0.11751 × AGAP2-AS1) + (0.28073 × HOXD-AS2) + (0.46398 × PAXIP1-AS2)—(0.45368 × WAC-AS1)- (1.23692 × SNAI3-AS1). According to the hazard ratio (HR) score obtained by the multivariate Cox regression analysis, CRNDE, AGAP2-AS1, HOXD-AS2, and PAXIP1-AS2 were risk factors, whereas SNAI3-AS1 and WAC-AS1 were protective factors. These six lncRNAs were used to construct an optimal prognostic risk model and a prognostic visual coexpression network of cellular SRlncRNAs-mRNAs (Figures 1C, D).

To assess the sensitivity and specificity of the riskscore in predicting the survival of glioma patients, the training cohort was categorized into low-risk group and high-risk group. The K-M curve results indicated that patients in the high-risk group had a significantly lower survival rate than those in the low-risk group (p < 0.001) (Figure 2A). The area under the ROC curve of the risk score (AUC) was 0.946, which was higher than the AUC values of other clinical parameters (Figure 2B). The ROC curves also gave AUC values of 0.904, 0.947, and 0.897 for overall survival (OS) at 1, 3, and 5 years, respectively (Figure 2C).

According to the heatmap, the expression of the 6 SRlncRNAs differed significantly between the high-risk group and the low-risk group (Figure 2D). Patients with glioma who had high risk scores had significantly worse survival rates than those who had low risk scores, as shown in a scatter plot (Figure 2E). Furthermore, the distribution plot of the riskscore corresponded to the patient population categorization (Figure 2F). Further prognostic results of the 6 SRlncRNAs examined by K-M curves revealed that greater expression of CRNDE, AGAP2-AS1, HOXD-AS2, and PAXIP1-AS2 and lower expression of SNAI3-AS1 and WAC-AS1 were substantially linked with shorter OS survival (p < 0.0001) (Figures 3A–F). Overall, the results confirmed that the 6 SRlncRNAs formed a highly accurate glioma prognostic risk model.

To validate the predictive power of the SRlncRNA prognostic signature, the risk scores were calculated for patients in both the testing cohort and the whole cohort, and patients were categorized into the low-risk group and high-risk group based on the median risk score. Analysis of OS by K-M curves for the testing cohort (p < 0.001) (Figure 4A) and the whole cohort (p < 0.001) (Figure 4C) demonstrated that these results were consistent with the training cohort. The ROC curves for the testing cohort (AUC = 0.917) (Figure 4B) and whole cohort (AUC = 0.914) (Figure 4D) demonstrated the accuracy of the predictions of the SRlncRNA signature for OS in glioma patients, which was verified further by the ROC survival curves and AUC values (Figures 4F, I). The heatmap illustrated the consistent expression profile of the training cohort (Figures 4E, G). Scatter plots indicated worse survival in the high-risk group than in the low-risk group in both the testing cohort and whole cohort, and risk score distribution plots verified a higher risk score in the high-risk group (Figures 4H, K, J, L). Furthermore, the constructed signature was evaluated using the TCGA-LGG, TCGA-GBM cohorts and the CGGA cohorts, the results indicated that the signature was an excellent predictor of survival of glioma patients (Supplementary Figures S1, S2). In conclusion, the SRlncRNA prognostic signature we constructed exhibited good predictive value.

To further assess the role of SRlncRNAs in glioma prognosis, we analyzed the correlation between SRlncRNAs and clinicopathological factors. As shown in Figure 5, the heatmap suggested that IDH mutation status, 1p19q codeletion status, grade, age, and survival status were significantly different in different groups (Figure 5A). The survival status was negatively correlated with a higher risk score (Figure 5B), followed by a positive correlation between glioma grade and a higher risk score (Figure 5C), and there was a lower risk and a better prognosis when an IDH mutation occurred (Figure 5D), as well as a lower risk with 1p19q codeletion (Figure 5E). In conclusion, there is a significant correlation between our constructed glioma risk score and the clinical factors of glioma.

Univariate and multivariate Cox regression analyses were used to explore whether the aforementioned 6 SRlncRNA prognostic signature was an independent prognostic signature for glioma. In univariate Cox regression analysis (Figure 6A), the risk score hazard ratio (HR) was 1.147 (95% CI 1.126–1.169) (p < 0.001), and in multivariate Cox regression analysis (Figure 6B), it was 1.064 (95% CI 1.034–1.096) (p < 0.001). Using age, sex, stage, 1p19q codeletion status, and the risk score, a nomogram plot was constructed to predict 1-year, 3-year, and 5-year survival in glioma patients (Figure 6C). Before the construction of nomogram, Schoenfeld test was performed to examine the quality of factors (Supplementary Figure S4). The nomogram’s prediction capacity was demonstrated by the calibration curve, and the C-index was 0.8547 (Figures 6D–F).

In the GSEA-GO analysis (Figure 7A), the highly expressed SRlncRNAs were mainly concentrated in the Immune response to tumor cell, immunoglobulin production, synaptic transmission, NF-kB signaling pathway, whereas the SRlncRNAs with low expression were mainly concentrated in ligand gated ion channel signaling pathway, regulation of trans-synaptic signaling. Furthermore, in the GSEA-KEGG analysis (Figure 7B), mismatch repair, immunodeficiency, p53 signaling pathway, autoimmune thyroid disease was enriched in the high expression group, whereas mTOR signaling pathway, WNT signaling pathway, ERBB signaling pathway, phosphatidylinositol signaling system were enriched in the low expression group.

To further explore the correlation between riskscore and immune cells, functions and checkpoint, we quantified the enrichment scores of ssGSEA. We found that all items except three immune cells (DCs, NK_cells, Th1_cells) were significantly different (p < 0.01, Figure 8A), indicating a significant change in the immunophenotype in the high-risk and low-risk groups. We further explored the expression of immune function and checkpoint-related markers in two groups and found all markers were also significantly different (Figures 8B, C). By assessing the immune subtypes, we discovered that the high-risk group was mostly concentrated in the C4 subtype (Lymphocyte-depleted), was statistically significant, and that the low-risk group was primarily concentrated in the C5 subtype (immunologically quiet) (Figure 8D).

The CIBERSORT algorithm was utilized to further evaluate the differences in infiltration of 22 subtypes of immune cells. The relative information proportion of the 22 immune subtypes in the whole cohort and the correlation between the 22 immune subtypes are shown (Figure 9A). Consequently, we investigated the relationship between tumor-infiltrating immune cells and SRlncRNAs. Out of our results, we found that CD8_Tcells, CD4_Tcell memory activated, T-cell follicular helper, T-cell gamma delta, MacrophagesM0-M1-M2 and neutrophils were more abundant in the high-risk group (p < 0.05). In contrast, the levels of NK cells activated, Monocytes, and mast cells activated were lower in the high-risk group (p < 0.05) (Figure 9B). We also explored the correlation between the SRlncRNA riskscore and tumor-infiltrating immune cells using the ESTIMATE Databases. From the results of ESTIMATE, all items were significant between high-risk and low-risk groups (Figures 9C–F). Based on the aforementioned analysis, we discovered that two groups had a significant and varied immune infiltration pattern, which may have different survival benefits.

To further identify the expression levels of SRlncRNAs, we selected SNAI3-AS1 and its target gene for the follow-up analysis. We found that SNAI3-AS1 was significantly under expressed in LGG and GBM by using the Gene Expression Profiling Interactive Analysis (GEPIA) database (Figure 10A). The expression levels of SNAI3-AS1 were confirmed by qRT–PCR in normal human astrocytes and two glioma cell lines, as well as in tumors and adjacent normal tissues from 24 glioma patients (Supplementary Table S4). The results revealed that SNAI3-AS1 was significantly downregulated in glioma cell lines and glioma samples (Figures 10B, C, paired t-test). The low expression of its target genes in LGG and GBM was also validated by the GEPIA database (Figure 10D), and the quantity and intensity of the target genes’ immunohistochemical staining in pathological specimens (Figures 10E–G) agreed with the findings above.