Ferroptosis-Related Gene Signature Predicts Glioma Cell Death and Glioma Patient Progression

Glioma is a fatal brain tumor characterized by rapid proliferation and treatment resistance. Ferroptosis is a newly discovered programmed cell death and plays a crucial role in the occurrence and progression of tumors. In this study, we identified ferroptosis specific markers to reveal the relationship between ferroptosis-related genes and glioma by analyzing whole transcriptome data from Chinese Glioma Genome Atlas, The Cancer Genome Atlas dataset, GSE16011 dataset, and the Repository of Molecular Brain Neoplasia Data dataset. Nineteen ferroptosis-related genes with clinical and pathological features of glioma were identified as highly correlated. Functional assays in glioma cell lines indicated the association of ferroptosis with temozolomide resistance, autophagy, and glioma cell migration. Therefore, the identified ferroptosis-related genes were significantly correlated with glioma progression.


INTRODUCTION
Glioma is the most common primary malignant tumor in the central nervous system (Rasmussen et al., 2017). It can be graded from I to IV in accordance with the 2016 World Health Organization (WHO) classification (Louis et al., 2016). Owing to high proliferative rate, heterogeneity of tumor cells and diffuse infiltrating property, high-grade glioma is difficult to completely removed by surgery (Lara-Velazquez et al., 2017;Manrique-Guzman et al., 2017;Ferguson and McCutcheon, 2018). Moreover, chemotherapy resistance often occurs and leads to treatment failure and tumor recurrence (Osuka and Van Meir, 2017). Glioblastoma (WHO grade IV glioma, GBM) has a median overall survival (OS) of only 14.6 months (Lara-Velazquez et al., 2017). Molecular markers, such as O6-methylguanine-DNA methyltransferase (MGMT) methylation, codeletion of the short arm of chromosome 1 and the long arm of chromosome 19 (1p/19q), mutations in isocitrate dehydrogenase (IDH), mutations in alpha thalassemia/mental retardation syndrome X-linked, and epidermal growth factor receptor are being used in molecular pathological diagnosis, treatment options, and prognostic evaluation with glioma patients (Zeng et al., 2015). These molecular markers play a central role in regulating tumor cell proliferation and death (Ceccarelli et al., 2016;Stupp et al., 2019). Numerous glioma therapies targeting these molecular markers are used in clinical trials, but few have ultimately succeeded. Thus, novel targets for glioma therapy need to be identified.
In the present study, we analyzed the differential expression of ferroptosis-related genes in glioma samples to identify the enriched pathways and their biological functions. We determined that ferroptosis-related genes were associated with the prognosis of glioma. Using the 19 identified ferroptosisrelated genes, we accurately predicted the outcome for glioma patients. The underlying mechanisms were ultimately confirmed in in vitro studies. Overall, our data suggest that ferroptosis-related genes play pivotal roles in glioma progression and are potential prognostic markers and therapeutic targets for glioma.

Datasets
All datasets used in this study were available to the public. Expression RNA-seq data and clinical and molecular information of patients were obtained from the Chinese Glioma Genome Atlas (CGGA) dataset 1 as a training set. The Cancer Genome Atlas (TCGA) dataset, the GSE16011 dataset, the Repository of Molecular Brain Neoplasia Data (REMBRANDT) dataset, and associated clinical information were obtained from the TCGA official website 2 , GEO website 3 , and GlioVis website 4 as validation sets. The characteristics of the glioma patients in these 4 datasets-OS, age, gender, WHO grade, TCGA subtypes, and molecular pathological features-have been recorded in previous studies (Gravendeel et al., 2009;Bao et al., 2014;Tomczak et al., 2015;Gusev et al., 2018) (Table 1).

RNA Sequencing of Glioma Cell Lines
Total RNA from U251MG, U87MG, U251TR, and U87TR cell lines were extracted and quantified. Up to 3 µg RNA per sample was used as the input material for library preparation. The library preparation was conducted using the standard Illumina HiSeq platform, and 150 bp paired-end reads from transcriptome sequencing were generated. After quality control procedures, the data were mapped to the human hg19 reference by STAR (Dobin et al., 2013). The expression level (Fragments per Kilobase Million, FPKM) of each gene was calculated based on the length of the gene and the read count mapped to this gene The RNA seq dataset was submitted to https://figshare.com/s/ 7af369b94634a51b11f1 for reference.

