Research Paper Volume 13, Issue 13 pp 17655—17672

Characterization of ferroptosis signature to evaluate the predict prognosis and immunotherapy in glioblastoma

Xiaopeng Zhu1, , Yuxiang Zhou1, , Yangqian Ou1, , Zebo Cheng1, , Deqing Han1, , Zhou Chu2, , Sian Pan3, ,

  • 1 Department of Neurosurgery, Zhuzhou Central Hospital, Zhuzhou 412000, Hunan Province, PR China
  • 2 Department of Operating Theatre, Zhuzhou Central Hospital, Zhuzhou 412000, Hunan Province, PR China
  • 3 Department of Rehabilitation Medicine, Zhuzhou Central Hospital, Zhuzhou 412000, Hunan Province, PR China

Received: March 29, 2021       Accepted: June 19, 2021       Published: July 9, 2021      

https://doi.org/10.18632/aging.203257
How to Cite

Copyright: © 2021 Zhu et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

Background: Glioblastoma (GBM) is the most common type of brain cancer with poor survival outcomes and unsatisfactory response to current therapeutic strategies. Recent studies have demonstrated that ferroptosis-related genes (FRGs) are linked with the occurrence and development of GBM and may become promising biological indicators in GBM therapy.

Methods: We systematically assessed the relationship between FRGs expression profiles and prognosis in glioma patients based on the Cancer Genome Atlas (TCGA) and Chinese Glioma Genome Atlas (CGGA) datasets to establish a risk score model according to the gene signature of multiple survival-associated DEGs. Further, the differences between the tumor microenvironment score, immune cell infiltration, immune checkpoint expression levels, and drug sensitivity in the high- and low-risk group are analyzed through a variety of algorithms in R software.

Results: GBM patients were divided into two subgroups (high- and low-risk) according to the established risk score model. Patients in the high-risk group showed significantly reduced overall survival compared with those in the low-risk group. Also, we found that the high-risk group showed higher ImmuneScore and StromalScore, while different subgroups have significant differences in immune cell infiltration, immune checkpoint expression levels, and drug sensitivity. In summary, we developed and validated an FRGs risk model, which served as an independent prognostic indicator for GBM. Besides, the two subgroups divided by the model have significant differences, which provides novel insights for further studies as well as the personalized treatment of patients.

Introduction

Clinically, glioblastoma (GBM) ranks as most prevalent central nervous system cancer comprising 55.4% of all and 15% of all the central nervous systems, with an occurring percentage of 3.2% out of 100,000 people [14]. World Health Organization (WHO) has classified glioblastoma as the grade IV of astrocytic tumors, characterized by mitotic activity, anaplasia, cytological atypia, microvascular proliferation, and necrosis [57]. In recent decades, although there are a variety of therapeutic strategies including surgery, radiotherapy, and chemotherapy, it remains incurable with recurrence and poor prognosis, and infected patients after diagnosis survive for 15 months with a 5.5% 5-year survival rate [810]. Meanwhile, the increasing incidence of glioma raises the importance and urgency of its diagnosis and therapies. To further improve the survival period and quality of life of patients, accurate molecular biomarkers are urgently needed.

Recently, ferroptosis is an emerging approach that has fetched much attention from researchers. It is considered an iron-dependent cell death, which is triggered by high levels of lipid hydroperoxide [1114]. Ferroptosis contribute to occurrence and development of diverse disorders e.g., blood diseases, kidney damage, nervous system diseases, and cancers [1517]. For instance, non-thermal plasma (NTP) splits ferritin and reduces Fe3+ to Fe2+, eliminating oral squamous cell carcinoma cells [18]. In lung cancer, blockade of NFS1, induces iron-starvation response and ferroptosis [19]. Besides, withaferin A can kill tumors in high-risk neuroblastoma via Kelch-like ECH-associated protein 1 (KEAP1) -Nrf2 axis, a noncanonical ferroptosis pathway [20]. Furthermore, it is reported that cisplatin together with inhibition of GPX4 can initiate ferroptosis and synergistically improved chemotherapeutic efficacy in GBM [21]. However, the effect of ferroptosis and ferroptosis-related genes in GBM is still not well studied, and they are great treatment value in GBM that is worthy of further researches.

In the current investigation, the features of ferroptosis-related genes (FRGs) in GBM were characterized using data from the TCGA and CGGA data sets. An individualized signature of FRGs was constructed and validated for GBM patients, which holds promising prospects for diagnosis and prognosis application in the future. In addition, the FRGs model classified GBM patients into two subgroups. Different subgroups have significant differences in tumor microenvironment score, immune cell infiltration, and immune checkpoint expression levels, which may provide help for the development of GBM novel immunotherapy. Also, we investigated candidate drugs targeting this FRGs-related signature via the publicly available drug sensitivity database. At last, we preliminarily performed the in vitro protein expression level of the genes in our model through western blot.

Materials and Methods

Data collection

