Identification and validation of an autophagy-related signature for predicting survival in lower-grade glioma

ABSTRACT

Abnormal levels of autophagy have been implicated in the pathogenesis of multiple diseases, including cancer. However, little is known about the role of autophagy-related genes (ARGs) in low-grade gliomas (LGG). Accordingly, the aims of this study were to assess the prognostic values of ARGs and to establish a genetic signature for LGG prognosis. Expression profile data from patients with and without primary LGG were obtained from The Cancer Genome Atlas (TCGA) and Genome Tissue Expression databases, respectively, and consensus clustering was used to identify clusters of patients with distinct prognoses. Nineteen differentially expressed ARGs were selected with threshold values of FDR < 0.05 and |log2 fold change (FC)| ≥ 2, and functional analysis revealed that these genes were associated with autophagy processes as expected. An autophagy-related signature was established using a Cox regression model of six ARGs that separated patients from TCGA training cohort into high- and low-risk groups. Univariate and multivariate Cox regression analysis indicated that the signature-based risk score was an independent prognostic factor. The signature was successfully validated using the TCGA testing, TCGA entire, and Chinese Glioma Genome Atlas cohorts. Stratified analyses demonstrated that the signature was associated with clinical features and prognosis, and gene set enrichment analysis revealed that autophagy- and cancer-related pathways were more enriched in high-risk patients than in low-risk patients. The prognostic value and expression of the six signature-related genes were also investigated. Thus, the present study constructed and validated an autophagy-related prognostic signature that could optimize individualized survival prediction in LGG patients.

KEYWORDS:

• Lower-grade glioma
• the Cancer Genome Atlas
• Chinese Glioma Genome Atlas
• autophagy-related gene
• prognosis signature

1. Introduction

Gliomas are the most frequent primary tumors of the central nervous system [1], and low-grade gliomas (LGGs), which originate from neuroepithelial tissue [2], are notably fatal in young adults and exhibit highly variable clinical manifestations that make them more difficult to diagnose [3]. Furthermore, despite improvements in the treatment of LGG, disease prognosis remains unsatisfactory owing to the conditions’ complex pathogenesis [4]. Owing to recent technological advances in next-generation sequencing, numerous studies have investigated and identified molecular biomarkers for LGG, including the O-6-methylguanine-DNA methyltransferase (MGMT) methylation and isocitrate dehydrogenase (IDH) mutation [5,6], thereby providing valuable insight into LGG pathogenesis and advancing the development of anti-cancer therapies [7]. However, the most widely used markers fail to adequately account for inter-individual variation or properly assess LGG risk [8,9]. Consequently, there is a pressing need to elucidate the underlying mechanisms of LGG and to develop effective biomarkers that facilitate the early identification and preventative treatment of high-risk patients.

Various relationships between cancer and autophagy have been described. Autophagy (i.e., type II cell death) is a highly conserved and important biological process that maintains intracellular homeostasis by removing aggregation-prone proteins and dysfunctional organelles [10] and plays a crucial role in a variety of diseases, including metabolic and cardiovascular disorders, as well as cancer [11]. However, increasing evidence shows that autophagy plays a bilateral role in cancer pathogenesis, depending on tumor type and stage [12,13]. Indeed, in some cases, autophagy protects cells from oxidative stress and harmful genomic mutations at the level of tumor initiation, thereby suppressing tumor formation [14,15], but in other cases, unbalanced autophagy can promote tumor growth and cancer cell fitness under certain stress conditions [16]. As such, recent studies have investigated the association between autophagy and gliomas and have reported that, in some cases, tumor growth can be inhibited by targeting autophagy [17]. In fact, it is even possible that temozolomide (TMZ), which is commonly used to treat pediatric and adult gliomas, extends the duration of patient survival by stimulating autophagy [18,19].

Previous studies have reported that the abnormal regulation of autophagy is directly related to tumorigenesis. We reasoned that autophagy plays an important role in its LGG progression and development. The expression and prognostic significance of genome-wide expression profiles of ARGs are largely unexplored, yet, in LGG. Accordingly, the aims of the present study were to identify differentially expressed ARGs in patients with and without primary LGG, to evaluate the function of the identified genes, and to develop and evaluate a prognostic ARG signature for LGG. In the present study, we applied bioinformatic methods on LGG and normal cases from the Cancer Genome Atlas (TCGA) database and the Genotype-Tissue Expression (GTEx) database to identify 19 differentially expressed ARGs. The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways analyses were applied for function analysis. Then, we constructed a 6-gene autophagy prognostic signature, which had good prognostic value and could be used as an independent prognostic indicator within the TCGA training group. In addition, the signature was further confirmed in the Chinese Glioma Genome Atlas (CGGA) database. In conclusion, the findings might contribute to further explore the functional role of ARGs in LGG pathogenesis and to improve clinical practice.

