Autophagy-Dependent Ferroptosis-Related Signature is Closely Associated with the Prognosis and Tumor Immune Escape of Patients with Glioma.

BACKGROUND: Ferroptosis is an autophagy-dependent form of cell death, sometimes called "ferritinophagy". Its related pathway has been proven to regulate the programmed death of glioma stem cells. Mining autophagy-dependent ferroptosis-related gene (AD-FRG) signature could facilitate the discovery of mechanisms and therapeutic targets showing drug resistance to chemotherapeutic drugs.
METHODS: We exhaustively searched HADB, MSigDB and FerrDb datasets and obtained 25 genes confirmed to exist in autophagy and ferroptosis death pathways. Glioma gene expression and clinicopathological data were collected from TCGA and CGGA datasets.
RESULTS: Lasso regression and Cox regression analysis were carried out to construct a nine AD-FRGs signature (SIRT1, MTDH, HSPB1, CISD2, HMOX1, ATG7, MTOR, PRKAA2 and EIF2AK4). ROC curve showed that nine genes signature could effectively predict 1- (AUC = 0.869), 3- (AUC = 0.922) and 5-year (AUC = 0.870) survival rates. Immunohistochemical images confirmed the protein expression level of the gene model. The prognostic nomogram of risk score, age, WHO grade, isocitrate dehydrogenase (IDH) wild-type condition, 1p/19q co-deletion state was built. The calibration curve demonstrated that the prediction of the nomogram is highly consistent with the actual results. Moreover, tumor microenvironment analysis showed that the high-risk group was associated with high immune infiltration status and high tumor purity. Correlation analysis showed that the expression of SIRT1, CISD2 and HSPB1 might be related to macrophage infiltration and immunotolerance in glioma tissues.
CONCLUSION: Based on autophagy-dependent ferroptosis-related genes, we established gene signature and nomogram that maybe effectively predict the overall survival rate of glioma and correlate with the immunosuppressive tumor microenvironment (TME).

Introduction

Autophagy is a conserved intracellular catabolism process that delivers cytoplasmic components to lysosomes of degradation.1 The role of autophagy is controversial, contextual in cancers.2 Autophagy is a conserved protein degradation process that fuels cell metabolism and energy pathways to maintain cell homeostasis. However, in oncology, autophagy has been described as a mechanism of cells developing resistance to chemotherapy, targeted therapy, or immunotherapy.3,4 Ferroptosis is not only a chemical element that regulates enzyme activity but also oxidizes a variety of substrates and causes biological damage.5 Ferroptosis-mediated oxidative stress is the central mechanism of Ferroptosis death.6,7 Ferroptosis, the newly detected form of cells death, depends on ferroptosis accumulation within cells and lipid peroxidation.8 In addition to inducing tissue damage and protective effects in neurodegenerative diseases, ferroptosis activation also indicates a tremendous anticancer activity.9 Studies have discovered that autophagia plays a pivotal role in ferroptosis induction.10 Although increased autophagy promotes ferroptosis accumulation and subsequent lipid peroxidation through selective degradation of anti-ferroptosis death regulators such as NCOA4 [nuclear receptor coactivator 4], GPX4 [glutathione peroxidase 4] and ARNTL/BMAL1 [aryl hydrocarbon receptor nuclear transport protein, etc],11,12 the underlying molecular mechanisms remain largely unknown.

There is a complicated feedback system between ferroptosis and autophagy. Increasing studies focus on the crosstalk between ferroptosis and autophagy, and show that ferroptosis depends on autophagy13,14. Evidence supported that excessive activation of autophagy cells can promote ferroptosis cell death through ferroptosis activating factor15. However, overexpression of autophagy inducers also leads to resistance to cancer treatment.16 Autophagy and ferroptosis play a critical role in cell fate, and a comprehensive analysis of shared genes between them could help understand feedback mechanism.16,17 Thus, we sought to identify critical genes that cause crosstalk through mining ferroptosis genes in autophagy pathways.

