A Novel Ferroptosis-related Gene Signature for Prognosis Prediction in Glioma

Background: Glioma is a fatal brain tumor characterized by invasive nature, rapidly proliferation and tumor recurrence. Despite aggressive surgical resection followed by concurrent radiotherapy and chemotherapy, the overall survival (OS) of Glioma patients remains poor. Ferroptosis is a unique modality to regulate programmed cell death and associated with multiple steps of tumorigenesis of a variety of tumors. Methods: In this study, ferroptosis-related genes model was identied by differential analysis and Cox regression analysis. GO, KEGG and GSVA analysis were used to detect the potential biological functions and signaling pathway. The inltration of immune cells was quantied by Cibersort. Results: The patients’ samples are stratied into two risk groups based on 4-gene signature. High-risk group has poorer overall survival. The results of functional analysis indicated that the extracellular matrix-related biologic functions and pathways were enriched in high-risk group, and that the inltration of immunocytes is different in two groups. Conclusion: In summary, a novel ferroptosis-related gene signature can be used for prognostic prediction in glioma. The ltered genes related to ferroptosis in clinical could be a potential extra method to assess glioma patients’ prognosis and therapeutic.


Background
Glioma is the most common primary malignant brain tumor, representing 81% of all malignant brain tumors [1]. According to World Health Organization, glioma is graded into 4 levels (I ~ IV). The 10-year survival rate of patients with low-grade glioma (WHO grade II ~ III glioma) is 47% with a median survival time of 11.6 years [2]. Glioblastoma (WHO grade IV glioma) patients who received post-operative radiation and concurrent temozolomide, have a median survival of 14.6 months [3]. Since patients with gliomas have such a poor prognosis, and currently the predictive prognosis ability of clinical testing indicators is weak and de cient, additional predictive models as novel biomarkers are needed to ll the vacancy. 2007 classi cation relied on histological analysis while 2016 classi cation incorporated molecular markers [4,5]. Some molecular markers of glioma: deletion of short arm Ch.1 and long arm of Ch. 19 (1p/19q)[6], methylation of O6methylguanine-DNA methyltransferase (MGMT) [7], mutation of isocitrate dehydrogenase 1 (IDH1)[8] were involved, but numerous clinical trails of glioma therapies that targeting these biomakers failed to demonstrate a robust therapeutic effect. Therefore, additional indicators need to be supplemented to comprehensively predict the prognosis of glioma patients.
Ferroptosis, characterized by iron-dependent lipid peroxidation, is a newly non-apoptotic form of cell death that is widely implicated in human pathological conditions [9]. It has been reported that ferroptosis plays a vital role in the pathogenesis of a growing list of disease, such as cancer, neurological disorders, including neurodegenerative diseases and brain damage, as a way to promote cell death [10][11][12][13][14], liver and lung brosis [15,16] and autoimmune disease [17,18]. Furthermore, some ferroptotic regulatory genes, such as GPX4 [19], ACSL4 [20], FANCD2 [21], SLC7A11 [22], have been reported to be closely related to tumor proliferation and migration. As such, the development of ferroptosis inductionbased cancer therapies is being actively explored, such as hepatocellular carcinoma [23], pancreatic cancer [24], gastric cancer [25] and non-small cell lung carcinoma [26]. Thus, modulating ferroptosis activities could hold promise for therapeutic development against malignant diseases.
Given this, our study aimed to screen out genes associated ferroptosis with glioma from multiple gene expression datasets, to explore novel biomakers that can help to supplementally and comprehensively predict the prognosis of glioma. Furthermore, we also explore ferroptosis-related genes and potential molecular mechanism of prognostic differences in glioma. Our study suggest that ferroptosis-related genes may be potential prognostic markers and therapeutic targets for glioma.

