A Novel Ferroptosis-Related Prognostic Signature Reveals Macrophage Infiltration and EMT Status in Bladder Cancer

Bladder cancer (BC) belongs to one of the most common and highly heterogeneous malignancies. Ferroptosis is a newly discovered regulated cell death (RCD), characterized by accumulation of toxic lipid peroxides, and plays a crucial role in tumor progression. Here, we conducted a comprehensive analysis on the transcriptomics data of ferroptosis-related genes in BC based on The Cancer Genome Atlas (TCGA) and three Gene Expression Omnibus (GEO) datasets. In our study, a 6-gene signature was identified based on the potential prognostic ferroptotic regulatory genes. Furthermore, our signature revealed a good independent prognostic ability in BC. Patients with low-risk score exhibited higher FGFR3 mutation rates while high risk score had a positive association with higher RB1 mutation rates. Meanwhile, higher proportions of macrophages were observed in high BC risk group simultaneously with four methods. Unexpectedly, the risk score showed a significant positive correlation with epithelial-mesenchymal transition (EMT) status. Functional assays indicated that CRYAB and SQLE knockdown was associated with attenuated invasion capacity. Our study revealed a ferroptosis-related risk model for predicting prognostic and BC progression. Our results indicate that targeting ferroptosis may be a therapeutic strategy for BC.


INTRODUCTION
Bladder cancer (BC) belongs to one of most common malignancies in genitourinary system which can been further classified into two subtypes: muscle-invasive bladder cancer (MIBC) and non-muscle-invasive bladder cancer (NMIBC) (Siegel et al., 2020). Conventional treatment options for BC include surgery, radiation, and cisplatin-based chemotherapy. Yet, despite considerable advances in diagnosis and treatment, BC still exhibits high rates of recurrence and metastasis due to the high-level heterogeneity and genomic instability of BC (Alfred Witjes et al., 2017;Dy et al., 2017). BC responds differently to treatment resulting from different driver events. For example, BC patients with a mesenchymal related signature appears resistant to platinum-based chemotherapy and sensitive to immunotherapy (Wang et al., 2018a). Therefore, there is an urgent need for the development of novel therapeutic strategies so as to improve outcome.
Conventional therapy aims to suppress tumor cells through activating a certain kind of regulated cell death (RCD) (Stockwell et al., 2017). Ferroptosis is a unique cell death pathway driven by iron-dependent lipid peroxidation. Increasing evidence has determined the pivotal role of ferroptotic regulatory genes, such as TP53 (Xie et al., 2017), CISD1 (Yuan et al., 2016), GPX4 (Ingold et al., 2018), FINO2 (Gaschler et al., 2018), in tumor progression of numerous cancer types. In addition, chemotherapeutic drugs can improve efficacy on various tumor cells when combined with ferroptosis inducer Belavgeni et al., 2019). Notably, previous studies indicated the significant role of ferroptosis in diagnosis and treatment management of BC (Drayton et al., 2014;Guo et al., 2020). Therefore, ferroptosis can serve as a potential target for intervention in BC patients. However, there are still few scientific studies on the correlation between BC and ferroptosis. Despite significant progress in BC, few have considered the use of ferroptosis-related gene characteristics to construct a prognostic signature in BC. The precise underlying molecular mechanism and critical molecules of ferroptosis in BC progression remain to be illuminated.
Epithelial-mesenchymal transition (EMT), a plastic process in which epithelial cells gain mesenchymal characteristics, plays an important role in embryonic development (Sha et al., 2020;Yang et al., 2020a). Growing evidences showed that this physiological process was closely related to enhanced capacity for BC invasion (Schulz et al., 2019). Activation of EMT was reported to promote cell growth and metastasis in the BC . And we also proved this assumption in our previous research (Yan et al., 2020b). Moreover, evidence has shown that EMT-related gene and immune cell infiltration could impact outcomes in BC patients treated with immunotherapy. Infiltrated immune cells play important roles in iron homeostasis and ferroptosis (Wang et al., 2018b). Many immune cells, such as Th1 cells and macrophages, are reported to be involved in the maintenance of iron metabolism (Ganz and Nemeth, 2015). In addition, immunoregulation was found to efficiently inhibit tumor progression by synergistic ferroptosis. Reports claimed that increased ferroptotic level enhances the anti-tumor therapeutic effect of immunotherapy. Meanwhile, immunotherapy can activate T cells infiltration to promote the lipid-ROS formation and ferroptosis in tumor cells (Wang et al., 2019).
Herein, we evaluated publicly available BC datasets and identified differentially expressed ferroptosis-related genes strongly correlated to the prognosis of BC. Then we constructed a predictive risk model and validated its prognostic accuracy. Alterations of mutation profile and immune cell infiltration were also explored. Furthermore, functional enrichment analysis and the underlying mechanisms were ultimately confirmed via in vitro experiments. Our findings may help lead to a deeper understanding of BC progression and further provide novel therapeutic targets for BC.

