A risk signature of four aging-related genes has clinical prognostic value and is associated with a tumor immune microenvironment in glioma

Selection the AGs association with the prognostic of patients with glioma

To find differentially expressed AGs, we initially selected 676 differentially expressed genes based on TCGA database (Figure 1A, 1C). Then we identified 4 differentially expressed AGs from 676 differentially expressed genes: LEP and TERT were upregulated while PON1 and SSTR3 were downregulated (Figure 1B, 1D). We then found that these AGs were correlated with the prognosis of patients with glioma (P < 0.01, Figure 1E). We also identified that the frequency of these survival-associated AG genetic alterations (< 1.8%, Figure 1F).

Figure 1. Selection of differentially expressed and survival-associated AGs based on TCGA datasets. (A) Heatmap of all genes with significant differences between lower-grade glioma (LGG) and glioblastoma (GBM) samples. (C) Volcano plot of all genes based on TCGA datasets, light blue represented downregulated of genes, and pink represented upregulated of genes. (B) Differentially expressed AGs are showed in heatmap between LGG and GBM samples. (D) Volcano plot of the four selected differentially expressed AGs, light blue represented downregulated of AGs, and pink represented upregulated of AGs. (E) Forest plot of the four differentially expressed AGs. (F) Genetic changes of the four survival-associated AGs.

Consensus clustering for four survival-associated AGs and with the prognoses of patients with glioma

To explore the association of four survival-associated AGs, we performed correlation analysis according to their mRNA expression level in the TCGA datasets. Our results revealed the expression of PON1 was crucially positive associated with SSTR3 in glioma, while there were a crucial negative association between the expression of PON1 and TERT, and the expression of LEP was negatively associated with PON1 and SSTR3 (Figure 2A). Consensus clustering analysis was used to sort samples into subtypes based on the expression profiles of the above-identified four survival-associated AGs in the TCGA datasets. The resulting cumulative distribution function (CDF) curves and SigClust analysis indicated a K value of = 2 (Figure 2B, 2C and Supplementary Figure 1), categorization of two subtypes (cluster1 (n = 317) and cluster2 (n = 311)) on the basis of different expression levels of the four survival associated AGs. The expression levels of the four selected AGs were statistically different between the two subtypes, which showed cluster2 with the upregulated expression levels of the risk factors (TERT and LEP) and clsuter1 with the upregulated expression of protective factors (SSTR3 and PON1, Figure 2D). Then we further revealed the gene expression profiles between the two subtypes were differentiated well by using principal component analysis (PCA, Figure 2E). The Kaplan-Meier (KM) curves indicated a poor prognosis for the samples in cluster2 (P < 0.001, Figure 2F). Furthermore, high grade, old age, mutant-type isocitrate dehydrogenase (IDH) status,1p19q non-codeletion status were presented in cluster2 than cluster1 (Figure 2D and Supplementary Table 1).

Figure 2. Two categories of patients based on distinct clinical characteristics and the OS according to the expression levels of the four survival-associated AGs in the TCGA datasets. (A) Spearman’s correlation analysis of the four survival-associated AGs (LEP, TERT, PON1, and SSTR3) in the TCGA datasets. (B) The relative change of the area under the cumulative distribution function (CDF) for K =2-9. (C) The consistency clustering CDF curve for K = 2-9. (D) Clinical characteristics of two clusters identified based on the expression levels of the four survival-associated AGs in the TCGA datasets. (E) Principal component analysis (PCA) for the total mRNA expression profile based on TCGA datasets. (F) Kaplan-Meier (KM) curves for 628 cluster1 and cluster2 glioma patients based on the TCGA datasets. * P < 0.05, ** P < 0.01, and *** P < 0.001.

Consensus clustering analysis revealed the potential cellular biological effects of four survival-associated AGs

Since cluster2 presented the low OS, the malignancy-related mechanisms were further explored in this subtype. We selected differentially expressed genes between cluster2 and cluster1, and analyzed some biological processes significantly correlated with cluster2. Compared with cluster1, genes whose products are involved in neutrophil degranulation, neutrophil activation, ras protein signal transduction and regulation of cell morphogenesis were positively enriched in cluster2 (Figure 3A). A Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis also revealed a crucial positive enrichment of genes involved in human papillomavirus infection, the wnt signaling pathway, cellular senescence, and the AMPK signaling pathway (Figure 3B). In addition, we further showed malignant hallmarks by carrying out a gene set enrichment analysis (GSEA), which indicated that the IL6 JAK STATS signaling (NES = 1.53, normalized P = 0.039), apoptosis (NES = 1.60, normalized P = 0.013), DNA repair (NES = 1.89, normalized P = 0.004), and G2M checkpoint (NES = 1.68, normalized P = 0.03) were significantly positively associated with cluster2 (Figure 3C). In conclusion, our results might provide novel insights for cellular biological function related to four survival-associated AGs.

