Prognostic Correlation of Autophagy-Related Gene Expression-Based Risk Signature in Patients with Glioblastoma

Introduction

Autophagy is an important selective recycling mechanism by which cell components are degraded in lysosomes to provide basic materials and energy for cells.¹ The autophagy process has a vital role in cancer initiation, malignant progression, and resistance to treatment.² Further, regulation of cancer stem cell homeostasis by autophagy is an important mechanism of tumor resistance, recurrence, and metastasis.^3,4 The discovery of numerous of autophagy-related genes has facilitated the discovery of the pathophysiological roles of autophagy.⁵ To illustrate the importance of autophagy, there is increasing evidence that abnormalities in autophagy-related genes play pathogenic roles in human disease, including cancer.¹ For example, Takamura et al found that various liver tumors developed in a mouse model with deletion of autophagy-related gene 5 (Atg5), which exhibited accumulation of p62, due to deficiency in autophagy, promoting tumor progression.⁶

Various relationships between cancer and autophagy have been described. FIP200 interacts with its autophagy partner, ATG13, to regulate autophagy activity in cancer, which is essential for tumor cell growth.^7,8 Further, deletion of the two autophagy regulators, PARKIN and BNIP3, leads to tumor formation and accelerated malignant progression.^9,10 Hence, understanding the role of autophagy and related genes in cancer is of great significance for the development of new therapeutic strategies.

Glioblastoma multiforme (GBM) is the most common malignant primary tumor of the central nervous system.¹¹ The current standard treatment for glioblastoma is surgical resection combined with chemoradiotherapy; however, this approach only extends median survival to 14.6 months.^11,12 According to the 2016 WHO guidelines for classification of central nervous system tumors, glioblastoma classification combines molecular phenotype with traditional pathology.¹³ Glioblastomas are divided into two groups based on IDH mutation status: IDH-wildtype and IDH-mutant GBM, with significant differences in clinical features and prognosis between the two groups.

In this study, we comprehensively analyzed autophagy-related genes in GBM. First, we collected autophagy-related gene data and identified 14 genes related to patient prognosis. These 14 genes could stratify patients into two groups with significant molecular and prognostic differences, suggesting a correlation with glioblastoma malignancy. Hence, we constructed 14 gene autophagy-related signature with good prognostic value and which can serve as an independent prognostic indicator in GBM. Moreover, functional analysis revealed the elevation of immune and inflammatory responses in GBM samples classified as high-risk using our 14 gene signature, which was directly related to tumor malignancy and poor prognosis. In conclusion, our study confirms the correlation between autophagy status and glioblastoma malignancy and prognosis, providing new insights relevant to individualized treatment.

Methods

Samples and Data Collected for This Study

