A Novel TGF-β Risk Score Predicts the Clinical Outcomes and Tumour Microenvironment Phenotypes in Bladder Cancer

Background The TGF-β pathway plays critical roles in numerous malignancies. Nevertheless, its potential role in prognosis prediction and regulating tumour microenvironment (TME) characteristics require further elucidation in bladder cancer (BLCA). Methods TGF-β-related genes were comprehensively summarized from several databases. The TCGA-BLCA cohort (training cohort) was downloaded from the Cancer Genome Atlas, and the independent validation cohorts were gathered from Xiangya Hospital (Xinagya cohort) and Gene Expression Omnibus. Initially, we identified differentially expressed TGF-β genes (DEGs) between cancer and normal tissues. Subsequently, univariate Cox analysis was applied to identify prognostic DEGs, which were further used to develop the TGF-β risk score by performing LASSO and multivariate Cox analyses. Then, we studied the role of the TGF-β risk score in predicting prognosis and the TME phenotypes. In addition, the role of the TGF-β risk score in guiding precision treatments for BLCA has also been assessed. Results We successfully constructed a TGF-β risk score with an independent prognostic prediction value. A high TGF-β risk score indicated an inflamed TME, which was supported by the positive relationships between the risk score, enrichment scores of anticancer immunity steps, and the infiltration levels of tumour-infiltrating immune cells. In addition, the risk score positively correlated with the expression of several immune checkpoints and the T cell inflamed score. Consistently, the risk score was positively related to the enrichment scores of most immunotherapy-positive pathways. In addition, the sensitivities of six common chemotherapeutic drugs were positively associated with the risk score. Furthermore, higher risk score indicated higher sensitivity to radiotherapy and EGFR-targeted therapy. On the contrary, patients with low-risk scores were more sensitive to targeted therapies, including the blockade of FGFR3 and WNT-β-catenin networks. Conclusions We first constructed and validated a TGF-β signature that could predict the prognosis and TME phenotypes for BLCA. More importantly, the TGF-β risk score could aid in individual precision treatment for BLCA.