Figure 3. The potential biological functions of four selected associated-survival genes. (A, B) Functional annotations of differentially expressed genes in clsuter2 compared with cluster1 based on TCGA datasets determined from Geno Ontology (GO) and Kyoto Encyclopedia of Gene and Genomes (KEGG) pathway analyses. (C) Malignancy hallmarks positively enriched in cluster2 determined using gene set enrichment analysis (GSEA) in the TCGA datasets.

Identification of the prognostic value of selected AGs and risk model derived from four select AGs

The least absolute shrinkage and selection operator (Lasso) Cox regression algorithm was used to build a prognostic risk model according to the expression levels of the four selected AGs in the TCGA datasets, and coefficients were obtained to calculate the risk scores for each patient with glioma (Supplementary Figure 2). We sorted glioma samples into two subtypes by the median risk scores. The KM curves revealed that the samples in the high-risk categories had a poorer OS than did those in the low-risk categories based on training and validation databases (P < 0.001, Figure 4A, 4B). The area under the time-dependent receiver operating characteristic (ROC) curve (AUC) values were calculated to assess the value of our four-AG risk model which was 0.805 for the TCGA dataset, 0.801 for the Chinese Glioma Genome Atlas (CGGA) datasets (Figure 4C, 4D). These results showed the accuracy of our four-AGs risk model for glioma prognosis. The risk plot distribution, survival status of patients with glioma, and heatmap of the expression of included genes were determined based on the TCGA and CGGA databases (Figure 4E, 4F).

Figure 4. Construction of a risk model and the association of the risk model with clinical characteristics of patients with glioma. (A, B) Kaplan-Meier (KM) curves for overall survival (OS) prediction based on the training (TCGA) and validation (CGGA) datasets. (C, D) Receiver operating characteristic (ROC) curves for the risk model both in the training and validation datasets. (E, F) Risk plot distribution, survival status of patients, and heatmap of expression of included genes in the training and validation datasets. (G–K) Relationships between the risk score and clinical characteristics (grade, age, 1p19q status, IDH status, and gender) of patients with glioma. Non-significant (ns) P > 0.05, * P < 0.05, ** P < 0.01, and *** P < 0.001.

Positive correlations of the risk model with the OS and clinical characteristics of patients with glioma

The correlations between the risk scores and clinical characteristics of patients with glioma were examined using the Wilcoxon test, which showed significant differences between risk scores in groups stratified by the Word Health Organization (WHO) tumor grade (P < 0.001), age (P < 0.001), 1p/19q status (P < 0.001) and IDH status (P < 0.001, Figure 4G–4K). Univariate Cox regression analyses revealed that the WHO grade, age, IDH status, 1p/19q status, and risk score were significantly correlated with prognosis of patients in the TCGA datasets. Multivariate Cox regression analyses showed that the WHO grade (P < 0.01), IDH status (P < 0.05), and risk score (P < 0.05) remained crucially associated with the OS of patients (Figure 5A). Similar conclusions were made based on the results of the multivariate Cox regression analysis using the validation dataset (Figure 5B), which revealed positive associations of the WHO grade (P < 0.001), IDH status (P < 0.01), 1p19q status (P < 0.001), and risk scores (P < 0.05) with prognosis. Therefore, the risk scores derived from the four selected AGs indicated a good prognostic performance for this set of AGs and that they could represent an independent prognostic index for glioma. A nomogram was then constructed for providing a prognosis of glioma that integrated tumor grade, IDH status, and risk scores based on the TCGA database (Figure 5C), and the C-index for survival prediction was 0.828. The calibration plots for providing prognoses of patients at 2-, 3-, and 5- years revealed an optimal agreement between the nomogram prediction and the actual observed outcomes (Figure 5D–5F). The results showed this four-AGs risk model to have accurate predictive value for prognosis and clinicopathological features.