It has been shown that autophagy-dependent ferroptosis was a new form of apoptosis that is dependent on iron-mediated excess lipid peroxidation.15,18 Studies have reported that SIRT3/AMPK/mTOR/GPx4 and COPZ1/NCOA4/FTH1-related pathways modulate autophagy-dependent ferroptosis cancer cell death.19,20 The studies confirmed that the increase of SIRT3 expression could stimulate AMPK mTOR to release iron from ferritin, which leads to free iron accumulation and lipid peroxidation, thereby promoting autophagy activation. The decrease of SIRT3 can inhibit ferroptosis through regulating the GPx4 level.19 NCOA4 is an autophagy receptor. Under iron deficiency, ferritin is degraded in lysosomes and releases iron to maintain intracellular iron level. The increase of free iron level can mediate the autophagic degradation of ferritin, that is, the iron storage complex of cytosol.11,21 When cells respond to oxidative stress, autophagy is activated, and level of NCOA4 or ATG7 feedback increases, and activates ferritinophagy.22,23 Traumatic iron death is detected in several human diseases, such as neurodegenerative diseases, diabetes mellitus, infectious diseases, and cancer.24,25 Imbalance between autophagy and iron death signal allows cancer cells to escape the programmed cell death program resulting in chemoradiotherapy resistance and seriously affecting patients’ prognosis.26 The search of AD-FRGs associated with prognosis facilitates the discovery of novel therapeutic targets to overcome drug resistance in cancer.

As a fatal brain tumor showing rapid proliferation and antitumor therapy characteristics,27 glioma is often classified into lower-grade gliomas (LGGs; WHO grade II–III) and glioblastoma (GBM; WHO grade IV).28,29 Among them, glioma that belongs to WHO grade IV, which is extremely aggressive, highly lethal, extremely malignant, with a high recurrence rate and a very poor prognosis, is individually named glioblastoma multiforme (GBM), which accounts for up to 56% of all gliomas and approximately 80% of all primary intracranial malignancies.30,31 Glioma cells not only show an uninhibited tumor growth pattern but also have an aggregation effect on cells at various levels of differentiation, thus forming a “tumor microenvironment”, and thus presenting the characteristics of tumor heterogeneity, which increases the resistance of glioma cells to chemotherapy and radiotherapy.32,33 Therefore, there is a need to better understand the mechanisms of glioma occurrence and development at the genetic and molecular levels, which in turn will provide new theoretical basis and methods for the diagnosis and treatment of glioma.

More and more studies have shown that prognostic markers help improve the survival rate of patients with glioma, and predict tumor progression earlier, so as to intervene in advance. For example, Tan et al34 identified six immune-related gene models that evaluate the prognosis of glioma, and Zhang et al35 reported that seven gene features could be used as prognostic factors of glioma patients. Although some prognostic gene signatures have been reported, none of them have been used in clinic. This indicates that more genetic signature studies are needed as candidate signatures. In this work, we identified a nine AD-FRGs signature strongly associated with prognostic risk of glioma. Previous studies have confirmed the existence of these genes in autophagy-dependent ferroptosis-related pathways. The prognostic nomogram model indicated that the biomarker had the excellent prognostic ability.

Materials and Methods

Obtaining Glioma Datasets

The TCGA RNA-seq database containing 697 glioma cases from the UCSC Xena browser was taken as a training cohort. We collected gene expression matrix from normal brain tissues from the Genome Tissue Expression (GTEx) dataset (n = 475). Normalized gene expression was measured by FPKM and further performed with log2-based transformation. Similarly, RNA sequence data and related clinicopathological information (such as gender, age, and grade, subtype, IDH condition, 1p/19q state, and survival information) were downloaded from the CGGA database (). CGGA cohort (mRNAseq_693, mRNAseq_325) was used as a validation cohort. After that, batch effect of data merged from different batches were eliminated using SVA package. Only cases with both gene expression profile and clinicopathological data were enrolled for subsequent study, otherwise removed.

Human Autophagy-Dependent Ferroptosis-Related Genes

Autophagy-related genes were obtained from the Human Autophagy Database (HADb, ) and the Molecular Signatures Database (MSigDB, ), respectively. After eliminating overlapping genes, we included a total of 531 autophagy-related genes used for analysis. In addition, we also collected 259 ferroptosis-related genes from the FerrDb dataset (), and distinguished differentially expressed genes (DEGs) between tumors (TCGA dataset) and adjacent nontumorous tissues (GTEx dataset). The false discovery rate (FDR) <0.05 was considered a significant difference. Then, the Venn diagram was used to display the DEGs, which also involved both autophagy and ferroptosis pathways.

Construction and Validation of an AD-FRG Signature