Data Acquisition and Processing
All datasets used in this study were available to the public. Data of gene expression and clinical information were obtained from the Cancer Genome Atlas (TCGA) data portal, 1 the GTEx database, 2 GSE13507 dataset, GSE31684 dataset, GSE48075 datasets. Data from the GTEx database were selected to expand the subset of data from TCGA data portal. Then, Robust Multiarray Average was used to normalize the raw expression data (Irizarry et al., 2003). Sixty ferroptosis-related genes were obtained from the previous literature (Stockwell et al., 2017;Bersuker et al., 2019;Doll et al., 2019;Hassannia et al., 2019). HPA 3 is a platform that contains representative immunohistochemistry images expression data for common kinds of cancers (Thul et al., 2017). In this study, images of protein expression of SQLE and CRYAB between normal and BC samples were directly visualized by HPA.

Construction and Validation of the Prognostic Model
Differentially expressed genes (DEGs) between normal and tumor tissues were identified with the "limma" R package with a false discovery rate (FDR) <0.05 as cutoff threshold (Ritchie et al., 2015). BC patients in TCGA database were randomized into two groups at a ratio of 3:1 using the R package "caret." Univariate Cox regression analysis of overall survival (OS) was used to identify prognostic ferroptosis-related genes by "survival" R package filtered by p < 0.05. Then, "Venn" R package was implemented to get the intersect genes between ferroptosis-related DEGs and prognostic genes. For the selection of predictor variables, Least Absolute Shrinkage and Selection Operator (LASSO) regression analysis was applied to construct a prognostic model using "glmnet" on R (Simon et al., 2011). At last, a ferroptosis-related prognostic model was identified by extracting the coefficients, and the risk score for each case was determined by multiplying the level of each selected gene with its corresponding coefficient. The BC tissues were then divided into low-and high-risk group according to the median risk score. Univariate and multivariate Cox regression analyses were utilized to explore whether the risk score calculated from our model could play as an independent prognostic factor for BC patients after considering other clinical factors including gender, age, stage, T and N stage. The results were acquired by application of the "forestplot" package. To assess the predictive power of the signature, survival analysis and area under the curve (AUC) was calculated with "timeROC" R package. Somatic mutation data, which stored in the form of Mutation Annotation Format (MAF), were visualized and analyzed using "maftools." For subgroup analysis, patients were divided into different groups based on features as follows: gender (Male or Female), age (≤70 years old or >70 years old), stage status (Stage I, II or III, IV) and N status (N0 or >N0).

Pathway Enrichment Analysis
To understand the underlying mechanisms affecting the risk signature, BC samples were divided into two cohorts depending on risk score of each patient, and then gene set enrichment analysis (GSEA) enrichment analysis was performed via "javaGSEA" to show the result (Subramanian et al., 2005). Nominal p < 0.05 and FDR < 0.05 were considered significant. Enrichment of hallmark pathways within modules on a single sample to calculate a signaling pathway variation score was implemented through gene set variation analysis via "gsva" R package (Hanzelmann et al., 2013).