Figure 5. The risk score could serve as an independent prognostic index and predict the prognoses of patients with glioma. (A, B) Univariate and multivariate Cox analysis of clinical characteristics and molecular features based on the training (TCGA) and validation (CGGA) datasets. (C) Nomogram based on the WHO grade, IDH status and risk score for providing prognoses of patients with glioma in the TCGA datasets. (D–F) Calibration plots used to indicate that the nomogram could effectively predicted the 2-,3-, and 5- years prognoses of patients with glioma.

The association between the risk model and the immune cell infiltration level in glioma

To explore the occurrence of any association between the risk score, immune score, and estimate score, the estimate formula was performed to calculate the immune score and estimate score of patients with glioma based on the TCGA dataset. The high-risk subtype was associated with a higher immune score and estimate score than was the low-risk subtype (P < 0.001, Figure 6A, 6B). Furthermore, our results also showed crucially positive associations of risk score with immune score (R = 0.43, P < 0.001, Figure 6C) and estimate score (R = 0.46, P < 0.001, Figure 6D) in glioma samples.

Figure 6. Relationships between the risk score and immune cell infiltration levels, and PD-L1 expression levels based on TCGA datasets. (A, B) Comparison of immune and estimate scores according to two subtypes grouped by the median risk score. (C, D) Relationship between the risk score and immune scores and the estimate score. (E) The violin plot showed distribution of the ratio differentiation of 22 kinds of immune cells in glioma samples according to two subtypes grouped by the median risk score. (F–J) Distribution of PD-L1 expression levels based on the risk score (F) and expression levels of the four selected AGs (LEP (G), TERT (H), PON1 (I), and SSTR3 (J)) of patients with glioma in the TCGA datasets. Non-significant (ns) P > 0.05, * P < 0.05, ** P < 0.01, and *** P < 0.001.

We evaluated the composition of 22 important immune fractions of samples based on the TCGA dataset using the CIBERSORT formula to explore the compositional fraction of 22 types of immune cells between high-risk and low-risk subtypes (Figure 6E). The results revealed that patients with glioma in the high-risk subtype included a higher fraction of memory B cells (P < 0.05), CD4 memory resting T cells (P < 0.001), CD4 memory activated T cells (P < 0.001), regulatory T (Tregs) cells (P < 0.001), resting NK cells (P < 0.001), M0 macrophages (P < 0.001), M1 macrophages (P < 0.001), M2 macrophages (P < 0.001), activated dendritic cells (P < 0.01), resting mast cells (P < 0.001), activated mast cells (P < 0.001), eosinophils (P < 0.01), and neutrophils (P < 0.001) than did those in the low-risk subtype. While naïve B cells, plasma cells, naïve CD4 T cells, activated NK cells, monocytes, and activated mast cells showed the opposite result (P < 0.001). Our results showed that the risk score was statistically positively correlated with immune cell infiltration levels in glioma.

Association of PD-L1 expression with the risk score and the four selected AGs

To explore the relationship between PD-L1 and the four selected AGs, expression levels of PD-L1 were calculated for the high-risk and low-risk subtypes based on the TCGA database. The PD-L1 expression was found to be upregulated in the high-risk subtypes compared to that in the low-risk subtypes (Figure 6F). Furthermore, the PD-L1 expression levels showed a positive correlation with risk factors (LEP and TERT, Figure 6G, 6H), while a crucially negative relationship was noted with protective factors (PON1 and SSTRS, Figure 6I, 6J).

Association of the six types of immune cells with the four selected AGs in the risk model

According to the associations of the risk score and the six types of immune cells (Figure 7A), the relationships between these immune cells and the four selected AGs were investigated. Glioma patients of the subtype expressing high levels of LEP had a lower fraction of naïve B cells, plasma cells, NK activated cells, and monocytes cells than did those of subtype expressing low levels of LEP, while the opposite was the case for CD4 memory resting T cells and M0 macrophages (P < 0.05, Figure 7B), in accord with the results grouped by the risk score. Similarly, the results for TERT, except for the CD4 memory resting T cells and NK activated cells (P > 0.05), were also in accord with the risk score results (P < 0.01, Figure 7C). But the results for PON1 and SSTR3, as protective factors, were not in accord with the risk score results; glioma patients of the high-LEP-expression subtype included a higher fraction of naïve B cells, plasma cells, NK activated cells, and monocytes than did those of low-LEP-expression subtype, while CD4 memory resting T cells and M0 macrophages showed the opposite result (P < 0.05, Figure 7D). Similarly, the results for SSTR3, except for M0 macrophage and monocytes (P > 0.05), were in contrast with the risk score results (P < 0.01, Figure 7E). The four selected AGs (LEP, TERT, PON1, and SSTR3) were concluded to have important functions in immune cell infiltration of glioma.