We used univariate Cox analysis to screen AD-FRGs with prognostic significance, and a forest plot to visualize the results. The LASSO-penalized Cox regression analysis was employed to establish a predictive model to reduce the overfitting risk via the “glmnet” R package. The risk score of each patient was calculated according to the risk score formula, and training cohort patients were further divided into low- or high-risk groups based on median risk scores. The survival difference between the two groups was shown by the Kaplan–Meier curve and compared by Log rank test using the “survival” and “survminer” R package. Time-dependent ROC curve and AUC were applied to evaluate the performance of the gene signature and clinical factors in predicting 1-, 3- and 5-year survival chances.

Verification of Gene Signature Expression and Mutation

Immunohistochemical images were retrieved from the Human Protein Atlas (HPA) database () for comparison of protein expression levels associated with gene signatures. Mutation levels of gene signatures were obtained from the cBioPortal website.

Development and Assessment of a Predictive Nomogram

Univariate and multivariate Cox regression analyses were first carried out for the construction of a nomogram. We used a forest plot to display P-value, HR, and 95% CI of each variable via the “forestplot” R package. The independent prognostic factors filtered by the above steps were used to construct a nomogram using the “rms” and “survival” R-package. Next, we evaluated whether the predicted survival result was close to the actual output in the calibration curve. Nomogram‐predicted survival and the autual survival rate were indicated on the X- and Y-axes, respectively, and the 45° lines showed the optimal predictions.

Gene Set Enrichment Analysis and Tumor Immune Microenvironment Analysis

We used GSEA software for enrichment analysis of biological function pathways to analyze gene function of signature. The enrichment results of the high-risk group were above the X-axis, and the results of the low-risk group were below the X-axis.36 Pathway with a nominal p-value less than 0.05 and FDR<0.1 was selected.

To study the role of AD-FRG signature in immune system response, we used ssGSEA, CIBERSORT method and ESTIMATE algorithm to compare the proportion of immune components and tumor purity in high- and low-risk glioma samples. SsGSEA was performed using the “GSEAbase” R package to quantify the enrichment activities of 16 immune cells and 13 immune functions.37 The proportion of 22 human hematopoietic cell phenotypes in the immune microenvironment was further analyzed using CIBERSORT algorithm.38 In addition, to assess the degree of immune infiltration and tumor purity, we performed the ESTIMATE algorithm to calculate the matrix, immunity, and estimated scores.39

Correlation Analysis Between AD-FRG Signature and Immune Cells

To further investigate the relationship between AD-FRG signature and immune cells associated with glioma, we calculated the correlation coefficient between gene expression and relative proportion of immune cells. The data were visualized using “corrplot” package in R studio. R > 0.2 was defined as relevant. Finally, co-expression analysis of gene expression and immune cell content was performed. A p < 0.001 indicated a significant correlation using Pearson correlation coefficient.

Statistical Analysis

All statistical analyses and charting were performed in R v. 3.6.2. Lasso-Cox regression analysis was used to construct a prognosis prediction model based on AD-FRGs. Nomograms and calibration curves were drawn in the “rms” R package. Pearson correlation coefficient and Spearman correlation coefficient were calculated to test the correlation between variables. P < 0.05 was considered as statistically significant.

Results

AD-FRGs Identification

The research was conducted following the flowchart is carried out (Figure 1). We collected and integrated 55,175 RNA expression profiles of 697 glioma patients from the TCGA dataset and identified a total of 199 autophagy-related and 303 ferroptosis death-related genes as differentially expressed gene in tumor tissues compared with normal tissues. Twenty-five AD-FRGs, the ferroptosis death function genes in autophagy pathways, were identified from the intersection of the identified DEGs (Figure 2A).

Figure 1

The study flow diagram. In this study, the AD-FRG signature was screened based on TCGA-GBM/LGG dataset. Then, a prognostic nomogram was established in the training and validation cohorts to predict the survival time, respectively. Finally, CIBERSORT and ESTIMATE algorithms are instigated further to explore the underlying mechanism of immune-related prognostic signature.

Figure 2

