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

An accumulation of studies has indicated aging to be a significant hazard factor for the development of tumors. Cellular senescence is positively associated with aging progress and aging-related genes (AGs) can regulate cellular senescence and tumor malignancy. While the association between AGs and the prognosis of patients with glioma is still unclear. In our study, we initially selected four survival-associated AGs and performed consensus clustering for these AGs based on The Cancer Genome Atlas (TCGA) database. We then explored the potential biological effects of four selected AGs. A prognostic risk model was constructed according to four selected AGs (LEP, TERT, PON1, and SSTR3) in the TCGA dataset and Chinese Glioma Genome Atlas (CGGA) database. Then we indicated the risk score was an independent prognostic index, and was also positively correlated with immune scores, estimate score, immune cell infiltration level, programmed death ligand 1 (PD-L1) expression, and expression of proinflammatory factors in patients with glioma. Finally, we performed the RT-qPCR and immunohistochemistry assay to validate our bioinformatics results. Thus, this study indicated the risk model was concluded to possibly have potential function as an immune checkpoint inhibitor and to provide promising targets for developing individualized immunotherapies for patients with glioma.


INTRODUCTION
Glioma is an aggressive brain tumor with high recurrence rates [1][2][3]. Universal treatment strategies for glioma involves surgical resection with postoperative chemotherapy and radiotherapy, but the clinical prognoses for patients suffering from glioma still remain poor due to its lethal malignancy [4][5][6][7]. Nowadays, a need still exists to find more advanced treatments for glioma.
Aging is an essentially universal characteristic of living organisms, is considered to involve a progressive decline of internal cellular functions and is a hot spot in tumor research [8][9][10]. Tumor, like the other diseases of aging, become much more prevalent beginning at around the midpoint of life. Cellular senescence plays significant role in contributing the aging progress and developing of tumor, while the mechanisms of it on tumor are extremely complex, which can both stimulate and suppress tumor malignancy [11][12][13]. For example, some in vivo experiments indicated cellular senescence could restrict tumorigenesis in early-stage prostate cancer and Braig et al. revealed H3K9me-mediated cellular senescence as a novel mechanism to suppress the AGING formation of lymphomas [14,15]. While some studies indicated an obvious ability of injected senescent fibroblasts to stimulate the proliferation of human epithelial tumor cells in immunocompromised mice [16,17], which was closely associated with the senescence-associated secretory phenotype (SASP) [18]. Aging-related genes (AGs) can regulate cellular senescence and play a key role in tumor malignancy [11,19]. There is, however, limited knowledge about the relationships between AGs and the prognosis of patients with glioma.
In our study, we selected four survival-associated AGs and performed consensus clustering for these AGs based on The Cancer Genome Atlas (TCGA) database. We then explored the potential biological effects of four selected AGs. Then a prognostic risk model was constructed, we indicated the risk score was an independent prognostic index, and was positively correlated with immune scores, estimate score, immune cell infiltration level, programmed death ligand 1 (PD-L1) expression, and proinflammatory factors expression in patients with glioma. Finally, we performed some laboratory experiment to validate our bioinformatics results. Our research revealed their underlying implication as biomarkers for predicting clinical prognosis of patients with 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).

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).

Consensus clustering analysis revealed the potential cellular biological effects of four survival-associated AGs
Since cluster2 presented the low OS, the malignancyrelated 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.

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  AGING 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).

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 AGING 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 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.

AGING
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.

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 AGING 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.
We evaluated the composition of 22 important immune fractions of samples based on the TCGA dataset using AGING 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

AGING
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.

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-LEPexpression 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.