Figure 7. Association of six types of immune cells with the risk score and four selected AGs in the TCGA datasets. (A) Distributions of the six types of immune cells (naïve B cells, plasma cells, CD4 memory resting T cells, NK activated cells, monocytes, and M0 macrophages) in the two subtypes grouped by the median risk score. (B–E) Comparison of the six types of immune cells (naïve B cells, plasma cells, CD4 memory resting T cells, NK activated cells, monocytes, and M0 macrophages) according to LEP (B), TERT (C), PON1 (D), and SSTR3 (E) expression levels. Non-significant (ns) P > 0.05, * P < 0.05, ** P < 0.01, and *** P < 0.001.

Association of proinflammatory factors with the risk score and the four selected AGs

An accumulation of research has indicated chronic inflammation associate with cellular senescence to have an important function in immune cell infiltration and major proinflammatory factors, including interleukin-1α (IL-1α), interleukin-1β (IL-1β), interleukin-6 (IL-6), interleukin-8 (IL-8) and interleukin-18 (IL-18) [20, 21]. In our study, we explored the association of major proinflammatory factors with the risk score and the four selected AGs, and showed the expressing features of eight major proinflammatory factors based on the TCGA datasets (Supplementary Figure 3). Our results revealed that the levels of expression of IL-1α, IL-1β, IL-6, IL-8, and IL-18 in the high-risk groups were statistically higher than those in the low-risk groups (P < 0.001, Table 1), in accord with the risk score results. This study also revealed that IL-1α, IL-1β, IL-6, IL-8, and IL-18 were expressed at higher levels in the high-LEP-expression group than in the low-LEP-expression group (P < 0.05, Table 1). And the results for TERT showed that, except for IL-1α (t = 1.640, P = 0.104), IL-1β (t = -2.646, P < 0.001), and IL-18 (t = -2.827, P < 0.001), the other interleukins were expressed at high levels in the high-TERT-expression group (P < 0.05, Table 1). The PON1 and SSTR3, as protective factors, the results showed that IL-1α, IL-1β, IL-6, IL-8, and IL-18 were expressed at lower levels in the high expression group than in the low expression group (P < 0.05, Table 1).

The mRNA and protein expression levels of four selected AGs in glioma

We performed the Real-time quantitative polymerase chain reaction (RT-qPCR) and immunohistochemistry assay to validate the bioinformatics results. The RT- qPCR assay showed that the four AGs (LEP, TERT, PON1, and SSTR3) were expressed to different extents in normal brain tissue (NBT), lower-grade glioma (LGG) and Glioblastoma (GBM) tissue in mRNA expression level (Figure 8A–8D), and the immunohistochemistry assay revealed that the protein expression of four AGs were also expressed to different extents in NBT and glioma tissue (Figure 8E–8H), in accord with the bioinformatics results.

Figure 8. Validation of the bioinformatics results using RT-qPCR and immunohistochemistry assay. (A–D) Comparison of LEP (A), TERT (B), PON1 (C), and SSTR3 (D) mRNA expression levels in normal brain tissue (NBT), lower-grade glioma (LGG), and glioblastoma (GBM) tissues by RT-qPCR assay. (E–H) Comparison of the protein expression of LEP (E), TERT (F), PON1(G) and SSTR3 (H) in NBT and glioma tissue by immunohistochemistry assay. Non-significant (ns) P > 0.05, * P < 0.05, ** P < 0.01, and *** P < 0.001.

Data collection and acquisition