Evaluation of Immune Cell Infiltration
Transcriptome profiles from TCGA-BLCA cohort were used to estimate tumor immune infiltrate populations by TIMER, quanTIseq, TIP, and ssGSEA algorithm. TIMER 4 is a method that it only makes estimations of six immune cell types, taking the tissue specificity into consideration . QuanTIseq performs an absolute quantification of cell types in the samples (Finotello et al., 2019). Tumor ImmunoPhenotype (TIP) pipeline was applied for immune activity estimation . Single sample Gene Set Enrichment analysis (ssGSEA) was applied to evaluate the enrichment scores.

Construction and Evaluation of the Nomogram
The clinical features and risk signature were extracted to construct a prognostic nomogram to assess the OS probability as a single numerical value. The "rms" R package was used to draw the nomogram (Feng et al., 2017). In addition, Calibration plots were generated to explore the performance characteristics of the nomograms. The clinical usefulness of the risk model was evaluated using decision curve analysis (DCA).

Western Blotting
Western Blotting was performed using the routine methods as described before (Yan et al., 2020b). Briefly, cells were lysed by RIPA buffer, then 20 µg proteins were separated by SDS-PAGE and transferred onto PVDF membrane. After blocking with 5% non-fat milk in PBST, the membrane was incubated with primary antibodies overnight at 4 • C. β-actin was used as a loading control. The membranes were then incubated with secondary anti-body respectively. The Western blots were visualized using the ECL substrate kit. The antibodies against CRYAB, SQLE were purchased from Proteintech (Rosemont, IL, United States). β-actin were purchased from Sigma-Aldrich (Saint Louis, MO, United States).

Transwell Invasion Assay
8.0-µm-pore polycarbonate membrane were used to assess cell invasion ability of BC cells as previously described (Yan et al., 2020b). Briefly, after the Matrigel was coagulated, 1 × 10 5 cells were plated on the Matrigel. After incubation, cells were fixed in 4% formaldehyde and stained with crystal violet. Cells were visualized using light microscopy. Figure 1 shows the study flow chart. A total of 403 BC patients from the TCGA-BLCA, 165 patients from the GSE13507, 93 patients from the GSE31684 and 73 patients from the GSE 48075 were finally included. After data collection, 60 ferroptosis-related genes were identified. DEGs between normal and tumor tissues within these genes were identified and were subsequently analyzed by Cox regression and LASSO regression to construct a prognostic ferroptosis-related risk signature. This signature was further evaluated using Cox regression analysis, ROC curve analysis, clinical stratification analysis, mutation correlation analysis, immune cell infiltration comparison, and external dataset validation. Finally, an optimized model and nomogram were established. The three GEO datasets were subjected to validation of all results. To determine the prognostic value of the above ferroptosis-related genes in BC, univariate Cox regression analysis identified 16 genes as significantly correlated with BC survival (Figure 2A), and 9 of them were differentially expressed between tumor tissues and adjacent carcinoma tissues ( Figure 2B). Expression correlations between these genes is shown in Figure 2C. LASSO and multivariable Cox regression analysis were performed to identify a risk  Figure 1). A 6-gene signature was identified and risk score was calculated by following formula:

Identification of Bladder Cancer-Specific Ferroptosis Related Genes and Construction of the Prognostic Model
Next, risk score of each BC tissue was calculated and patients were assigned into a low-or a high-risk group according to the median cut-off value ( Figure 2D). Survival status of the BC patients was visualized in Figure 2E. The proportion of mortality in high-risk group was 48%, much higher than 27% in the low-risk group. PCA analysis showed clear separation of patients in different risk groups ( Figure 2F). Kaplan-Meier survival curve revealed significant difference in OS between the groups. As demonstrated, more deaths happened in the highrisk groups than those in low-risk groups, suggesting that the risk model could accurately distinguish BC patients with poor prognosis. We then evaluated the predictive power and accuracy of the risk signature with ROC curve analysis. A time-dependent ROC curve was performed and the AUC was 0.709 at 3-years and 0.725 at 5-years (Figures 2G,H).

