Multiomics analysis of ferroptosis-related molecular subtypes in muscle-invasive bladder cancer immunotherapy

Introduction

Bladder cancer is the most common malignancy of the genitourinary system; the incidence of bladder cancer is ninth, and the mortality rate is thirteenth of all malignant tumors worldwide (1). Muscle invasive bladder cancer (MIBC) accounts for approximately 20% of newly diagnosed bladder cancers. Although radical cystectomy (RC) and pelvic lymph node dissection are the standard treatments for MIBC and are effective treatments for improving the survival rate and avoiding local recurrence as well as distant metastasis (2,3), approximately 50% of patients still have disseminated micrometastases. Therefore, the combination of systematic treatment and local treatment plays a key role in reducing the recurrence rate (2,4). Since the 1980s, cisplatin-based combination chemotherapy has been considered the first-line treatment for MIBC, but its use is limited by its toxic effects. In recent years, immunotherapy has been applied in the clinic and has shown strong antitumor activity in a variety of tumor treatments, including melanoma (5), non-small cell lung cancer (6), renal cell carcinoma (7) and so on, helping some patients with advanced cancer achieve long-term survival and even achieve the goal of clinical treatment. However, only some patients can benefit from treatment with immune checkpoint inhibitors. PD-1/PD-L1 inhibitors play an important role in regulating the function of CD8+ T cells in the tumor microenvironment (TME), and their efficacy is also influenced by the TME. Therefore, research on the TME in MIBC is attracting an increasing amount of attention.

A typical hallmark of tumors is immune dysregulation, and the interaction between ferroptosis and lipid metabolism is critical in the regulation of tumor immunity. Ferroptosis is a nonapoptotic, peroxidation-driven regulatory form of cell death that was proposed by Dixon in 2012. It is characterized by iron-dependent lipid metabolism and the accumulation of reactive oxygen species (ROS) (8). A large number of studies have confirmed that ferroptosis plays a key role in killing tumor cells and inhibiting tumor growth in cancers such as non-small cell lung cancer (9), breast cancer (10), leukemia (11), pancreatic cancer (12) and hepatocellular carcinoma (13). Therefore, inducing ferroptosis may be a new treatment strategy for cancer. However, there is no clear evidence that ferroptosis can inhibit the progression of bladder cancer. Mazdak et al. detected serum iron levels in 51 patients with bladder cancer and 58 healthy controls and found that the serum iron level in patients with bladder cancer was lower than that in healthy controls (14). In another study, the combination of gallium and transferrin resulted in an increase in free iron levels in bladder cancer cells, and a higher concentration of free iron was beneficial to inhibiting the proliferation of bladder cancer cells (15).

Different immune cells in the TME show different sensitivities to ferroptosis. Ferroptosis is a metabolic vulnerability of activated CD8+ T cells, while inhibition of ferroptosis can promote the survival and antitumor effects of CD8+ T cells in tumors (16). In tumors with different immunophenotypes, ferroptosis inducers might have distinct impacts on cancer immunity, which may determine the efficacy of immune checkpoint inhibitor immunotherapy (17). Hence, a better understanding of ferroptosis-related molecular characteristics may contribute to the progression of cancer immunotherapy research and provide a theoretical basis for clinical trials to help improve treatment effects. In this study, we aimed to comprehensively characterize the ferroptosis-related molecular subtypes of the TME and immunotherapy for MIBC via multiomics data. We present the following article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-22-1653/rc).

Methods

Ferroptosis-related genes