The mRNA expression files and related clinicopathological data were acquired from the TCGA (n = 628) (https://portal.gdc.cancer.gov/) and CGGA (https://www.cgga.org.cn/) (n = 620)) datasets (Supplementary Table 2). A list of AGs was acquired from the human aging genome resource dataset (http://genomics.senescence.info/genes/) [45], and copy number alterations data were acquired from cBioPortal websites (http://www.cbioportal.org/) [46].

Selected survival-associated AGs and investigation of their potential cellular biological effects

We identified differentially expressed genes between LGG and GBM samples, and selected the differentially expressed AGs from differentially expressed genes based on TCGA datasets, according to the standards of | log2(Fold change) | > 1 and P < 0.05. We selected survival-associated AGs according to the standard of P < 0.01 and | hazard ratio | > 1 by univariate Cox analysis using the R package “survival”. Then we sorted patients with glioma into two clusters basing on the four selected AGs expression profiles of samples from TCGA datasets using the R package “ConsensusClusterPlus”. The Euclidean distance was used to compute the similarity distance between patients with glioma, and the k-means method was utilized for clustering based on 50 iterations, which each iteration included 80% of patients. The optimal number of clusters was identified by CDF and consensus matrices. Then we validated the classification results by performing PCA. The molecular mechanisms of selected survival-associated AGs were explored by performing GO pathway analyses, KEGG pathway analyses, and GSEA.

Establishing a risk model and exploring the prognostic value of risk score

To further investigate the biological functions of the four selected survival-associated AGs, we performed Lasso Cox regression algorithm based on the expression levels of survival-associated AGs, and four selected AGs were identified based on the minimum criteria to establish the risk score (Supplementary Figure 2). Each risk score was defined by the coefficients derived from Lasso algorithm, and the formula was as followed: risk score = [TERT expression * (0.226371)] + [LEP expression * (0.30102)] + [PON1 expression * (-0.34145)] + [SSTR3 expression * (-0.08439)]. We calculated the risk score for glioma sample both in training and validation datasets. We then sorted the glioma samples into two high-risk and low-risk subtypes. A cluster heat map was constructed to display the correlations between candidate genes and risk scores. And ROC curves were constructed to evaluate the prediction efficiency of the risk model. We also constructed nomograms using the R package “rms” and evaluated their performance. The relationships between clinical characteristics and risk score were explored both in training and validation datasets.

Determination of associations between the level of immune infiltrates, PD-L1 expression, proinflammatory factors, and the prognostic model

We evaluated immune scores of glioma samples using the estimate formula [47]. We calculated the immune and estimate scores according to gene expression profiles and determined their associations with the risk score and the four selected AGs. An estimation of the compositional fraction of 22 types of immune cells of each sample was calculated using CIBERSORT (http://cibersort.stanford.edu/), a tool developed to calculate cellular compositions of tumors according to gene expression data in the TCGA datasets [48]. The relationship between the prognostic model, PD-L1 expression and some proinflammatory factors were investigated based on TCGA datasets, and we normalized the gene expression in the TCGA datasets using formula: log₂ (N+1).

Use of RT-qPCR and immunohistochemistry assay to validate bioinformatics results

We collected NBT and glioma tissues specifically 6 NBTs, 10 LGG tissues and 9 GBM tissues from the Second Affiliated Hospital of Nanchang University from May 2017 to June 2020 (Supplementary Table 3). Our study was approved by the Ethics Committee of this hospital. RT-qPCR was conducted using a LightCycler® 480 real-time PCR System based on the manufacturer’s instructions. The expression levels of the four selected AGs were calculated using the 2-ΔΔCt method. Primer sequences for the four AGs were as follows: LEP forward 5’-GCAGTTGCGCAAGTTGTGAT-3’ and reverse 5’-GATGGGCTTCTTGGGCCTTG-3’; TERT forward 5’- CTTGCGGAAGACAGTGGTGA -3’ and reverse 5’-TCCGGGCATAGCTGGAGTA-3’; PON1 forward 5’-CACCCGATGGCAAGTATGTCT-3’ and reverse 5’-GGCATCCAACCCAAAGGTCT-3’; SSTR3 forward 5’-CCCCATGGGCAGGCAAATA-3’ and reverse 5’-CGAGGAGGCATTCTCAGGTT-3’. We performed the immunohistochemistry assay on human tissues by methods described previously [49].

Statistical analyses

KM curves were used to contrast OS between pairs of subtypes. The Lasso Cox regression formula was used to establish a risk model. Wilcoxon test was used to compare the risk score and four AGs between pairs of subtypes in the different clinical characteristics and PD-L1 expression. Univariate and multivariate Cox regression analyses were used to determine the independent prognostic value of the risk score and the nomogram. The statistical analyses were carried out using R programming language v3.6.3, SPSS Statistics software 26.0, and Prism 8.0.

Data availability statement

The data for the study were acquired from the TCGA and CGGA datasets and the human aging genome resource.