Glioblastoma samples (n = 138) with mRNA sequencing data and clinicopathological information were collected from the Chinese Glioma Genome Atlas (CGGA) dataset (http://www.cgga.org.cn). Sample pathology was diagnosed by two experienced neuropathologists, according to the 2016 WHO Classification of CNS tumors. Overall survival (OS) was defined as the interval from diagnosis to death or last follow-up. In addition, mRNA expression and survival data from 689 (152 sequencing and 537 microarray included) glioblastoma samples were downloaded from The Cancer Genome Atlas (TCGA) dataset (https://cancergenome.nih.gov/), along with data from 159 glioblastoma samples in the GSE16011 dataset.¹⁴ Thus, mRNA expression and survival data from a total of 986 glioblastoma samples were included in this study. Our CGGA dataset was approved by the Beijing Tiantan Hospital Capital Medical University Institutional Review Board (IRB) and was conducted according to the principles of the Helsinki Declaration. Written informed consent was obtained from all patients. The specimens were collected under IRB KY2013-017-01 and were frozen in liquid nitrogen within 5 min of resection.

Human Autophagy-Related Gene Set

Genes (n = 232) involved directly or indirectly in autophagy assembled in the Human Autophagy Database (HADb, http://autophagy.lu/clustering/index.html), as previously described,¹⁵ were combined with 394 autophagy-related genes acquired from the GO_AUTOPHAGY gene set in the Molecular Signatures Database v6.2 (MSigDB, http://software.broadinstitute.org/gsea/msigdb); by eliminating overlapping genes in the two databases, an autophagy-related gene set comprising 531 genes was finally constructed.

Gene Screening and Gene Signature Development

We performed univariate Cox regression with 531 autophagy-related genes in two datasets (138 RNA-seq data from CGGA and 152 RNA-seq data from TCGA). Fourteen genes were significantly associated with overall survival in two datasets (p value < 0.05). We checked collinearity between the genes by calculating the variance inflation factor (VIF) and tolerance statistic. All VIF values were below the commonly suggested cut-off of 10 and all tolerance values were greater than 0.1. And the regression coefficients for each gene calculated using the CGGA dataset were used to build a gene signature, as follows:

where “expr” indicates gene expression level and “β” indicates the regression coefficient for each gene.

A risk score was developed by a linear combination of the expression levels of the 14 genes, weighted by regression coefficients calculated, by univariate Cox regression analysis for CGGA data.

Patients were assigned into high- or low-risk groups using the median calculated risk score as a cut-off value.

Gene Ontology and Gene Set Enrichment Analysis (GSEA)

Genes significantly correlated with risk score (Pearson |R| > 0.4, p < 0.05) in the CGGA dataset were filtered by Pearson correlation analysis. Then, risk score positively or negatively associated genes were chosen to perform gene ontology (GO) analysis in DAVID (https://david.ncifcrf.gov/).¹⁶ Gene set enrichment analysis (GSEA) was conducted to identify gene sets that differed significantly between the high- and low-risk patient groups, using GSEA 3.0 software (http://www.broadinstitute.org/gsea/index.jsp).¹⁷ Max gene set size was 500 and min size was 15. Number of permutations was 1000. In the analysis results, it was generally believed that the pathways were significantly enriched when |NES| > 1, NOM p-value < 0.05, and FDR q-value < 0.25.

Estimate

“Estimation of STromal and Immune cells in MAlignant Tumours using Expression data” (ESTIMATE) is a method to predict the ratio of stromal and immune cells in the bulk tumor with gene expression signatures.¹⁸ Using the ESTIMATE R package (https://sourceforge.net/projects/estimateproject/), we calculated the immune score, stromal score, ESTIMATE score as well as tumor purity of all samples from CGGA and TCGA datasets.

Statistical Analysis

R language (version 3.5.1, https://www.r-project.org/) was the main data analysis and mapping tool. The 14 genes were clustered using the “cluster” function provided in the R package pheatmap. The “Coxph” function in the survival package was used for univariate and multivariate Cox regression analysis. Other figures were drawn using R packages, including ggpubr, pROC,¹⁹ gglpot2, Hmisc, estimate¹⁸ and circlize. Kaplan-Meier estimates were used for survival analysis, with a two-sided log-rank test. Differences in pathological features between the two patient groups were evaluated using the Chi-square test. Values of p < 0.05 were considered significant.

Results

Gene Screening in GBM

We analyzed 531 autophagy-related genes, including 232 genes from the Human Autophagy Database (HADb) and 394 genes from the GO_AUTOPHAGY gene set from MSigDB. To screen for genes associated with clinical prognosis, univariate Cox regression analysis was performed in two datasets (138 RNA-seq data from CGGA and 152 RNA-seq data from TCGA, respectively). Finally, 14 genes were identified as significantly associated with OS (p < 0.05). Clustering of the 138 patients with GBM from the CGGA dataset with expression profiles for the 14 genes identified two robust groups (Figure 1A). Survival analysis indicated a significant difference in OS between the two groups (Figure 1B, p = 0.0023). Validation analysis using data from 152 patients with GBM from TCGA dataset generated similar results (Figure S1). These findings suggest that the 14 genes associated with autophagy have strong prognostic value; therefore, these genes were chosen for inclusion in a prognostic gene signature. Next, we generated a heatmap demonstrating the relationship between gene signature risk score and the 14 genes (Figure 1C and Figure S2). The trend was striking: the expression of the top nine genes (MTMR14, LENG9, P4HB, TCIRG1, HSPA5, DRAM1, CTSD, S100A8 and CCL2) increased with signature risk score and the bottom five genes (MSTN, UBQLN4, PHYHIP, RRAGB and ZKSCAN3) was opposite. Hazard ratio (HR) values for the 14 genes were calculated by univariate Cox regression from the CGGA dataset (Table 1). HR values for the top nine genes were >1, while those for the other five genes were <1. These suggested that the potential prognostic value of these 14 genes in GBM remained to be further investigated.

Table 1 The HR and p values of 14 Genes

Figure 1 Fourteen autophagy-related genes were identified. (A) COX univariate analysis screened out 14 genes and grouped CGGA glioblastoma patients into two groups. (B) The prognosis of the two groups was significantly different (p = 0.0023). (C) The expression pattern of 14 genes in CGGA database. The top nine genes increased with risk score, while next five genes were reversed. (D) The coefficients of 14 genes calculated in CGGA dataset.

Construction of a Signature Predicting Prognosis in GBM Using the 14 Autophagy-Related Genes

Next, the 14 selected genes were used to construct a signature. We calculated the coefficients of the 14 genes using univariate Cox regression analysis of the CGGA dataset (Figure 1D). The algorithm was as follows: signature risk score = (MTMR14*1.61) + (LENG9*0.63) + (P4HB*0.44) + (TCIRG1*0.39) + (HSPA5*0.37) + (DRAM1*0.30) + (CTSD*0.24) + (S100A8*0.17) + (CCL2*0.17) + (MSTN*–0.19) + (UBQLN4*–0.52) + (PHYHIP*–0.72) + (RRAGB*–0.85) + (ZKSCAN3*–6.64). To evaluate the classification efficacy of the signature, we divided CGGA patients into high- and low-risk groups, based on the median risk scores, and observed significant clinical and molecular differences between the two groups (Figure 2A and Table 2). In the high-risk group, patients tended to be older than those in the low-risk group (p < 0.001). Further, more classical and mesenchymal tumor subtypes were found in the high-risk (94%) than in the low-risk (46%) group (p < 0.001). In addition, we found that IDH wild-type patients accounted for a large proportion in the high-risk group (97%), while IDH wild-type patients comprised only 55% in the low-risk group (p < 0.001). Moreover, 58% and 29% of samples carried MGMT promoter methylation in the low- and high-risk groups, respectively (p < 0.01), while 74% and 46% of low-risk and high-risk patients, respectively, had received chemotherapy (p < 0.01).

Table 2 Clinicopathological Factors of Glioblastoma Patients in CGGA Cohorts

Figure 2 Correlation between signature risk score and clinicopathology in CGGA dataset. (A) The clinicopathologic information of 138 glioblastoma patients in CGGA. Patients were arranged by the increasing risk score and divided into two groups based on median risk score. (B) Distribution of the risk score in patients stratified by IDH1 status, MGMT promoter methylation, and TCGA subtypes. (C) ROC curve was utilized to evaluate the predictive efficacy of risk score for IDH1 mutation, MGMT promoter methylation, and TCGA mesenchymal subtypes. **P < 0.01, ***P < 0.001, ****P < 0.0001. Abbreviations: CGGA, Chinese Glioma Genome Atlas; IDH1, isocitrate dehydrogenase 1; MGMT, methylguanine methyltransferase; TCGA, The Cancer Genome Atlas.

Next, we divided patients into groups according to their molecular characteristics and compared the risk scores between the groups (Figure 2B). We observed significantly higher risk scores in IDH wild-type patients and our signature could accurately predict the IDH mutation status (AUC = 0.894, specificity = 0.838, sensitivity = 0.909, Figure 2C). In MGMT promoter unmethylated patients and mesenchymal subtype patients, we observed the same situation.

Using coefficients derived from the CGGA dataset, we formulated a risk signature score for each patient in TCGA dataset as a validation cohort. Consistent with the above results, there were significant differences between high- and low-risk groups (Figure S3 and Table S1). In summary, patients classified as high risk using our 14-gene signature had more malignant characteristics than those classified as low risk.

Prognostic Value of Autophagy-Related Signature in GBM

Given the close relationship between the 14 autophagy-related gene signature and patient clinical characteristics, we sought to assess the clinical prognostic efficacy of the signature. We calculated risk scores for each patient using the 14 gene expression signature and divided them into high-risk and low-risk groups, based on the median risk score cut-off value. OS was shorter in the high-risk than the low-risk group (p = 0.00037, Figure 3A). According to the WHO CNS tumor classification (2016), glioblastoma (GBM, WHO IV) can be further divided into IDH-mutant and IDH-wild-type. Due to clinical and prognostic differences between the two subgroups, we next verified the prognostic value of our signature in the two subgroups. We found patients in high-risk group lived shorter than low-risk group both in IDH-wildtype GBM (p = 0.01) and IDH-mutant GBM (p = 0.037).

Figure 3 Prognostic significance of the autophagy-related gene signature in GBM. Patients were from CGGA (A), TCGA with RNA-seq data (B), TCGA with microarray data (C) and GSE16011 (D). Kaplan–Meier survival analysis showed that high-risk group conferred a significantly worse prognosis in GBM, especially in GBM with IDH wild-type. P-value shown in survival was determined by a log-rank test between the two groups.

In addition to CGGA as training cohort, TCGA sequencing data, microarray data and GSE16011 microarray data were included as the verification cohort (Figure 3B–D). Among IDH-wildtype patients, the prognosis remained significantly different (p < 0.05). And among IDH-mutant patients, the high-risk group had a worse survival trend. Similar survival curves in three verification cohorts further demonstrated the excellent prognostic value of our signature, especially in IDH-wildtype GBM.

The Autophagy-Related Signature Is an Independent Prognostic Indicator in Patients with GBM

To verify the independence of signature, we carried out univariate and multivariate COX regression analysis. Our result indicated that autophagy-related signature was significantly associated with overall survival (HR=1.921, 95% CI=1.013–3.644, p value=0.045), independent of clinical and pathological indicators, including age, gender, KPS, IDH1 status, MGMT promoter methylation, chemotherapy and radiotherapy (Table 3). Consistently, autophagy-related signature was verified as independent indicator in TCGA verification cohort (HR=1.718, 95% CI=1.122–2.629, p value=0.013).

Table 3 Variables Related to OS in Glioblastomas: Univariate and Multivariate Analysis

Functional Characteristics of the GBM Autophagy-Related Signature

To study alterations in signature-related functions, gene ontology (GO) analysis of signature-related genes was conducted using the CGGA dataset. Pearson correlation analysis was conducted to filter 3504 signature-related genes, including 1669 positively (Pearson R > 0.4, p < 0.05) and 1835 negatively (Pearson R < −0.4, p < 0.05) correlated genes. Positively correlated genes were mainly involved in the tumor immune environment, including “inflammatory response”, “immune response”, “extracellular matrix organization”, “chemotaxis”, and “leukocyte migration”, among others (Figure 4A, Table S2); however, negatively correlated genes were involved in normal physiological processes, such as “transcription, DNA-templated”, “nervous system development”, “neurotransmitter secretion”, and “regulation of exocytosis” (Figure 4B, Table S3).

Figure 4 Functional characteristics of autophagy-related signature in GBM. Gene Ontology (GO) analysis of biological processes was performed to analyze signature positively related genes (A) and negatively related genes (B). Gene set enrichment analysis (GSEA) revealed Go terms enriched in high-risk group (C) and low-risk group (D).

Further, we performed GSEA in high-risk and low-risk patients (Table S4, enrichment in high-risk group; Table S5, enrichment in low-risk group). Consistently, we found GO terms related to immunity were enriched in high-risk group, covering “leukocyte chemotaxis” (NES=2.161, FDR q-value=0.131), “inflammatory response” (NES=2.142, FDR q-value=0.053) and “positive regulation of cytokine production” (NES=2.109, FDR q-value=0.047, Figure 4C). Whereas Go terms including “regulation of synapse assembly” (NES= −2.049, FDR q-value=0.238), “regulation of synapse organization” (NES= −2.044, FDR q-value=0.168) and “metencephalon development” (NES= −1.980, FDR q-value=0.233) were enriched in low-risk group, implying normal neurodevelopment (Figure 4D).

To investigate the immune status in high-risk tumors, we carried out ESTIMATE analysis (Table S6, CGGA; Table S7, TCGA). In both datasets, we found that immune, stromal, and ESTIMATE scores were significantly elevated in the high-risk groups, while tumor purity clearly decreased (Figure 5A and B, p value < 0.001). Next, we selected immune checkpoints (PD1, PD-L1, TIM-3, B7-H3, B7-H4, LAG3)²⁰ and inflammatory genes (INF-α, INF-β, TNF-α, IL-6, IL-17, CCL2, CXCL2 and HLA-A)²¹ for Pearson correlation analysis with risk score (RS, Figure 5C and D). We found immune checkpoints were positively related to risk score (RS), indicating immune suppression in high-risk tumors. Among inflammatory genes, the expression of IL-6, CCL2, CCXCL2 and HLA-A were obviously correlated with risk score, implying enrichment of T cells and macrophage-associated immunity in the high-risk group.

Figure 5 Comparison of difference in immune status between high-risk and low-risk groups. With ESTIMATE tumor purity algorithm, we calculated immune score, stromal score and ESTIMATE score of two group patients in CGGA (A) and TCGA (B) dataset. Circos plots showed correlation analysis between risk score (RS) and immune checkpoints/inflammatory genes in CGGA (C) and TCGA (D) dataset. In circos plots, red/green scale bar represented Pearson’s correlation coefficient, which had a value between +1 and −1, where 1 is total positive linear correlation, 0 is no linear correlation, and −1 is total negative linear correlation.

Discussion

Numerous studies demonstrate that disorders of autophagy are closely related to the occurrence of various cancers. Autophagy deficiency induced by allele deletion of Beclin1 increases tumor cell apoptosis and delays the progression of breast cancer in the presence of P53.²² Common human cancer cell lines with RAS activating mutations generally exhibit higher autophagy activity, and the downregulation of important autophagy-related proteins significantly inhibits the growth of these cells.²³ In lung cancer, knockout of the autophagy gene, Atg7, led to tumor progression to oncocytomas, rather than adenomas and carcinomas, and tumor proliferation was markedly inhibited.²⁴ Numerous malignancies are associated with autophagy activation, with proteins related to autophagosome formation and maturation detected;^25,26 however, due to its clear heterogeneity, autophagy-related gene expression in glioblastoma remains controversial. Previous studies demonstrated that downregulation of LC3B-II and Beclin leads to decreased autophagy, inducing glioblastoma progression.²⁷ In contrast, more recent studies show that high levels of LC3/beclin expression are associated with poor prognosis in patients with GBM.²⁸ In some cell stress situations, such as during activation of malignant phenotypes, changes in the expression of a subset of genes may regulate autophagy. Therefore, study of the expression of autophagy-related gene sets has clinical implications.

In the present study, we screened 14 autophagy-related genes that could be used to cluster patients with glioblastoma into two groups, with significant differences in clinicopathology and prognosis. Next, we constructed a signature based on 14 genes, reflecting autophagy status in glioblastoma. We found that most high-risk patients were IDH wild-type, MGMT promoter non-methylated, and had more malignant tumor subtypes (including classical and mesenchymal), predicting poor prognosis. Cox regression analysis further confirmed the independent prognostic value of our signature.

Both functional analysis using DAVID and GSEA revealed an increase in immune and inflammatory responses in the high-risk group, indicating an interaction between autophagy and the immune microenvironment. A remarkable recent study suggested that autophagy regulates tumor immunity and thereby modulates malignancy.²⁹ In breast cancer with up-regulated autophagy, CD8+ cytotoxic T lymphocytes were enriched and immunosuppressive FOXP3+ regulatory T cells were reduced.³⁰ By targeting BECN1 to inhibit autophagy, Baginska et al restored granzyme B levels in hypoxic cells in vitro, and induced NK-mediated tumor cell killing in vivo.³¹ The expression of the immune checkpoint regulator, PD-L1, in lung cancer cells was down-regulated by activation of autophagy using an mTOR inhibitor.³² To further investigate the relationship between our autophagy-related signature and immunity in glioblastoma, we conducted ESTIMATE analysis, and found an increase in immune/stromal components, and a decrease in tumor purity, in high-risk patients. We also conducted further analysis of immune checkpoint and inflammatory genes and found that an immunosuppressive state and T cell- and macrophage-related immune factors were enriched in the high-risk group.

Much preclinical evidence suggests that inhibition of autophagy improves clinical outcomes in patients with cancer.^33,34 Meanwhile, it has been widely reported that inhibition of autophagy increases the sensitivity of cancer to anti-cancer drugs.^35,36 Further, concurrent regulation of autophagy can significantly improve the therapeutic effects of cancer immunotherapy³⁷ For example, alpha-TEA can stimulate autophagy to enhance antigen presentation to CD8 T cells as a strategy for improving tumor immunotherapy.³⁸ Spleen tyrosine kinase (Syk)-mediated autophagy stimulates cell surface major histocompatibility protein expression and CD4 T cell activation to enhance anti-tumor immunotherapy.³⁹ In the future, the relationship between autophagy and immunity needs to be further studied, and alterations in glioma cells autophagy may play an important role in immune regulation.

Conclusions

In summary, here we identify 14 autophagy-related genes that can be used to stratify patients with glioblastoma by clinical and pathological features. Further, we built a 14-autophagy-related gene expression-based signature, which could accurately predict important molecular markers and survival in patients with glioblastoma. Our study provided a new understanding of autophagy status in glioblastoma.