RNA sequencing data (TCGA-GBM) of 169 GBM samples and 5 normal tissue samples with their clinical details were retrieved from TCGA (http://cancergenome.nih.gov/). However, some samples were excluded which were without survival time and finally 165 GBM samples from TCGA were considered in this study (Supplementary Table 1). Using DESeq2 package, we carried out preprocessing of raw data. Subsequently, the |log2 fold change (FC)| ≥ 1 and false discovery rate (FDR) < 0.05 were employed to identify differentially expressed FRGs. The clinical information and RNA-seq transcriptome data of the samples (dataset ID: mRNAseq_325, dataset ID: mRNAseq_693) were retrieved from the CGGA (http://www.cgga.org.cn). In this study, TCGA data was used as the training set, and the CGGA data was used as the validation set. Next, 60 genes involved in ferroptosis were obtained from published studies [2225] which are outlined in Supplementary Table 2.

Design of a PPI network

With the aim of analyzing the protein-protein interaction information between different genes, the STRING (https://string-db.org/) online platform was utilized to evaluate the interactions among those FRGs and hided disconnected nodes in the network. The Cytoscape_v3.7.0 software was applied in the construction of a visual PPI network [2629].

KEGG pathway and GO enrichment analysis

The “clusterProfiler” R package [30] was applied in the enrichment analysis of differentially expressed FRGs based on Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) and analyses. P < 0.05 represented statistical significance.

Prognostic model construction and evaluation

Based on the FRGs in TCGA first selected prognosis-related FRGs by univariate Cox regression, a multivariate Cox proportional hazards regression model was designed for predicting GBM patients’ prognosis. In this model, we determined risk scores for each sample according to the following formula:

RiskScore=i=1nExpiβi

Where Exp is the gene expression value, β stands for the regression coefficient [3133].

Further, the GBM patients were assigned into low- and high-risk groups based on risk score cutoff values. The groups were subjected to survival analysis using the Kaplan-Meier method with the log-rank test. In addition, an ROC curve was developed using the “SurvivalROC” R package which was then applied in evaluation of the prognosis prediction ability of the designed model. A nomogram was drawn to assess the OS based on the “rms” R package. The performance of the nomogram was further clarified using CGGA as the validation cohort. P < 0.05 represented statistical significance.

Generation of ImmuneScore and StromalScore

The ratio of immune to stromal components for every sample in the tumor microenvironment was determined using the ESTIMATE package [34], which was shown as two scores: StromalScore and ImmuneScore. Higher scores correlated with higher ratios of in the tumor microenvironment.

Calculation of immune cell type fractions

The CIBERSORT was applied for estimating 22 immune cell types fractions (neutrophils, eosinophils, activated mast cells, resting mast cells, activated dendritic cells, resting dendritic cells, macrophages M2, macrophages M1, macrophages M0, monocytes, activated NK cells, resting NK cells, gamma delta T cells, T cells regulatory, follicular helper T cells, activated memory CD4 T cells, resting memory CD4 T cells, naive CD4 T cells, CD8 T cells, plasma cells, memory B cells, and naive B cells) between subjected with high -and low-risk scores.

GSEA functional analysis

To characterize pathways associated with low- and high-risk subgroups, the Gene Set Enrichment Analysis (GSEA software, version 4.0.1) was applied. We set number of random sample permutations at 1000, and significance cut-off values at P < 0.05.

Cell culture

The normal human glial cell line (HEB) was bought from Otwo Biotech (China). HS 683 cell line, H4 cell line, and U251 cell line were all purchased from the Cell Resource Center, Peking Union Medical College. These cells were passaged with DMEM (Hyclone), enriched with 10% FBS (Hyclone) in 37 Celsius and 5% CO2 incubator.

Western blot

The protein expression levels were determined as previously reported with some modifications [3537]. Antibodies used were: Anti-STEAP3 (1/1000 dilution, Abcam, ab151566, UK), Anti-CRYAB (1/1000 dilution, Abcam, ab76467, UK), Anti-MT1G (1/1000 dilution, Omnimab, #OM263051, USA).

Statistical analysis

The “survminer” package in R was applied in univariate and multivariate Cox regression analyses. The ROC curve and AUC were obtained using the “timeROC” package in R. The “ggplotify” “VennDiagram,” “maftools,” “plot3D,” “cowplot,” “ggforest,” and “ggplot2,” packages in R were used for visualization. The results were expressed as mean ± SEM. Statistical analysis was done using SPSS.22 and R software version 4.0.3. P < 0.05 was chosen as the threshold of statistically significance.

Ethics approval and consent to participate

Animal and human experiments were not conducted in this study.

Availability of data and material

Technical appendix, statistical code, and dataset are available from the corresponding author at [090102080@m.fafu.edu.cn].

Results

Determination of differently expressed FRGs in GBM patients

Firstly, we extracted FRGs from the TCGA-GBM as presented in (Figure 1A). The differential gene expression of FRGs between the two groups identified 17 upregulated and 15 downregulated FRGs (Figure 1B). The STRING and Cytoscape were used to visualize the interactions among FRGs. Moreover, FRGs disconnected nodes in the network were not shown (Figure 1C).

Identification of FRGs in GBM. (A) A Venn diagram indicating that 60 FRGs were identified in the TCGA-GBM cohorts. (B) Volcano plot showing DEGs among FRGs in GBM. (C) A PPI network on the relationship between up-regulated and down-regulated DEGs. Blue or yellow are up-or down-regulated DEGs, respectively. (D) GO analysis of up-regulated DEGs. (E) KEGG analysis of up-regulated DEGs. (F) GO analysis of down-regulated DEGs. (G) KEGG analysis of down-regulated DEGs.

Figure 1. Identification of FRGs in GBM. (A) A Venn diagram indicating that 60 FRGs were identified in the TCGA-GBM cohorts. (B) Volcano plot showing DEGs among FRGs in GBM. (C) A PPI network on the relationship between up-regulated and down-regulated DEGs. Blue or yellow are up-or down-regulated DEGs, respectively. (D) GO analysis of up-regulated DEGs. (E) KEGG analysis of up-regulated DEGs. (F) GO analysis of down-regulated DEGs. (G) KEGG analysis of down-regulated DEGs.

Results of GO and KEGG analysis

Based on GO analysis, upregulated differently expressed FRGs were strongly linked to cellular response to oxidative stress, NADP metabolic process, regulation of oxidative stress−induced intrinsic apoptotic signaling pathway, cellular response to oxidative stress, response to oxidative stress, intrinsic apoptotic signaling pathway, apoptotic signaling pathway, cellular response to chemical stress, cell−substrate junction, focal adhesion, invadopodium, lamellipodium membrane, and intrinsic apoptotic signaling pathway (Figure 1D).

The GO results exhibited that downregulated differently expressed FRGs showed strongly enrichment in sulfur compound biosynthetic process, organic acid transport, carboxylic acid transport, several metabolic processes including purine nucleoside bisphosphate, fatty−acyl−CoA, ribonucleoside bisphosphate, nucleoside bisphosphate, and glutamate, fatty−acyl−CoA biosynthesis, organic hydroxy compound biosysnthesis, ER membrane protein complex, autolysosome, lipid droplet, microbody, peroxisome, peroxisomal membrane, medium−chain fatty acid−CoA ligase activity, arachidonate−CoA ligase activity, long−chain fatty acid−CoA ligase activity, and decanoate−CoA ligase activity (Figure 1F).

Moreover, the upregulated differently expressed FRGs were mainly enriched in Ferroptosis, Fluid shear stress, and atherosclerosis (Figure 1E), while downregulated FRGs were significantly enriched for Ferroptosis, Proximal tubule bicarbonate reclamation, Terpenoid backbone biosynthesis, Steroid biosynthesis, 2−Oxocarboxylic acid metabolism, Phenylalanine metabolism, AMPK signaling pathway, Peroxisome, PPAR signaling pathway, Adipocytokine signaling pathway, Fatty acid degradation, Alanine, aspartate and glutamate metabolism, Arginine biosynthesis, Fatty acid metabolism, Fatty acid biosynthesis (Figure 1G).

Prognosis-related FRGs selection

A total of 6 prognostic-related hub FRGs were found (Figure 2A). Three of the six hub FRGs could independently predict the outcomes of GBM patients (Figure 2B, Table 1).

Prognostic significance of the FRGs signature derived risk scores. (A, B) Univariate and multivariate Cox analysis evaluating the prognostic-related genes in the TCGA (A) and CGGA cohort (B). The Kaplan-Meier survival curves for the high- and low-risk groups in TCGA (C) and CGGA cohort (D). (E, I) The predictive efficiency of the FRGs risk signature on the 1-, 3-, and 5-years survival rate in TCGA (E) and CGGA cohort (I) via ROC curve. (F, J) Heat maps of these three FRGs (CRYAB, MT1G, STEAP3) expression profiles in TCGA (F) and CGGA cohort (J). (G–I, K–L) Distribution of risk score and patient survival time, and status of GBM in TCGA (G, H) and CGGA cohort (K, L). The black dotted line is the optimal cut-off value for dividing patients into low-risk and high-risk groups. (M–P) Univariate and multivariate Cox analyses for evaluating the independent prognostic value of the FRGs signature in TCGA (M, N) and CGGA cohort (O, P).

Figure 2. Prognostic significance of the FRGs signature derived risk scores. (A, B) Univariate and multivariate Cox analysis evaluating the prognostic-related genes in the TCGA (A) and CGGA cohort (B). The Kaplan-Meier survival curves for the high- and low-risk groups in TCGA (C) and CGGA cohort (D). (E, I) The predictive efficiency of the FRGs risk signature on the 1-, 3-, and 5-years survival rate in TCGA (E) and CGGA cohort (I) via ROC curve. (F, J) Heat maps of these three FRGs (CRYAB, MT1G, STEAP3) expression profiles in TCGA (F) and CGGA cohort (J). (GI, KL) Distribution of risk score and patient survival time, and status of GBM in TCGA (G, H) and CGGA cohort (K, L). The black dotted line is the optimal cut-off value for dividing patients into low-risk and high-risk groups. (MP) Univariate and multivariate Cox analyses for evaluating the independent prognostic value of the FRGs signature in TCGA (M, N) and CGGA cohort (O, P).

Table 1. Univariate and multivariate Cox regression analysis.

GeneUnivariate analysisMultivariate analysis
HR (95% CI)pHR (95% CI)p
HSPB11.002 (1.000–1.004)0.029
CRYAB0.999 (0.999–1.000)0.0170.852 (0.711–1.022)0.085
MT1G1.027 (1.011–1.043)0.0001.269 (1.032–1.561)0.024
NCOA40.981 (0.963–0.998)0.032
PTGS21.084 (1.005–1.169)0.035
STEAP31.017 (1.005–1.029)0.0061.226 (1.003–1.499)0.047

Construction and validation of prognosis risk score model

A prognosis prediction model was designed from the 3 hub FRGs. Using the model, a risk score for a given patient could be determined as:

RiskScore=i=1nExpiβi

The model assigned the TCGA-GBM patients into low- and high-risk subgroups. The findings revealed that those in high-risk subgroup had worse OS relative to those in low-risk subgroup as displayed in Figure 2C. Analysis of a time-dependent ROC and AUC) of the FRGs model found a prognostic ability of 0.708 (one-year), 0.707 (two-year), and 0.683 (three-year), indicating a good performance (Figure 2E). Figure 2F2H displays the heat map, risk score, and survival status of subjects in the 3 FRGs in both subgroups. We applied the formula to validate the prognostic significance of the model in other GBM patient cohorts using CGGA dataset as validation. Similarly, those with high-risk scores showed poorer OS than subjects with low-risk scores in the CGGA dataset (Figure 2D). Figure 2J2L shows the survival status, risk score, and expression heat map of CGGA cohorts both risk groups, the FRGs signature built in this study has stable prognostic ability and a similar tendency of gene expression in both groups. The results presented in Figure 2M2P depicted that risk score can could independently assess the progress of GBM patients. Hence the model had satisfactory specificity and sensitivity.

Design of a nomogram and drug relevance using 3 hub FRGs

Briefly, the 3 FRGs were applied to design a nomogram (Figure 3A). The calibration curve results as presented in Figure 3B3D exposed that the survival rate obtained by the model was nearly equal to the actual survival rate. Additionally, the relationship between level of the three FRGs and drugs was explored (Figure 3E). The drug data is obtained from the cellminer database (https://discover.nci.nih.gov/cellminer/loadDownload.do).

Construction of a nomogram and drug relevance based on the 3 hub FRGs. (A) Validation of the nomogram in the TCGA cohort. (B–D) Calibration maps used to predict the 1-year (B), 3-year (C), and 5-year survival (D). (E) The correlation between gene expression levels and drugs. The top 16 most relevant were visualized.

Figure 3. Construction of a nomogram and drug relevance based on the 3 hub FRGs. (A) Validation of the nomogram in the TCGA cohort. (BD) Calibration maps used to predict the 1-year (B), 3-year (C), and 5-year survival (D). (E) The correlation between gene expression levels and drugs. The top 16 most relevant were visualized.

Significant differences between low- and high-risk patients

Next, the patients were scored by the prognostic FRGs model, and grouped into low-risk and high-risk groups on the basis of median score. The results of PCA and t-SNE supported the classification of GBM patients into two subgroups by our FRGs signature (Figure 4A4B). Additionally, verification in the CGGA cohort was also carried out and the outcomes are shown in Figure 4C4D.

Analysis of differences between high- and low-risk subgroups (tumor microenvironment, immune cell infiltration, and immune checkpoint regulators). PCA (A) and t-SNE (B) analysis supported the stratification into two GBM subclasses the TCGA cohort. PCA (C) and t-SNE (D) analysis supported the stratification into two GBM subclasses the CGGA cohort. The comparison of stromal scores (E) and immune scores (F) in high- and low-risk subgroups. (G) The comparison of immune cell fractions between high- and low-risk subgroups. (H) The pathways enriched in high-risk GBM through GSEA analysis by enrichment map. (I–P) The key Immune checkpoint regulators with significant differential expression in the high- and low-risk subgroups.

Figure 4. Analysis of differences between high- and low-risk subgroups (tumor microenvironment, immune cell infiltration, and immune checkpoint regulators). PCA (A) and t-SNE (B) analysis supported the stratification into two GBM subclasses the TCGA cohort. PCA (C) and t-SNE (D) analysis supported the stratification into two GBM subclasses the CGGA cohort. The comparison of stromal scores (E) and immune scores (F) in high- and low-risk subgroups. (G) The comparison of immune cell fractions between high- and low-risk subgroups. (H) The pathways enriched in high-risk GBM through GSEA analysis by enrichment map. (IP) The key Immune checkpoint regulators with significant differential expression in the high- and low-risk subgroups.

The ESTIMATE algorithm was adopted to assess the TCGA-GBM tumor microenvironment, and the results showed higher ImmuneScore and StromalScore in high-risk patients (Figure 4E4F). The Supplementary Table 3 outlines the features of patients in both risk groups. The CIBERSORT algorithm analysis results found 21 immune cell types in TCGA-GBM samples (Figure 4G). Through GSEA enrichment analysis revealed that high-risk patient had enrichment in ferroptosis-related pathways, such as AMINO SUGAR AND NUCLEOTIDE SUGAR METABOLISM, APOPTOSIS, CELL CYCLE, ETHER LIPID METABOLISM, GALACTOSE METABOLISM, and GLYCOLYSIS GLUCONEOGENESIS (Figure 4H, Table 2).

Table 2. Gene set enrichment analysis in high-risk group.

NAMEESNESNOM p-valFDR q-val
KEGG_CELL_CYCLE–0.590–1.7520.0210.494
KEGG_RNA_DEGRADATION–0.548–1.6010.0360.433
KEGG_AMINO_SUGAR_AND_NUCLEOTIDE_SUGARMETABOLISM0.6141.8120.0050.052
KEGG_APOPTOSIS0.5931.99400.011
KEGG_ETHER_LIPID_METABOLISM0.4481.4630.0400.187
KEGG_GLYCOLYSIS_GLUCONEOGENESIS0.4831.6010.0180.139
KEGG_GALACTOSE_METABOLISM0.6091.7570.0010.062
KEGG_GLYCOSAMINOGLYCAN_DEGRADATION0.7921.98000.010
KEGG_GLYCOSPHINGOLIPID_BIOSYNTHESIS_GANGLIO_SERIES0.6761.7070.0090.087
KEGG_JAK_STAT_SIGNALING_PATHWAY0.4621.5810.0200.148
KEGG_LYSOSOME0.6712.05000.011
KEGG_STARCH_AND_SUCROSE_METABOLISM0.4631.5990.0220.135
NOM P value < 0.01 was statistically significant. Abbreviations: FDR, False discovery rate.

In previous researches, immune checkpoint inhibitors e.g., PD-L1, PD-L2, and LAG3, were proposed as treatments for cancer. In the present study, CD274 (PD-L1), STEAP3, IL1B, PDCD1LG2 (PD-L2), and MT1G expressions were enriched in patients with high risk scores whereas the expression of CRYAB, LAG3, and IL12A was markedly enriched in those with low risk (Figure 4I4P).

To further study characterize drug responses in patients, the R package “pRRophetic” [35] was adopted to determine the half-maximal inhibitory concentration (IC50) of every GBM patient using the Genomics of Drug Sensitivity in Cancer (GDSC) website. Consequently, 24 drugs showed distinct estimated IC50 between low and high-risk GBM patients (Figure 5A5X).

Drug sensitivity analysis to drugs of high- and low-risk subgroups. Differential chemotherapeutic responses in high- and low-risk patients (A–X).

Figure 5. Drug sensitivity analysis to drugs of high- and low-risk subgroups. Differential chemotherapeutic responses in high- and low-risk patients (AX).

Validation of the three genes in GBM cells

Western blot assay was conducted to quantify the protein level of the three genes in three GBM cell lines (HS 683 cells, H4 cells, and U251 cells), and in an HEB cell line (Figure 6A). Results presented in Figure 6 indicate that CRYAB (Figure 6B) and STEAP3 (Figure 6D) were upregulated in GBM cells than in HEB cell line, whereas MT1G (Figure 6C) was downregulated.

Validation of the differential expression of the three genes in GBM cells. (A) Western blot images and the relevant quantification (B–D) of CRYAB, MT1G, and STEAP3. Data are shown as mean ± SEM from three independent experiments, *P **P ***P

Figure 6. Validation of the differential expression of the three genes in GBM cells. (A) Western blot images and the relevant quantification (BD) of CRYAB, MT1G, and STEAP3. Data are shown as mean ± SEM from three independent experiments, *P < 0.05, **P < 0.01, ***P < 0.001.

Discussion

GBM is one of the deadliest cancers with only a few approved therapies at present. Despite extensive studies on it over the last several decades, it remains incurable and fatal disease with little improvement in the survival rates for patients [3639]. Thus, it boosted researchers to further explore on the mechanisms of GBM for novel prognostic biomarkers and therapeutic targets.

Ferroptosis is a form of cell death driven by iron dependent phospholipid peroxidation, which completely differs from other cell death such as necrosis, apoptosis, and autophagy for unique cell morphology, gene expression, and metabolic pathways. Cells undergoing ferroptosis gain a round shape before death with its distinctive feature, the smaller mitochondria, and increased mitochondrial membrane density. Previously, many pieces of researches have revealed that cancer cells in several cancers e.g., ovarian, lung, and breast cancers have high levels of iron relative to normal cells [4045]. Moreover, it is reported that salinomycin can kill cancer stem cells by sequestering iron in lysosomes causing lysosomal membrane permeabilization and ferroptosis, suggesting that therapies targeting ferroptosis might be a possible strategy for cancer [46]. In another research, authors concluded that withaferin A can kill tumors in high-risk neuroblastoma via KEAP1 -Nrf2 axis, a noncanonical ferroptosis pathway [20]. Also, compelling evidence reported the lethality of GPX4 inhibitors in drug-resistant cells via ferroptosis, supporting that targeting GPX4 could act as a potential therapy to prevent drug resistance [4749]. And in GBM, cisplatin co-delivered with small interference RNAs of GPX4 showed a significant superior therapeutic effect through a mechanism related to ferroptosis compared with cisplatin only in vitro and in vivo [21]. These evidences support the perspective applications of ferroptosis in GBM therapies. Nevertheless, the mechanisms of ferroptosis in GBM and its potential ability to prognosis are still not clear.

Therefore, we carried out current study to systematically evaluate the significance of FRGs in GBM with advanced computational tools. We assessed the expression profile of FRGs in GBM to get 32 DEGs. The relationship within these differentially expressed FRGs visualized by the PPI network suggesting TP53 with the most connected lines might play an important part in ferroptosis of cancer cells, while TP53 has already been one of the crucial dysregulated genes in most cancer types, indicating ferroptosis could also be one fundamental mechanism in GBM development. Next, we conducted GO and KEGG in the upregulated FRGs and downregulated FRGs. Results displayed that the upregulated FRGs are enriched in responses to oxidative stress except for ferroptosis itself. The close connection of FRGs and oxidative stress in GBM is readily comprehensible for the definition of ferroptosis cell death as an iron-dependent accumulation of lipid peroxidation. Some radical-trapping antioxidants can prevent ferroptosis. Meanwhile, the enrichment of the downregulated FRGs is associated with metabolism and biosynthesis as well as metabolism-related signaling pathways like PPAR and AMPK pathways. At present, metabolic changes are essential evidence for ferroptosis research. For example, iron abundance and lipid peroxidation level are two of the most critical indicators of ferroptosis. These findings revealed that ferroptosis, as well as FEGs, may have a significant impact on GBM development and give direction for in-depth studies on ferroptosis in GBM.

Based on the 32 differentially expressed FRGs, we observed 3 prognosis-associated candidate hub FRGs (CRYAB, MT1G, and STEAP3) as independent predictors to design a risk score model for prognostic prediction. The strong correlation between the low-risk subgroup of GBM and high survival rate suggested that those patients with low risk predicted by the FRGs signature are prone to a better prognosis with decent efficacy in both TCGA and CGGA cohorts. Furthermore, we performed various analyses to assess the correlation of FRGs signature and tumor microenvironment, which displayed higher ImmuneScore and StromalScore in the high-risk subgroup. More specifically, patients with high risk have a higher fraction of B cell memory (p = 0.024), T cell CD4 memory activated (p = 0.009), NK cells resting (p = 0.032), monocytes (p = 0.014), dendritic cells activated (p = 0.027), and neutrophils (p < 0.001), while those with low risk generate more T cells follicular helper (p < 0.001), and macrophages M1 (p = 0.013). For gene expressions of key immune checkpoints, there were significant differences in CD274, STEAP3, CRYAB, IL1B, PDCD1LG2, LAG3, IL12A, and MT1G in high- and low-risk subgroups. Our results uncovered that the FRGs signature may participate in tumor immunity and guide stratification and therapeutic strategies of immunotherapies in the future. In recent times, immunotherapies have made great progress in several solid tumors but could not improve the survival rate of GBM patients. Integrating immunotherapies with ferroptosis can provide a new insight to solve the challenges of current immunotherapies. Previous research found that during the anti-PD-L1 treatment, the dramatically elevated lipid peroxidation specific to ferroptosis along with inhibition of ferroptosis signaling cascades contributed thereby reducing tumor cells sensitivity. Above all, the functional mechanism of CRYAB, MT1G, and STEAP3 in ferroptosis and immunity should be further explored for potential treatments of GBM.

Thus, we analyzed the correlation between gene expression of CRYAB, MT1G, and STEAP3 and drugs, and also investigated the drug sensitivity in low- and high-risk subgroups, introducing possible drugs for certain subgroups divided by FRGs signature. Combination therapy with checkpoint inhibitors and drugs has proven to be a reliable therapy and is associated with a better prognosis [5053].

At last, we preliminarily measured the protein expression of CRYAB, MT1G, and STEAP3 in GBM cell lines (HS 683, H4, and U251) compared with that in HEB cells. The rise of CRYAB and STEAP3 expression and fall of MT1G expression in vitro are consistent with their RNA level in TCGA sequencing data, supporting our bioinformatics analysis.

In summary, our novel FRGs-associated prognostic model for GBM could greatly provide access to the ferroptosis and pathogenesis in GBM and guide new prognostic biomarkers as well as therapeutic strategies for GBM.

Conclusions

Overall, our study comprehensively analyzed the FRGs in GBM and builds an FRGs model (CRYAB, MT1G, and STEAP3) for prognosis and stratification of GBM patients. Subgroups with relative low or high-risk classified by the model have differences in ImmunoScore and StromalScore in tumor microenvironment and fraction of multiple immune cells, expression of immune checkpoint genes, and drug sensitivity. These findings can provide new insights for the development of new immunotherapy for GBM.

Supplementary Materials

Supplementary Tables

Author Contributions

Zhou Chu and Sian Pan designed the study, analyzed data, and wrote the manuscript. Yuxiang Zhou, Yangqian Ou, Zebo Cheng, and Deqing Han analyzed data and contributed in writing the manuscript. Xiaopeng Zhu performed the experiments, analyzed data, and wrote the manuscript.

Acknowledgments

We would like to thank everyone who contributed to this article. We thank Dr. Li Ming from Beijing Foreign Studies University for his assistance in language modification.

Conflicts of Interest

The authors declare no conflicts of interest related to this study.

Funding

This work was supported by the Scientific research project (2019) of health commission of Hunan (B2019200) and the Science and technology innovation project of Hunan (2018SK52802).

References

  • 1. Grimm SA, Chamberlain MC. Brainstem glioma: a review. Curr Neurol Neurosci Rep. 2013; 13:346. https://doi.org/10.1007/s11910-013-0346-3 [PubMed]
  • 2. Ostrom QT, Bauchet L, Davis FG, Deltour I, Fisher JL, Langer CE, Pekmezci M, Schwartzbaum JA, Turner MC, Walsh KM, Wrensch MR, Barnholtz-Sloan JS. The epidemiology of glioma in adults: a "state of the science" review. Neuro Oncol. 2014; 16:896–913. https://doi.org/10.1093/neuonc/nou087 [PubMed]
  • 3. Chen R, Smith-Cohn M, Cohen AL, Colman H. Glioma Subclassifications and Their Clinical Significance. Neurotherapeutics. 2017; 14:284–97. https://doi.org/10.1007/s13311-017-0519-x [PubMed]
  • 4. Miller JJ, Shih HA, Andronesi OC, Cahill DP. Isocitrate dehydrogenase-mutant glioma: Evolving clinical and therapeutic implications. Cancer. 2017; 123:4535–46. https://doi.org/10.1002/cncr.31039 [PubMed]
  • 5. Xu S, Tang L, Li X, Fan F, Liu Z. Immunotherapy for glioma: Current management and future application. Cancer Lett. 2020; 476:1–12. https://doi.org/10.1016/j.canlet.2020.02.002 [PubMed]
  • 6. GLASS Consortium. Glioma through the looking GLASS: molecular evolution of diffuse gliomas and the Glioma Longitudinal Analysis Consortium. Neuro Oncol. 2018; 20:873–84. https://doi.org/10.1093/neuonc/noy020 [PubMed]
  • 7. Oberndorfer S, Hutterer M. Palliative care in glioma management. Curr Opin Oncol. 2019; 31:548–53. https://doi.org/10.1097/CCO.0000000000000584 [PubMed]
  • 8. Poff A, Koutnik AP, Egan KM, Sahebjam S, D'Agostino D, Kumar NB. Targeting the Warburg effect for cancer treatment: Ketogenic diets for management of glioma. Semin Cancer Biol. 2019; 56:135–48. https://doi.org/10.1016/j.semcancer.2017.12.011 [PubMed]
  • 9. Alarcón S, Niechi I, Toledo F, Sobrevia L, Quezada C. Glioma progression in diabesity. Mol Aspects Med. 2019; 66:62–70. https://doi.org/10.1016/j.mam.2019.02.002 [PubMed]
  • 10. Jung E, Alfonso J, Osswald M, Monyer H, Wick W, Winkler F. Emerging intersections between neuroscience and glioma biology. Nat Neurosci. 2019; 22:1951–60. https://doi.org/10.1038/s41593-019-0540-y [PubMed]
  • 11. Liang C, Zhang X, Yang M, Dong X. Recent Progress in Ferroptosis Inducers for Cancer Therapy. Adv Mater. 2019; 31:e1904197. https://doi.org/10.1002/adma.201904197 [PubMed]
  • 12. Li J, Cao F, Yin HL, Huang ZJ, Lin ZT, Mao N, Sun B, Wang G. Ferroptosis: past, present and future. Cell Death Dis. 2020; 11:88. https://doi.org/10.1038/s41419-020-2298-2 [PubMed]
  • 13. Hassannia B, Vandenabeele P, Vanden Berghe T. Targeting Ferroptosis to Iron Out Cancer. Cancer Cell. 2019; 35:830–49. https://doi.org/10.1016/j.ccell.2019.04.002 [PubMed]
  • 14. Zhuo S, Chen Z, Yang Y, Zhang J, Tang J, Yang K. Clinical and Biological Significances of a Ferroptosis-Related Gene Signature in Glioma. Front Oncol. 2020; 10:590861. https://doi.org/10.3389/fonc.2020.590861 [PubMed]
  • 15. Tang R, Hua J, Xu J, Liang C, Meng Q, Liu J, Zhang B, Yu X, Shi S. The role of ferroptosis regulators in the prognosis, immune activity and gemcitabine resistance of pancreatic cancer. Ann Transl Med. 2020; 8:1347. https://doi.org/10.21037/atm-20-2554a [PubMed]
  • 16. Wang Y, Zhao G, Condello S, Huang H, Cardenas H, Tanner EJ, Wei J, Ji Y, Li J, Tan Y, Davuluri RV, Peter ME, Cheng JX, Matei D. Frizzled-7 Identifies Platinum-Tolerant Ovarian Cancer Cells Susceptible to Ferroptosis. Cancer Res. 2021; 81:384–99. https://doi.org/10.1158/0008-5472.CAN-20-1488 [PubMed]
  • 17. Huang HX, Yang G, Yang Y, Yan J, Tang XY, Pan Q. TFAP2A is a novel regulator that modulates ferroptosis in gallbladder carcinoma cells via the Nrf2 signalling axis. Eur Rev Med Pharmacol Sci. 2020; 24:4745–55. https://doi.org/10.26355/eurrev_202005_21163 [PubMed]
  • 18. Sato K, Shi L, Ito F, Ohara Y, Motooka Y, Tanaka H, Mizuno M, Hori M, Hirayama T, Hibi H, Toyokuni S. Non-thermal plasma specifically kills oral squamous cell carcinoma cells in a catalytic Fe(II)-dependent manner. J Clin Biochem Nutr. 2019; 65:8–15. https://doi.org/10.3164/jcbn.18-91 [PubMed]
  • 19. Alvarez SW, Sviderskiy VO, Terzi EM, Papagiannakopoulos T, Moreira AL, Adams S, Sabatini DM, Birsoy K, Possemato R. NFS1 undergoes positive selection in lung tumours and protects cells from ferroptosis. Nature. 2017; 551:639–43. https://doi.org/10.1038/nature24637 [PubMed]
  • 20. Hassannia B, Wiernicki B, Ingold I, Qu F, Van Herck S, Tyurina YY, Bayır H, Abhari BA, Angeli JPF, Choi SM, Meul E, Heyninck K, Declerck K, et al. Nano-targeted induction of dual ferroptotic mechanisms eradicates high-risk neuroblastoma. J Clin Invest. 2018; 128:3341–55. https://doi.org/10.1172/JCI99032 [PubMed]
  • 21. Zhang Y, Fu X, Jia J, Wikerholmen T, Xi K, Kong Y, Wang J, Chen H, Ma Y, Li Z, Wang C, Qi Q, Thorsen F, et al. Glioblastoma Therapy Using Codelivery of Cisplatin and Glutathione Peroxidase Targeting siRNA from Iron Oxide Nanoparticles. ACS Appl Mater Interfaces. 2020; 12:43408–21. https://doi.org/10.1021/acsami.0c12042 [PubMed]
  • 22. Liang JY, Wang DS, Lin HC, Chen XX, Yang H, Zheng Y, Li YH. A Novel Ferroptosis-related Gene Signature for Overall Survival Prediction in Patients with Hepatocellular Carcinoma. Int J Biol Sci. 2020; 16:2430–41. https://doi.org/10.7150/ijbs.45050 [PubMed]
  • 23. Stockwell BR, Friedmann Angeli JP, Bayir H, Bush AI, Conrad M, Dixon SJ, Fulda S, Gascón S, Hatzios SK, Kagan VE, Noel K, Jiang X, Linkermann A, et al. Ferroptosis: A Regulated Cell Death Nexus Linking Metabolism, Redox Biology, and Disease. Cell. 2017; 171:273–85. https://doi.org/10.1016/j.cell.2017.09.021 [PubMed]
  • 24. Bersuker K, Hendricks JM, Li Z, Magtanong L, Ford B, Tang PH, Roberts MA, Tong B, Maimone TJ, Zoncu R, Bassik MC, Nomura DK, Dixon SJ, Olzmann JA. The CoQ oxidoreductase FSP1 acts parallel to GPX4 to inhibit ferroptosis. Nature. 2019; 575:688–92. https://doi.org/10.1038/s41586-019-1705-2 [PubMed]
  • 25. Doll S, Freitas FP, Shah R, Aldrovandi M, da Silva MC, Ingold I, Goya Grocin A, Xavier da Silva TN, Panzilius E, Scheel CH, Mourão A, Buday K, Sato M, et al. FSP1 is a glutathione-independent ferroptosis suppressor. Nature. 2019; 575:693–98. https://doi.org/10.1038/s41586-019-1707-0 [PubMed]
  • 26. Xu N, Wu YP, Ke ZB, Liang YC, Cai H, Su WT, Tao X, Chen SH, Zheng QS, Wei Y, Xue XY. Identification of key DNA methylation-driven genes in prostate adenocarcinoma: an integrative analysis of TCGA methylation data. J Transl Med. 2019; 17:311. https://doi.org/10.1186/s12967-019-2065-2 [PubMed]
  • 27. Xu M, Li Y, Li W, Zhao Q, Zhang Q, Le K, Huang Z, Yi P. Immune and Stroma Related Genes in Breast Cancer: A Comprehensive Analysis of Tumor Microenvironment Based on the Cancer Genome Atlas (TCGA) Database. Front Med (Lausanne). 2020; 7:64. https://doi.org/10.3389/fmed.2020.00064 [PubMed]
  • 28. Chen L, Lu D, Sun K, Xu Y, Hu P, Li X, Xu F. Identification of biomarkers associated with diagnosis and prognosis of colorectal cancer patients based on integrated bioinformatics analysis. Gene. 2019; 692:119–25. https://doi.org/10.1016/j.gene.2019.01.001 [PubMed]
  • 29. Li L, Chen X, Hao L, Chen Q, Liu H, Zhou Q. Exploration of immune-related genes in high and low tumor mutation burden groups of chromophobe renal cell carcinoma. Biosci Rep. 2020; 40:BSR20201491. https://doi.org/10.1042/BSR20201491 [PubMed]
  • 30. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012; 16:284–87. https://doi.org/10.1089/omi.2011.0118 [PubMed]
  • 31. Cao M, Cai J, Yuan Y, Shi Y, Wu H, Liu Q, Yao Y, Chen L, Dang W, Zhang X, Xiao J, Yang K, He Z, et al. A four-gene signature-derived risk score for glioblastoma: prospects for prognostic and response predictive analyses. Cancer Biol Med. 2019; 16:595–605. [PubMed]
  • 32. Sheng Y, Yanping C, Tong L, Ning L, Yufeng L, Geyu L. Predicting the Risk of Melanoma Metastasis Using an Immune Risk Score in the Melanoma Cohort. Front Bioeng Biotechnol. 2020; 8:206. https://doi.org/10.3389/fbioe.2020.00206 [PubMed]
  • 33. Huang R, Mao M, Lu Y, Yu Q, Liao L. A novel immune-related genes prognosis biomarker for melanoma: associated with tumor microenvironment. Aging (Albany NY). 2020; 12:6966–80. https://doi.org/10.18632/aging.103054 [PubMed]
  • 34. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, Treviño V, Shen H, Laird PW, Levine DA, Carter SL, Getz G, Stemke-Hale K, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013; 4:2612. https://doi.org/10.1038/ncomms3612 [PubMed]
  • 35. Geeleher P, Cox N, Huang RS. pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS One. 2014; 9:e107468. https://doi.org/10.1371/journal.pone.0107468 [PubMed]
  • 36. Hellmark T, Segelmark M. Diagnosis and classification of Goodpasture's disease (anti-GBM). J Autoimmun. 2014; 48–49:108–12. https://doi.org/10.1016/j.jaut.2014.01.024 [PubMed]
  • 37. Wang S, Yao F, Lu X, Li Q, Su Z, Lee JH, Wang C, Du L. Temozolomide promotes immune escape of GBM cells via upregulating PD-L1. Am J Cancer Res. 2019; 9:1161–71. [PubMed]
  • 38. Mukherjee S, Fried A, Hussaini R, White R, Baidoo J, Yalamanchi S, Banerjee P. Phytosomal curcumin causes natural killer cell-dependent repolarization of glioblastoma (GBM) tumor-associated microglia/macrophages and elimination of GBM and GBM stem cells. J Exp Clin Cancer Res. 2018; 37:168. https://doi.org/10.1186/s13046-018-0792-5 [PubMed]
  • 39. Martikainen M, Essand M. Virus-Based Immunotherapy of Glioblastoma. Cancers (Basel). 2019; 11:186. https://doi.org/10.3390/cancers11020186 [PubMed]
  • 40. Shoja Z, Chenari M, Jafarpour A, Jalilvand S. Role of iron in cancer development by viruses. Rev Med Virol. 2019; 29:e2045. https://doi.org/10.1002/rmv.2045 [PubMed]
  • 41. Torti SV, Torti FM. Iron and cancer: more ore to be mined. Nat Rev Cancer. 2013; 13:342–55. https://doi.org/10.1038/nrc3495 [PubMed]
  • 42. Wei J, Gao X, Qin Y, Liu T, Kang Y. An Iron Metabolism-Related SLC22A17 for the Prognostic Value of Gastric Cancer. Onco Targets Ther. 2020; 13:12763–75. https://doi.org/10.2147/OTT.S287811 [PubMed]
  • 43. Shen Y, Li X, Zhao B, Xue Y, Wang S, Chen X, Yang J, Lv H, Shang P. Iron metabolism gene expression and prognostic features of hepatocellular carcinoma. J Cell Biochem. 2018; 119:9178–204. https://doi.org/10.1002/jcb.27184 [PubMed]
  • 44. Torti SV, Manz DH, Paul BT, Blanchette-Farra N, Torti FM. Iron and Cancer. Annu Rev Nutr. 2018; 38:97–125. https://doi.org/10.1146/annurev-nutr-082117-051732 [PubMed]
  • 45. Liu Y, Zhang X, Zhang J, Tan J, Li J, Song Z. Development and Validation of a Combined Ferroptosis and Immune Prognostic Classifier for Hepatocellular Carcinoma. Front Cell Dev Biol. 2020; 8:596679. https://doi.org/10.3389/fcell.2020.596679 [PubMed]
  • 46. Mai TT, Hamaï A, Hienzsch A, Cañeque T, Müller S, Wicinski J, Cabaud O, Leroy C, David A, Acevedo V, Ryo A, Ginestier C, Birnbaum D, et al. Salinomycin kills cancer stem cells by sequestering iron in lysosomes. Nat Chem. 2017; 9:1025–33. https://doi.org/10.1038/nchem.2778 [PubMed]
  • 47. Yang WS, SriRamaratnam R, Welsch ME, Shimada K, Skouta R, Viswanathan VS, Cheah JH, Clemons PA, Shamji AF, Clish CB, Brown LM, Girotti AW, Cornish VW, et al. Regulation of ferroptotic cancer cell death by GPX4. Cell. 2014; 156:317–31. https://doi.org/10.1016/j.cell.2013.12.010 [PubMed]
  • 48. Seibt TM, Proneth B, Conrad M. Role of GPX4 in ferroptosis and its pharmacological implication. Free Radic Biol Med. 2019; 133:144–52. https://doi.org/10.1016/j.freeradbiomed.2018.09.014 [PubMed]
  • 49. Gaschler MM, Andia AA, Liu H, Csuka JM, Hurlocker B, Vaiana CA, Heindel DW, Zuckerman DS, Bos PH, Reznik E, Ye LF, Tyurina YY, Lin AJ, et al. FINO2 initiates ferroptosis through GPX4 inactivation and iron oxidation. Nat Chem Biol. 2018; 14:507–15. https://doi.org/10.1038/s41589-018-0031-6 [PubMed]
  • 50. Lynes JP, Nwankwo AK, Sur HP, Sanchez VE, Sarpong KA, Ariyo OI, Dominah GA, Nduom EK. Biomarkers for immunotherapy for treatment of glioblastoma. J Immunother Cancer. 2020; 8:e000348. https://doi.org/10.1136/jitc-2019-000348 [PubMed]
  • 51. McGranahan T, Therkelsen KE, Ahmad S, Nagpal S. Current State of Immunotherapy for Treatment of Glioblastoma. Curr Treat Options Oncol. 2019; 20:24. https://doi.org/10.1007/s11864-019-0619-4 [PubMed]
  • 52. Li L, Zhu X, Qian Y, Yuan X, Ding Y, Hu D, He X, Wu Y. Chimeric Antigen Receptor T-Cell Therapy in Glioblastoma: Current and Future. Front Immunol. 2020; 11:594271. https://doi.org/10.3389/fimmu.2020.594271 [PubMed]
  • 53. Zhao J, Chen AX, Gartrell RD, Silverman AM, Aparicio L, Chu T, Bordbar D, Shan D, Samanamud J, Mahajan A, Filip I, Orenbuch R, Goetz M, et al. Immune and genomic correlates of response to anti-PD-1 immunotherapy in glioblastoma. Nat Med. 2019; 25:462–69. https://doi.org/10.1038/s41591-019-0349-y [PubMed]