INTRODUCTION
Bladder cancer (BLCA) is the second frequently diagnosed urinary cancers with high heterogeneity, causing more than 17000 deaths in the US in the past year alone (1). A variety of treatment options, including surgery, neoadjuvant or adjuvant chemotherapy, immunotherapy, targeted therapy, and radiotherapy, provide competent tools for urologists to cure patients with BLCA (2). Nonetheless, most BLCA patients are not sensitive to these current treatments, which leads to high overall mortality (3)(4)(5)(6). Hence, it is imperative to explore new prognostic biomarkers and treatment response predictors, which promoting individual precision treatment for BLCA.
Currently, the revolution of immune checkpoint blockades (ICB), including anti-PD-L1 or anti-PD-1 therapies, brings hope to patients with advanced BLCA (7). The response rate of ICB mainly depends on the tumour microenvironment (TME) (8). The TME is a comprehensive system consisting of tumour cells, tumour-infiltrating immune cells (TIICs), and an extracellular matrix. The categories and infiltration levels of TIICs, including CD8 T cells and Tregs, influence the activities of the anticancer immune response (9)(10)(11). The status and components of the TME are regulated by several mechanisms. One of the most important of these mechanisms is TGF-b (12)(13)(14)(15)(16).
The TGF-b pathway significantly regulates cancer biology, especially in anticancer immune responses (17). On the one hand, TGF-b suppresses tumour progression and strengthens TME homeostasis by inducing cancer cell apoptosis (18,19). On the other hand, autocrine and paracrine TGF-b stimulates EMT, commonly mediated by the transcription factors SNAIL and SLUG, leading to resistance to anticancer treatments (17)(18)(19). TGF-b promotes immune evasion of cancer cells by suppressing the proliferation, differentiation, and immunocompetence of many immune cells, such as cytotoxic cells, DCs, and NK cells (16,(20)(21)(22). Consistently, activation of the TGF-b pathway indicated poor survival and resulted in resistance to ICB. Despite the prominent role of TGF-b in regulating multiple cancer processes, therapies targeting the TGF-b pathway have not been well explored (10). With the improvement of genomic sequencing methods, more and more genomic signatures have been created to predict the prognosis and therapeutic opportunities for cancers (8,23,24). In BLCA, Stojnev et al. explored the roles of three key components (TGF-b1, Smad2, and Smad4) of the canonical TGFb pathway in predicting the prognosis (25). However, few studies have been conducted to systematically correlate the TGF-b signature with the tumour microenvironment phenotypes for BLCA.
This study, first established a novel TGF-b signature by integrating several independent BLCA cohorts, including the TCGA-BLCA and Xiangya cohorts. The TGF-b signature was closely related to TME characteristics and could predict several therapeutic opportunities for BLCA.

Training Cohort
We downloaded TCGA-BLCA mRNA sequencing data (FPKM) and clinical data from The Cancer Genome Atlas (TCGA). The FPKM values were then switched into TPM values. The TCGA-BLCA cohort included 403 BLCA samples and 19 paired normal tissues.

Validation Cohorts
Xiangya cohort (internal validation set) was collected in our hospital. The Xinagya cohort included 57 BLCA samples sequenced by the BGISEQ-500 platform (BGI-Shenzhen, China) (GSE188715) (8,23). The TPM values of the mRNA expression matrix were used to do the analysis. Another two external GEO cohorts, GSE32894 (224 samples) and GSE48075 (73 samples), with detailed survival information, were gathered from the Gene Expression Omnibus (GEO). The original RNA expression matrices of these two GEO cohorts were directly downloaded and eligible to be directly analysed. Three immunotherapy associated cohorts were collected, including IMvigor210 cohort, GSE78220, and GSE100797.
Clinicopathologic information of these data sets is provided in Supplementary Table 1.

Screening Differentially Expressed TGF-b Genes
We identified TGF-b DEGs between bladder cancer and normal tissues by using the empirical Bayesian approach of the limma R package. |log 2 (fold change) | > 1 and the adjusted P-value < 0.05 were defined as standardization to determine significant TGF-b DEGs (26). Then, we used the ClusterProfiler R package to perform Gene Ontology (GO) analyses and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses based on the TGF-b DEGs (26). Furthermore, the protein-protein interaction network of those TGF-b DEGs was plotted by using the STRING database.

Development and Validation of the TGF-b Risk Score
First, the prognostic TGF-b DEGs were screened by using univariate Cox analysis and survival R package in the TCGA-BLCA cohort. Second, we performed the least absolute shrinkage and selector operation (LASSO) regression on the prognostic TGF-b DEGs to reduce the dimensionality of high-dimensional data in the TCGA-BLCA cohort. The variables with relatively small contributions to the outcomes were given zero coefficients in the LASSO regression analysis. Then, the genes with nonzero coefficients were ultimately selected for multivariable Cox regression analysis and further construct the TGF-b risk score. The computation formula of the risk score is as follow: Risk score = S bi * RNAi, where bi is the coefficient of the i-th gene in multivariable Cox regression analysis.
Then, all patients were classified into two groups (high-risk and low-risk score groups) based on the median of the TGF-b risk scores (23,27,28). To illustrate the survival difference between the high-risk and low-risk groups, Kaplan-Meier survival analysis was used to generate the survival curves. The prognostic significance of the TGF-b risk score was estimated by using the log-rank test. In addition, we plotted receiver operating characteristic (ROC) curves by using the 'survival ROC'R package. We calculated the area under the curve (AUC) to assess the accuracy of the TGF-b risk score in predicting prognosis (26). We subsequently correlated the TGF-b risk score with the tumour stage and grade. Finally, the role of the TGF-b risk score was further validated in several independent BLCA validation sets, including the GSE32894, GSE48075, and Xiangya cohorts.

Development and Validation of a Comprehensive Nomogram
We performed univariate and multivariate Cox regression analyses to screen the independent prognostic characters based on clinicopathologic characters and TGF-b risk score. Only 19 patients diagnosed with high tumour grade existed in the TCGA-BLCA cohort. As a result, the tumour grade was not included in the Cox survival analysis. Furthermore, we integrated the prognostic factors in univariate Cox analysis to construct a comprehensive nomogram to predict the prognosis of BLCA patients. The statistical performance of the comprehensive nomogram was evaluated by plotting ROC and calibration curves in the TCGA-BLCA cohort and several external validation cohorts.
Defining the Immunological Characters in the BLCA TME The anticancer immune response in the BLCA TME comprises several steps: tumour cells releasing cancer cell antigens (Step 1); tumour cells presenting cancer antigens to immune cells in the TME (such as dendritic cells) (Step 2); antigen-presenting cells (APCs) carrying antigens priming and activating the immune system in the TME (Step 3); the initiation of the immune system results in the release of related chemokines and cytokines, hence recruiting immune cells infiltrating the TME (Steps 4-5); and cytotoxic immune cells, such as natural killer T cells and CD8+ T cells, existing in the TME recognizing cancer cells (Step 6) and killing them (Step 7) (29). The vitality of seven steps decides on the fate of the tumour cells (8,23). In addition, we calculated the levels of tumour-infiltrating immune cells (TIICs) by adopting several independent methods, including TIMER, Quan TIseq, TIP, XCell, Cibersort-ABS, and MCP-Counter, based on the RNA-seq data (30)(31)(32)(33)(34)(35)(36). Pre-existing anticancer immunity in the TME can be reflected by the T cell inflamed score (TIS), which can predict the clinical response to immune checkpoint block (ICB) (37). Then, we used ssGSEA to calculate the enrichment scores of several immunotherapy response-related pathways (8,38). Finally, we screened and collected 20 inhibitory immune checkpoints, including PD-1, PD-L1, and CTLA-4. These immunological features in the TME have been well described in our previous studies (8,23).

Predicting the Molecular Subtypes in BLCA
BLCA is a highly heterogeneous tumour. The confirmation of the individual molecular subtype contributes to personalized precision medicine for BLCA patients. Thus, we have lucubrated molecular subtype systems in our previous study, including CIT, Lund, MDA, TCGA, Baylor, UNC, and Consensus subtypes (8,(39)(40)(41)(42)(43). We adopted ConsensusMIBC and BLCAsubtyping R packages to define individual subtypes. In addition, we collected 12 heterogeneous signatures of BLCA. We further correlated the TGF-b risk score with these molecular subtypes and BLCA-specific signatures. After all samples were preliminarily assigned into basal or luminal subtypes, we plotted ROC curves to evaluate the predictive accuracy of the TGF-b signature in evaluating the molecular subtypes.

Prediction of the Clinical Sensitivities to Several Treatment Options for BLCA
Since chemotherapy is critical for terminal BLCA patients, we used the pRRophetic package to evaluate the clinical response to six common chemotherapy drugs (docetaxel, bleomycin, paclitaxel, camptothecin, vinblastine, and cisplatin) based on data of the Genomics of Drug Sensitivity in Cancer (GDSC) (https://www. cancerrxgene.org/) (44). We compared the difference in the 50% inhibitory concentration (IC50) of the six chemotherapeutics mentioned above between the high-risk and low-risk score groups. Targeted therapies, radiotherapy, and optional therapeutic schedules are also crucial. Therefore, several potential predictors related to the response of radiotherapy and targeted therapy were gathered from our previous studies (8,23).

Statistical Analysis
Data visualization and statistical analyses were conducted by using R software (Version: 4.0.5). We adopted Pearson or Spearman correlation analyses to ascertain the relationships between continuous variables. As for continuous variables fitting normal distribution between binary groups, t test was used. Otherwise, we used the Mann-Whitney U test. Univariate and multivariate Cox regression analysis were used to evaluating the correlation between the risk score and prognosis. We calculated the accuracy of the TGF-b risk score in predicting prognosis and molecular subtypes by plotting the ROC curves. Analyses with two-sided P less than 0.05 were considered statistically significant.

Screening of TGF-b DEGs and Corresponding Functional Enrichment Analysis
Sixty TGF-b DEGs were screened between BLCA and normal tissues (Supplementary Table 2 Figure 1C). For the GO_CC analysis, these TGF-b DEGs were enriched in transcription regulator complex, RNA polymerase II transcription regulator complex, and phosphatase complex. For the GO_MF analysis, DEGs were enriched in transforming growth factor-beta binding, SMAD binding, and cytokine binding. Moreover, KEGG analyses demonstrated that TGF-b DEGs were also enriched in the TGF-b signalling pathway, Hippo signalling pathway, and bladder cancer (Supplementary Figure 1D). As expected, the TGF-b signalling pathway was distinctly the most common enrichment pathway of these DEGs. The PPI network analysis results indicated that these TGF-b DEGs were closely related to each other (Supplementary Figure 1E).

Development and Validation of TGF-b Risk
Score for BLCA First, we identified 12 prognostic genes from the 60 TGF-b DEGs in TCGA-BLCA cohort (Supplementary Table 2). Second, we performed a LASSO regression based on these 12 prognostic genes and screened eight genes with nonzero coefficients ( Figures 1A,  B). The minimum lamda was 0.019. Eventually, we identified five optimal TGF-b genes (MAP2K1, PPP2CB, LRRC32, ID2, and FNTA) based on the Akaike information criterion (AIC). The results of univariate Cox regression analyses based on these five optimized genes included in the risk score model is shown in Supplementary Table 3A. We then calculated the TGF-b risk score by performing multivariate Cox regression on these five genes. Figure 1C and Supplementary Table 3B displayed the coefficients of these five genes. In TCGA-BLCA cohort, patients were assigned into high-risk and low-risk score groups based on the median risk score. Patients with high-risk score had poorer OS than patients with low-risk score ( Figure 1D). Meanwhile, the accuracy of the TGF-b risk score in predicting OS was acceptable ( Figure 1E). To adjust the bias caused by clinicopathological heterogeneity, we performed further subgroup analyses based on different clinicopathological features. As shown in Supplementary Figure 2, patients with high-risk score had poorer prognosis regardless of the subgroups. In addition, we successfully validated the role of the TGF-b risk score in three independent validation cohorts, including the Xiangya cohort, GSE32894, and GSE48075 ( Figures 1F-K). In these validation cohorts, patients with high-risk score had poorer OS. Especially in the Xiangya cohort, the accuracy of the TGF-b risk score in predicting the three-year OS was more than 0.9 ( Figure 1G).

Development and Validation of a Comprehensive Nomogram Based on Clinicopathological Characteristics and TGF-b Risk Score
We analysed the correlations between TGF-b risk score and tumour grade, tumour stage. As expected, patients with advanced grades and stages had higher risk scores, which was in accord with the adverse prognostic role of the TGF-b risk score (Figures 2A, B). We successfully validated these results in several external cohorts (Supplementary Figure 3). In addition, we analysed the relationships between the risk score and histological subtypes, muscle invasiveness, and metastasis (Supplementary Figure 4). We found a higher risk score in cancers with nonpapillary histology and muscle invasiveness (Supplementary Figures 4A, B, D, F). Patients with metastasis had higher risk score compared with patients without metastasis (Supplementary Figures 4C, E). The significant difference was not significant, which may be attributed to the small number of patients with metastasis. The univariate Cox analysis results based on several common clinicopathological characteristics showed that tumour stage, age, and TGF-b risk score were prognostic factors ( Figure 2C). We subsequently performed multivariate Cox analysis and further confirmed the independent prognostic role of the TGF-b risk score ( Figure 2D). As a result, the TGF-b risk score was a robust prognostic predictor for BLCA. Then, we integrated the TGF-b risk score with two other prognostic factors (stage and age) to exploit a comprehensive nomogram ( Figure 2E). The accuracy of the nomogram in predicting 1-, 3-, and 5-year overall survival was 0.729, 0.746, and 0.762, separately, which was obviously higher than the accuracy of the TGF-b risk score ( Figure 2F). The nomogram-predicted OS was consistent with the real OS as shown in the calibration curves ( Figure 2G). We appraised and validated the ability of the nomogram in predicting prognosis among three external BLCA cohorts, including the Xiangya cohort, GSE32894, and GSE48075. In the Xiangya cohort, the accuracy in predicting 1-, 2-, and 2.5year OS was 0.834, 0.851, and 0.830, separately ( Figure 3A). The nomogram-predicted OS was highly similar to the real OS ( Figure 3B). Similarly, we observed the same results in GSE32894 and GSE48075 ( Figures 3C-F TGF-b Risk Score Predicted the Immune Phenotypes of the TME and the Therapeutic Response to ICB The immune status of the TME influences the fate of cancer cells and predicts sensitivity to ICB. First, we analysed the relationships between the risk score and the activities of the anticancer immunity cycles. Activities of several steps, such as the release activities of cancer cell antigens, priming, and activation, were significantly upregulated in patients with highrisk score ( Figure 4A). This indicated that the anticancer  immune response of the TME in patients with higher risk scores was more activated and immunogenic. Subsequently, the activities of several steps that recruited TIICs, such as CD8 + T cells, Th1 cells, NK cells, and DCs, were also significantly upregulated in patients with higher risk score. Therefore, the anticancer killing activity was consistently higher in patients with high risk score. Besides, we verified that the infiltration levels of TIICs, including CD8 + T cells, macrophages, NK cells, DCs, and Th1 cells, were positively correlated with the risk score in six independent algorithms ( Figure 4B). Relationships between the  Table 4. These data revealed that a high-risk score predicted a hot TME status in BLCA. In line with the common sense that a hot TME was more sensitive to ICB, the TGF-b risk score positively correlated with the TIS (Spearman R = 0.34) ( Figure 4C). Furthermore, the risk score was positively related to the expression of 20 inhibitory immune checkpoints and the scores of several immunotherapy response related pathways, including the IFN-g signature, DNA repair, and base excision repair ( Figure 4D). Correlations between the risk score and immune checkpoints are shown in Supplementary Table 5. Notably, the risk score positively correlated with CD274 (PD-L1) (P < 0.05). Overall, a high TGF-b risk score indicated an inflamed TME phenotype and may be more sensitive to ICB.
To adjust the bias caused by clinicopathological heterogeneity, we performed further subgroup analyses based on different clinicopathological features. As expected, the correlations between TGF-b risk score and TME immune status were highly consistent in all the subgroups (Supplementary Figures 5-10). Therefore, the TGF-b risk score was a robust predictive biomarker for TME immune phenotypes.  In addition, we analysed the correlations between the risk score and TME stromal components, such as tumour fibroblasts and several stromal signatures. Tumour fibroblasts, as essential stromal components of the TME, were reported to have a close connection with the TGF-b pathway. Consistently, our results displayed a remarkable positive correlation between the TGF-b risk score and cancer-associated fibroblasts in several independent algorithms (Supplementary Table 4). In addition, the risk score positively correlated with many stromal signatures, such as pan-F-TBRS (panfibroblast TGFb response signature) and EMT (epithelial-mesenchymal transition) signatures (Supplementary Figure 11).

TGF-b Risk Score Predicted Molecular Subtypes and Guided Precision Medicine for BLCA
The molecular subtypes could guide individual precision treatment for BLCA. Therefore, we first correlated the TGF-b risk score with several traditional molecular subtype classification systems, including the UNC, TCGA, MDA, Lund, CIT, Consensus, and Baylor subtypes. As shown in Figure 5A, cancers with higher risk score were more likely to be basal subtype, featured by higher enrichment scores of basal differentiation and immune activation. However, cancers with lower risk score were more likely to be luminal subtype, featured by higher enrichment scores of luminal differentiation and urothelial differentiation. Furthermore, the results of ROC curves demonstrated that the accuracy of the TGF-b risk score in predicting the molecular subtypes of UNC, TCGA, MDA, Lund, CIT, Consensus, and Baylor systems were 0.79, 0.76, 0.80, 0.85, 0.81, 0.80, and 0.66, respectively ( Figure 5B).
The IC50 of six common chemotherapy drugs was lower in patients with high risk score ( Figure 5C). Similarly, the enrichment scores of EGFR ligands and radiotherapy responserelated pathways were higher in patients with high risk score ( Figure 5D). These data revealed that patients in the high risk score groups may benefit more from chemotherapy, EGFRtargeted therapy, and radiotherapy. However, several oncogenic signatures which supressed the anticancer immunity were enriched in patients with low-risk score ( Figure 5D). Therefore, blocking these oncogenic pathways may improve the prognosis of patients with low TGF-b risk scores.

Roles of TGF-b Risk Score in the Xiangya Cohort and External Validation Cohorts
We further explored the associations between the TGF-b risk score and TME immune phenotypes, molecular subgroups, and therapeutic response in our Xiangya cohort. First, we proved that the TGF-b risk score was positively related to many anticancer immunity steps ( Figure 6A left panel). In addition, the risk score positively correlated with the levels of several TIICs, including CD8 + T cells, macrophages, and DCs, in six independent algorithms ( Figure 6D). The TGF-b risk score was also positively related to the enrichment of ICB response-related signatures, the expression of 20 immune checkpoints, and TIS ( Figures 6A right panel, B, C). The results of the Xiangya cohort proved that a high risk score indicated an inflamed phenotype. Then, we further validated the role of TGF-b risk score in differentiating the molecular subtypes in the Xiangya cohort ( Figures 7A, B). Finally, we confirmed the closed associations between the risk score and clinical responses to several treatment options, including targeted therapies, EGFR therapies, and radiotherapy ( Figure 7C). It was apparent that patients with high-risk scores may benefit more from radiotherapy and EGFR targeted treatment. However, patients with low-risk scores may benefit more from targeted treatments, such as blocking out the WNT-b-catenin networks or FGFR3 pathways ( Figure 7C). All of these results were well-validated in two other cohorts, GSE32894 and GSE48075 (Supplementary Figures 12, 13). Finally, we correlated the risk score with therapeutic response and prognosis of patients in three immunotherapy associated cohorts (Supplementary Figure 14). In the IMvigor210 cohort, we did not observe relationship between the risk score and prognosis in the whole cohort or in three immune phenotypes. Most of patients in IMvigor210 cohort received chemotherapy before immunotherapy, which may diminish the correlations between the TGF-b risk score and therapeutic response and prognosis. In the GSE100797 cohort, patients in the high risk score group had better prognosis and higher ICB response rates, though the statistical difference was not significant. Similar results were observed in the GSE78220 cohort.

DISCUSSION
It is widely known that the TGF-b pathway is a critical factor for tumour development. However, TGF-b had two-sided effects (45). From one point of view, TGF-b protected tumour cells from malignant evolution. The inverse aspect stated that TGF-b modulated processes such as cell invasion and immune suppression, which cancer cells may take advantage of (17). In addition, the expression and activation of TGF-b could be considered potential targets for antitumor therapy (46). With revealing the mechanisms of TGF-b in regulating tumour biology, more and more attention has been paid to the therapeutic ability of TGF-b (17).
Currently, the TGF-b signalling pathway is extensively used for prognosis prediction in various tumours, such as ovarian cancer, hepatocellular carcinoma, non-small-cell lung cancer (NSCLC), and colorectal cancer (47)(48)(49)(50). Previous works have shown that TGF-b is critical for development and progression of bladder cancer. For instance, TGF-b-induced transgelins could promote bladder cancer metastasis (51), and the TGF-b-related miR-1305/TGF-b2/smad3 pathway was shown to help circRIP2 accelerate bladder cancer development (52). These studies indicated a close correlation between TGF-b and bladder cancer. However, little work has systematically contributed to exploring the role of TGF-b-related signatures in regulating TME characteristics and predicting prognosis in BLCA. Thus, we thoroughly explored TGF-b in BLCA in this study to compensate for this gap. We screened various differentially expressed TGF-b-related genes. Then, we selected five vital  genes among them to construct a TGF-b signature according to the standard of our method. These genes included FNTA, PPP2CB, LRRC32, ID2, and MAP2K1. These five genes were correlated with tumour progression. FNTA was regulated by scaffold association factor B (SAFB) to promote RAS membrane association, which was buttressed to exploit a rational approach to anticancer therapy (53). It was remarkable that our team once intensively studied FNTA and found its association with BLCA. The androgen receptor promoted BLCA by regulating the expression level of circRNA-FNTA (54); PPP2CB was reported as protein phosphatase 2A. This enzyme regulated various of cellular biological processes, including tumour development (55). It is one of the indicators for 8p suppressor gene (56). LRRC32 was upregulated on CD4 + CD25 + FOXP3 + Treg cells; thus, it was always considered a Treg-specific activation marker (57). GARP is regarded as a receptor for LAP/latent TGF-b (58). The GARP and GARPligand signalling pathways was demonstrated to be a target for cancer immunotherapy (58). ID2 plays a critical role in the survival of aggressive cancer cells (59) and positively modulates the activity of HIF2a in glioblastoma (60). MAP2K1, also known as MEK1, plays a critical role in inhibiting the progression of histiocytic neoplasms (61). Based on these five genes, we integrated multiple BLCA datasets (TCGA-BLCA, GSE32894, GSE48075), in addition to the Xiangya cohort, to develop and validate a novel TGF-b risk score. TGF-b signalling was proven to be a dominant suppressor of adaptive and innate immune responses during tumour progression (62). We regarded the risk score as reflecting the characteristics of the TME from several aspects. First, we verified that the TGF-b risk score predicted the survival and immune characteristics in the TME. As displayed in our plots, patients with higher risk scores always had lower OS rates, higher tumour grade and stage. Second, we screened the difference in the enrichment scores of immunotherapy response-related signatures between the high-risk and low-risk score groups. The IFN-g signature, APM signal, and base excision repair were obviously enriched in patients with high risk score. Chemotherapy significantly improved the survival for patients with advanced BLCA (63). Therefore, the development of accurate indicators for chemotherapy is extremely urgent. Interestingly, we found that patients with high-risk scores might benefit more from chemotherapy.
Here, the TGF-b risk score had positive correlations with the enrichment scores of several anticancer immunity cycles, TIS, and TIICs. As a result, the activity of pre-existing anticancer immunity in the TME was higher in patients with high risk score (37). Based on the theory of immune checkpoints inhibiting anticancer immunity in the TME (64), we consistently correlated the risk score with the expression of representative immune checkpoints, including CD274, CTLA4, and CD276. The correlation was positive and significant between the risk score and ICI. The results showed that a tumour TME with a higher risk score indicated a better effect of immunotargeted therapy and more sensitivity to immunotherapy. However, inadequacy and doubt still existed. The results also showed that the risk score positively correlated with macrophages (M2), which were considered a signal to inhibit anticancer immunity ( Figure 4B). Hence, the pre-existing anticancer activities of patients with high risk scores may be restrained by the higher level of tumour infiltrated macrophages (M2) and the over-expression of many inhibitory immune checkpoints (such as CD274 and LAG3). Treatments, especially ICB, were suitable for patients with highrisk scores whose anticancer immunity was suppressed to a certain extent (8). In contrast, patients with low-risk scores always had lower TIS and lower expression of immune checkpoints, which indicated that they might not be suitable for ICB.
Anil Korkut et al. presented a pan-cancer analysis of genomic alterations that regulated TGF-b signalling in a large sample set. The results identified mutation hotspots in the TGF-b superfamily, which indicated potential biomarkers for further therapeutic and cancer diagnostic studies. Their work demonstrated the significance and prospect of a thorough investigation of TGF-b (65). Several studies focusing on the prognostic prediction of TGF-b have been published. Slavica Stojnev et al. evaluated the associations of TGF-b1, Smad4, and Smad2 with clinicopathologic characteristics and prognosis in urothelial bladder cancer (25). We compared our work with these previous works. First, our key genes were screened from numerous TGF-b-related genes with more comprehensive algorithms, such as the lasso algorithm and Akaike information criterion. Second, we constructed a five-gene signature and correlated it with molecular subtypes, TME immune characteristics and therapeutic opportunities. Third, we validated our model in several external cohorts. Of note, we did not perform immunohistochemical staining, which could be considered a flaw. The construction and validation of scoring models has been a core of ceaseless studies in oncology research. Numerous risk models related to EMT, immune TME, and hypoxia have been lucubrated and reported with the potential to predict prognosis and guide precision medicine for BLCA (23,25,(65)(66)(67)(68). It was remarkable that our previous work clarified thought in screening key genes, constructing a gene model, associating the model with immune characteristics, and evaluating the accuracy for predicting prognosis and treatment response. In the current study, we adopted a similar thought and obtained results as expected.
However, drawbacks inevitably existed in several aspects of our study. First, all of our outcomes, including the construction of the signature and the nomogram, were based on public data and clinical cohorts. Although we also validated these results in the Xiangya cohort, molecular mechanisms and experiments in vivo or in vitro may be necessary. Second, only two external cohorts and the Xiangya cohort (only 57 samples included) contributed to the verification of our outcomes. More prospective data are urgently needed for further validation of the TGF-b signature due to our clinical data limitations. Third, the clinical appliance of the TGF-b signature needs more exploration. The risk model was developed based on the RNAseq data. As a result, qPCR may be a more suitable method to calculate the risk score during clinical practice. In addition, we will design a corresponding genetic testing kit based on these five genes of the risk model. Fourth, we dogmatically consider the median of the TGF-b risk score as the cut-off value in all analyses, which was needed to be further validated.
In summary, we took one small step forwards in the sphere of exploring rational prognostic predictors for BLCA. We developed and validated a well-grounded TGF-b risk score based on realistic clinical data. The TGF-b risk score was proven to stratify the clinical outcomes and TME phenotypes of BLCA. It also indicated the sensitivity of several treatment options for BLCA. We found that patients with high-risk scores may benefit more from ICB, radiotherapy, chemotherapies, and EGFR-targeted therapy. However, patients with low-risk scores may be more suitable to receive other targeted therapies, including blocking out the FGFR3 or WNT-b-catenin networks.

ETHICS STATEMENT
The collection of the Xiangya cohort was approved by the Ethics Committee of the Xiangya Hospital of Central South University.

ACKNOWLEDGMENTS
We sincerely thank all participants in the study. Supplementary Table 1 | TGF-b genes and detailed information of four BLCA cohorts applied in our analysis.
Supplementary Table 2 | Sixty differentially expressed TGF-b genes between BLCA and normal tissues.
Supplementary Table 3 | (A) Coefficients of 5 optimal candidate genes included in the TGF-b risk score with univariate Cox regression analysis. (B) Coefficients of 5 optimal candidate genes included in the TGF-b risk score algorithm.