Validation of the 6-Gene Risk Model
To evaluate whether the risk model exhibits similar predictive performance and accuracy in other BC patient cohorts, robustness of the 6-gene signature was then tested in TCGA internal cohort and another three independent GEO BC cohorts. Survival analysis results confirmed that patients of high-risk correlated with poor OS (Figures 3A-D). Figures 3E-H, the AUC was 0.677, 0.684, 0.625, 0.690 for 5-year survival prediction and 0.648, 0.689, 0.603, 0.636 for 3-year survival prediction in TCGA, GSE48075, GSE31684, and GSE13507. Heat maps were shown to present the expression levels of the six genes and clinical features ordered by risk score in the above four datasets. Clinical and pathological features, such as TNM stage and grade were enriched BC cases with highrisk score (Figures 3I-L). Therefore, we believe that the ferroptosis-based risk model indicates good prediction stability and performance.

Prognosis Analysis of the 6-Gene Risk Signature With Clinicopathological Features
Relationship between OS, clinicopathological parameters and risk score were further analyzed to assess whether the risk score was an independent predictor of outcome. Univariate and multivariate Cox analysis were used based on OS of BC patients, using the co-variables including risk score, age, gender, stage, T and N stage to validate the independence of the risk model among other clinic-pathologic characteristics. Results from univariate and corresponding multivariable Cox regression analysis demonstrated that our risk model could serve as an independent prognostic factor (p < 0.05; Figures 4A, B). Then we further evaluated the effects of clinical characteristics, such as different stage and lymph node metastatic status on outcome among BC patients, and found that high risk score was highly correlated with poor outcome (Figures 4C-J). These results demonstrated that the ferroptosis-related prognostic model

Construction and Evaluation of a Predictive Nomogram
Furthermore, we utilized a quantitative method by integrating the risk score and clinical features to construct a nomogram ( Figure 5A). Female patients had a higher risk of a poor prognosis while the higher age, tumor stage and risk score indicated a lower survival rate in patients. Combined with our risk model and clinical feature, net benefits were presented in our DCA curve, and we found that compared to a single factor, the combined showed the optimal net benefit ( Figure 5B). Meanwhile, the calibration plot for the probability of 3-, 5-, and 10-year survival showed an optimal consistency between observation and predictive curves (Figures 5C-E). These results indicate that the nomogram has proper clinical applicability.

Correlation of Mutation and Immune Cell Infiltration Landscape With the Risk Model
To determine whether the risk model was correlated with tumor mutation of BC, relationship between the model and gene mutations was analyzed. The profile of somatic mutation is visualized in Supplementary Figure 2. The mutation profile distribution based on patients' risk score were shown in Figures 6A,B. Notably, the low-risk tumor groups have significantly lower rates of FGFR3 mutation while the overall mutation rates of RB1 appeared to be higher in the high-risk subgroup. FGFR mutation was reported to be significantly enriched in the luminal papillary subtype, characterized by lower stage, while high RB1 mutation frequencies are significantly high in basal-squamous and neuronal subtypes, indicating lower survival rate. Then we evaluated the TIICs proportions with the aid of TIMER, quanTIseq, TIP and ssGSEA, four well-accepted methods for comprehensive analysis of TIICs. Two distinct patterns of immune infiltrate were observed in high and low-risk cohorts. Surprisingly, tissues with high risk generally had higher density of macrophage (Figures 6C-F). Increased proportions of tumor-infiltrating macrophages have been reported to be involved in tumor progression and metastasis of BC. These evidences demonstrated that high RB1 mutation and tumorinfiltrating macrophages may account for the poor prognosis in high-risk group.

Identification of the 6-Gene Risk Model Correlated Biological Pathways
Gene set enrichment analysis (GSEA) analysis of BC patients with a different risk score was performed to investigate the possible biological function of the risk model in the carcinogenesis of BC. The results showed that, in four BC groups, high-risk score was positively correlated to important tumor related pathways, especially "HALLMARK_EMT" and "HALLMARK_HYPOXIA" (Figures 7A-D). Meanwhile, the GSVA scores of such two pathways were significantly higher in high-risk score cohort, respectively ( Figure 7E). Survival analysis showed a positive correlation between EMT score and the OS in patients with BC (Figure 7F), while the difference was not significant for HYPOXIA score (Figure 7G). In addition, expression of EMTrelated markers was elevated in high-risk group (Figure 7H). These results indicated that the risk model played a critical role in predicting EMT-related function in BC.

