A Novel Survival Model Based on Ferroptosis-related Gene Signature for Overall Survival Prediction in Bladder Cancer

Background: The effective treatment and prognosis prediction of bladder cancer(BLCA) remains a medical problem. Ferroptosis is an iron-dependent form of programmed cell death. Ferroptosis are closely related to tumor occurrence and progression, but the prognostic value of ferroptosis-related genes (FRGs) in BLCA remains to be further clarified. In this study, we identified a FRGs signature with potential prognostic value for patients with BLCA. Methods: The corresponding clinical data and the mRNA expression profile of BLCA patients were downloaded from The Cancer Genome Atlas (TCGA). Univariate Cox regression was used to extract FRGs related to survival time, Cox regression model was applied to construct a multigene signature. Both principal component analysis (PCA) and single-sample gene set enrichment analysis (ssGSEA) were performed for functional annotation. Results: Clinical traits were combined with FRGs, so that 15 prognostic-related FRGs were identified by Cox regression. High expression of CISD1, GCLM, CRYAB, SLC7A11, TFRC, ACACA, ZEB1, SQLE, FADS2, ABCC1, G6PD and PGD are related to poor survival rates of BLCA patients. Multivariate Cox regression constructed a prognostic model with 7 FRGs and divided patients into two risk groups. Compared with the low-risk group, the overall survival(OS) of patients in the high-risk group was significantly lower (P <0.001). In multivariate regression analysis, the risk score was shown to be an independent predictor of OS (HR> 1, P <0.01). ROC curve analysis verified the predictive ability of the model. In addition, the two risk groups displayed different immune statuses in the ssGSEA and different distributed patterns in PCA. Conclusion: Our research suggests that a new gene model related to ferroptosis can be applied for the prognosis prediction of BLCA. Targeting FRGs may be a treatment option for BLCA.


Itroduction
Bladder cancer(BLCA) is a global problem and has been reported as the ninth most common tumor in the world [1,2]. More than 200,000 people died from the disease, and over 549,000 new cases were diagnosed in 2018 [1,2]. In the past 20 years, the number of BLCA incidents has been increasing worldwide. Due to the population aging, environmental pollution and smoking [3], the number of BLCA may increase in the future [4,5]. Although many treatments including surgery, chemotherapy, radiotherapy have been performed on BLCA patients, BLCA has a high risk of progression, metastasis, recurrence and poor prognosis [6,7]. Therefore, it is very important that we should explore and discover reliable prognostic biomarkers to guide clinical treatment and to improve the prognosis of BLCA, so as to improve the survival and prognosis of BLCA patients.
Ferroptosis is an iron-dependent form of programmed cell death, driven by the accumulation of lipid peroxides [8,9], which is different from autophagic cell death or traditional apoptosis or necrosis [8]. In recent years, promoting ferroptosis has become a promising treatment option that causes cancer cell death, especially for malignant tumors that are resistant to conventional treatment [10,11]. In addition to drugs that cause ferroptosis, many genes have also been identified as regulators or markers of ferroptosis. Ferroptosis regulatory genes such as GPX4 [12], P53 [13], DPP4 [14], HSPB1 [15] and FANCD2 [16] are closely related to tumor occurrence and progression. More and more literature reports that a variety of tumor cells include hepatocellular carcinoma cells [17], ovarian cancer cells [18] and adrenocortical carcinomas cells [19] are sensitive to ferroptosis. In addition, the combination of chemotherapeutic drugs and erastin, a ferroptosis inducer can improve its therapeutic effect on lung cancer cells [20], acute ovarian cancer cells [21], gastric cancer cells [22] and myeloid leukemia cells [23]. Therefore, ferroptosis may become a potential target for cancer treatment. However, the relationship between the prognosis of BLCA patients and the expression of ferroptosis-related genes (FRGs) has not been studied in detail.
In present study, the corresponding clinical data and the mRNA expression profile of BLCA patients was downloaded from the The Cancer Genome Atlas (TCGA) database. Then a prognostic multigene signature was constructed, based on 7 FRGs to predict survival in BLCA patients. Overall, our data indicate that FRGs play a key role in the pathogenesis of BLCA and are potential therapeutic targets and prognostic markers for BLCA.

