Biological Functions and Prognostic Value of Ferroptosis-Related Genes in Bladder Cancer

Background: Every year, nearly 170,000 people die from bladder cancer worldwide. A major problem after transurethral resection of bladder tumor is that 40–80% of the tumors recur. Ferroptosis is a type of regulatory necrosis mediated by iron-catalyzed, excessive oxidation of polyunsaturated fatty acids. Increasing the sensitivity of tumor cells to ferroptosis is a potential treatment option for cancer. Establishing a diagnostic and prognostic model based on ferroptosis-related genes may provide guidance for the precise treatment of bladder cancer. Methods: We downloaded mRNA data in Bladder Cancer from The Cancer Genome Atlas and analyzed differentially expressed genes based on and extract ferroptosis-related genes. We identified relevant pathways and annotate the functions of ferroptosis-related DEGs using Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis and Gene Ontology functions. On the website of Search Tool for Retrieving Interacting Genes database (STRING), we downloaded the protein-protein interactions of DEGs, which were drawn by the Cytoscape software. Then the Cox regression analysis were performed so that the prognostic value of ferroptosis-related genes and survival time are combined to identify survival- and ferroptosis-related genes and establish a prognostic formula. Survival analysis and receiver operating characteristic curvevalidation were then performed. Risk curves and nomograms were generated for both groups to predict survival. Finally, RT-qPCR was applied to analyze gene expression. Results: Eight ferroptosis-related genes with prognostic value (ISCU, NFE2L2, MAFG, ZEB1, VDAC2, TXNIP, SCD, and JDP2) were identified. With clinical data, we established a prognostic model to provide promising diagnostic and prognostic information of bladder cancer based on the eight ferroptosis-related genes. RT-qPCR revealed the genes that were differentially expressed between normal and cancer tissues. Conclusion: This study found that the ferroptosis-related genes is associated with bladder cancer, which may serve as new target for the treatment of bladder cancer.


INTRODUCTION
The mortality of Bladder cancer ranks 13th among cancers (Antoni et al., 2017). To make a diagnosis and give treatment of bladder cancer, prognosis and risk assessment are important. The TNM staging system has been widely applied to predict the prognosis in cancer patients (Kluth et al., 2015). In recent years, the diagnostic or prognostic prediction model established by combining molecular markers and staging systems has improved the diagnosis and treatment ability of bladder cancer (Chen et al., 2020a;Yang et al., 2020). However, most molecular markers are derived from differential gene analysis, and the function of genes has not been further studied. Therefore, it is necessary to explore new functional biomarkers.
Ferroptosis is a new type of programmed cell death that is iron-dependent (Latunde-Dada, 2017;Wang et al., 2020b;Riegman et al., 2020). Ferroptosis is associated with many important physiological processes such as mitochondrial function, lipid metabolism and oxidative stress, so that induction of ferroptosis play a pivotal role in inhibiting the progression of tumor (Torti and Torti, 2020). There are increasing indications that ferroptosis can produce a marked effect in tumor treatment, and several studies have focused on cancer therapies that rely on ferroptosis (Conrad et al., 2020;Wu et al., 2020;Zuo et al., 2020). Wang et al. identified branched-chain amino acid aminotransferase 2 (BCAT2) as a novel inhibitor of ferroptosis using CRISPR/ Cas9-based assays in hepatoma cells . Li et al. (2019) constructed CaO 2 and Fe 3 O 4 nanoparticles and utilized cascade reaction-mediated ferroptosis to synergize with immune regulation in cancer treatment (Li and Rong, 2020). In a related bladder cancer study, Chen et al. (2020b) found that the quinazolinyl-arylurea derivatives can induce ferroptosis based on the structural modification of sorafenib. Therefore, genes associated with ferroptosis may predict the survival of patients with bladder cancer and new personalized therapeutic targets.
We collected mRNA expression data and corresponding clinical data from BLCA-TCGA and analyzed the biological functions of the ferroptosis-related DEGs. We then constructed a prognostic 8-gene signature associated with ferroptosis, which was validated by survival analysis and risk assessment. In addition, we performed biological functions analysis to investigate the potential mechanism of ferroptosis-related genes. These results will help to afford new insights for the prognosis and treatment of bladder cancer.