Identification of the candidate genes. (A) Venn diagram to identify autophagy-dependent ferroptosis-related DEGs between TCGA dataset and GTEx dataset. (B and C) GO and KEGG enrichment analysis of the 23 prognosis-related DEGs. (D) Forest plots showing nine-gene signature expression. (E) Kaplan–Meier survival analyses in training cohort. (F) ROC analysis in training cohort. (G) Kaplan–Meier survival analyses in verification cohort. (H) ROC analysis in verification cohort. *P<0.05, ***P<0.0001.

Functional enrichment analysis revealed that the most relevant signaling pathways were “autophagy-animal”, “longevity regulatory pathway”, and “FOXO signaling pathway” (Figure 2B). Meanwhile, the most enriched terms of Gene Ontology were “regulation of autophagy”, “positive regulation of autophagy”, “macroautophagy”, “regulation of macroautophagy,” and “cellular response to chemical stress” (Figure 2C).

Establishment and Validation of AD-FRGs Prognostic Model

Overall, 666-glioma cases from the TCGA dataset were included after the survival analysis as a training set, and 23 prognostically relevant AD-FRGs were identified (Table 1). We used LASSO regression to remove CO-expressed genes to avoid overfitting signature ( and ). The PPI network downloaded from the STRING database indicated the interactions among the nine AD-FRGs (). Finally, multivariate Cox regression analysis was employed to identify nine AD-FRGs (SIRT1, MTDH, HSPB1, CISD2, HMOX1, ATG7, MTOR, PRKAA2, EIF2AK4) to construct gene signature. Specifically, six AD-FRGs were of HR > 1 and three AD-FRGs were of HR < 1 (Figure 2D).

Table 1

Twenty-three Prognosis-Related Genes Based on Autophagy-Dependent Ferroptosis-Related Genes

ID	HR	HR.95L	HR.95H	p value
BECN1	0.0735	0.0242	0.2229	0
BNIP3	0.1607	0.0704	0.3668	0
EIF2AK4	16.5955	6.5337	42.1526	0
HMOX1	10.4122	6.8787	15.7609	0
IFNG	4.7681	2.8673	7.9293	0
MAPK3	0.0085	0.0026	0.0281	0
MTDH	4.4488	1.3093	15.1168	0.0168
PRKAA1	54.0037	18.7485	155.554	0
PRKAA2	0.1594	0.0966	0.2631	0
SESN2	5.176	2.604	10.2882	0
SIRT1	0.0095	0.0053	0.0172	0
ATG7	4.7031	1.8873	11.7203	0.0009
MT3	4.6779	2.634	8.3077	0
MTOR	0.4541	0.2221	0.9286	0.0305
PIK3CA	0.127	0.0731	0.2206	0
ATG13	0.0049	0.0021	0.0114	0
ATP6V1G2	0.1224	0.0938	0.1597	0
CISD2	488.8903	171.7366	1391.746	0
DDIT3	7.8979	4.7319	13.1821	0
FBXW7	0.4947	0.313	0.7819	0.0026
HSPB1	78.1609	41.6145	146.8026	0
KEAP1	0.005	0.001	0.0254	0
MAPK8	0.0205	0.0127	0.0333	0

The risk score formula for each patient is as follows: Risk Score= (0.993×EIF2AK4 Expression) + (0.525×HMOX1 Expression) + (2.846×MTDH Expression) + (−0.515×PRKAA2 Expression) + (−2.905×SIRT1 Expression) + (−1.682×ATG7 Expression) + (1.175×MTOR Expression) + (3.361×CISD2 Expression) + (2.315×HSPB1 Expression).

Cases were assigned to high-risk (n = 333) and low-risk (n = 333) groups according to the median value of the risk score (0.681) ( and ). Kaplan–Meier survival analysis showed that the OS of the high-risk group was significantly lower than that of the low-risk group (Figure 2E; P < 0.001). The AUC values at 1, 3, and 5 years were 0.870, 0.922, and 0.869, respectively, highlighting that the AD-FRG signature had a strong prediction ability (Figure 2F).

Validation of AD-FRG Signature in External Datasets

To further validate the reliability of the AD-FRGs prognostic signature model, we divided 657 patients from the CGGA cohort into high-risk (n = 328) and low-risk (n = 329) groups using the corresponding median risk as the cutoff ( and ). The survival time of the high-risk population was significantly lower than that of the low-risk population (Figure 2G), showing that AD-FRG signatures had high prognostic prediction ability (Figure 2H).

Expression and Mutation of Nine Model Genes