Datasets and experimental process
All datasets we used are available to public. We selected TCGA dataset (Genotype-Tissue Expression (GTEx) dataset as control dataset) as a derivation cohort. CGGA , Gravendeel and REMBRANDT datasets were chosen as validation cohorts.
The RNA-seq data, clinical information and molecular information were downloaded from UCSC-XENA (http://xena.ucsc.edu), CGGA (http://www.cgga.org.cn) and GlioVis (http://gliov is.bioin fo.cnio. es/). 687 patients from TCGA (5 normal samples and 2 samples without clinical messages were excluded), 105 controls from GTEx, 693 patients from CGGA, 276 patients from Gravendeel (8 normal samples were excluded) and 383 patients from REMBRANDT (92 sample with normal or missing histological classi cation were excluded) were nally enrolled in this study. The details of clinical and molecular information of those samples are shown in Table 1. The process of the study is shown in Figure 1.
The construction and validation of gene-associated prognostic model Univariate and multivariate Cox regression analyses were applied to analyze the relationships between patient' overall survival (OS) and the expression level of each ferroptosis-related genes. Under threshold of p<0.01 in the univariate Cox regression analysis, the gene was selected signi cant. We nally selected 4 genes by using the multivariate Cox regression to assess the contribution of gene as independent prognostic factor for patient survival and created a model. A 4-gene signature model of risk score was established based on a linear combination of the regression coe cient derived from the multivariate Cox regression model (coe cient) multiplied by its expression level. The prognosis index (PI) = (coe cient*expression level of FANCD2) + (coe cient*expression level of IGFBP3) + (coe cient*expression level of LPCAT3) + (coe cient*expression level of HMOX1). Cut-off value is the limit of the maximum value obtained by the Yoden index from receiver operating characteristic (ROC) curve by "survivalROC" [30] package. According to the cut-off value, the patients were grouped under low-risk and high-risk. The "survminer" (https://cran.rproject.org/web/packages/survminer/survminer.pdf) package in R was utilized to conduct Kaplan-Meier (KM) survival curves for the samples with low-and high-risk group. Univariate and multivariate Cox regression analyses were employed to certify whether the predictive power of prognostic model could be independent of other clinical variables (containing Age, Gender, tumor Grade, IDH1 Status and histology) for glioma patients. The hazard ratio (HR) with 95% con dence intervals (CI) were calculated.

Functional Enrichment Analysis
Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) are used to detect the potential biological functions and signaling pathway by the "clusterPro ler" [31] package. GeneMANIA (http://www.genemania.org) is used to conduct the gene-gene interaction network and potential pathways based on the DEGs between high-risk and low-risk group. Gene Set Variation Analysis (GSVA) [32] is used to estimate variation of gene set enrichment in four datasets, thereby allowing the evaluation of pathway enrichment for each sample.

Immune in ltration
To explore the immune in ltration in high risk and low risk group related to ferroptosis in glioma, "CIBERSORT" [33] package was used to assess the proportions of 22 immune cell subtypes based on expression le. The perm was set at 1000. Samples with P < 0.05 in CIBERSORT analysis result were used in further analysis. Mann-Whitney U test was used to compare differences in immune cell subtypes in the high-risk score and low-risk score groups.
Statistical Analysis R software (version 4.0.2) and SPSS software (version 23.0) were used to complete all the statistical analysis. The logrank test, and the differences between categorical variables were compared by using the Chi-square test. Independent prognostic factors were assessed by univariate and multivariate Cox regression. GO and KEGG analyses were calculated by Mann-Whitney U test with P values adjusted by the BH method. P < 0.05 was considered statistically signi cant.

Prognostic model construction in TCGA cohort
The risk score was calculated according to the expression of these ferroptosis-related genes in TCGA dataset, as formula: (0.504*expression level of FANCD2 + 0.264*expression level of HMOX1 + 0.235*expression level of IGFBP3 + 0.267*expression level of LPCAT3) (Supplementary table 2). Patients were ranked by the score, and they were strati ed into a low-risk group (n=436) and a high-risk group (n=251) based on the optimal cut-off expression value. Patients in the low-risk group had prolonged OS than high-risk group ( Figures 3A-C). The predictive performance of the risk model for OS was estimated by time-dependent ROC curves. The results showed that the area under the curve (AUC) reached 0.843 at 1 year, 0.846 at 2 years, and 0.849 at 3 years respectively ( Figure 3D).

External validation of 4-gene prognosis model
From the analysis above, the model was able to effectively predict glioma patients' prognosis in TCGA dataset. Subsequently, CGGA, Gravendeel and REMBRANDT datasets were performed to verify the e cacy and availability of the model. The patients in CGGA, Gravendeel and REMBRANDT datasets were also grouped into low-risk and high-risk respectively. The results in three datasets showed that high-risk group had a markedly poor prognosis than low-risk group Clinical characteristics of the patients based on risk score The heatmaps ( Figure 4A) showing the differential expression of the 4 selected genes, and baseline characteristics of the patients in different risk groups (Table 3) suggested that clinical and molecular features, such as WHO grade, age, classical subtypes, mesenchymal subtypes, and IDH wild types were enriched in high-risk group. Furthermore, with an increase in glioma grade, the risk score increased. The highest increase in risk score was found in the WHO grade IV patients, whereas the lowest increase in risk score was observed in the WHO grade II patients in the TCGA, CGGA, Gravendeel and REMBRANDT datasets ( Figure 4B). Meanwhile, Kaplan-Meier plots of overall survival based on different clinical characters presented apparent differences between the high-risk group and low risk group for survival ( Figures 5A-H). The survival-predictor scores from these models were highly predictive of survival in datasets (p < 0.001). The 4-gene mod was able to distinguish patients with obviously distinct outcomes across subsets, demonstrating the risk score based on four gene was greater prognostic power as compared with only used of clinical subgroups (WHO grade, IDH1, MGMT, 1p/19q).

Independent Prognostic Value of the 4-gene Signature
The independence of the clinical prognostic signi cance of the signature in glioma was veri ed. The risk score was signi cantly associated with OS in univariate Cox regression analysis in TCGA, CGGA, Gravendeel and REMBRANDT datasets, and the Hazard ratio (HR) of four datasets were: 9.57, 3.64, 2.80, 3.18 (95%CI, Under threshold of p<0.001) respectively (Table 4) The Clinical Predictive Performance of the Nomogram Based on the four ferroptosis-related gene signature (risk score) and clinical factor, a nomogram to predict patients' prognosis in the TCGA database was created ( Figure 7A). In the study, we used a bootstrap method to verify the developed nomogram with the C-index of 0.869 95%CI), which suggested that the predictive model had good predictive performance. Furthermore, the calibration curves in 1-, 2-, 3-year survival of patients also showed good consistency compared with the ideal model, further indicating that the nomogram was stable in predicting the prognosis of glioma patients ( Figures 7B-D). This suggests that basing therapy strategy on our nomogram will improve clinical outcome ( Figure 7E).

Functional Annotation of the 4-gene Signature
To clarify the potential functional characteristics associated with risk score, the DEGs between the low-risk and high-risk groups were used to perform GO enrichment and KEGG pathway analyses. The results of GO enrichment were concentrated in three areas: binding with skeleton-associated, changing of receptor channel and extracellular matrix constituent, extracellular matrix could be a pro-ferroptosis stress, which may mediate cell adhesion, vesicular bodies and exosomes resulted in poor prognosis for high-risk groups. The results indicated that the biological processes associated with extracellular matrix and cell adhesion were enriched in high-risk group (P. adjust < 0.05, Figure 8A). KEGG pathway analysis also showed that the ECM−receptor interaction pathway and the Focal adhesion pathway were signi cantly enriched in high-risk group (P. adjust < 0.05, Figure 8B). Meanwhile, the results of GeneMANIA also con rmed that the high-risk group was obviously gathered in the structure, organization and disassembly of extracellular matrix ( Figure  8C). These enriched biologic functions thus provided a deeper understanding of the patients' poor prognosis.

Relationship between ferroptosis and Immune in ltration
Analysis of hallmark pathway gene signatures indicated that signaling pathways gathering at various biological progress were signi cantly different between high-and low-risk patients ( Figure 9A reported that due to the multifaceted nature and the involvement of numerous metabolic pathways directly impinging on ferroptosis, it is highly likely that certain key nodes of the ferroptosis process are targeted by the immune system [43]. To further explore the correlation between the model and immune status, enrichment scores of the diverse immune cell subpopulations in TCGA and CGGA were quanti ed by Cibersort. Various immune cells (CD4+T cell, CD8+T cell, Tregs, monocytes, macrophages M0/M2, NK cells) were obviously different between the high risk and low risk ( Figure 9G-H).
The immune in ltrations of the two databases perfectly shows that the four ferroptosis-gene model are closely relevant to the immune.

Discussion
Glioma is the main malignant tumor in primary brain tumors. At present, the prognostic prognosis signature in brain tumors include IDH1, MGMT, etc. subtype classi cation is helpful to better predict the survival time of subgroups based on the biological detection indexes [44]. The combination of IDH1, MGMT and other biological indicators helped to easily achieve this expectation. Otherwise, the molecular pathway mechanism of glioma prognosis in other elds needed to be expanded. Ferroptosis is a regulatory form of non-apoptotic cell death, which is characterized by the accumulation of iron dependent lethal lipid reactive oxygen species (ROS). There had been a break-through in tumor and ferroptosis. Some studies had shown that ferroptosis and glioma was closely correlative, but there were few data analysis studies to clearly explain the relationship between iron death and glioma prognosis. A prognosis prediction model was created by using four public databases based on ferroptosis-related genes to further explore the prognosis correlated ferroptosis to glioma as well as the potential biological relevance.
It had been reported that the ferroptosis-related genes had diverse functions in uencing tumor's progress, such as migration, invasion, angiogenesis and immune in ltration. Based on GO and KEGG function analysis, cell adhesion, extracellular matrix, exosomes were obviously enriched, which was correlated with immune microenvironment in glioma reported by previous studies. gliomas were characterized by profound immunosuppression mediated by secreted (TGF-β, IL-10) and cell surface (CD95L, PD-L1) immunosuppressive factors, by in ltration of immune inhibitory cells in the tumor microenvironment and by anatomical peculiarities of the brain [45]. Immune in ltration in our study also presented that ferroptosis prognosis genes were closely related to the immune. Interestingly, previous studies showed that HMOX1 and IGFBP-3 were signi cantly associated with immune. The enzyme heme oxygenase-1 (HOMX1) was part of an endogenous defense system implicated in the homeostatic response, data were available from a wide spectrum of physiopathologic conditions that link HOMX1 to modulation of angiogenesis and the immune function. The immunodulatory function of HOMX1 is well documented and its immunosuppressive role is widely accepted [46]. It had reported that Hemin treatment resulted in enhanced CTL effector function both in normal and tumor microenvironments.
Insulin-like growth factor 1 (IGF-I) and IGF-II were members of the insulin superfamily of growth-promoting peptides and were among the most abundant and ubiquitous polypeptide growth factors [47]. The IGFBPs were secreted proteins that were also found intracellularly and that interact with many ligands other than the IGFs. IGFBP-3 was determined to activate the type 1 TGFβ receptor and its downstream effectors SMAD2 and SMAD3 to inhibit breast cancer cell growth [48]. Phospholipid-remodeling enzyme, lysophosphatidylcholine acyltransferase 3 (Lpcat3), as a critical determinant of membrane phospholipid composition, was important drivers of ferroptosis. Loss of Lpcat3 or activation of SREBP-2 in Apcmin mice markedly promote intestinal tumor formation [49]. In other word, Lipid metabolism marked as LPCAT3 might be potential correlated to immune as HMOX1 and IGFBP3.but it needed to be veri ed by further studies.
Moreover, IGFBP-3 had been acting as a marker of prognosis in various tumor. Four genes prognosis model with FANCD2, HMOX1, LPCAT3 and IGFBP3 were identi ed in 4 datasets. We attested that the 4-gene signature is a clinical promising biomarker, which could classify glioma patients into subgroups. Our function analysis presented that ferroptosis was closely associated with the immune and further predicted the existence of potential relationship between the lipid metabolism and the immune in glioma. Previous studies showed that blocking key enzymes, such as fatty acid synthase (FASN), which is required for the synthesis of palmitate from acetyl-CoA and malonyl-CoA48, and ELOVL fatty acid elongase 2 (ELOVL2), which catalyses the elongation of fatty acids, suppressed glioma tumor growth in GBM xenograft models, indicating that fatty acids are required for glioma growth and survival [50,51]. But it's needed to further understand the mechanisms of ferroptosis as well as its effect to glioma.

Conclusions
In conclusion, a novel ferroptosis-related gene signature, which could classify glioma patients into subgroups, can be used for prognostic prediction in glioma. The 4-gene model in clinical could be a potential extra method to assess glioma patients' prognosis and therapeutic.

Availability of data and materials
The data sets used and/or analyzed during the current study are available from the corresponding author on reasonable request.
Ethics approval and consent to participate Not applicable

Consent for publication
We assure that the material is original and it has not been published elsewhere yet.  Table 1 The clinical and molecular characteristic of samples in GTEx, TCGA, CGGA, Gravendeel and REMBRANDT  Table 3 Baseline characteristics of the patients in different risk groups    Identi cation of 4 ferroptosis-related candidate genes in TCGA cohort. A. The volcano gure shows the differential expression of 23 ferroptosis-related genes in 4352 DEGs between normal and tumor tissues. B. Venn diagram to identify differentially expressed ferroptosis-related genes that were correlated with OS. C. Forest plot shows the coe cient of the 6 genes. TP53 and SLC40A1 were not signi cant. D. The correlation plot shows the co-linearity between the 4 genes. E.
The heatmap shows the expression of 4 genes in normal and glioma tissues.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.