Cell Viability Assay
Cell Counting Kit-8 Assay Kits (CCK-8; Dojindo, Kumamoto, Japan) were used to determine cell viability. The cells were seeded at a density of 2 × 103 cells/well in 96-well plates and were incubated with serum-containing media for 24 h. The cells were treated with erastin at varying concentrations (0, 10, 20, 30, 40, and 50 µM) on U251MG, U251TR, U87MG, and U87TR cells for 24 h. The medium was replaced with 100 µL of the fresh medium, and 10 µL of the CCK-8 solution was added to each well. Blank controls without cells were prepared. Subsequently, 100 µL of each sample was transferred to a new 96-well plate and was analyzed with the plate reader Infinite 200 PRO (Tecan Life Sciences). The absorbance values were determined at 450 nm.

Western Blot Assay
The procedures for Western blot were described in detail in a previous study (Liu H. J. et al., 2018). In the current study, after incubation with primary antibodies, the HRP labeled anti-rabbit or anti-mouse secondary antibody was used for incubation at room temperature for another period of 1 h. Specific protein bands were detected using an ECL Western blotting kit, following the recommended procedure. β-actin was used as an internal control for sample loading and standardization.

Migration Assay
The Transwell system (24 wells, 8 µm pore size with a polycarbonate membrane) was used for in vitro migration assays. The U251MG, U251TR, U87MG, and U87TR cells were pretreated with erastin (50 µM) or without erastin. A total of 1 × 10 5 cells were suspended in 100 µL serum-free medium and then added to the upper chamber. The medium containing 10% FBS and erastin at different concentrations was added to the lower chamber. The cells in the upper chamber were carefully wiped using a cotton swab after 2 h for the U87MG and U87TR cells or 7 h for the U251MG and U251TR cells.
The cells that were attached to the lower surface of the filter were fixed with 4% paraformaldehyde and then stained with crystal violet. The migrated cells on the lower surface of the membrane filter were photographed using an Axio Observer3 microscope (Carl Zeiss), and 5 randomized fields were counted using the Image J software. All experiments were repeated three times.

Gene Signature Building
We chose 40 genes, which were verified to be involved in ferroptosis (Stockwell et al., 2017). Univariate Cox analysis was first performed, and genes with P values less than or equal to 0.0001 were retained. We then developed a risk score model for each patient on the basis of the expression level of the 19 genes, and their regression coefficients were derived from the univariate Cox regression analysis. The risk score was expressed as (expr gene1 × coefficient gene1 ) + (expr gene2 × coefficient gene2 ) + · · · + (expr gene19 × coefficient gene19 ). The regression coefficients derived from the training set were then applied into the three other validation sets (TCGA, GSE16011, and REMBRANDT datasets) to calculate the risk scores.

Statistical Analysis
The receiver operating characteristic (ROC) curve and a nomogram were used to show the predictive accuracy of the gene signature with the R package "survival ROC" (Heagerty et al., 2000) and "rms" (Feng et al., 2017). The area under the curve (AUC) of the ROC curves, concordance index (C-index), and calibration curves of nomograms were calculated. Gene ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were applied in DAVID 5 (Dennis et al., 2003) for functional annotation of the genes positively correlated with risk score in the 4 cohorts. Gene Set Variation Analysis (GSVA) was conducted to detect the difference in expression with risk score by using the R package GSVA (Hanzelmann et al., 2013). We classified the patients in each dataset into two groups (high-risk-score and lowrisk-score groups) in accordance with the median risk score of the ferroptosis-related gene signature. Kaplan-Meier curves and the log-rank test were used to assess the differences in OS between two groups. Independent prognostic factors were identified by univariate and multivariate Cox regression analysis. All statistical analyses were conducted using SPSS, GraphPad Prism 7, or the R software. P < 0.05 was considered statistically significant.