In TCGA GBM/LGG cohort, compared with adjacent non-tumor brain tissues, MTDH, HSPB1, CISD2, HMOX1 and mTOR were high-expressed in glioma, while the expression of SIRT1, ATG7, PRKAA2 and EIF2AK4 were low-expressed. To confirm these nine gene expressions in clinical practice, we further validated the expression of proteins encoded by nine genes signature using clinical samples from HPA. Figure 3 shows immunohistochemical images of nine genes signature in glioma and the normal cerebral cortex. The mutation status of nine gene signatures was identified in the cBioPortal database, indicating that MTOR, EIF2AK4, HMOX1, and MTDH have significant mutations. The mutations of ATG7, SIRT1, EIF2AK4, HMOX1, MTDH and MTOR are shown in Figure 4.

Figure 3

Protein expression levels of gene signatures in human tissues, SIRT1 (A), ATG7 (B), PRKAA2 (C), EIF2AK4 (D), HMOX1 (E), MTDH (F), MTOR (G), CISD2 (H), and HSPB1 (I).

Figure 4

Mutation characteristics of model genes (A), ATG7 (B), SIRT1 (C), PRKAA2 (D), EIF2AK4 (E), HMOX1 (F), MTDH (G), MTOR (H), and HSPB1 (I). *There is currently information about targeted drugs associated with this mutation.

Independent OS Prediction by the AD-FRG Signature

To further confirm whether the risk score of generated AD-FRG signature was an independent predictor of patients with glioma, we performed univariate and multivariate Cox regression analysis to evaluate the ability of risk score and other clinical parameters to predict OS. Univariate and multivariate Cox analysis has shown that age, WHO grade, IDH condition and 1p/19q state and risk score were independent risk factors for OS (Figure 5; P < 0.05).

Figure 5

Independent predictive power of the risk score and individual clinicopathological covariates. Univariate (A and C) and multivariate (B and D) analyses grounding the nine-gene signature and clinical covariates in TCGA training (A and B) and CGGA validation (C and D) cohorts.

To make the results easier to interpret, subgroup analysis was carried out according to WHO grade, age, gender, IDH condition, 1p/19q co-deletion state, and treatment or non-treatment. For glioma patients, the risk scores of the age > 50 subgroup, WHO stage IV subgroup, IDH wild-type subgroup and 1p19q non-codel subgroup were much higher than the relative group (; P < 0.001). Finally, we analyzed the relativity between individual genes and clinicopathological features. The outcome indicated that EIF2AK4, HMOX1, SIRT1, CISD2 and HSPB1 could effectively distinguish clinical grades ().

Constructing Prediction Nomogram

To develop a clinical strategy for predicting OS in glioma patients, we integrated the AD-FRG signature and multiple clinical factors to build a nomogram to evaluate the OS at 1, 3, 5 and 10 years. A nomogram was constructed based on the independent predictors of the CGGA cohort (age, WHO grade, IDH1 wild-type and 1p/19q codeletion) (Figure 6A). The calibration curve verifies the accuracy of the nomogram (Figure 6B). The AUC value was 0.830, which was higher than that of clinical predictors alone (Figure 6C).

Figure 6

Construction of a nomogram. Construction of the nomogram based on prognostic markers in the validation (A) cohorts. Calibration plot assessed the predicted power of OS in the validation (B) cohorts. ROC analysis verified the prognostic effects of risk score and clinical covariates in the validation (C) cohorts.

The AD-FRGs Signature Was Related to Immune Cell Infiltration

Through GSEA analysis, the molecular biological functions of prognostic models were mainly found to be related to some carcinogenic pathways, such as cell cycle, ECM–receptor interaction and ERBB signaling pathway (Figure 7A and B). In addition, AD-FRGs signature was associated with multiple immune disease-related pathways, for instance, autoimmune thyroid disease, intestinal immune network for IgA production, Leishmania infection and systemic lupus erythematosus (Figure 7A and B).

Figure 7

Gene Set Enrichment Analysis (GSEA) and tumor immune microenvironment analysis. (A and B) Enriched gene sets annotated by the C7 collection between the high- and low-risk groups in the TCGA cohort. The ssGSEA was used to estimate the 16 immune cell scores (C) and 13 immune-related functions (D) in the TCGA cohort. The 22 immune-related terms were used to assess immune cells in the tumor microenvironment (E). Stromal (F), immune (G) and ESTIMATE (H) scores between high- and low-risk groups was estimated by the ESTIMATE algorithm. *P<0.05, **P<0.01, ***P<0.0001.