Data Source
The mRNA data and clinical data of bladder cancer patients were obtained from the Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). Clinical features of bladder urothelial carcinoma patients (BLCA) from the TCGA database are shown in Table 1. The ferroptosis-related genes were then downloaded on the FerrDb database (http://www. zhounan.org/ferrdb) (Zhou and Bao, 2020). A list of related genes is provided in the Supplementary Material (Supplementary Table S1). Combined with the transcriptome sequencing map of bladder cancer, ferroptosis-related genes of bladder cancer were obtained. All data were publicly available online.

Data Processing of Ferroptosis-Related Differentially Expressed Genes
Using the R software, we analyzed ferroptosis-related differentially expressed genes (DEGs) to identify differences between tumor and sample groups. The difference is significant when | log2Foldchange | > 1 and false discovery rate (FDR) < 0.05.

Functional and Pathway Enrichment
Kyoto Encyclopedia of Genes and Genomes (KEGG) is a pathway-related database that systematically analyzes gene   Frontiers in Molecular Biosciences | www.frontiersin.org November 2021 | Volume 8 | Article 631152 3 function and links genomic information and functional information. Gene Ontology (GO) is an international standardized classification system for gene functions that divides the functions of genes into three parts: cellular component (CC), molecular function (MF), and biological process (BP) (Yu et al., 2012). Bar graphs and bubbles were

Protein-Protein Interaction Network
We analyzed the PPI network of ferroptosis-related DEGs in STRING (https://string-db.org/). Subsequently, the PPI network was analyzed and visualized by Cytoscape 3.7.1, and the most relevant subnetwork module was obtained using Cytoscape plug-in MCODE.

Identification of Prognostic Ferroptosis-Related Differentially Expressed Genes
We performed Kaplan-Meier survival analysis to assess the effect of ferroptosis-related genes on overall survival (OS) in bladder cancer patients. Then we conducted Cox regression. Prognostic ferroptosis-related genes (p < 0.05) were screened out. We also retained the ferroptosis-related DEGs if | log2Foldchange | >1 and

Construction and Analysis of Prognostic Models
We analyzed the prognostic value and survival time of ferroptosis-related DEGs using Cox regression analysis to identify ferroptosis-related genes associated with survival and generated forest maps. Univariate Cox regression analysis was conducted to identify the OS-related core genes, and genes with p < 0.1 were utilized for the subsequent multivariate Cox regression analysis. Thereafter, samples from TCGA were randomly divided into two groups in a 1:1 ratio to construct the best prognostic model, including a training group (n 197) and a testing group (n 196). We identified eight survival-related ferroptosis-related genes were by this model and obtained the correlation coefficients of each gene. Exp(GENE n) *coef (GENE n), where "Exp" represents the gene expression, and "coef" represents the coefficient score from multivariate Cox analysis. The risk score of each patient was calculated based on this model. Then, we analyzed the survival analysis, receiver operating characteristic curve (ROC) and risk curves of each group. Furthermore, by univariate and multivariate analysis, a nomogram was generated to provide patient survival information.

Preparation of Bladder Cancer Clinical Tissue Samples and Adjacent Tissue Samples
Throughout the study, we obtained informed consent from each patient for the collection and analysis of tissue samples, which was approved by the Union Hospital, Tongji Medical College, and Huazhong University of Science and Technology ethics review boards. We immediately frozen and stored the tissue at −80°C after extracting it. Total Ribo Nucleic Acid Extraction, Reverse Transcription, and Real-Time

Quantitative-Polymerace Chain Reaction Analysis
We extracted total RNA from tissue samples using Trizol reagent and used the nanodrop2000 to assess RNA integrity. We used reverse transcription reagents to react with extracted RNA to produce cDNA complementary to the single strand of RNA, which was then quantified in real time using the SYBR Green PCR kit. Each gene's cycle threshold (CT) was recorded. The relative mean value of each gene's expression was calculated, and the 2 -ΔΔCT method was used to determine whether the differences between groups were statistically significant, and the mean value was used as the final experimental result for replicate wells. All procedures were carried out in accordance with the manufacturer's instructions.

Statistical Analysis
The predictive significance of ferroptosis associated genes and other clinical characteristics was determined by using univariate and multivariate Cox regression. The effect of this model was assessed by using Kaplan Meier curves and the log rank test method to validate the difference in survival time between different risk groups based on the ferroptosis prognostic model. The effectiveness of this ferroptosis prognostic model to predict the prognosis of bladder cancer patients was assessed using a receiver operating characteristic curve (ROC). p values less than 0.05 were deemed statistically significant in two-sided statistical testing. The graphing software graphpad prism version 8.0.1 and the R language version 4.0.2 were used in this investigation. The R packages mainly used in this paper include limma, pheatmap, survival, iGraph, reshape2, glmnet, clusterprofiler, org. hs.eg.db, enrichplot, ggplot2.

Identification of Differentially Expressed Ferroptosis-Related Genes
The mRNA expression levels of all ferroptosis-related genes were compared between 19 control bladder tissues and 412 bladder cancer tissues in the TCGA-BLCA dataset. 103 ferroptosis-related DEGs were obtained, including 57 upregulated ferroptosis-related genes and 46 downregulated ferroptosis-related genes, and plotted volcano and heat maps ( Figure 1).

Functional Enrichment Analyses of Differentially Expressed Ferroptosis-Related Genes
To study the biological functions of 103 differentially expressed ferroptosis-related genes, GO function and KEGG pathway enrichment were performed, and bar plots and bubbles were plotted. As shown in Figure 2, the top 10 most relevant genes were selected. GO results indicate that differentially expressed ferroptosis-related genes are involved in some important biological processes, such as response to stimuli, maintenance of cell homeostasis, cell autophagy, apoptosis, and a variety of iron-related functions. KEGG pathway enrichment analysis indicated that 103 genes were mainly involved in PD-L1/PD-1 pathway in cancer, multiple signaling pathways related to cell differentiation, tumorigenesis and development, and autophagyrelated pathways.

Protein-Protein Interaction Network Construction
We used STRING tools to predict the protein interactions among the differentially expressed ferroptosis-related genes. We found 101 nodes and 348 edges in the PPI network. Low degree values were associated with smaller node sizes, while lower co-expression values were associated with smaller edge sizes. Thereafter, a network diagram of 102 genes was drawn using the Cytoscape software, as shown in Figure 3A. In addition, two key subnets were extracted and visualized using the MCODE plug-in, as shown in Figure 3B. Two most significant MCODE modules visualization and two most significant MCODE components form the PPI network as shown in Figures 3C-E. The results of KEGG pathway enrichment and GO function on the genes of the two subnetworks were shown in Table 2 and Table 3.

Construction and Analysis of Prognostic Signature
103 differentially expressed ferroptosis-related genes interacting with survival time were analyzed to identify the prognostic value of by Cox regression; 18 survival-related genes were obtained, and forest maps were generated ( Figure 4A). The samples from TCGA were randomly divided into a training set (n 197) and test set (n 196) according to a 1:1 ratio. Eight survivalrelated genes were screened out through the training set and are used to construct the best prognostic signature, and the correlation coefficients for each gene were obtained. These eight genes were: iron-sulfur cluster assembly enzymes (ISCU), nuclear factor, erythroid 2 like 2 (NFE2L2), MAF bZIP transcription factor G (MAFG), zinc finger E-box-binding protein 1 (ZEB1), voltage-dependent anion channel 2 (VDAC2), thioredoxin-interacting protein (TXNIP), stearoyl-CoA desaturase (SCD), and Jun dimerization protein 2 (JDP2) ( Figure 4B). Among them, ISCU, NFE2L2, and TXNIP are classified as low-risk genes, while ISCU, MAFG, ZEB1, VDAC2, SCD, JDP2 are classified as high-risk genes. In addition, details of the prognostic signature are listed in Table 4.
We obtained a risk score for each patient and obtained the median of the risk scores for all patients according to the risk formula. We classified patients below the median as a low-risk group and those above the median as a high-risk group. The patients were divided into high-risk and low-risk groups in the training and test groups, respectively, according to the risk score ( Figure 5A). In the training group, there were 98 patients in the high-risk group and 99 patients in the low-risk group. In the test set, there were 93 patients in the high-risk group and 103 patients in the low-risk group. As shown in the results, patients with high-risk scores had shorter survival times than those with low-risk scores. The ROC curve of survival prediction showed satisfactory performance ( Figure 5B). The AUC (Area Under Curve) value in the training set was 0.766, and that in the test set was 0.674 ( Figure 6). The closer the AUC value is to 1.0, the more authentic the detection method is. The same analysis was also performed in the test group with consistent results ( Figure 5C,D).
We performed univariate and multivariate independent prognostic analysis. Age, tumor stage, and the risk score could be used as independent prognostic factors for survival of bladder cancer patients in the training and validation sets (p < 0.05) in the univariate independent prognostic analysis. Sex, age, tumor stage, and the risk score were independent prognostic factors of bladder cancer (p < 0.05) in the multivariate independent prognostic analysis. For the test set, age was not found to be a prognostic factor (p 0.889) (Figures 7A-D). All the prognostic factors were included in the nomograms ( Figure 7E).
Finally, the nomogram of these eight prognostic genes was drawn in the training set. The RNA expression of eight genes was used as a parameter to draw the dotted lines in the nomogram and the total score can predict the survival rates.

Validation of Ferroptosis-Related Genes
We further validated the expression of these eight genes in 12 pairs of clinical sample tissues (Figure 8). The results showed that the expression levels of ISCU and JDP2 were lower in tumor samples and higher in paracancerous tissues. In clinical sample tissues, the expression of SCD was exactly the opposite. Although the p-value was not less than 0.05, we found that MAFG with VDAC2 presented higher expression in paracancerous tissues than in cancer tissues, and TXNIP presented lower expression in paracancerous tissues than in cancer tissues. The low abundance of nfe2l2 expression in tissues precludes its detection by qPCR. Because the number of clinical samples was small and the plotting of survival curves was difficult, we divided the patients into high-risk and low-risk patients by the risk model derived in the previous study and found a statistical difference in survival time between these two groups. Low-risk patients had a longer survival time than high-risk patients ( Figure 8H).

DISCUSSION
Bladder cancer has a certain recurrence rate after operation. Therefore, it is meaningful to find new diagnostic biomarkers or therapeutic targets for bladder cancer. Ferroptosis is associated with many physiological and pathological processes, such as tumor suppression, neuronal degeneration, antiviral immune response, and ischemia-reperfusion injury. From the perspective of cancer therapy, we hope to eliminate cancer cells by promoting ferroptosis or inhibiting ferroptosis in healthy cells. However, the specific mechanism of ferroptosis remains to be investigated in bladder cancer.
We integrated mRNA data from TCGA database and identified 103 ferroptosis-related DEGs, including 57 upregulated and 46 downregulated ferroptosis-related genes. In addition, through functional analysis, we found that ferroptosis-related DEGs were involved in apoptosis and immune checkpoint pathways. We then established a prognostic model based on eight ferroptosis-related genes to predict the accurate survival for patients with bladder cancer based on these results. By combining important clinical data and risk scores, a nomogram was constructed to predict survival rates.
Ferroptosis-related genes in this study included ISCU, NFE2L2, MAFG, ZEB1, VDAC2, TXNIP, SCD, and JDP2. These genes have been studied in metabolic processes or tumor development. ISCU is an important scaffold protein and can assemble iron and sulfur into iron-sulfur clusters These clusters participate in the maturation of [2Fe-2S] and [4Fe-4S] proteins in the mitochondria and cytoplasm and the regulation of iron metabolism (Duan et al., 2009). As an important redox center, iron-sulfur clusters participated in many physiological functions such especially in metabolism (Xu and Møller, 2011). NFE2L2 is transcription factor, which regulate genes that contain antioxidant response elements in their promoters. NFE2L2 plays a vital role in the treatment of neurodegenerative diseases and regulation of ferroptosis (Liu et al., 2020;Song and Long, 2020). MAFG is a small MAF protein (Liu et al., 2018). MAFG-mediated cisplatin resistance lies in lower reactive oxygen species production; therefore, MAFG may be a potential therapeutic target (Vera-Puente et al., 2018). ZEB1 is a transcriptional repressor. Studies have shown that ZEB1 is abnormally expressed in many liver diseases, including hepatocellular carcinoma . Studies have found that ZEB1 overexpression increases sensitivity to ferroptosis (Lee et al., 2020). VDAC2 is a member of the voltage-dependent anion channel pore-forming protein family, which is thought to play an important role in the metabolite diffusion process on the outer mitochondrial membrane. Moreover, it can participate in the mitochondrial apoptotic pathway by regulating the activity of BCL2-antagonist/killer 1 protein. TXNIP is a thioredoxinbinding protein, which participates in cellular redox signaling and plays an important role in oxidative stress. TXNIP has been found to be involved in the autophagy and ferroptosis of a class of protein disulfide isomerase inhibitors in the treatment of glioblastoma (Kyani et al., 2018). SCD encodes enzymes involved in fatty acid biosynthesis. In one study, SCD expression in bladder cancer patients was found to be significantly associated with poor prognosis. Loss-of-function experiments further revealed that inhibiting SCD expression can reduce cell proliferation and invasion (Piao et al., 2019). JDP2, a member of the bZIP transcription factor family and acts as a broad AP-1 inhibitory protein. JDP2 plays an important role in inhibiting AP-1-driven biological processes, such as cell cycle regulation, cell differentiation, apoptosis, and tumorigenesis Frontiers in Molecular Biosciences | www.frontiersin.org November 2021 | Volume 8 | Article 631152 (Chen et al., 2017). In hepatocellular carcinoma, miR-501 can regulate JDP2 to promote tumor cell development (Yu et al., 2019). In conclusion, these eight genes are mainly involved in iron metabolism, lipid metabolism, and amino acid metabolism.
In general, we investigated the biological function and prognostic value of ferroptosis-related genes and established a new prognostic model based that are closely related to the overall survival of bladder cancer patients.
The main limitation of this study is the small number of clinical samples used for validation, which requires further investigation. And we need to experimentally validate the hypothesis. Also, more validation research is needed to explore the mechanism of ferroptosis-related genes in bladder cancer. In spite of all its limitations, our results suggest that prognostic markers based on ferroptosisrelated genes may be a reliable predictive tool for the prognosis of bladder cancer patients.

DATA AVAILABILITY STATEMENT
The data analyzed in this study were obtained from publicly available datasets. The mRNA data and clinical can be found here: The Cancer Genome Atlas (https://portal.gdc.cancer.gov/). The ferroptosis-related genes are downloaded from the FerrDb database (http://www.zhounan.org/ferrdb).