Identification of 19 Ferroptosis-Related Genes
To characterize the ferroptosis-related gene expression pattern in gliomas, we examined the RNA-seq data of glioma patients from 4 datasets (CGGA, TCGA, REMBRANDT, and GSE16011). First, univariate Cox regression analysis was performed to identify the genes related to patient survival in the CGGA dataset. Subsequently, 19 of the 40 ferroptosis-related genes were selected in glioma patients (SAT1, ATP5G3, FANCD2, HSPB1, HMGCR, CBS, GCLC, GCLM, CD44, ALOX12B, ALOX5AP, CISD1, NFE2L2, EMC2, ALOX5, DPP4, AKR1C2, LPCAT3, and NCOA4, P < 0.0001). A gene-based prognostic model was then established to evaluate the risk of each patient as described in the methods. Details on the 19 genes are presented in Table 2. The risk scores of each patient in these four datasets (CGGA, TCGA, REMBRANDT, and GSE16011) were calculated. Heat maps were shown to present the different expression levels of the 19 selected genes and clinical information ordered by risk score in the CGGA ( Figure 1A), TCGA ( Figure 1B), GSE16011 ( Figure 1C), and REMBRANDT ( Figure 1D) datasets. In the CGGA cohort, with an increase in risk score, the expression levels of EMC2, AKR1C2, ALOX12B, CISD1, CBS, HMGCR, NCOA4, and GCLC were distinctly downregulated. Meanwhile, the expression levels of FANCD2, LPCAT3, ATP5G3, HSPB1, ALOX5, ALOX5AP, CD44, GCLM, SAT1, DPP4, and NFE2L2 were upregulated. Clinical and molecular features, such as WHO grade, age, classical subtypes, mesenchymal subtypes, and IDH wild types were enriched in high-risk-score gliomas. The other three cohorts showed similar patterns of genes and clinical information. Notably, the MGMT gene with promoter unmethylation was enriched in high-risk-score glioma in the TCGA dataset. These results suggest that the risk score of signatures of the 19 ferroptosis-related genes positively correlated with glioma malignancy, and the expression levels of the 19 genes in the gliomas exhibit a pattern similar to those of the other three datasets.

Association of Risk Scores With Patient Prognosis and Glioma Grade
The distribution of risk scores in the CGGA, TCGA, GSE16011, and REMBRANDT datasets is presented (Figures 2A-D). Glioma patients were divided into low-risk and high-riskscore groups on the basis of their median risk scores. The OS of each patient was shown in the CGGA, TCGA, GSE16011, and REMBRANDT (Figures 2E-H) datasets.
The patients with a low-risk-score had a markedly lower mortality rate than those with a high-risk-score in these four datasets. Meanwhile, with an increase in glioma grade, the risk score increased. The highest increase in risk score