To further examine the relation between AD-FRGs signature and immune microenvironment, we quantified infiltration scores of immune cells and immune-related functions with ssGSEA. Notably, the low-risk group was related to a higher immune score of the tumor immune microenvironment (Figure 7C and D). In addition, the above differences were also confirmed in the CGGA cohort ().

Subsequently, 22 immune cells in glioma tissue between high-risk and low-risk groups were studied by the CIBERSORT algorithm (the data were visualized by violin map). We found high-risk groups showed significantly higher levels of T cells CD8, T cells CD4 memory resetting, T cells CD4 memory activated, T cells gamma delta, macrophages M0, macrophage M1, mass cells resting, mass cells activated, eosinophils and neutrophils, and lower levels of macrophage M2 (Figure 7E; P < 0.001). Kaplan–Meier survival analysis of immune cells showed that macrophage M0, macrophage M1, neutrophils, eosinophils, mast cells activated, monocytes, T cells CD4 memory activated, T cells CD4 memory resting, and T cells gamma delta had good ability to predict OS (Figure 8; P < 0.001).

Figure 8

Immune microenvironment and prognosis. K-M survival analysis indicated that macrophases M0 (A), macrophage M1 (B), neutrophils (C), eosinophils (D), mast cells activated (E), monocytes (F), T cells CD4 memory activated (G), T cells CD4 memory resting (H), and T cells gamma delta (I) were significantly correlated with survival duration.

After processing by the ESTIMATE algorithm, higher estimated scores were discovered in the high-risk group. Similarly, the proportion of immune and stromal cells was related to the high-risk group (Figure 7F-H).

In addition, the correlation analysis between tumor clinical grade and immune cell infiltration showed that a higher proportion of macrophage M0, macrophage M1 and T cells CD4 memory resting was related to higher tumor stage, while a lower proportion of mast cells activated was associated with higher tumor stage ().

Co-Expression Analysis of AD-FRGs Signature and Immune Cell Signature

Finally, to further investigate the correlation between AD-FRG and immunocyte signatures (Figure 9A), we calculated the correlation coefficients between them. SIRT1 was negatively correlated with the immune scores of T cells CD4 memory resting (R = −0.33, P < 0.001), macrophage M0 (R = −0.35, P < 0.001), macrophage M1 (R = −0.42, P < 0.001) and neutrophils (R = −0.32, P < 0.001), but positively correlated with the mast cells activated. CISD2 was positively correlated with the immune scores of macrophage M0 (R = −0.43, P < 0.001), macrophage M1 (R = −0.3, P < 0.001) and neutrophils (R = 0.38, P < 0.001), but negatively correlated with the risk score. HMOX1 was significantly correlated with macrophage M0 (R = 0.34, P < 0.001) and neutrophils (R = 0.35, P < 0.001). HSPB1 was significantly correlated with T cells CD4 memory resting (R = 0.3, P < 0.001) and macrophage M0 (R = 0.32, P < 0.001) (Figure 9B-M).

Figure 9

Correlation between AD-FRG signature and immune cells. (A) Co-expression model of immunocyte components and AD-FRG signature. The proportion of T cells CD4 memory resting (B), macrophage M0 (C), macrophage M1 (D) and neutrophils (E) were apparently negatively correlated with the expression level of SIRT1, mast cells activated was just the opposite (F). The proportion of macrophage M0 (G), macrophage M1 (H) and Neutrophils (I) were positively related to the expression level of CISD2. The proportion of macrophage M0 (J) and neutrophils (K) were significantly positively correlated with HMOX1 expression. The proportion of T cells CD4 memory resting (L) and macrophage M0 (M) was positively correlated with the HSPB1 expression.

AD-FRGs Were Associated with Chemotherapy and Immunotherapy Responses

The relationship between drug sensitivity of glioma cell lines and relative expression levels of AD-FRGs was explored using the GDSC database. We further analyzed the correlation between the expression of AD-FRGs and the IC50 of multiple chemotherapeutic agents (Figure 10A). Glioma with high expression of AD-FRGs had higher drug resistance to temozolomide (TMZ), methotrexate, vincristine, and vinorelbine. The potential sensitivity of the two groups to immune checkpoint (CD274, CTLA4, LAG3, and PDCD1) inhibitors was also analyzed. Pods plots showed the relationship between the risk score and four immune checkpoint targets (Figure 10B). The expression of immune checkpoints in high-risk patients increased significantly.