Data collection
The public available RNA sequencing (RNA-seq) data for BLCA with corresponding clinical information of 433 patients contained BLCA sample(n=414) and normal bladder samples(n=19) were downloaded from TCGA website (https://gdc.cancer.gov/). We used "limma" R package to identify the gene expression profiles in TCGA BLCA dataset. The data from TCGA are available publicly. Therefore, this study does not need to be approved by the local ethics committee.

Construction a prognostic-related FRGs risk score model
FRGs expression were conbined with clinical features. Univariate Cox regression was used to extract FRGs related to survival time, using p<0.05 as the threshold. The STRING database (version 11.0) was used to generate an interactive network of prognostic-related FRGs [26]. Multivariate Cox regression was used to construct a multigene signature and calculate the risk score for each patient [27,28].

Independent prognostic analysis
The clinical characteristics and risk scores were compared with survival time and survival status to determine whether the risk score can be used as an independent prognostic factor. The receiver operating characteristic curve (ROC) that changes over time was used to evaluate the accuracy of predicting the prognosis of BLCA, and the package "survivalROC" in R was used to measure the area under the curve (AUC) of risk score, age, gender, stage, T, N, M.

Clinical relevance analysis
Box charts were used to analyze whether clinical features were associated with FRGs in the risk score model.

Bioinformatics analysis
Principal component analysis (PCA) was performed by using the "scatterplot3d" R package to profile the expression patterns of grouped samples. The single-sample gene set enrichment analysis (ssGSEA) in the R package "gsva" was utilized to calculated the immune infiltrating score of 16 immune cells and the expression of 13 immune-related pathways [29]. The relevant annotated gene set file was provided in Supplementary Table S2. In order to present the potential functions of differentially expressed genes (DEGs) between the low-risk and high-risk groups, the "clusterProfiler" R package was used to perform Gene Ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses.

Cell culture and siRNA knockdown
The human BLCA cell lines T24 cells were obtained from the Type Culture Collection in Chinese Academy of Sciences. T24 cells were cultured in DMEM medium (Gibco, Grand Island, NY, USA) containing 10% foetal bovine serum at 37 °C with 5% CO2. T24 cells were transfected with 20 nmol/L siRNAs using Lipofectamine 2000 RNAiMAX reagent (Invitrogen).

Cell proliferation assay
We seeded 2000 T24 cells per well into 96-well culture plates for 5 days in triplicate wells. Cell Counting Kit 8 (Gibco) was used according to the manufacturer's instructions. Then, 10 µL CCK-8 reagent was added to each well and incubated cells for 1-2 h. The optical density (OD) values of each well were determined at 450 nm using a microplate reader.

Statistical analysis
All statistical analyses were conducted using R software (Version 3.5.3) and SPSS (Version 23.0). The overall survival(OS) between the high-risk and low-risk group was compared by K-M analysis with the log-rank test. Univariate and multivariate Cox regression was employed to analysis independent predictors of OS.
Mann-Whitney test was used to compare the ssGSEA scores of immune cells or immune-related pathways between the different grops. If not specified instructions, Pvalue < 0.05 was considered statistically significant.

Identification of prognostic and FRGs in TCGA-BLCA cohort
The flowchart was a brief summary of the present study in Figure 1 and PGD are related to poor survival rates of BLCA patients. In contrast, the high expression of ACSL4, ALOX5 and NOX1 is associated with a higher survival rate in BLCA patients ( Fig.2a-b). The protein-protein interaction(PPI) network between 15 ferroptosis-related proteins indicated that PGD, G6PD and ABCC1 were the hub genes between these FRGs (Figure 2c). The correlation coefficient between these FRGs were presented in Figure 2d.

Construction a prognostic FRGs risk score model
According to the above prognostic-related FRGs, 7 prognostic-related FRGs were screened out through Multivariate Cox regression analyses ( Table 1). The prognostic model was used to establish and calculate the risk scores of all patients.
The patients were divided into a high-risk group (n=201) or a low-risk group (n=202).
The survival of the low-risk group was significantly higher than that of the high-risk

Independent Prognostic Analysis and Clinical Relevance Analysis
Among all the clinical features, univariate Cox regression analysis revealed that age, clinical stage, T stage, N stage, and risk score were responsible for OS in TCGA-BLCA (p < 0.05) (Figure 4a). Among them, risk score was an independent prognostic factor in multivariate Cox regression for BLCA patients (p < 0.05), but age, clinical stage, T stage, N stage were not (Figure 4b). Furthermore, we performed ROC curve to represent the prognostic prediction of this clinical features. The AUC of gender was 0.429; all other factors including age, clinical stage, T stage, M stage, N stage and risk score were over 0.500. Among them, the AUC score of risk score was 0.729 ( Figure 4c). Therefore, it is proved that the risk score calculated by the model can exactly forcast the 5-year survival rate of BLCA patients. Box charts were used to analyze the correlation between FRGs and clinical traits. The results were presented here: T stage was related to ALOX5, FADS2, and ZEB1; N stage was associated with FADS2, GCLM, and ZEB1; and clinical stage was associated with ALOX5, FADS2, GCLM, NOX1 and ZEB1 (Figure 4d-f).

The analysis of the high-and low-risk population by PCA and ssGSEA
We used PCA to detect the distribution patterns of different risk states between the high-and the low-risk group through genome-wide expression, FRGs and 7-FRGs signature sets. In the 7-FRGs signature set, the high-and low-risk group were clearly separated ( Figure 5a). While, the genome-wide expression and FRGs sets cannot distinguish high-and low-risk groups well. In other words, 7 FRGs were used to divide BLCA patients into two parts, which indicates that there is a big difference between BLCA patients in the high-and low-risk group (Figure 5b-c).
In order to explore the correlation between immune status and risk score, we

Functional analyses in the TCGA cohort
In order to clarify the biological functions and pathways related to risk score, DEGs between the high-and low-risk groups were used for GO enrichment and KEGG pathway analysis. The molecular functions and biological processes of these genes were obviously enriched in epidermal cell differentiation, skin development, keratinocyte differentiation, collagen-containing extracellular matrix and glycosaminoglycan binding. KEGG pathway analyses also demonstrated that PPAR signaling pathway, Retinol metabolism and Vascular smooth muscle contraction was the most paramount pathways in the TCGA cohorts. (P < 0.05, Figure. 6a-b).

Validation of gene expression and function
Next, we used ZEB1and FADS2 gene to verify our model. The Human Protein Atlas database analysis showed that ZEB1and FADS2 expression were significantly higher in BLCA tissues compared to normal bladder tissues (Figure 7a, b, d, e).
CCK8 proliferation analysis demonstrated that compared with control T24 cells, ZEB1and FADS2 knockdown T24 cells showed a significant reduction in proliferation (P < 0.05, Figure 7c, f).

Discussion
The conventional treatment of BLCA is surgery and chemotherapy, but this is The XC-system transfers cystine from the cell and converts it into cysteine for the synthesis of glutathione(GSH). The production and maintenance of GSH is essential to protect cells from oxidative stress. So GSH scavenge lipid peroxidation to inhibit ferroptosis [30]. The risk score prognostic model proposed in this study was composed of 7 FRGs (X1). These genes were roughly divided into two categories, ACSL4, ALOX5, ACACA, ZEB1 and FADS2 were lipid metabolism, NOX1 and GCLM were (anti)oxidant metabolism [9,10].
ACSL4 is one of the importent enzymes in the pro-ferroptotic lipid metabolism.
Knockdown of ACSL4 can reduce the ferroptosis of epithelial cells in lung ischemiareperfusion injury, and can protect lung epithelial cells [31]. ALOX5 is an ironcontaining non-heme dioxygenase that can limit leukotriene biosynthesis and catalyze the peroxidation of arachidonic acid, thereby mediating lipid peroxidation. Knockout ACACA can inhibit ferroptosis induced by pharmacological agents [34].
ZEB1 plays an important role in cell lipid metabolism. More and more evidences show that ZEB1 can regulates lipid uptake, accumulation and mobilization, and also affects plasma membrane remodeling associated with EMT to change the process of cell ferroptosis [35]. The knockout of it can inhibit ferroptosis caused by GPX4 depletion [35]. FADS2 is related to many chronic diseases, including obesity, type 2 diabetes and metabolic abnormalities. FADS2 is essential for maintaining the body's long-chain polyunsaturated fatty acid homeostasis. By knocking down FADS2 in lung cancer cells, cell growth rate is significantly reduced, intracellular iron and lipid ROS levels increased, and protein kinase-induced cell death increased [36,37]. NOX1 is a type of NADPH oxidases, which is the main source of reactive oxygen species (ROS) that regulates redox sensitive signaling pathways. The use of NOX1 inhibitors can significantly reduce the ROS, lipid ROS and cell death of non-small cell lung cancer induced by Erastin [38]. Because there are few researches on these genes, it remains to be clarified whether they play a role in the prognosis of BLCA patients by affecting the ferroptosis process.
By investigating the specific functions of 7 FRDs, previous studies have shown that most of these genes played a key role in cancer cells including BLCA. The abnormal expression of ZEB1 in BLCA is related to the differentiation and metastasis of bladder tumors. We used the Human Protein Atlas to verify the expression of ZEB1 in BLCA. ZEB1 was found to be highly expressed in the BLCA tissue when compared with normal bladder tissue. This protein can be used as a candidate target for prognosis prediction and early diagnosis [39]. In BLCA cells, knocking down NOX1 can reduce cell ROS production, leading to BLCA cell apoptosis [40].
Knockdown of ACSL4 expression can significantly attenuate the lipid peroxidation and ferroptosis induced by sorafenib in Huh7 cells, and save the sorafenib-induced growth inhibition of xenograft tumors in vivo. It is shown that ACSL4 is a key factor in sorafenib-induced ferroptosis and is useful for predicting the sensitivity of sorafenib in HCC [41]. Overexpression of ALOX5 in a PTC cell line (BCPAP) can increase cell invasion across the ECM barrier and the secretion of MMP-9. ALOX5 can be used as a new mediator for cell invasion induced by MMP-9 [42]. The mRNA and protein expression levels of FADS2 in human mesenteric tissues is reduced. The increased expression of FADS2 convert n-3 fatty acids into decomposable lipid mediators, resulted in a significant reduction in the infiltration of pro-inflammatory macrophages and weakened expression of inflammatory cytokines or adipokines.
FADS2 may improve Crohn's Disease treatment [43]. The Human Protein Atlas database analysis shows that FADS2 expression is significantly higher in BLCA tissues compared to normal bladder tissues. These results were consistent with previous bioinformatics analysis results (Figures 2a, b). Therefore, genes have different manifestations in diseases, and their specific roles have to be clarified.
Although the mechanism of ferroptosis has been the focus of research in the past, the potential relationship between ferroptosis and tumor immunity is still elusive. Our ssGSEA results showed that the scores of Macrophages, Mast_cells, pDCs, Tfh, Treg, APC_co_stimulation, CCR, Check-point and T_cell_co-inhibition were significantly higher in the high-risk group. It can be effectively assumed that ferroptosis may be closely related to tumor immunity. Studies have shown that the increase in tumorrelated Treg cells [44] or macrophages [45] are related to the poor prognosis of HCC patients. In this study, the antigen presentation process is different between the highrisk group and the low-risk group obviously. One speculation is that ferroptosis cells may release lipid mediators, to recruit antigen presenting cells (APCs) to the location of ferroptotically dying cells [46]. In addition, a higher risk score is associated with decreased anti-tumor immunity, including the fractions of pDCs, Tfh and Mast_cells, and the activity of Check-point and T_cell_co-inhibition. Therefore, the weakened anti-tumor immunity of high-risk patients may lead to poor prognosis. CCK8 proliferation assay shows that the proliferation ability of ZEB1and FADS2 knockdown T24 cells were significantly reduced compared with control T24 cells.
This suggests that ZEB1and FADS2 may play an oncogenic role in BLCA, but the specific mechanism needs further study.
This study also has several limitations. First, our prognostic model is constructed with retrospective data from TCGA databases. Some cellular studies and animal experiments should be conducted on the 7 FRGs alone or in combination to check the predictive accuracy of model and discover potential mechanisms. Second, it is necessary to use prospective real-world data to verify its clinical applicability. Third, using only a single hallmark to establish a prognostic model is flawed, as many important prognostic genes in BLCA may have been excluded. At last, the links between immune status and risk score have not yet been experimentally verified.

Conclusion
In conclusion, we construct a new prognostic model with 7-FRGs, which has been proved to be independently related to OS and can accurately predict the prognosis of BLCA. Understanding the underlying mechanism and significance of these FRGs for BLCA can provide insights for determining the therapeutic target of BLCA.

Acknowledgements
Not applicable.

Authors' contributions
Yingchun Liang collected data, analyzed data, and wrote the manuscript; Fangdie Ye analyzed data and edited the manuscript; Chenyang Xu collected and analyzed data; Lujia Zou analysed data; Yun Hu analysed data; Haowen Jiang designed research, edited the manuscript, and coordinated the project. All authors gave final approval of the manuscript.

Funding
Not applicable.

Availability of data and materials
All data generated or analyzed during the present study were downloaded from TCGA database (https://portal.gdc.cancer.gov/repository?facetTab=cases).