The “WP_FERROPTOSIS” gene sets (65 genes) were downloaded from The Molecular Signatures Database (MSigDB; https://www.gsea-msigdb.org/gsea/msigdb). The 125 ferroptosis-related genes were obtained from FerrDb database (http://www.datjar.com:40013/bt2104/), which provides ferroptosis-related gene regulators, including drivers, suppressors, and markers.

Data collection and preprocessing

Two MIBC cohorts (TCGA-BLCA, GSE13507), which had detailed follow-up information, were included in this study by searching The Cancer Genome Atlas (TCGA) database (https://tcga-data.nci.nih.gov/) and Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/). For the TCGA-BLCA cohort, high-throughput data, including RNA-sequencing data [fragments per kilobase of exon per million reads mapped (FPKM)], somatic mutation data of MIBC samples with detailed clinicopathological information (age, sex, grade, TNM staging, survival status and follow-up time) were downloaded. The inclusion criteria included the following: (I) the survival time of patients was more than 30 days; (II) the patient’s clinical information must include T (T2-T4) and N stages; and (III) the pathological type must be urothelial carcinoma. Consequently, 367 MIBC samples from TCGA were enrolled as the training set, while 61 samples from the GEO database were included as the external validation set. Patient clinical information is shown in Table S1. The batch effects between two datasets were corrected with the “ComBat” algorithm of the sva package (18).

Clustering analysis

Before clustering, univariate Cox regression analysis was used to evaluate the correlation between ferroptosis-related genes and overall survival (OS) in the TCGA cohort. Genes with P<0.05 were retained for nonnegative matrix factorization (NMF) analysis clustering conducted via the NMF package in R (19). When the cophenetic correlation coefficient declined sharply and the colors in the groups in the heatmaps were more uniform, the k-value was chosen as the optimal number of clusters. According to the clusters, Kaplan–Meier overall survival curves were drawn using the survival package in R, followed by the log-rank test. Principal components analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE) were presented to verify the classification performance on the basis of NMF analysis.

Differential expression and functional annotation analysis

Differentially expressed genes (DEGs) were filtered between two molecular subtypes via the limma package with cutoffs of false discovery rate (FDR) <0.05 and |log2 fold change |≥1. Their underlying functions were predicted through GO enrichment analysis and Gene Set Enrichment Analysis (GSEA) via the clusterProfiler package in R (20). The P value was adjusted by the Benjamini–Hochberg method. Adjusted P value <0.05 was considered significant.

Quantitative estimation of immune cell infiltration in MIBC

A group of immune cell gene markers, consisting of 782 genes, which represent 28 immune cell types related to innate and adaptive immunity, were obtained from previous study to estimate the infiltration level of different immune cell types in the TME (21). The assessed immune cell types included B cells, T cells, and natural killer (NK) cells and so on. Subsequently, the expression profiles of each sample were used to estimate the infiltration level of each immune cell type in MIBC using the single-sample gene set enrichment analysis (ssGSEA) algorithm with the R package “GSVA” (22). The differences in the immune infiltration levels between subtypes were calculated via the Wilcoxon rank-sum test and presented by heatmap.

Estimation of stromal and immune cells in MIBC

The levels of infiltrating stromal and immune cells in MIBC were estimated for each sample based on the gene expression profiles utilizing the estimation of stromal and immune cells in malignant tumour tissues using expression data (ESTIMATE) algorithm (23). By combining stromal and immune scores, ESTIMATE scores were determined. The tumor purity of each sample was then calculated according to the ESTIMATE scores.

Response to immune therapy and tumor mutation burden (TMB) between subtypes

The response to immunotherapy was assessed by the Tumor Immune Dysfunction and Exclusion (TIDE) website (http://tide.dfci.harvard.edu/login/). The TIDE algorithm was used to estimate the response of each sample to anti-PD-1/PD-L1 and anti-CTLA4 immunotherapy based on the gene expression profiles of the TCGA cohort. TMB was defined as the ratio of the total count of variants to the whole length of exons. The differences in the expression levels of TIDE scores and TMB levels were compared by the Wilcoxon rank-sum test.

Establishment of a signature based on DEGs in two molecular subtypes

Prognosis-related DEGs with P<0.05 were screened by univariate Cox regression analysis and then used for the next calculation by least absolute shrinkage and selector operation (LASSO) analysis. A prognostic signature was developed for multivariate Cox regression analysis. The risk score of each sample was calculated using the formula Risk score = ExpGene1 × CoefGene1 + ExpGene2 × CoefGene2 + … ExpGene(n) × CoefGene(n). In this equation, “ExpGene” represents gene expression, and “CoefGene” is the regression coefficient. MIBC patients from TCGA and GEO databases were divided into high- and low-risk groups in line with the cutoff value of risk scores. Kaplan-Meier curves were used to compare the differences in OS between the two groups via the survival package. Time-dependent receiver operating characteristic curves (ROCs) for one-, three- and five-year OS were conducted to assess the predictive power of the risk score using the survival ROC package. Furthermore, the stability and reliability of the signature was verified in the external validation group by the same methods.

Development and evaluation of a nomogram

With an attempt to predict the 1-, 3- and 5-year survival probability of each MIBC patient, we constructed a nomogram based on the five-gene signature and independent clinicopathological characteristics by univariate and multivariate Cox regression analysis. The performance of the nomogram was validated using time-dependent ROC curves and calibration plots.

Statistical analysis

All statistical analysis and figures exhibition were implemented using R 4.1.0 (R Foundation, Vienna, Austria). Log-rank test was performed to assess the difference between the survival curves. Variables of two subtypes or two risk groups were analyzed with Student’s t-test or Mann–Whitney U-test. A hazard ratio (HR) and a 95% confidence interval (CI) were evaluated by univariable and multivariate Cox regression models. If not specified above, P<0.05 was considered statistically significant.

Ethical consideration

The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Results

Characterization of two ferroptosis-related molecular subtypes for MIBC

A total of 190 ferroptosis-related genes were retrieved from the list “WP_FERROPOTOSIS” gene set and from FerrDb database. In the TCGA cohort, a total of 41 genes were associated with MIBC prognosis (all P<0.05). Based on the expression profiles of prognostic and ferroptosis-related genes, MIBC samples from TCGA were clustered via the NMF package. When starting from k=2, the cophenetic correlation coefficient started to decrease (Figure 1A). When k=2, the cophenetic correlation coefficient declines sharply, and the colors in the groups in the heatmap are more uniform (Figure 1B) compared to other k value (Figure 1C). Therefore, 367 MIBC samples were clustered into two molecular subtypes, C1 (n=200) and C2 (n=167). A prognostic difference based on OS (Figure 1D) was observed according to this subtype classification. Patients with the C1 subtype had a better prognosis than those with the C2 subtypes (P<0.001). Also C1 and C2 can be distinguished by dimensional reduction methods such as PCA (Figure 1E) and t-SNE (Figure 1F).

Figure 1 NMF identifies two distinct ferroptosis-related molecular subtypes for MIBC in TCGA-BLCA dataset. (A) Factorization rank for k=2–10. (B) The single-heatmap of the consensus matrix when the consensus clustering k=2. The columns and rows are sorted through hierarchical clustering according to the Euclidean distance of the average link. (C) The heatmap of the consensus matrix when the consensus clustering k=2–10. (D) Kaplan-Meier OS curves for the two clusters in TCGA cohort. The assessment of difference was achieved by log-rank test. (E) The PCA and (F) t-SNE scatter plots are in support of the classification into two MIBC molecular subtypes based on the ferroptosis-related gene expression profiles. The colors are indicative of samples from two molecular subtypes. NMF, nonnegative matrix factorization; MIBC, muscle invasive bladder cancer; TCGA-BLCA, The Cancer Genome Atlas-Bladder Cancer; OS, overall survival; PCA, principal components analysis; t-SNE, t-distributed stochastic neighbor embedding.

DEGS between subtypes and enrichment analysis

A total of 2,237 DEGs from the whole expression profile were identified between C1 and C2. Among them, there were 1,558 up- and 679 downregulated genes in C2 compared to C1 (Figure 2A,2B). The results of GO functional annotation analysis of the DEGs (Figure 2C) showed that the most significantly enriched biological processes (BPs) included external encapsulating structure organization, extracellular matrix organization, and extracellular structure organization. The most significantly enriched cellular components (CCs) included collagen-containing extracellular matrix, collagen trimer, and cornified envelope, and the molecular functions (MFs) included receptor ligand activity, signaling receptor activator activity, and cytokine activity. GSEA pathway enrichment analysis (Figure 2D) showed that the significant pathways included cytokine-cytokine receptor interactions, ECM receptor interactions, and cell adhesion molecules. Therefore, from the enrichment analysis of DEGs, we found that there may be a relationship between cytokines and ferroptosis-related subtypes. In our study, the expression of CXCL1/2/5/16 cytokines in the C1 subtype was significantly lower than that in the C2 subtype (Figure 2E-2H).

Figure 2 The difference between C1 and C2. (A) Heatmap depicting all DEGs between C1 (green) and C2 (red). (B) Volcano plots up-regulated genes (red dots) and down-regulated genes (green dots) in C1 compared to C2. (C) Bubble plots of the top ten GO function annotation analysis results including BP, CC and MF categories. (D) Top five enriched KEGG signaling pathways in C2 cluster. Differences in expression patterns of CXC cytokines between two MIBC molecular subtypes. As depicted in the box plots, the expression levels of (E) CXCL1 (F) CXCL2 (G) CXCL5 (H) CXCL16 are visualized in MIBC samples between C1 and C2. FDR, false discovery rate; DEGs, differentially expressed genes; GO, gene ontology; BP, biological processes; CC, cellular component; MF, molecular function; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Immune cell infiltration analysis between subtypes

Here, we first assessed the differential sensitivity to immunotherapy between the two ferroptosis-related molecular subtypes in MIBC. Specifically, C2 displayed higher infiltration levels of neutrophils (P<0.001) (Figure 3A), activated B cells (P<0.001) (Figure 3B), regulatory T cells (P<0.001) (Figure 3C), and macrophages (P<0.001) (Figure 3D) than C1. Additionally, C1 exhibited distinctly higher levels of CD56dim natural killer cells (P=0.045) (Figure 3E) and monocytes (P=0.029) (Figure 3F) than C2. Overall, the degree of immune infiltration of C2 was higher than that of C1 (Figure 3G).

Figure 3 Immune cells infiltration between subtypes. Patterns of (A) neutrophil, (B) activated B cell, (C) regulatory T cell, (D) macrophage, (E) CD56dim natural killer cell and (F) monocyte between subtypes. (G) The heatmap visualizing the expression patterns of 28 immune cells in the two subtypes.

Another important aspect of the TME is stromal cells, which are also believed to play an important role in tumor growth, disease progression, and drug resistance (24). The patterns of stromal scores (Figure 4A), immune scores (Figure 4B), ESTIMATE scores (Figure 4C) and tumor purity (Figure 4D) in each MIBC sample from TCGA were evaluated via the ESTIMATE algorithm. The patients in the C2 subtype had relatively high levels of stromal, immune and ESTIMATE scores in comparison to C1 (all P<0.001), while we found lower levels of tumor purity in C2 than in C1 (P<0.001).

Figure 4 ESTIMATE, TIDE and TMB. Patterns of (A) stromal cell scores, (B) immune cell scores, (C) ESTIMATE scores and (D) tumor purity between subtypes. (E,F) Box plots showing the correlation between TIDE/TMB levels and molecular subtypes. ESTIMATE, Estimation of Stromal and Immune Cells in Malignant Tumour Tissues using expression data; TIDE, Tumor Immune Dysfunction and Exclusion; TMB, tumor mutation burden.

In addition, we used the TIDE algorithm to predict the response of the C1/C2 groups to immunotherapy and found that C2 displayed higher TIDE levels than C1 (P<0.001), suggesting that C2 was less likely to respond to immunotherapy than C1 (Figure 4E). TMB represents the number of mutations per megabase (Mut/Mb) of DNA that were sequenced in a specific cancer. We know that somatic cell mutations can be expressed at the RNA or protein level, producing new antigens, protein fragments or polypeptides. These new proteins are recognized by the autoimmune system as non-self-antigens, activate T cells, and cause immune responses. Therefore, as the number of accumulated genetic variations per megabase increases, many new antigens can be created (25). However, there was no significant difference in TMB between subtypes (Figure 4F).

Development of a subtype-specific prognostic five-gene signature for MIBC

To explore the relationships between the DEGs and the prognosis of MIBC patients, univariate Cox regression analysis was performed with the 2,237 DEGs using the survival R package. A total of 641 DEGs were found to be associated with OS in MIBC patients (P<0.05). Then, 53 prognostic genes were selected from the 641 prognostic DEGs using LASSO Cox regression (Figure 5A,5B). Finally, to identify the optimal prognosis-related DEGs, we used multivariate Cox proportional hazards regression to further screen 5 pivotal genes to evaluate independent prognostic value (Table 1). Among these 5 pivotal genes, the high expression of 3 genes (SLC1A6, UPK3A and SLC19A3) is closely related to the shorter survival of patients, so their high expression may be a risk factor, while the remaining 2 genes (CCL17 and UGT2B4) are often prognostic protective factors, and their high expression is associated with longer survival. The expression level of each pivotal gene was weighted by the regression coefficient of the multivariate Cox regression model, and an individual risk scoring model was established. Through this risk scoring model, we calculated the risk score of each sample in the TCGA cohort and GEO cohort. All TCGA patients were separated into high- and low-risk groups in accordance with the cutoff values of risk scores. The higher the risk score, the greater the number of patients who died. The high-risk group had higher expression levels of 3 genes (SLC1A6, UPK3A and SLC19A3) than the low-risk group (Figure 5C). As shown in the ROC curves, the AUCs for one-, three- and five-year OS were 0.73, 0.69, and 0.71, respectively, suggesting the good performance of the risk prediction model for OS in the TCGA cohort (Figure 5D). Likewise, in the GEO cohort, the AUCs for one-, three- and five-year OS were 0.67, 0.76, and 0.72, respectively (Figure 5E). Patients with high-risk scores exhibited poorer OS (P<0.001) in the two cohorts (Figure 5F,5G).

Figure 5 Development of a prognostic five-gene signature for MIBC. (A) 10-time cross-validation for tuning parameter selection in the LASSO Cox model. (B) Plots of the LASSO coefficients. (C) The risk score rank (up), distribution of survival status (alive or dead; middle) and expression patterns of five genes in high- and low-risk groups. (D) Time-dependent ROC curves of five-gene signature for one-, three- and five-year OS time in TCGA cohort. (E) Time-dependent ROC curves of five-gene signature for one-, three- and five-year OS time in GEO cohort. (F,G) Kaplan-Meier OS curve for high- and low-risk groups in TCGA cohort and GEO cohort. MIBC, muscle invasive bladder cancer; LASSO, Least Absolute Shrinkage and Selector Operation; ROC, receiver operating characteristic; GEO, Gene Expression Omnibus.

A nomogram integrating the five-gene signature and clinical factors improves predictive power for MIBC prognosis

Considering that the sample size of N3 in the N stage is too small, forcible modeling will reduce the accuracy of the model. Therefore, we constructed a nomogram by combining the five-gene signature and clinical factors, including age, T stage (T2-4), and N stage (N0-2), for predicting MIBC patient OS (Figure 6A). From the nomogram, by adding up the corresponding scores for each independent risk factor, we obtain a total score that yields the corresponding predicted patient survival rates at 1, 3, and 5 years. We further evaluated whether the integration of the five-gene signature and clinical factors could boost the predictive efficiency for MIBC prognosis in the TCGA cohort. The AUCs for one-, three- and five-year OS were 0.77, 0.71, and 0.72, respectively, suggesting the good performance of the nomogram (Figure 6B). We can use the information of each patient to infer its prognosis based on the nomogram (Figure 6C). Calibration plots confirmed that the nomogram-predicted probabilities of one- (Figure 6D), three- (Figure 6E) and five-year (Figure 6F) OS had high consistency with the actual survival. Collectively, the nomogram integrating the five-gene signature, age, T stage, and N stage could enhance the predictive power of MIBC patient prognosis.

Figure 6 A Nomogram Integrating the five-gene signature and clinical factors improves predictive power for MIBC Prognosis. (A) Forest plots showing the multivariate Cox regression analyses results of the five-gene signature and clinical factors with OS. (B) Time dependent ROC curves of five-gene signature and clinical factors for one-, three- and five-year OS time in TCGA cohort. (C) A nomogram incorporating five-gene signature and clinical factors improves predictive efficacy for MIBC prognosis. (D-F) Calibration plots displayed the actual and nomogram-predicted probability of one-, three- and five-year OS in TCGA cohort. ROC, receiver operator characteristic; AUC, area under the curves; MIBC, muscle invasive bladder cancer; OS, overall survival; TCGA, The Cancer Genome Atlas.

Discussion

Although RC is the main treatment for MIBC at present, Herr et al. showed that for patients treated with RC alone, the postoperative recurrence rate of T1 and T2 stage is 20–30% and that of T3 and T4 stage is 50–90% (26). For patients with high risk or recurrence, adjuvant chemotherapy or radiation therapy is required after surgery. However, its practical clinical application is somewhat limited considering the adverse effects of chemotherapy and radiotherapy. The presence of immunological checkpoint inhibitors has radically changed the therapeutic pattern of bladder cancer and expanded the therapeutic strategy (27). Compared with traditional chemotherapy and radiotherapy, immunotherapy has the advantages of good targeting, fewer adverse reactions and better curative effects. Currently, major authoritative national and international guidelines recommend immune checkpoint inhibitors as second-line treatment for patients with bladder cancer who have lost the opportunity for resection and metastases and as first-line treatment for PD-L1-positive patients who are not suitable for platinum-based chemotherapy (28,29).

Notably, ferroptosis is an iron-dependent form of cell death characterized by an accumulation of lipid peroxides, playing a significant role in mediating immune cells and immunotherapies. Recent study has shown that the combination of immunotherapy and ferroptosis inducers can lead to a strong immune response, which in turn promotes ferroptosis in tumor cells to enhance the antitumor effect. In tumor patients treated with combination therapy, the level of ferroptosis in tumor cells correlated significantly with the prognosis of survival and treatment efficacy (30). Therefore, ferroptosis is closely related to immunotherapy of tumors. Considering that only 30% of tumor patients currently respond to immunotherapy, it is necessary for us to identify sensitive patients for precise treatment. In this study, we attempted to differentiate the sensitivity of MIBC patients to immunotherapy by different subtypes of ferroptosis-related molecules.

First, based on 41 ferroptosis-related genes associated with prognosis, we identified two ferroptosis-associated molecular subtypes in the TCGA cohort by NMF analysis: C1 and C2. Furthermore, two data dimensionality reduction algorithms, PCA and t-SNE, were used to visualize the differentiation of the two subtypes. Additionally, we further performed prognostic survival analysis for both subtypes, and the results suggested that patients in the C2 subtype group had a poorer prognosis (P<0.001). What exactly are the differences between the two subtypes?

To further investigate the two ferroptosis-associated molecular subtypes in MIBC, we explored the differences from the perspective of gene expression profiles. Genes between the two subtypes were subjected to differential analysis and pathway enrichment analysis, and the results suggested that DEGs between the two subtypes were associated with cytokines, especially chemokines that play important roles in angiogenesis (31), hematopoiesis (32), tumor cell metastasis (24) and so on. Among these cytokines, we found that some of them were significantly associated with prognosis. Zangouei et al. previously found that CXCL1/2/5/16 cytokines were associated with prognosis or infiltration of bladder cancer (33). With increased levels of CXCL1, patients face poor survival and advanced grade and stage (34). CXCL2/16 upregulation was associated with advanced stage and poor survival in bladder cancer patients (35,36). There was a direct correlation between CXCL5 downregulation and decreased bladder tumor cell migration and invasion, which was due to the reduced and increased expression of Snail and E-cadherin, respectively (37). Therefore, the expression of these cytokines in the C1 subtype was significantly lower than that in the C2 subtype, suggesting that the C2 had a poor prognosis in another aspect (Figure 2E-2H).

It is well known that the TME is closely related to tumor prognosis, and the sensitivity of patients to immunotherapy depends largely on the TME. In addition to tumor cells, the TME is also associated with normal stromal cells, immune cells, vascular endothelial cells and blood cells in blood vessels. The most important are immune cells, different types of which play different roles in antitumor and tumor immune escape processes. Therefore, to explore the abundance of immune infiltration in the two subtypes, we utilized the ssGSEA algorithm and found a higher abundance of immune cell infiltration in the C2 subtype, mainly in neutrophils, activated B cells, regulatory T cells, and macrophages. Recent study targeting the TME have found that several immune-infiltrating cells exhibit important immunosuppressive functions, including myeloid suppressor cells (MDSCs) (38), tumor-associated macrophages (TAMs), regulatory T cells (Tregs) (39), and tumors associated neutrophils (TANs) (40). Some TAMs show M1-like macrophage phenotype, especially in the early stage of tumorigenesis, expressing anti-tumor effect. By reprogramming TAMs, the vast majority of TAMs show an M2-like macrophage phenotype, which plays an important role in promoting tumorigenesis, progression and metastasis (41). In this study, all four types of cells were highly infiltrated in C2 cluster, which may affect its prognosis. To further explore the relationship between ferroptosis and immune cell infiltration in this study, we used correlation analysis to study the relationship between ferroptosis prognosis-related genes and T cells in this study, and expressed it in the form of heatmap (Figure S1). We found that IFNG, ZEB, SLC39A14, and AGPS were all positively correlated with the infiltration degree of activated CD8 T cells, and all were ferroptosis-promoting genes. The two genes, PLA2G6 and SRC, showed a negative correlation with activated CD8 cells, and they were both ferroptosis suppressor genes. Another important aspect of the TME is stromal cells, which are also thought to play an important role in tumor growth, disease progression, and drug resistance (24). Similarly, using the Estimate algorithm, we found higher stromal cell scores and lower tumor purity for the C2 subtype, which matches the immune infiltration results we obtained previously using the ssGSEA algorithm.

Since the two subtypes have different immune microenvironments, are there differences in their effects on immunotherapy? Based on the TIDE algorithm, immune checkpoint blockade (ICB) was found to have a low therapeutic potential for C2. Through the above study, we found that patients with the C2 subtype of MIBC have a higher degree of tumor immune infiltration, lower tumor purity, and are less sensitive to immunotherapy, which reminds us of the need for alternative treatment options for this subtype of patients.

However, how can MIBC patients belong to the C1 or C2 subtype be distinguished? We previously performed NMF analysis using 48 ferroptosis-associated prognostic genes to derive two subtypes. Then, these 48 genes were analyzed for differences between subtypes, and genes with P values less than 0.05 and large differences (|log2 fold change>2|) were identified. SLC7A11, SLC16A1, and CAV1 ranked in the top three highly expressed genes in the C2 subtype. The cystine/glutamate antiporter SLC7A11 functions to import cystine for glutathione biosynthesis and antioxidant defense, and its overexpression promotes tumor growth partly by suppressing ferroptosis (42). Li et al. found that silencing lncRNA SLC16A1 can induce ferroptosis through miR-1433p/SLC7A11 signaling in renal cancer (43). In addition, overexpression of CAV1 in HNSCC inhibited the process of ferroptosis, leading to aggressive phenotypes and worse prognosis (44). SLC7A11, SLC16A1 and CAV1 are ferroptosis inhibitor-related genes and are closely related to the occurrence, development and prognosis of various malignant tumors, while there is a lack of research on these three genes in MIBC, which provides a new direction for subsequent basic research on MIBC. Whether these three genes can be used to distinguish MIBC as the C1 or C2 subtype will depend on more studies in the future.

Currently, it is a hot topic to predict the prognosis of tumor patients by signature, and it has great significance for clinical treatment and basic research guidance. This study constructed a subtype-specific signature of 5 pivotal genes (SLC1A6, UPK3A, SLC19A3, CCL17 and UGT2B4) from DEGs between different subtypes after using univariate Cox regression analysis, LASSO Cox regression analysis and multivariate Cox regression analysis. Among 5 pivotal genes, SLC1A6, UPK3A, and SLC19A3 were risk factors for the prognosis of MIBC in this study, while CCL17 and UGT2B4 were protective factors. Through the difference analysis between C1 and C2, CCL17 and UGT2B4 were relatively low expressed in C2. Li et al. found that bladder cancer patients with high expression levels of CCL17 were associated with a significantly better prognosis (45), which is consistent with our study. Zeng et al. found that high expression of UGT2B4 is associated with low-grade prostate cancer, and UGT2B4 upregulation in tumors is associated with upregulation of metabolic pathways, such as novel IMP biosynthetic processes, glutamine and monocarboxylic acid metabolism (46). UP3KA is not only a prognostic risk factor, but Lai et al. found that the measurement of UPK3A in urine is a sensitive new marker with good performance for detecting bladder cancer (47). SLC1A6, a member of the SLC1A family, encodes a transmembrane transporter that mediates L-glutamate and L/D-aspartate uptake, and its overexpression reduces the response of nasopharyngeal carcinoma radioresistant cells to cisplatin and radiation sensitivity. Zou et al. concluded that the expression level of SLC1A6 was an independent predictor of poor prognosis in bladder cancer patients through multivariate COX regression analysis (48). SLC19A3, which was relatively high expression in C2, encodes hTHTR2 and is mainly expressed in the intestine, liver, kidney and placenta. Mutations in it cause alkaloid-responsive basal ganglia disease and thiamine-responsive encephalopathy (49). However, research on UGT2B4 and SLC19A3 in bladder cancer is still lacking. Following validation, this signature could robustly and independently predict the OS of MIBC patients. To more accurately predict the prognosis of MIBC patients, we also constructed a nomogram combining the signature and other clinical factors. The prediction system can guide the establishment of personalized examination procedures for MIBC patients and boost the effective use of medical resources.

There are also some limitations in our study. First, the multigene prognostic models were constructed and validated on public databases, which need further bioinformatic and experimental verification of our own data. Second, we identified the difference between the two subtypes, but the underlying mechanisms and new markers for distinct subtypes still warrant further investigation.

Conclusions

In this study, we obtained two ferroptosis-related molecular subtypes using NMF analysis of ferroptosis-related genes and found that the C2 subtype had more abundant immune infiltration, lower tumor purity, higher cytokine expression, less sensitivity to immmunotherapy and poorer prognosis. Following the prognostic correlation analysis by DEGs of the two subtypes, five genes, SLC1A6, UPK3A, SLC19A3, and CCL17 UGT2B4, were constructed for prognostic models that had better predictive performance and credibility both at the genetic level alone and in combination with clinical data.

Acknowledgments

The authors gratefully acknowledge the data generated by TCGA database and GEO database used in this study.

Funding: None.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-22-1653/rc

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-22-1653/prf

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-22-1653/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Antoni S, Ferlay J, Soerjomataram I, et al. Bladder Cancer Incidence and Mortality: A Global Overview and Recent Trends. Eur Urol 2017;71:96-108. [Crossref] [PubMed]
2. Stein JP, Lieskovsky G, Cote R, et al. Radical cystectomy in the treatment of invasive bladder cancer: long-term results in 1,054 patients. J Clin Oncol 2001;19:666-75. [Crossref] [PubMed]
3. Stein JP, Quek ML, Skinner DG. Lymphadenectomy for invasive bladder cancer: I. historical perspective and contemporary rationale. BJU Int 2006;97:227-31. [Crossref] [PubMed]
4. Yafi FA, Aprikian AG, Chin JL, et al. Contemporary outcomes of 2287 patients with bladder cancer who were treated with radical cystectomy: a Canadian multicentre experience. BJU Int 2011;108:539-45. [Crossref] [PubMed]
5. Larkin J, Minor D, D'Angelo S, et al. Overall Survival in Patients With Advanced Melanoma Who Received Nivolumab Versus Investigator's Choice Chemotherapy in CheckMate 037: A Randomized, Controlled, Open-Label Phase III Trial. J Clin Oncol 2018;36:383-90. [Crossref] [PubMed]
6. Broderick SR. Adjuvant and Neoadjuvant Immunotherapy in Non-small Cell Lung Cancer. Thorac Surg Clin 2020;30:215-20. [Crossref] [PubMed]
7. Motzer RJ, Escudier B, McDermott DF, et al. Nivolumab versus Everolimus in Advanced Renal-Cell Carcinoma. N Engl J Med 2015;373:1803-13. [Crossref] [PubMed]
8. Dixon SJ, Stockwell BR. The role of iron and reactive oxygen species in cell death. Nat Chem Biol 2014;10:9-17. [Crossref] [PubMed]
9. Dixon SJ, Lemberg KM, Lamprecht MR, et al. Ferroptosis: an iron-dependent form of nonapoptotic cell death. Cell 2012;149:1060-72. [Crossref] [PubMed]
10. Ma S, Henson ES, Chen Y, et al. Ferroptosis is induced following siramesine and lapatinib treatment of breast cancer cells. Cell Death Dis 2016;7:e2307. [Crossref] [PubMed]
11. Trujillo-Alonso V, Pratt EC, Zong H, et al. FDA-approved ferumoxytol displays anti-leukaemia efficacy against cells with low ferroportin levels. Nat Nanotechnol 2019;14:616-22. [Crossref] [PubMed]
12. Yamaguchi Y, Kasukabe T, Kumakura S. Piperlongumine rapidly induces the death of human pancreatic cancer cells mainly through the induction of ferroptosis. Int J Oncol 2018;52:1011-22. [Crossref] [PubMed]
13. Sun X, Ou Z, Chen R, et al. Activation of the p62-Keap1-NRF2 pathway protects against ferroptosis in hepatocellular carcinoma cells. Hepatology 2016;63:173-84. [Crossref] [PubMed]
14. Mazdak H, Yazdekhasti F, Movahedian A, et al. The comparative study of serum iron, copper, and zinc levels between bladder cancer patients and a control group. Int Urol Nephrol 2010;42:89-93. [Crossref] [PubMed]
15. Martin-Sanchez D, Fontecha-Barriuso M, Sanchez-Niño MD, et al. Cell death-based approaches in treatment of the urinary tract-associated diseases: a fight for survival in the killing fields. Cell Death Dis 2018;9:118. [Crossref] [PubMed]
16. Tang R, Xu J, Zhang B, et al. Ferroptosis, necroptosis, and pyroptosis in anticancer immunity. J Hematol Oncol 2020;13:110. [Crossref] [PubMed]
17. Xu H, Ye D, Ren M, et al. Ferroptosis in the tumor microenvironment: perspectives for immunotherapy. Trends Mol Med 2021;27:856-67. [Crossref] [PubMed]
18. Leek JT, Johnson WE, Parker HS, et al. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics 2012;28:882-3. [Crossref] [PubMed]
19. Gaujoux R, Seoighe C. A flexible R package for nonnegative matrix factorization. BMC Bioinformatics 2010;11:367. [Crossref] [PubMed]
20. Yu G, Wang LG, Han Y, et al. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 2012;16:284-7. [Crossref] [PubMed]
21. Charoentong P, Finotello F, Angelova M, et al. Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cell Rep 2017;18:248-62. [Crossref] [PubMed]
22. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 2013;14:7. [Crossref] [PubMed]
23. Yoshihara K, Shahmoradgoli M, Martínez E, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun 2013;4:2612. [Crossref] [PubMed]
24. Quail DF, Joyce JA. Microenvironmental regulation of tumor progression and metastasis. Nat Med 2013;19:1423-37. [Crossref] [PubMed]
25. Addeo A, Friedlaender A, Banna GL, et al. TMB or not TMB as a biomarker: That is the question. Crit Rev Oncol Hematol 2021;163:103374. [Crossref] [PubMed]
26. Herr HW, Dotan Z, Donat SM, et al. Defining optimal therapy for muscle invasive bladder cancer. J Urol 2007;177:437-43. [Crossref] [PubMed]
27. Song YP, McWilliam A, Hoskin PJ, et al. Organ preservation in bladder cancer: an opportunity for truly personalized treatment. Nat Rev Urol 2019;16:511-22. [Crossref] [PubMed]
28. Alfred Witjes J, Lebret T, Compérat EM, et al. Updated 2016 EAU Guidelines on Muscle-invasive and Metastatic Bladder Cancer. Eur Urol 2017;71:462-75. [Crossref] [PubMed]
29. Flaig TW, Spiess PE, Agarwal N, et al. Bladder Cancer, Version 3.2020, NCCN Clinical Practice Guidelines in Oncology. J Natl Compr Canc Netw 2020;18:329-54. [Crossref] [PubMed]
30. Wang W, Green M, Choi JE, et al. CD8+ T cells regulate tumour ferroptosis during cancer immunotherapy. Nature 2019;569:270-4. [Crossref] [PubMed]
31. Strieter RM, Burdick MD, Gomperts BN, et al. CXC chemokines in angiogenesis. Cytokine Growth Factor Rev 2005;16:593-609. [Crossref] [PubMed]
32. Youn BS, Mantel C, Broxmeyer HE. Chemokines, chemokine receptors and hematopoiesis. Immunol Rev 2000;177:150-74. [Crossref] [PubMed]
33. Zangouei AS, Hamidi AA, Rahimi HR, et al. Chemokines as the critical factors during bladder cancer progression: an overview. Int Rev Immunol 2021;40:344-58. [Crossref] [PubMed]
34. Miyake M, Lawton A, Goodison S, et al. Chemokine (C-X-C) ligand 1 (CXCL1) protein expression is increased in aggressive bladder cancers. BMC Cancer 2013;13:322. [Crossref] [PubMed]
35. Lee JT, Lee SD, Lee JZ, et al. Expression analysis and clinical significance of CXCL16/CXCR6 in patients with bladder cancer. Oncol Lett 2013;5:229-35. [Crossref] [PubMed]
36. Zhang H, Ye YL, Li MX, et al. CXCL2/MIF-CXCR2 signaling promotes the recruitment of myeloid-derived suppressor cells and is correlated with prognosis in bladder cancer. Oncogene 2017;36:2095-104. [Crossref] [PubMed]
37. Zheng J, Zhu X, Zhang J. CXCL5 knockdown expression inhibits human bladder cancer T24 cells proliferation and migration. Biochem Biophys Res Commun 2014;446:18-24. [Crossref] [PubMed]
38. Kumar V, Patel S, Tcyganov E, et al. The Nature of Myeloid-Derived Suppressor Cells in the Tumor Microenvironment. Trends Immunol 2016;37:208-20. [Crossref] [PubMed]
39. Corthay A. How do regulatory T cells work? Scand J Immunol 2009;70:326-36. [Crossref] [PubMed]
40. Shaul ME, Fridlender ZG. Tumour-associated neutrophils in patients with cancer. Nat Rev Clin Oncol 2019;16:601-20. [Crossref] [PubMed]
41. Salmaninejad A, Valilou SF, Soltani A, et al. Tumor-associated macrophages: role in cancer development and therapeutic implications. Cell Oncol (Dordr) 2019;42:591-608. [Crossref] [PubMed]
42. Koppula P, Zhuang L, Gan B. Cystine transporter SLC7A11/xCT in cancer: ferroptosis, nutrient dependency, and cancer therapy. Protein Cell 2021;12:599-620. [Crossref] [PubMed]
43. Li YZ, Zhu HC, Du Y, et al. Silencing lncRNA SLC16A1-AS1 Induced Ferroptosis in Renal Cell Carcinoma Through miR-143-3p/SLC7A11 Signaling. Technol Cancer Res Treat 2022;21:15330338221077803. [Crossref] [PubMed]
44. Lu T, Zhang Z, Pan X, et al. Caveolin-1 promotes cancer progression via inhibiting ferroptosis in head and neck squamous cell carcinoma. J Oral Pathol Med 2022;51:52-62. [Crossref] [PubMed]
45. Li Y, Chen X, Li D, et al. Identification of prognostic and therapeutic value of CC chemokines in Urothelial bladder cancer: evidence from comprehensive bioinformatic analysis. BMC Urol 2021;21:173. [Crossref] [PubMed]
46. Zeng T, Fedeli MA, Tanda F, et al. Whole-exome Sequencing of Prostate Cancer in Sardinian Identify Recurrent UDP-glucuronosyltransferase Amplifications. J Cancer 2021;12:438-50. [Crossref] [PubMed]
47. Lai Y, Ye J, Chen J, et al. UPK3A: a promising novel urinary marker for the detection of bladder cancer. Urology 2010;76:514.e6-11. [Crossref] [PubMed]
48. Zou X, Wei Y, Qi T, et al. A novel 6-gene signature derived from tumor-infiltrating T cells and neutrophils predicts survival of bladder urothelial carcinoma. Aging 2021;13:25496-517. [Crossref] [PubMed]
49. Porta M, Toppila I, Sandholm N, et al. Variation in SLC19A3 and Protection From Microvascular Damage in Type 1 Diabetes. Diabetes 2016;65:1022-30. [Crossref] [PubMed]