Figure 10

Prediction of response to chemotherapy and immunotherapy. (A) IC50 values of six cytotoxic chemotherapeutic drugs (temozolomide, methotrexate, vincristine, vinorelbine, bleomycin, and etoposide) selected by the pRRophetic algorithm in TCGA cohort, Spearman correlation coefficient was calculated to test the correlation between variables. (B) The expression of immune checkpoints (CD274, CTLA4, LAG3, and PDCD1) in both groups in TCGA cohort. ***P<0.0001.

Discussion

Due to the aggressiveness and recurrence of brain glioma, the current treatment effect is far from satisfactory.40 Genes of the autophagy-dependent ferroptosis-related pathway have been increasingly discovered in glioma tissues, but none has integrated these genes to establish a gene prognosis model. Studies have described ferritinophagy as a form of autophagy-dependent programmed cell death.41 The link between autophagy and ferroptosis may imply a complex interplay between metabolic dysfunction and oxidative stress.13,42 Therefore, it may be helpful to build a ferritinophagy-related RNA signature to predict the glioma patients’ survival and improve therapeutic approaches.

We were the first to screen ferroptosis genes in autophagy pathway. The current analysis and validation demonstrate that the prognosis model constructed by nine autophagy-dependent ferritinophagy-related genes (EIF2AK4, HMOX1, MTDH, PRKAA2, SIRT1, ATG7, MTOR, CISD2, and HSPB1) can effectively predict the survival time. Previous studies have confirmed the critical role of these genes in the autophagy-dependent iron death pathway. Our study confirmed that the knockdown of these genes is more conducive to prognosis of glioma, except EIF2AK4.

Our results showed that the down-regulation of the three genes (ATG7, SIRT1, PRKAA2) played a protective role in glioma. Autophagy-related protein 7(Atg7) can trigger regulatory cell death (RCD) in glioma and often serve as a target for inducing autophagic-dependent death.43 Studies have shown that SIRT1 inhibition can inhibit tumor growth.44 PRKAA2 is also known as AMP-activated protein kinase (AMPK) again. Inhibition of AMPK has been found to reduce the activity of GBM tumor cells.45

Next, of the six risk genes, only general control nonderepressible 2 (GCN2; EIF2AK4) was downregulated. EIF2AK4 is an environmental sensor, activated EIF2AK4 up-regulates ATF4 to induce the expression of Sestrin2, which is vital for tumor cell survival under oxidative stress.46 However, there is no study reporting the effect of activating EIF2AK4 on the treatment of gliomas.

In addition, five risk genes (HMOX1, MTDH, MTOR, CISD2 and HSPB1) were up-regulated compared with in other tumors. In many reports, mechanistic target of rapamycin (mTOR) signaling pathway is closely related to tumor growth, metabolism, apoptosis and autophagy. MTOR inhibitors have been widely used in tumor therapy in recent years. It has been demonstrated that the deficiency of SIRT3 is resistant to autophagy-dependent ferritinophagy through restraining the AMPK/mTOR pathway and elevating GPX4 expression levels.19 Inhibition of PI3K/Akt/mTOR (mechanistic target of rapamycin) pathway can be used as a therapeutic target for GBM.47 MTOR-related pathways are complex, and their roles in autophagy and ferroptosis in glioma are widely studied. Upregulation of CDGSH iron sulfur domain 2 (CISD2) can inhibit beclin-1-mediated autophagy and promote glioma proliferation.48 Previous studies indicated that CISD2 overexpression maintains the efficacy of sulfasalazine by regulating iron death.49 However, no studies have confirmed whether CISD2 can overcome drug resistance in glioma. Elevated expression levels of heat shock protein beta-1 (HSPB1) could activate SIRT2-mediated G6PD and promote glioma cell proliferation,50 butHSPB1 can inhibit TfR1 recovery and reduce intracellular free iron concentration.51,52 When HSPB1 is inhibited, iron ions accumulate and activate elastin-induced ferroptosis,53 providing a research direction for the development of new anticancer drugs. Metadherin (MTDH) is highly expressed in gliomas, which can induce epithelial-mesenchymal transition (EMT)-like changes,54 and MTDH is considered as a promising biomarker and therapeutic target for gliomas. However, the specific role of MTDH in glioma pathogenesis remains unclear. Heme oxygenase 1 (HMOX1) is high-expressed in the hypoxic region of GBM and triggers mitophagy cell death of glioma cells.54,55 On the other hand, silencing HMOX1 can suppress the invasion of GBM,56,57 but the mechanism remains to be further studied.