DISCUSSION
There are currently many pieces of evidences of agedependent changes being correlated with proinflammatory nature and involved in chronic inflammatory microenvironment [22][23][24]. Increasing in inflammation and immune cell infiltration levels with age-dependent changes may contribute to formation and malignancy of tumors [25], while the functions of them in the malignancy of glioma are still unclear. Exploring the molecular mechanisms of AGs is important for identifying the role of age-dependent changes in glioma. Few studies have to date systematically investigated the molecular mechanisms of AGs in glioma and the association between the expression profile of AGs and the OS of patients with glioma.
In our study, we performed a systematic analysis to select the differentially expressed and survival-associated AGs. Consensus clustering analysis was used to sort glioma samples into two clusters according to the expression profile of four survival-associated AGs (TERT, LEP, PON1, and SSTR3). The cluster1/2 subtypes could affect prognosis and clinical characteristics of patients with glioma, and revealed strong correlations with some significant biological processes and signaling pathways, including human papillomavirus infection, the wnt signaling pathway, cellular senescence and the AMPK signaling pathway. We then established a risk model with the four selected AGs and showed that the risk score could be regarded as an independent prognostic index for predicting the prognosis of patients with glioma based on both training (TCGA) and validation (CGGA) datasets. In our risk model, TERT and LEP acted as risk factors, while PON1 and SSTR3 acted as protective factors. For TERT, it had been shown the promoter mutations were correlated with poor prognosis and shorter survival of patients with glioma [26,27], and some researches also identified that TERT promoter mutations was an independent prognostic factor in the other human cancer [28,29]. Yeh et al. found LEP expressed at higher levels in glioma than in astrocytes, and indicated that high levels of expression of LEP could promote glioma cell migration and invasion by increasing the MMP-13 expression levels [30]. Zhang et al. identified low LEP expression level was correlated with short OS and low complete remission rate in acute myeloid leukemia [31]. The molecular mechanisms of AGING

AGING
Interestingly, our Gene Ontology (GO) analysis results showed that neutrophil activation involved in immune response, neutrophil activation and neutrophil mediated immunity were enriched in cluster2. Neutrophils have been shown to have important functions for the innate immune response and to be involved in immune cell functions [35][36][37]. Therefore, we explored the association of risk score with immune score, PD-L1 expression, and compositional fraction of 22 types of immune cells in patients with glioma based on TCGA datasets. The PD-L1 expression levels and immune score of the high-risk subtypes were crucially higher than those of the low-risk subtypes, and positive associations between the risk score and the levels of most of the immune cells were observed: high-risk subtypes contained higher fractions of CD8 T cells, CD4 memory resting T cells, CD4 memory activated T cells, regulatory T cells, resting NK cells, M0 macrophages, M1 macrophages, M2 macrophages, activated dendritic cells, resting mast cells, eosinophils and neutrophils than did the low-risk subtypes. Our study also showed a significant positive association of the risk score with immune cell infiltration, and its correlation with the growth and differentiation of macrophages and T cells. Some SASP components play a crucial role in the development of inflammation, which is relevant to the potential functions of senescent cells in tumor malignancy [25,38], and the proinflammatory factors of the SASP can influence the immune cell infiltration levels [39,40]. So we also substantiated a crucially positive correlation of the risk score with IL-1α, IL-1β, IL-6, IL-8, and IL-18 expression levels in patients with glioma. We substantiated the idea of TERT and LEP being the risk factors, and LEP and PON1 as protective factors, affecting the immune cell infiltration level and the levels of proinflammatory factors in glioma. Jiang et al. showed a correlation between TERT alterations and tumor immune microenvironment, and thus suggested that TERT can serve as potential biomarkers for individualized immunotherapy [41]. LEP, shown to play a crucial role in the initiation of the immune system, had also been indicated to be one of the mediators of inflammation responsible [42,43]. Aharoni et al. showed PON1 to have an important function as an anti-inflammatory, which could reduce sustained pro-inflammatory responses [44]. While there are few studies reported the association between the SSTR3 expression and immune microenvironment in human cancer.
In our study, we developed a risk model with four selected AGs, and the results revealed that the risk score had a prognostic value and was associated with tumor immune microenvironment in glioma. Nevertheless, there were still some limitations of our study that needed to be overcome. First, the bioinformatics results were validated based on the TCGA and CGGA datasets The relationship between the risk signature of four survival-associated AGs and clinical prognostic value of patient with glioma were not subjected to external verification due to the lack of our own adequate available data; external validation should be performed based on our own data in the future. Second, we only performed the RT-qPCR and immunohistochemistry assay to substantiate our bioinformatics results; we should further perform more laboratory experiments to confirm this conclusion.
In conclusion, we established a risk model with four selected AGs, which showed potential function as immune checkpoint inhibitors and hence promise as targets for developing individualized immunotherapy for patient with glioma.

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.
AGING 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: log2 (N+1).

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.

AUTHOR CONTRIBUTIONS
Haitao Luo was involved in the conception and design of the experiment, Chuming Tao helped perform data analysis, and Xiaoyan Long performed the RT-qPCR and immunohistochemistry assays. Kai Huang drafted the manuscript, and Xingen Zhu contributed to the revision. All authors edited and approved the final version of manuscript.