Functions of the Identified Biomarkers in BC Invasion
EMT progress was previously reported to contribute to cancer cell invasion, resulting in rapid tumor development. Within the six genes, higher level of CRYAB, SQLE, and ZEB1 significantly correlated with advanced stage (Figures 8A-C). Considering the important role of ZEB1 in tumor metastasis, we thus focused on another two genes. Then we performed functional studies after transfection using specific small interfering RNAs (siRNAs) (Figures 8D, E). The results showed that knockdown of CRYAB and SQLE in BC cell lines inhibited cell migration and invasion ability (Figures 8F, G). Moreover, overexpression of these two genes promoted the invasive ability of bladder tumor cell (Figures 8D-G). Immunohistochemical analysis indicated the positive staining intensity for CRYAB and SQLE in BC tissues as significantly stronger than in the normal urothelium tissues. Moreover, BC tissues with high malignancy exhibited strong staining intensity than tissues with low malignancy by HPA ( Figure 8H). These results confirmed that CRYAB and SQLE were correlated well with the grade of tumorigenesis.

DISCUSSION
Bladder cancer (BC) has been reported among the most common malignancies, presenting a major health burden for society. How to optimize therapeutic protocols for mortality reduction of patients with advanced BC remains a challenge. Therefore, there is an urgent need to identify key biomarkers that affect the outcome of BC. With future expansion of the database and multi-omics data, improved data mining algorithms can have an essential impact on tumor biology (Angus et al., 2019;Nacev et al., 2019;Liu et al., 2020;Tabassum et al., 2020). Transcriptome profiling provided us with novel insights into assessing prognostic of the individual patient when combining the corresponding clinical information. Multiple risk signatures have been identified in various tumors for the past few years (Yu et al., 2019;Yan et al., 2020a). However, most studies have failed to be applied to clinical practice owning to lack of systematic evaluation and a broad roll-out. Therefore, these questions illustrate the urgent need to identify the prognostic factors of BC to predict to identify high-risk populations.
Increasing evidence has shown that induction of cell death is among the most effective anticancer strategy and ferroptosis belongs to a RCD in which iron metabolism plays a vital role (Badgley et al., 2020;Yang et al., 2020b). However, molecular changes and mechanism of ferroptosis in BC has yet to be elucidated. In the current research, we have identified a ferroptosis-related risk model for BC. We used BC samples from TCGA as training group, GSE13507, GSE48075, and GSE31684 as validation groups. Our risk model is composed of6 differentially expressed ferroptosis-related genes (CISD1, CRYAB, FTH1, ACACA, ZEB1, SQLE) that correlate with patient outcome. We found that risk scores of our model were highly associated with stage and metastatic status. We further demonstrated that the risk model serves as an accurate predictor of BC survival and an independent predictor for BC prognosis by validating it in another three independent BC datasets. Besides, the nomogram performed well with a good calibration, indicating that the model is an accurate prognostic tool. Our results suggested that the model can well distinguish BC patients and predict prognosis, thereby helping to develop optimal treatment options based on risk score.
We identified several ferroptosis-related genes that predict the outcome of BC patients. Most of the genes have been reported in previous studies to be closely related to ferroptosis and cancer development of malignancies. CDGSH iron sulfur domain 1 (CISD1) is an iron-containing protein and could negatively modulates ferroptosis. Increasing stabilization of CISD1 inhibits erastin-induced mitochondrial iron uptake and oxidative damage (Yuan et al., 2016). Crystallin alpha B (CRYAB, also known as HSP beta-5) belongs to the small heat shock protein (HSP20) family and is reported to regulate iron uptake and GPX4 abundance (Stockwell et al., 2017). Ferritin heavy chain 1 (FTH1) is one of the subunits of Ferritin and plays an essential role in cellular iron balance in ferroptosis. Knockdown of FTH1 could lead to iron overabsorption and promote ferroptosis in the intestines of mice (Sun et al., 2016). acetyl-CoA carboxylase alpha (ACACA) is implicated in catalyzing the committed steps in fatty acids' biosynthesis A and was reported to suppress FIN56-, but not erastin-or RSL3-induced ferroptosis (Shimada et al., 2016). Zinc finger E-box-binding homeobox 1 (ZEB1) function as a transcription factor that influences cell developmental and homeostasis cell fates. Knockout of ZEB1 suppressed GPX4-depletion-induced ferroptosis (Viswanathan et al., 2017).
Squalene monooxygenase (SQLE) is a key rate-limiting enzyme in the biosynthesis of cholesterol. Overexpression of SQLE sensitizes ALCL cells to ferroptosis (Garcia-Bermudez et al., 2019).
Studies have indicated that gene mutations play an important role in the development of carcinoma including BC (Lawson et al., 2020). Mutations of the several genes, such as p53, FGFR3 and RB1, have been reported to be involved in BC with a high incidence (Robertson et al., 2018). Therefore, we further attempt to study whether there were gene mutation alterations between different groups in BC patients. Interestingly, in our study, higher FGFR3 mutation rate was observed in low-risk BC patients. FGFR mutation was significantly enriched in the luminal papillary subtype, characterized by lower stage, lower risk for progression and papillary morphology. These results suggested specific inhibitors designed against FGFR3 as a treatment option for these patients. Conversely, patients of high-risk score accounted for higher RB1 mutation. High RB1 mutation frequencies are significantly high in basal-squamous and neuronal subtypes, indicating worse OS and resistance to platinum-based chemotherapy. Immune cells were reported to be involved in the process of iron metabolism and ferroptosis level could promote the anti-tumor efficacy of immunotherapy (Nairz and Weiss, 2020;Tang et al., 2020). Interestingly, the infiltrating immune cells exhibit significant disparities among different groups in our study. Additionally, with four wellaccepted methods, our signature identified higher proportion of macrophages in the high-risk group compared with the low-risk group. Previous researches have shown that increased proportions of tumor-infiltrating macrophages contributed to a propensity of metastases and poor prognosis for BC (Chen et al., 2018;Kobatake et al., 2020). Therefore, high infiltration of tumorassociated macrophages in patients with high risk may be an explanation for their poor outcome.
To investigate the underlying mechanisms, GSEA and GSVA were subsequently subjected to the pathway enrichment of each group. Interestingly, the enrichment analysis revealed that high risks were significantly enriched in EMT and hypoxia hallmark pathways. Of these two, only EMT was significantly associated with poor prognosis. EMT is an essential biological process, which plays a vital role in modulating tissue homeostasis and development (Bakir et al., 2020). Besides, EMT has been implicated in tumor progression and metastasis. In our study, high levels of EMT gene expression were observed in highrisk BC patients. Among genes of the risk model, we identified three genes (CRYAB, SQLE, and ZEB1) that were found to be was significantly and positively associated with clinical stage. Numerous studies have reported that ZEB1 contributes to cancer progression, while the role of CRYAB and SQLE in BC is rarely reported. We found that knockdown of CRYAB or SQLE significantly attenuated the invasive abilities of BC cells. Nevertheless, the potential mechanism of the two biomarkers might contribute to the carcinogenesis of BC remain unknown, and further investigation of potential mechanisms is needed.
Nonetheless, some limitations are inevitable. This risk model is highly dependent upon public databases. Further investigations need to be undertaken in future clinical researches. Moreover, protein level could differ from RNA expression, making it unavoidable to validate its clinical utility in more sets.
In summary, this research is the first to identify a novel signature of ferroptosis-related genes for predicting outcomes of BC patients. Results suggest that variation in gene mutation, immune response and EMT status might be several possible reasons for this model's prognostic ability. These findings may help unveil new targets for the prevention, diagnosis, and treatment of BC.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.

AUTHOR CONTRIBUTIONS
BS and SX: conception and design. YY and FZ: acquisition of the data. YY and JC: analysis of the data. YY, ZH, and XC: writing, review, and revision of the manuscript. PT and ZW: administrative, technical, or material support. BS: study supervision. All authors read and approved the final manuscript.