Validation of 19 Ferroptosis Related Gene Signature in Survival Using Kaplan-Meier Curves
Patients with different kinds of glioma, such as pan-glioma and WHO grade II-IV glioma, were divided into two groups based on their median risk scores. The Kaplan-Meier curve for the CGGA dataset showed that the high-risk patients had significantly shorter OS than low-risk patients in the glioma (Figure 3A), WHO grade II glioma (Figure 3B), WHO grade III (Figure 3C), and WHO grade IV glioma ( Figure 3D) groups. Consistency of results was validated for the TCGA (Figures 3E-H (Figures 4A-D). A 5-year survival nomogram prediction model was then built with independent prognostic parameters for the OS of patients in the CGGA dataset ( Figure 4E). The C-indices were 0.79 in the primary CGGA dataset, 0.849 in the TCGA dataset, 0.75 in the GSE16011 dataset, and 0.653 in the REMBRANDT dataset as validation. Meanwhile, the calibration plot for the probability of 5-year survival showed an optimal agreement between observation and prediction in these datasets ( Figure 4F). These results indicated that the signature of the 19 ferroptosis-related genes was a reliable prognostic indicator in glioma patients.

Functional Annotation of the Signature of 19 Ferroptosis-Related Genes
To clarify the potentially functional characteristics of the signature of 19 ferroptosis-related genes in glioma, GO analysis and KEGG analysis were conducted. We first created a list of genes that were positively correlated with risk score, as revealed by Pearson correlation analysis in the CGGA (Pearson R > 0.6, p < 0.05), TCGA (Pearson R > 0.5, p < 0.05), GSE16011 (Pearson R > 0.6, p < 0.05), and REMBRANDT (Pearson R > 0.6, p < 0.05) datasets (Supplementary Table S1). We then explored the biofunctions of the genes in each dataset by GO analysis and KEGG analysis using DAVID Bioinformatics Resources 6.8. Genes with the p value of gene functions > 0.05 were excluded. The GO terms of the remaining genes in the CGGA  Table S5) datasets were listed. These gene functions were excluded from the intersection of these four datasets, as reflected in the Venn diagram of the resulting GO terms in these datasets ( Figure 5A). Ten of the GO terms were included in all datasets ( Figure 5B): inflammatory response, innate immune response, antigen processing and presentation of peptide antigen via MHC class I, cellular defense response, antigen processing and presentation, immune response, response to wounding, cell death, and actin cytoskeleton organization and protein processing. These functions play a key role in tumorigenesis and progression. KEGG analysis was performed based on the aforementioned gene list. Gene pathways with a p-value > 0.05 were excluded. The remaining gene pathways in the CGGA (Supplementary  Table S9) datasets are listed. The results of the KEGG pathway analysis in these datasets are reflected in the Venn diagram in Figure 5C. Seven cancer-related terms in KEGG were included in at least three datasets (Figure 5D), which was consistent with the GO results. The terms were hsa04670: leukocyte transendothelial migration, hsa04210: apoptosis, hsa04610: complement and coagulation cascades, hsa04142: lysosome, hsa05416: viral myocarditis, hsa04510: focal adhesion, and hsa04810: regulation of actin cytoskeleton. Subsequently, 14 gene sets related to the GO terms were applied to GSVA analyses for validation. Consistent with the GO results, heat maps showed that the risk score was positively associated with these gene sets in the CGGA, TCGA, GSE16011, and REMBRANDT datasets ( Figure 5E). All aforementioned results suggest that the signature of the 19 ferroptosis-related genes is related to cancer biology, particularly cell death, migration, and immune-related functions.

Association of Ferroptosis With Glioma Drug Resistance
To determine the relationship between ferroptosis and glioma drug resistance, we periodically treated glioma cell lines U87 and U251 with TMZ and constructed their corresponding drug-resistant strains U87TR and U251TR, respectively. The mRNA expression profiling technique was used to detect the expression levels of ferroptosis-related genes in different glioma cell lines. Changes in the expression pattern of the 19 ferroptosisrelated genes were consistent (Figures 6A,B). Differences in the risk score of the characteristics of the ferroptosis-related genes were consistent as well. Meanwhile, the risk scores of the drug-resistant cells were lower than those of the normal tumor cells (Figures 6C,D). In the CGGA glioma dataset, the risk score of the signature of the 19 ferroptosis-related genes was negatively correlated to the expression level of the well-recognized gene MGMT for TMZ resistance (Figure 6E), verifying the results for the risk scores of resistant cells cultured in vitro. We used erastin to treat the normal and drug-resistant glioma cell lines at increasing doses and CCK-8 assay to evaluate cell proliferation activity. We found that the U87TR and U251TR cell lines were more sensitive to erastin, compared with the U87MG and U251MG cell lines (Figures 6F,G). These results indicate that compared with normal glioma cells, TMZ-resistant cells are more likely to induce ferroptosis. The lower their risk scores, the more likely the occurrence of ferroptosis.