2. Materials and methods

2.1 Data collection

RNA sequencing data for a cohort of 508 patients with primary LGG and corresponding clinical data were downloaded from The Cancer Genome Atlas (TCGA) database [20], and expression data for patients without primary LGG (n = 1152) were obtained from the Genome Tissue Expression (GTEx) portal (https://www.gtexportal.org/) [21]. The expression data (fragments per kilobase of exon per million fragments, FPKM) from both databases were then log2-transformed for normalization and combined using the ‘limma’ package in R, and the batch effect was processed [22]. Expression data and clinical information for patients with primary LGG were also obtained from the Chinese Glioma Genome Atlas (CGGA) database (http://www.cgga.org.cn/) to be used as an external validation cohort [23].

2.2 Acquisition of autophagy-related genes

A total of 232 ARGs that were either directly or indirectly related to autophagy were obtained from the Human Autophagy Database (HADb) database (http://autophagy.lu/clustering/index.html) [24] and an additional 396 non-duplicate ARGs were obtained from Molecular Signatures Database v7.2 (MSigDB, https://www.gsea-msigdb.org/gsea/index.jsp.) [25], resulting in a final set of 539 non-duplicate ARGs for subsequent analysis.

2.3 ARG consensus clustering

The TCGA-LGG cohort data were clustered by ARG expression profiles using the ‘ConsensusClusterPlus’ package in R to investigate whether ARG expression is associated with prognosis [26], and principal component analysis (PCA) was used to better visualize gene expression patterns in diverse clusters. Overall survival (OS) analysis between different clusters was estimated using Kaplan-Meier method.

2.4 Identification and analysis of differentially expressed ARGs

Differentially expressed genes (DEGs) were identified using TCGA and GTEx data for patients with and without LGG, using the ‘limma’ package in R, and threshold values of FDR < 0.05, and |log2 fold change (FC)| ≥ 2. P-values were adjusted to control for the false discovery rate (FDR) using the Benjamini-Hochberg method [27]. The intersection of the DEGs and ARGs were regarded as the significantly differential expressed ARGs and was then visualized via Venn diagrams. A protein-protein interaction (PPI) network was generated by STRING (https://string-db.org/) [28] and Pearson correlation analysis was used to evaluate the relationship among different expressed ARGs.

To identify the major biological attributes of the differentially expressed ARGs, Gene Ontogeny (GO) and pathway (Kyoto Encyclopedia of Genes and Genomes, KEGG) analyses were performed, and the annotation analysis results were visualized using the ‘ggplot2’ and ‘GOplot’ packages in R [29]. Statistical significance was maintained at an adjusted P-value of < 0.05.

2.5 Construction of autophagy-related prognostic signature

Strong correlations between ARG expression and OS were identified by using univariate Cox regression (P < 0.05) to analyze DEG expression data from TCGA-LGG patients. Patients with primary LGG were randomly separated into a training cohort (60%) and internal validation cohort (40%) using the ‘caret’ package [30]. Subsequently, stepwise multivariate Cox regression analysis was performed to identify optimal prognostic gene sets and their coefficients in the TCGA training cohort, and the gene sets were used to construct an autophagy-related prognostic risk signature. The risk score formula was calculated as follows: Risk score = $\sum _{i=1}^{n}Coef\left(i\right)\times x\left(i\right)$, where Coef(i) denotes the coefficient of a specific gene and x(i) indicates the expression level of the same gene. The prognostic signature was used to calculate a risk score for each patient with LGG patient, and then the patients were separated into low- and high-risk groups according to the median risk score. Survival analysis and area under the receiver operating characteristic (ROC) curve (AUC) were used to assess prediction performance of the constructed signature, and both univariate and multivariate Cox regression analyses were used to investigate whether the signature can be used as an independent prognostic factor. The performance of the ARG-based prognostic signature was further evaluated using the TCGA testing cohort, entire TCGA cohort, and CGGA cohort.

2.6 Clinical applications of the ARG signature

Relationships between clinicopathological variables of LGG patients and the predicted high- and low-risk groups were compared. To further investigate the prognostic value of the two groups in different clinical subgroups, hierarchical survival analysis was performed for age (g40 and >40 years), sex (female and male), grade (G2 and G3), IDH status (mutation and wild-type), radiotherapy (yes and no), and chemotherapy (yes and no), as indicated by information in the TCGA database. Based on this analysis, two additional clinical characteristics from the CGGA database, (1p19q codeletion status and MGMT methylation status) were added.

2.7 Gene set enrichment analysis

Gene set enrichment analysis (GSEA) was conducted using GSEA software (http://www.broadinstitute.org/gsea/index.jsp) [31], with a minimum and maximum gene set size of 15 and 500, respectively, and 10,000 permutations. Statistical significance was set at P < 0.05, and FDR values of <0.25 were regarded as to be statistically significant.

2.8 Analysis of signature-related genes

OS analysis was used to determine the prognostic value of each signature gene to better investigate the function of these genes in disease progression. The Human Protein Atlas database (HPA; https://www.proteinatlas.org/) is an interactive open-access database that provides protein expression profiles for a variety of human proteins [32]. Protein expression levels in tissues from patients with and without LGG were evaluated using immunohistochemistry analysis.

2.9 Quantification of ARGs in patient samples

Quantitative reverse transcription-polymerase chain reaction (qRT-PCR) was used to quantify the expression of signature-related genes in 45 tissues samples, including 17 from patients with LGG, 20 from patients with high-grade glioma (HGG), and eight from peritumoral brain regions. All samples were obtained from Shandong Provincial Hospital, which is affiliated with Shandong University (Jinan, Shandong, China), between 2018 and 2019, with written informed consent. All pathological features were reviewed by two pathologists, and none of the patients had received chemotherapy or radiotherapy before surgery. All procedures were approved by the ethics committee of our hospital.

Briefly, qRT-PCR was performed using SYBR GreenER qPCR SuperMix reagents (Invitrogen, Carlsbad, CA, USA), a Bio-Rad iCycler (Bio-Rad Laboratories, Inc., Hercules, CA, USA), and GAPDH as an endogenous control. Gene expression was calculated using the 2− ΔΔCt method, and all primers were synthesized by Guangzhou Rui Bo Biological Technology (Guangzhou, China) (Table 1). All data represent means (n = 3).

Table 1. Primer sequences used for qRT-PCR in this study

Gene Name	Forward Primer	Reverse Primer
EEF2	5ʹ-ATACCCTGTATGAAGGGAAGCC-3’	5ʹ-CTTACCCCGAAGTTACGTCTTTC-3’
GABARAP	5ʹ-AGAAGAGCATCCGTTCGAGAA-3’	5ʹ-CCAGGTCTCCTATCCGAGCTT-3’
HSPA5	5ʹ-CACGGTCTTTGACGCCAAG-3’	5ʹ-CCAAATAAGCCTCAGCGGTTT-3’
PEA15	5ʹ-GGAGAGCCACAACAAGCTG-3’	5ʹ-CCATAGTGAGTAGGTCAGGACG-3’
TP53	5′-GAGGTTGGCTCTGACTGTACC-3′	5′-TCCGTCCCAGTAGATTACCAC-3′
CHMP4A	5ʹ-ATTTCTCGGCCTATGGGCTTT-3’	5ʹ-CGGCAGATGAGTAGAAGGTACA-3’
GAPDH	5′-CTGGGCTACACTGAGCACC-3′	5′-AAGTGGTCGTTGAGGGCAATG-3′

3. Results

Abnormal levels of autophagy have been implicated in the pathogenesis of a variety of diseases, including cancer. However, ARGs have rarely been analyzed in patients with LGG. We reasoned that autophagy plays an important role in its LGG progression and development. Accordingly, the aims of the present study were to identify differentially expressed ARGs in patients with and without primary LGG, to evaluate the function of the identified genes, and to develop and evaluate a prognostic ARG signature for LGG.

3.1 ARG consensus clustering

To characterize the association between autophagy and the outcome of patients with LGG, consensus clustering was performed using the expression profiles of 539 ARGs from the TCGA database. Patients could be divided into two robust clusters (Clusters 1 and 2) using k = 2 as the optimal value (Figure 1a–d), and principal component analysis of the expression profiles of Clusters 1 and 2 revealed distinct distributions (Figure 1e), with patients in Clusters 1 and 2 distributed in two different directors. In addition, the samples in Cluster 1 exhibited a significantly shorter OS than those in Cluster 2 (Figure 1f). Together, these results indicated that patient prognosis is related to autophagy and that autophagy may play a significant role in LGG occurrence and progression.

Figure 1. Clustering analysis of autophagy-related genes (ARGs) in patients with low-grade glioma (LGG). (a) Consensus clustering of LGG samples from The Cancer Genome Atlas (TCGA) dataset for k = 2 to 5. Optimal clustering is represented by k = 2. (b) Consensus clustering cumulative distribution function (CDF) for k = 2 to 9. (c) Relative change in area under CDF curve for k = 2 to 9. (d) Tracking plot for k = 2 to 9. (e) Principal component analysis of ARG mRNA expression profiles from patients with LGG. (f) Survival analysis in Clusters 1 and 2

3.2 Differentially expressed ARGs

Sequence data from the TCGA and GTEx databases were compared to identify DEGs in tissues from patients with and without LGG. A total of 1085 significantly aberrantly expressed genes (345 upregulated and 740 downregulated in LGG) were identified, and the outcomes were visualized using volcano plots and heatmaps (Supplementary Figure 1). A Venn diagram was then used to determine the intersection of DEGs and ARGs. As a result, 19 differentially expressed ARGs, including 10 upregulated genes (EEF2, EGFR, EIF4EBP1, HIF1A, HSPA5, MYC, PEA15, TP53, GFAP, and TREM2) and nine downregulated genes (FKBP1B, GABARAP, MAP1LC3A, SPNS1, FTH1, SNCA, CHMP4A, PIK3R2, and XBP1) were selected for further analysis (Figure 2a–b). The expression patterns of the differentially expressed ARGs were visualized intuitively using a box plot (Figure 2c), and the results of the STRING analysis indicated that TP53, EGFR, and HSPA5 were the hub genes (Figure 2d). The correlation between these genes is presented in Figure 2e.

Figure 2. Differentially expressed autophagy-related genes (ARGs) in patients with and without low-grade glioma (LGG). (a) Venn diagram. (b) Heatmap. (c) Expression patterns. (d) Protein-protein interaction network. (e) Correlation

3.3 Functions of differentially expressed ARGs

GO enrichment and KEGG pathway analyses were performed to further investigate the potential biological mechanisms and signaling pathways of the nineteen differentially expressed ARGs. As expected, the differentially expressed ARGs were enriched in several autophagy-related biological processes (e.g., autophagy, processes utilizing autophagic mechanisms, and responses to nutrient levels), cellular components (e.g., autophagosomes, late endosomes, and secondary lysosomes), and molecular function (ubiquitin protein ligase binding, ubiquitin-like protein ligase binding, and protein phosphatase binding). Bubble plots of the relationships between ARGs and GO enrichment terms are presented in Figure 3a and b. KEGG pathway analysis suggested that these differentially expressed ARGs were significantly associated with central carbon metabolism in cancer, Kaposi sarcoma-associated herpesvirus infection, and proteoglycans in cancer. Interestingly, the genes were confirmed to be involved in several cancers, including LGG and a number of cancer-related biological processes, including PI3K-Akt, ErbB, and HIF-1 signaling (Figure 3c and d). Taken together, these results indicate that the differentially expressed ARGs identified here play a critical role in the regulation of LGG malignancy.

Figure 3.. GO term and KEGG pathway enrichment of differentially expressed autophagy-related genes (ARGs) in patients with and without low-grade glioma (LGG). (a) Enriched GO terms. (b) Enriched GO terms, with the x- and y-axes representing z-score and negative log P-value, respectively. (c) Enriched KEGG terms. (d) Enriched KEGG terms, with the x- and y-axes representing z-score and negative log P-value, respectively

3.4 Construction of LGG prognostic signature

Since autophagy is clearly associated with LGG prognosis and the identified differentially expressed ARGs are likely associated with LGG malignancy, an autophagy-related prognostic signature was constructed in order to establish a more accurate process for evaluating the prognosis of patients with LGG. Univariate Cox regression analysis revealed that nine of the 19 differentially expressed ARGs were significantly associated with OS (P < 0.05; Figure 4a). Of those nine genes, HSPA5, TP53, GFAP, and TREM2 were recognized as risk genes (hazard ratio, HR > 1), whereas EEF2, GABARAP, MYC, PEA15, and CHMP4A were recognized as protective genes (HR < 1).

Figure 4. Construction of autophagy-related signature for patients with low-grade glioma (LGG). (a) Univariate analysis of differentially expressed autophagy-related genes (ARGs) to identify genes correlated with overall survival. (b) Kaplan-Meier curve for high- and low-risk patients. (c) Time-dependent ROC curve for 1-, 3-, and 5-year survival rates. (d) Distribution of risk scores in patients with LGG. (e) Scatterplots of patients with LGG and different survival statuses. (f) Expression of genes in high- and low-risk patients with LGG. (g) Univariate Cox regression analysis of clinicopathological parameters (including risk score) and patient survival. (h) Multivariate Cox regression analysis of clinicopathological parameters (including risk score) and patient survival

Patients with LGG in the TCGA dataset (n = 508) were randomly assigned to a training group (n = 307) or testing group (n = 201). The nine prognostic ARGs were subjected to the multivariate Cox hazards model, resulting in six target genes (EEF2, GABARAP, HSPA5, PEA15, TP53, and CHMP4A), and the six genes were used to construct a prognostic LGG signature for the TCGA training cohort. The risk score for each patient was calculated as follows: Risk score = (-expression value of EEF2 × 0.758) + (-expression value of GABARAP × 0.649) + (-expression value of HSPA5 × 0.391) + (expression value of PEA15 × 1.039) + (expression value of TP53 × 0.459) + (-expression value of CHMP4A × 1.871).

Using the median risk score, patients were stratified into high- and low-risk groups within the TCGA training cohort. The results of Kaplan-Meier analysis suggested that low-risk patients survived significantly longer than high-risk patients (Figure 4b). Indeed, the 1-, 3-, and 5-year OS rates were 98.5, 91.7, and 79.8%, respectively, for the low-risk group and only 89.1, 55.2, and 38.9% for the high-risk group. Moreover, the constructed signature demonstrated better prediction ability with AUC values of 0.870 (1-year), 0.835 (3-year), and 0.781 (5-year), respectively (Figure 4c). Patients were also ranked based on risk scores in order to analyze the survival distribution (Figure 4d), and the vital status of individuals with different risk scores was assessed by scatterplot (Figure 4e). The heatmap illustrates the expression of the six signature-related genes in LGG patients at different risk levels (Figure 4f).

Univariate and multivariate analyses were also conducted to evaluate whether the signature could be used as an independent risk indicator for LGG outcome. Both IDH status (HR = 0.159, P = 0.003) and signature-based risk score (HR = 1.314, P < 0.001) exhibited prognostic value (Figure 4g), and after multivariate analysis, the signature-based risk score (HR = 1.359, P = 0.001) remained a significant prognostic factor. Univariate and multivariate analyses identified the signature-based risk score as an independent predictive factor for patients with LGG (Figure 4h).

3.5 Internal validation of LGG prognostic signature

Two internal validation cohorts (TCGA testing cohort and entire TCGA cohort) were used to confirm the predictive power of the signature. The results of Kaplan-Meier analysis indicated that high-risk patients had poor survival outcomes, regardless of cohort (P < 0.05; Figure 5a and b). The 1-, 3-, and 5-year survival rates for the testing cohort were 0.856, 0.787, and 0.67, respectively (Figure 5c), whereas those for the entire TCGA cohort were 0.868, 0.810, and 0.735 (Figure 5d). The distribution of gene risk score and survival status and expression values of the six signature genes in the testing cohort and entire TCGA cohort are shown in Figure 5e and f).

Figure 5. Evaluation of autophagy-related prognostic low-grade glioma (LGG) signature. Kaplan–Meier survival curves of autophagy-related signature in the testing cohort (a) and the entire TCGA cohort (b). Time-dependent ROC curve analysis of the testing cohort (c) and the entire TCGA cohort (d). Risk score distribution, survival status, and risk gene expression in the testing cohort (e) and entire TCGA cohort (f)

3.6 External validation of LGG prognostic signature

To confirm that the constructed signature could be applied to different populations, the same formula was applied to 408 LGG records from the CGGA database. The results of Kaplan-Meier analysis indicated that the OS of high-risk patients was significantly lower than that of low-risk patients (Supplementary Figure 2a), and the 1-, 3-, and 5-year survival rates for the independent external validation cohort were 0.748, 0.715, and 0.6464, respectively (Supplementary Figure 2b). As in the TCGA cohort, patient mortality rate increased increasing risk score (Supplementary Figure 2c). Therefore, the risk signature accurately predicted the outcomes of patients with LGG.

3.7 Relationship between prognostic signature and clinical factors

The potential prognostic value of the constructed signature in subgroups stratified by age, sex, grade, IDH status, and radiotherapy and chemotherapy history were investigated. In all cases (except in the IDH subgroups), the results of Kaplan-Meier analysis indicated that the OS of high-risk patients was significantly lower than that of low-risk patients (Figure 6). However, in the CGGA cohort, this pattern was only observed in a few subgroups, including age > 40 years, females, grade 3, receiving radiotherapy, receiving chemotherapy, and 1p19q non-codel (Supplementary Figure 3).

Figure 6. Outcome prediction of the autophagy-related signature in stratified patients with low-grade glioma (LGG). Survival analysis of prognostic LGG signature in patients stratified by age (≤ 40 and > 40), gender (female and male), grade (G2 and G3), IDH (mutant- and wild-type), radiotherapy (no and yes), and chemotherapy (no and yes)

The association between signature and clinical variables was also investigated (Figure 7a). Older, female, grade 3, radiotherapy, and IDH wild-type patients tended to yield higher scores than younger, male, grade 2, non- radiotherapy, and IDH mutant-type patients, respectively, in the TCGA cohort (Figure 7b–g), and in the CGGA cohort, IDH1 wild-type, 1p19q non-codel, and MGMT unmethylated patients tended to yield higher scores than IDH1 mutant-type, 1p19q codel, and MGMT methylated patients (Supplementary Figure 4).

Figure 7. Relationship between clinical factors and the autophagy-related prognostic signature for low-grade glioma (LGG). (a) Expression of autophagy-related genes (ARGs) and distribution of clinicopathological features in high- and low-risk patients. (b-g) Effects of age, gender, grade, IDH, radiotherapy, and chemotherapy, respectively, on risk score distribution. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001, ns: no significant

3.8 Gene set enrichment analysis

To identify pathways that could play an important role in cancer progression in high-risk patients, GSEA analysis was performed on the basis of risk signature. High-risk patients were significantly involved in multiple autophagy-and cancer-related signaling pathways (Figure 8a). The significantly enriched autophagy-related pathways included Fc gamma R-mediated phagocytosis, the JAK-STAT signaling pathway, and lysosomes, whereas the significantly enriched cancer-related pathways included cancer, small cell lung cancer, P53 signaling pathway, MAPK signaling pathway, PPAR signaling pathway, cell cycle, and focal adhesion. Meanwhile, low-risk patients were associated with nonmalignant pathways (Figure 8b).

Figure 8. Significantly enriched KEGG pathways in The Cancer Genome Atlas (TCGA) cohort. Representative KEGG pathways in the high‐risk patients (a) and low-risk patients (b)

3.9 Evaluation of signature genes

Kaplan-Meier analysis was also used to investigate the prognostic value of each signature gene. Both upregulation of the two risk genes (HSPA5 and TP53) and downregulation of the four protective genes (EEF2, GABARAP, PEA15, and CHMP4A) were strongly correlated with inferior OS in the TCGA-LGG cohort (Figure 9), as well as in the CGGA cohort (Supplementary Figure 5), which indicated that the genes possessed significant prognostic value.

Figure 9. Overall survival analysis of signature-related autophagy-related genes (ARGs) in patients with and without in low-grade glioma (LGG). (a-f) EEF2, GAGBRAP, HSPA5, PEA15, TP53, and CHMP4A, respectively

The protein expression levels of the six prognostic ARGs in tissues from patients with and without LGG were measured using clinical specimens from the Human Protein Atlas (HPA) database (Figure 10). As anticipated, expression of the two high-risk genes (HSPA5 and TP53) was greater in tumor tissues than in normal tissues, whereas the expression of two of the protective genes (EEF2 and PEA15) was lower. These results are compatible with our data. However, neither CHMP4A nor GABARAP were detected in the HPA database for either LGG or control tissues and were not investigated further.

Figure 10. Immunohistochemistry of prognostic signature-related genes in patients with and without in low-grade glioma (LGG). (a-e) EEF2, HSPA5, PEA15, TP53, and CHMP4A, respectively

The mRNA expression of the six signature-related genes in 37 LGG and eight normal tissue samples was further validated by qRT-PCR. HSPA5 and TP53 were significantly upregulated in glioma tissues, whereas EEF2, GABARAP, PEA15, and CHMP4A were significantly downregulated (Figure 11), as described above. In addition, HSPA5, TP53, and GABARAP expression was strongly correlated with tumor progression, with the highest and lowest HSPA5 and TP53 expression observed in HGG tissues and normal tissues, respectively, and the opposite trend observed for GABARAP expression.

Figure 11. Differential expression of signature-related genes in patients with and without primary glioma. (a-f) EEF2, GAGBRAP, HSPA5, PEA15, TP53, and CHMP4A, respectively. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001

4. Discussion

Low-grade gliomas (LGGs), which constitute ~29% of all brain and central nervous system tumors [33]. Despite tremendous progress in multimodal treatment, the treatment outcomes of patients with LGG remain unsatisfactory [34]. Increasing evidence suggests that autophagy plays a significant role in the development and progression of a variety of tumors, including gliomas [35,36], and even though autophagy functions as a protective factor, which can suppress tumorigenesis during the early stages of disease, it can provide energy to cells and promote cancer initiation and progression during later stages [16]. Thus, elucidating the role of ARG expression in LGG is of great clinical significance.

The present study identified a group of 19 differentially expressed ARGs that could be used to assign patients with LGG to one of two subgroups, which had significantly different prognoses. Functional analysis suggested that the differentially expressed ARGs were closely related to cancer and revealed that the dysregulation of autophagy is important in cancer initiation and progression. Using the results of multivariate Cox analysis and risk-scoring methods, six autophagy genes (EEF2, GABARAP, HSPA5, PEA15, TP53, and CHMP4A) were retained as components of a prognostic signature that could be used to separate high- and low-risk patients in the TCGA training cohort. Given the possibility of overtraining during construction of the signature, the prognostic value of the signature was validated using two internal validation cohorts (TCGA testing cohort and entire TCGA cohort) and a single external validation cohort (CGGA cohort). The results suggested that the autophagy-related signature was robust and reproducible in patients with LGG, and stratified analyses also demonstrated that the signature was associated with both clinical features and LGG prognosis.

Autophagy is a highly conserved catabolic process that generally occurs in cancer cells [11] and plays a complex role in cancer development, depending on tumor stage, type, and genetic background [37]. In addition, increasing evidence has implicated the prognostic value of autophagy in a variety of cancers. In the present study, six signature-related genes were identified, with two (HSPA5 and TP53) considered risky genes and four (EEF2, GABARAP, PEA15, and CHMP4A) considered protective genes.

HSPA5/GRP78/BiP (heat shock protein 5), which is a master regulator of the unfolded protein response, is upregulated by glucose starvation, and plays a role in tumor cell growth by facilitating proper protein folding [38,39]. The protein is expressed in the ER but is also upregulated in many tumors [40] and, interestingly, has been reported to contribute to survival from head and neck cancer, possibly by preventing non-thermal plasma-treated solution (NTS)-induced lysosome inhibition [41]. However, the protein has also been associated with tumor progression, tumor size, and adverse prognoses in melanoma patients [42], as well as with tumor promotion in gastric cancer, hepatocellular carcinoma, and colon cancer [43–45].

TP53 (tumor protein p53) is a well-known tumor suppressor that plays an essential role in the maintenance of genome stability under a variety of cellular stresses. The protein represents one of the most common alterations found in human cancers, and mutant TP53 proteins can actively contribute to tumorigenesis [46]. Many studies have demonstrated that TP53 plays a role at various levels of the autophagy process (including regulation of the AMPK and Akt/mTOR pathways, ARGs, and hypoxia-inducible factor-1 target genes), thereby exerting oncogenic effects in several human malignancies [47]. Yang et al. reported that autophagy appeared necessary for optimal pancreatic tumor growth and progression in the physiological setting of TP53 loss during tumor progression [48].

EEF2 (Eukaryotic elongation factor 2) is a well-studied gene that plays a key role in protein synthesis by mediating the translocation step in peptide chain elongation [49]. In addition, the associated protein EEF2 is considered a tumor-associated antigen that is upregulated in most types of cancers and that plays a tumorigenic role in the growth of cancer cells [50]. However, some studies have reported that the drug nelfinavir, by inhibiting EEF2, can reduce global protein synthesis ability and affect the translation of specific cancer-promoting mRNAs [51].

GABARAP (gamma-aminobutyric acid receptor-associated protein) has been identified as a potential linker between GABAA receptors and microtubules [52] and is generally localized to the pericentriolar material of the centrosome, where it contributes to the formation of autophagosomes during starvation [53]. However, despite previous reports that the protein plays a role in receptor transport, the biological function of GABARAP remains unclear. For example, Klebig et al. reported that GABARAP mRNA and protein expression were significantly lower in breast tumors and in normal tissues and suggested that GABARAP transfection could actually be used to suppress breast cancer [54].

PEA15 (phosphoprotein enriched in astrocytes 15) is an acidic, serine-phosphorylated, 15-kDa phosphoprotein that plays an important role in intracellular signaling and the regulation of a variety of cellular processes, including apoptosis and proliferation, dependent upon stimulation [55]. Indeed, the protein has been reported to suppress the proliferation and invasion of certain cancerous cells, and the antitumor activity of PEA-15 in ovarian cancer is mediated, at least in part, by the induction of autophagy, which involves activation of the ERK1/2 pathway [56]. However, increasing evidence indicates that PEA-15 plays an oncogenic role in a variety of cancers [55].

CHMP4A (Charged multivesicular body protein 4A) belongs to the chromatin-modifying protein/charged multivesicular body protein family and is a component of endosomal sorting complex required for transport III, which is involved in the formation of endocytic multivesicular bodies [57]. However, few studies have investigated the role of CHMP4A in cancers, and the specific mechanism needs further study.

GSEA revealed that these high-risk groups were mainly enriched in KEGG pathways related to the tumorigenesis and progression of LGG, including cell cycle, JAK/STAT signaling, and MAPK signaling. Autophagy regulates cell cycle progression by removing key cell cycle regulators and other signal transducers [58]. Previous studies have reported that continued activation of the JAK/STAT signal can promote tumor progression [59]. The MAPK pathway regulates the activity and expression of key inflammatory mediators, including cytokines and proteases, which may function as potent cancer promoters [60]. The present study reveals that ARGs are involved in a variety of biological processes and signaling pathways, which suggests their importance in LGG initiation and development.

Although the present study provided a comprehensive analysis of autophagy-related gene expression profiles based on the TCGA and GETx databases, the study was also limited by selection and recall bias, due to the study’s retrospective design, by the need for further verification, and by the need for functional experiments to elucidate the mechanisms underlying the signature-related genes.

5. Conclusion

The present study established a reliable signature of six ARGs (EEF2, GABARAP, HSPA5, PEA15, TP53, and CHMP4A) for the prognosis of patients with LGG. A signature-based risk score was identified as an independent predictor of LGG prognosis, and GSEA revealed that the autophagy- and cancer-related pathways are enriched in high-risk patients. More importantly, the present study might shed new light on the potential association between autophagy and LGG, and the constructed signature can be used in both clinical practice and in the design of individualized therapeutic regimens.

Notes on Contributions

S.F. defined the research theme, discussed analyses, and drafted the manuscript. H.L., X.D. and P.D. performed the literature search and collected data for the manuscript. Q.P. and H.G. were responsible for the integrity of the entire study and manuscript review. All authors read and approved the final manuscript.

Highlights

(1) Autophagy-related genes were screened for association to low-grade glioma prognosis.

(2) A genetic signature was established to predict prognosis from low-grade glioma.

(3) The genetic signature was an independent prognostic factor for lower-grade glioma.

Ethics approval and consent to participate

This study was approved by the Ethics Committee of Shandong Provincial Hospital Affiliated to Shandong University (Jinan, Shandong, China). All participants were informed consent. This study was conducted in accordance with the Declaration of Helsinki.

Acknowledgements

We express our gratitude to all of the medical staff who contributed to uploading and organizing the TCGA, GTEx, and CGGA databases.

Availability of data and materials

Data associated with this manuscript can be made available from the corresponding authors on reasonable request.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Supplementary material

Supplemental data for this article can be accessed here.

Additional information

Funding

This work was supported by Jinan Science and Technology Plan (Grant No. 201907049), China National Natural Science Foundation (Grant No. 81771270), Shandong Provincial Natural Science Foundation (Grant No. 202019088), and Shandong Province Science and Technology Development Program (Grant No. 2015GSF118164).