We found significant immune infiltration in tumor tissues of high-risk patients. Noticeably, high-risk patients had higher tumor purity and immune score, indicating that despite a high permeability of immune cells, immune cells are mainly concentrated around the stroma and fail to penetrate tumor tissue to clear tumor cells.58,59 In other words, AD-FRGs may be closely related to the immune escape of glioma. It should be mentioned that the degree of infiltration of tumor-associated macrophages (TAM), such as macrophage M0 and macrophage M1, is particularly significant in patients in the high-risk group.60 In recent years, TAM has been considered as a critical factor of tumor immune escape and can cause immunosuppression. Our study showed that macrophage M0 and macrophage M1 were significantly correlated with a shorter survival of high-risk patients and can effectively distinguish between different tumor stages. In addition, we found that, as a protective gene, SIRT1 was associated with reduced macrophage M0 and macrophage M1. High levels of SIRT1 may affect the prognosis by inhibiting immune infiltration, but the specific mechanism needs to be confirmed by further research. The high expression of HMOX1, CISD2, and HSPb1 were related to increased proportion of macrophage M0 and macrophage M1, which is indicative of a poor prognosis.

We found significant immune infiltration in tumor tissues of high-risk patients. Moreover, high-risk patients had higher tumor purity and immune score. It has been found that despite the high permeability of immune cells, the immune cells are mainly concentrated around the stroma and unable to penetrate tumor tissue to clear the tumor cells.58,59 In other words, AD-FRGs may be closely related to the immune escape of glioma. Interestingly, the degree of infiltration of tumor-associated macrophages (TAM), such as macrophages M0 and M1, is particularly significant in patients in the high-risk group.60 In recent years, TAM has been considered as a critical factor of tumor immune escape and can cause immunosuppression. Our study showed that macrophage M0 and macrophage M1 were significantly correlated with shorter survival of high-risk patients and could effectively distinguish different tumor stages. In addition, we found that down-regulation of SIRT1 was associated with reducing macrophage M0 and macrophage M1. The high expression of HMOX1, CISD2, and HSPb1 were related to increased proportion of macrophage M0 and macrophage M1, which is associated with a poor prognosis.

CD4+ T cells play crucial roles for host defense and immune-mediated disease by their ability to differentiate into specialized subpopulations, such as Th1 cells, Th2 cells, Th17 cells and regulatory T cells.61 Among them, Th1 cells are strongly associated with good clinical outcome for almost all cancer types.62 Regulatory T cells, characterized by the CD4+ CD25+ FOXP3+ phenotype, are believed to dampen T-cell immunity and to be the main obstacle tempering immunotherapy.63 In the present study, we found that all CD4+ T cells could be an adverse predictor for patients with glioma. A possible explanation for the detrimental effect of CD4+ T cells on patient prognoses is that protumoral T cell subpopulations (regulatory T cells) might be more prevalent than those with antitumoral phenotypes (eg, Th1 cells) in glioma tissues. Testing this hypothesis may be the subject of further investigations.

We first screened autophagy-dependent ferroptosis-related genes and established a prognostic prediction model with a nine-genes signature. Our prognostic model exhibited good predictive performance and effectively discriminated high-risk patients with glioma. However, our study has certain limitations. Firstly, we failed to obtain certain important clinicopathological information from the TCGA database. Moreover, the practicality and stability of the nomogram have not been further validated in clinical studies. Lastly, the mechanism by which AD-FRGs affect the immune microenvironment needs further investigation.

Conclusion

In conclusion, we developed an AD-FRG signature, which maybe help us understand the role of AD-FRG in glioma and facilitate decision-making in clinical practice. In addition, we explored the relationship between AD-FRG and immune infiltration, SIRT1, CISD2 and HSPB1 might be related to macrophage infiltration and immunotolerance in glioma tissues, providing a new target for immunotherapy.