Relationship Between Ferroptosis and Autophagy and Apoptosis
In the CGGA dataset, we analyzed the risk scores of the signature of the 19 ferroptosis-related genes and autophagy markers (SQSTM1, MAPLC3C, BECN1, and ATG9B), apoptosis markers (CASP3, BCL2, BAX, and BAK1), and immune checkpoint markers (PD-1, TIM3, LAG3, and B7-H3). We found that the highest risk scores were positively correlated with these markers, except BCL2 (Figures 7A-C). Moreover, we determined the relationship between erastin-induced ferroptosis and autophagy, apoptosis, and immune function in the GBM cell lines. Erastin at different concentrations (0-50 µM) was used to treat the U251 and U87 GBM cell lines. Erastin treatment increased the expression of LC3 (an autophagy marker) and PD-L1 (an immune checkpoint marker) ( Figure 7D). However, treatment with erastin failed to change the expression level of cleaved PARP, an apoptotic molecule. These data suggest that erastin-induced ferroptosis is closely related to autophagy and immune function in tumor cells but not to apoptosis.

Ferroptosis Is Positively Associated With Glioma Cell Migration
To determine the relationship between ferroptosis and glioma cell migration, we analyzed the relationship between the risk scores of the signature of 19 ferroptosis-related genes and migration-related genes. The risk score was positively correlated to the oncogenes S100A4, TWIST1, CDH2, and POSTN, which were critically involved in glioma migration and invasion ( Figure 8A). We then determined the migration ability of glioma cell lines (U251 and U87) treated with a vehicle or erastin by using Transwell migration assay. Compared with the vehicle-treated U251 and U87 cells, the erastin-treated U251 and U87 cells increased in migratory ability (Figures 8B-D). These data suggest that ferroptosis enhances cell migration in glioma cells.

DISCUSSION
Selective induction of cancer cell death is the most effective anticancer therapy. Increasing evidence has shown that ferroptosis, a recently discovered PCD, plays a crucial role in tumorigenesis and cancer therapeutics. However, profiling of ferroptosis in glioma has yet to be clarified. In this study, we used high-throughput expression analysis to investigate variations in expression profiling of ferroptosis-related genes in glioma.
On the basis of these analyses, we identified the signature of the 19 ferroptosis-related genes. We further explored tumor clinicopathological features and prognosis, which were closely related to tumorigenesis. Previous studies suggest that ferroptotic cell death results from fatal lipid peroxidation . In this regard, the accumulation of intracellular iron caused by the depletion of ferritin or iron transporters and subsequent peroxidation are fundamental mechanisms that lead to the accumulation of lipid peroxides and ferroptosis (Stockwell and Jiang, 2019). Among the 19 identified genes related to ferroptosis, CISD1 alleviated iron accumulation in the mitochondria, and NCOA4 promoted intracellular iron accumulation, whereas HSPB1 decreased intracellular iron accumulation. Polyunsaturated fatty acyl tail-containing phospholipids, the substrate of lipid peroxidation, is considered essential for ferroptosis. Lipid peroxidation is negatively regulated by AKRIC2 and NFE2L2 but positively regulated by DPP4, ALOXs, HMGCR, LPCAT3,  and SAT1. Malfunctioning of scavenging lipid peroxide also leads to ferroptosis. After being transported intracellularly by the system XC-, cystine may become cysteine, which is used to synthesize glutathione. GPX4 and GSH scavenge lipid peroxide to inhibit ferroptosis (Lewerenz et al., 2006). In the signature genes, CD44, CBS, and GCLC were positive regulators, whereas GCLM, ATP5G3, and FNDC2 were negative regulators for synthesizing GPX4. Thus, the coexistence of the aforementioned factors may trigger ferroptosis (Figure 9). In the current study, we found that higher risk scores were associated with worse prognosis in glioma patients. Moreover, high risk scores were insensitive to ferroptosis in glioma cell lines. To verify this finding, cell death pathways were shown to be frequently disturbed in human cancers (Bebber et al., 2020).
By focusing on the specific function of the 19 ferroptotic genes, previous studies have demonstrated that most of these genes played a pivotal role in cancer cells, including glioma. For instance, high expression of HSPB1 is correlated with a low survival rate in glioma patients. HSPB1 also enhances proliferation via SIRT2-mediated G6PD activation in glioma cells (Ye et al., 2016). The increased expression of SAT1 in GBM was related to the resistance of tumor cells to radiotherapy (Brett-Morris et al., 2014). During lung carcinogenesis, NFE2L2 accelerates progression via the KRAS signaling pathway (Satoh et al., 2013). Decreased expression of CBS promotes the formation of glioma tumors by increasing the HIF-2α protein levels and HIF-2 target gene expression (Takano et al., 2014). These genes are either positively or negatively related to the regulation of ferroptosis and are distinctly expressed in the highrisk-score or low-risk-score groups. This occurrence is consistent with their gene functions in cancers. Notably, not all 19 genes had expression levels consistent with their functions in previous reports. Busek et al. (2012) indicated that DPP4, which is highly expressed in gliomas with a high risk score, inhibits glioma cell growth independent of its enzymatic activity. Another study showed that DPP4 was related to the stemness of glioma stem cells (Sakamoto et al., 2019). Sohn et al. (2013) demonstrated that CISD1 that was highly expressed in human epithelial breast cancer cells suppressing CISD1 expression significantly inhibited cell proliferation and tumor growth. ARK1C2 was found to positively regulate proliferation in cancer cells (Ji et al., 2004). We found that ARK1C2 expression was low in high-risk-score groups. Therefore, the specific role of these genes in glioma has to be clarified.
We further demonstrated that the risk scores of the signature of the 19 ferroptosis-related genes were highly associated with the WHO tumor grades. The ROC curve generated using the risk scores of the signature of the 19 ferroptosis-related genes predicted patient OS. In addition, the signature of the 19 ferroptosis-related genes was independent of other clinical factors, including age, radiotherapy, grade, and chemotherapy. These results suggested that the activated process of ferroptosis in glioma cells were associated with improved survival of glioma patients.
Functional annotation of the signature of the 19 ferroptosisrelated genes showed that biological functions such as immune response, cellular defense response, actin cytoskeleton  organization, cell death, and protein processing could contribute to the poor clinical outcome of patients. KEGG analysis indicated that the signature of the 19 ferroptosis-related genes corresponding to biological functions were closely related to apoptosis, focal adhesion, regulation of actin cytoskeleton, leukocyte transendothelial migration, and lysosome pathways, which are strongly linked to tumorigenesis. Wang et al. (2019) recently reported that CD8-positive T cells induced ferroptosis in tumor cells. Chen et al. (2017) reported that ATF4, a gene that sensitizes tumor cells to ferroptosis, promotes glioma cell migration. Hong et al. (2017) suggested the involvement that the p53-independent CHOP/PUMA axis in response to ferroptosis inducers, which could play a key role in ferroptotic agent-mediated sensitization to TRAIL-induced apoptosis. Emerging evidence suggests that ferroptosis often shares common pathways with other types of biological functions, including cell death (Bock and Tait, 2019), immune response (Stockwell and Jiang, 2019), and migration (Chen et al., 2017). We found that the risk scores of the signature of the 19 ferroptosis-related genes was negatively correlated to the expression of the well-recognized gene MGMT for TMZ resistance, suggesting the association of ferroptosis with glioma drug resistance. We further found that in glioma cell lines, erastin-induced ferroptosis is closely related to autophagy but not apoptosis because erastin increased the autophagy marker LC3 level but not the apoptosis marker PARP levels.
Notably, erastin promoted cell migration in glioma cell lines. Increased cell migration likely occurred prior to the ferroptosis of the glioma cells.
In summary, this study is the first to investigate ferroptotic gene expression patterns in glioma patients and identify their relationship to patient outcome. The ferroptotic signature identified in our study exhibited potential as a biomarker of OS in glioma patients. Understanding the mechanisms underlying ferroptosis and its effect on OS, as well as its implications for the treatment of glioma, can provide insights into the identification of therapeutic targets for glioma.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/Supplementary Material.

ETHICS STATEMENT
This research protocol was approved by the ethical committee of the institutional review board of Capital Medical University and a consent form was signed by each participating patient.

AUTHOR CONTRIBUTIONS
HL, CZ, and TJ: conception and design, data analysis and interpretation (statistical analysis), and manuscript writing and revision. HL, HH, and GL: development of methodology. HL, YZ, FW, XL, and KW: data acquisition. All authors read and approved the final manuscript.