No ethnical approval nor informed consent was required in this study to the public availability of data in the public database. Based on the GeneCards database, 258 ERS-related genes were obtained (Supplementary Tables S1, S2), and the screening criteria were protein-coding genes with a correlation score greater than 12. Obtain BLCA gene expression matrix, clinical information and mutation information through the TCGA website (Blum et al., 2018). GSE cohorts were gathered from Gene Expression Omnibus (GEO), namely, GSE32894 and GSE13507 (Barrett et al., 2013). The training dataset was TCGA-BLCA, accepting samples from 400 patients with survival information. 224 samples of GSE32894 cohort and 165 samples of GSE13507 cohort with valid survival time was used for validation (Kim et al., 2010; Sjödahl et al., 2012). Obtain the expression matrix and patient information for the IMvigor210 cohort via the “IMvigor210CoreBiologies” R package (Necchi et al., 2017). TISCH database was used to characterize the expression of the signature genes in different kinds of cells of bladder cancer tissue (Sun et al., 2021).

To investigate the pathways associated with ERS-related genes in BLCA, we used the “clusterProfiler” package in R to perform GO and KEGG enrichment analysis on ERS related genes (Yu et al., 2012) and the pathways obtained from MsigDB website (Liberzon et al., 2015). ERS-related pathways (Supplementary Table S3), tumor-related pathways (Supplementary Table S4) and metabolism-related pathway were also available from the MsigDB website (https://www.gsea-msigdb.org/gsea/msigdb/). In addition, to explore the mechanisms of each molecular subtyping more precisely, we collected stromal-related pathways from a study conducted by Mariathasan et al. (2018) (Supplementary Table S5). The “GSVA” package was used to calculate the enrichment score for TCGA samples (Hänzelmann et al., 2013).

We performed consensus clustering analysis of ERS-related genes using the “ConsensusClusterPlus” package. The same clustering method was used for both the training and validation sets (Wilkerson and Hayes, 2010). GSE32894 cohort used for independent validation of our molecular subtypes.

The “maftools” package was used for analysis of mutations in BLCA subtypes (Mayakonda et al., 2018). The top ten genes with the highest degree of mutation in each subtype are shown by waterfall diagram. The “TCGAmutations” package was used for obtaining TMB scores of TCGA patients.

We used the “estimate” package for calculating immune, stromal and estimate scores. By using the TIMER method, immune infiltration analysis was performed on each sample of the four subtypes, thus comparing the immune microenvironment of each subtype (Li et al., 2017). To make the results more reliable, we still used CIBERSORT, CIBERSORT-ABS, QUANTISEQ, MCPCOUNTER, XCELL and EPIC methods to compare the immune infiltration of each subtype (Becht et al., 2016; Aran et al., 2017; Finotello et al., 2019; Racle and Gfeller, 2020; Le et al., 2021). Finally, we included the cancer immunity cycle (Chen and Mellman, 2013). The cycle reflects the 7 steps in which immune cells exert their anti-tumor effects (Supplementary Table S6).

In order to achieve precise treatment of patients with bladder malignancies, many experts have already performed molecular subtypes of BLCA. There are seven main molecular typing methods (Consensus, TCGA, MDA, Lund, CIT, UNC, Baylor) (Sjödahl et al., 2012; Choi et al., 2014; Damrauer et al., 2014; Rebouissou et al., 2014; Robertson et al., 2017; Mo et al., 2018; Kamoun et al., 2020)that are currently recognized. “ConsensusMIBC” and “BLCAsubtyping” R packages were used to divided samples of TCGA and GSE32894 to different subtypes in our study. Although these different methods could classify samples into a variety of types, the vast majority of patients can be categorized as luminal and basal subtypes. Additionally, we collected 12 specific classical signatures of BLCA. These results were summarized in a figure in the form of a heat map.

We used a public database, Genomics of Drug Sensitivity in Cancer (https://www.cancerrxgene.org/), to assess the response of four subtypes to chemotherapy with four drugs (Yang et al., 2013). The drug response was estimated by calculating the half-maximal inhibitory concentration (IC50). Additionally, the tumour immune dysfunction and exclusion (TIDE) algorithm (Jiang et al., 2018) and Immunophenoscore (IPS) were used for predicting the sensitivity of immunotherapy (Charoentong et al., 2017).

First, through univariate cox analysis and “survival” R package, we identified prognosis-related genes from 258 ERS-related genes. Then, to avoid overfitting of the model, Lasso regression analysis was performed on prognosis-related genes. Genes with non-zero coefficients at the best lambda value were included in multivariable Cox analysis. Finally, we obtain the prognostic risk score signature for ERS-associated genes. The following is the prognostic scoring formula: E R S   s c o r e = Ʃ β i * R N A i , w h e r e   β i   i s   t h e   c o e f f i c i e n t   o f   t h e   i − t h .

After that, we divided the patients into high ERS score and low ERS score groups according to their median ERS scores. The Kaplan-Meier method was applied to plot the survival curves, while the log-rank test was applied to calculate statistical significance. The feasibility and accuracy of ERS score for predicting 1-, 3-, and 5-year outcomes were further assessed by generating ROC curves. Meanwhile, the ERS scores were also highly effective in two independent validation sets GSE32894 and GSE13507.

For individualized assessment, we included stage, age, gender and ERS scores of the sample in univariate and multivariant Cox analysis. After that, we divided the age by 70 years into younger and older. The relationship between ERS score and age, gender, stage, survival, T-stage, N-stage and M-stage was further evaluated. We also evaluated the clinical outcomes of different types of patients in the high and low ERS score groups. Finally, stage, age and ERS score were used to construct comprehensive nomogram to enable prognostic assessment of individualized patients. In addition, we explored the correlation between the ERS score, the four molecular subtypes, and stage.

The data were analyzed through the use of R software (version 4.1.3). We analyzed the correlations between variables using Pearson or Spearman coefficients. Statistical tests were two-sided, and the level of significance was set at p < 0.05.

An RNeasy Mini Kit (Qiagen, Valencia, CA) was used to extract total RNA of specimens. According to the manufacturer’s instructions RNA-seq was performed using the QuantSeq kit FWD HT kit (Lexogen) using 500 ng input RNA. We used NEBNext UltraTM RNA Library Prep Kit for Illumina (NEB, USA) to generate Sequencing libraries. Consequently, 125- to 150-bp paired-end reads were generated using an Illumina HiSeq platform. Additionally, 6 μg of total RNA per sample was used as the input material for the small RNA library using the NEBNext Multiplex Small RNA Library Prep Set for Illumina (NEB, Ipswich, MA). 50-bp single-end sequencing was performed on an Illumina Hiseq 4000 platform. HISAT2 aligner was used to align the raw reads to the human reference genome GRCh38. Finally, we obtained the baseline phase RNA sequencing results for 29 patients (named the TMU-BLCA cohort). All patients included in this study signed an informed consent form and were approved by the Ethics Committee of the Second Hospital of Tianjin Medical University. Clinical information for patients can be found at Supplementary Table S7.

The BCa cell lines T24, UMUC3 and EJ were purchased at the Chinese Academy of Sciences Cell bank, and the 253J-BV was presented by Professor Li Lei of the First Affiliated Hospital of Xi’an Jiaotong University. All cells were cultured in RPMI1640 or MEM with 10% fetal bovine serum (FBS) at 37°C in a 5% CO2/95% air incubator.

For RT-qPCR, total RNA was extracted from 20 paired samples of bladder cancer and adjacent normal tissues by surgical excision in the Second Hospital of Tianjin Medical University. The specific situation of the patients is reflected in Supplementary Table S8. The total RNA isolation and subsequent RT-qPCR were conducted as previously described (Li et al., 2022). The primers sequences (Sangon Biotech) were as follow: ATP2A3: F CAT​CCT​GAC​GGG​TGA​ATC​TGT R TGC​CCG​ATG​TGA​TAT​TGG​TGC. STIM2: F AGA​CAA​CAA​TGT​CAA​AGG​AAC​GA R ACT​CCG​GTC​ACT​GAT​TTT​CAA​C. VWF: F CCG​ATG​CAG​CCT​TTT​CGG​A R TCC​CCA​AGA​TAC​ACG​GAG​AGG. P4HB: F GGC​TAT​CCC​ACC​ATC​AAG​TTC R TCA​CGA​TGT​CAT​CAG​CCT​CTC. GAPDH: F CGG​AGT​CAA​CGG​ATT​TGG​TC R TTC​CCG​TTC​TCA​GCC​TTG​AC.

ERS-related genes are mainly enriched in endoplasmic reticulum function-related GO and KEGG terms. In addition to this, they are also enriched in Human cytomegalovirus infection, Pathways of neurodegeneration-multiple diseases and response to topologically incorrect protein (Figure 1A). Overall, ERS related genes included in the study were significantly associated with endoplasmic reticulum stress, and they were also significantly enriched in the immune pathway. It confirmed a non-negligible relationship between ERS and the immune microenvironment.

Consensus clustering analysis of ERS genes classified BLCA into four subtypes (Figures 1B, C). These ERS-related genes were significantly different among the four subtypes, and more highly expressed genes were present in C1 and C2 (Figure 1D). Kaplan-Meier survival curves indicated that the C2 had the worst prognostic outcome, and the C3 had the best prognostic outcome (Figure 1E). With the same clustering approach, we observed the same phenomenon in the validation queue (Figures 3A–C). Combined with the clinical data, the C3 had younger age, better N-stage, pathological stage and T-stage (Supplementary Figures S1A–D). Similarly, the same phenomenon is observed in the verification cohort (Supplementary Figure S2). Thus, C3 may still be in the early stages of tumor development and patients have a better prognosis.

To further investigate the characteristics of ERS in the four subtypes, we enriched the ERS-related pathways (Figure 2A). The results showed that ERS activity was significantly elevated in C1 and C2. Hypoxia affects ERO1α-mediated protein post-translational folding and disulfide bond formation leading to ERS (Chen and Cubillos-Ruiz, 2021). Interestingly, the activity of hypoxia-related pathways was also significantly elevated in C1 and C2. Similarly, we enriched tumor-associated pathways to explored the characteristics of the four isoforms precisely (Figure 2B). C3 was mainly associated with fatty acid metabolism, while C2 was mainly associated with immune pathways. Tumor-specific pathways such as MTOR signaling pathway, p53 signaling pathway as well as ERBB signaling pathway were enriched in C1, while others (NOTCH signaling pathway, MAPK signaling pathway, JAK-STAT signaling pathway, TGF-β signaling pathway, WNT signaling pathway and Hedgehog signaling pathway) were enriched in C2. C4 was associated with cell cycle and DNA damage repair related pathways. In addition, C2 had higher Angiogenesis, EMT1, EMT2, EMT3, and Pan-F-TBRs scores which means a higher stromal level may contribute to its worse prognosis (Figure 2C). We also explored metabolism-related pathways. The results showed that most metabolism related pathways were enriched in C1 and C2, while pyruvate metabolism, fatty acid metabolism, and tyrosine metabolism related pathways were enriched in C3 subgroup. It implied that differences in the metabolic microenvironment of different subtypes may influence to some extent the clinical outcome of patients (Supplementary Figure S2E). Mutational analysis revealed that unlike C3, which exhibited high tumor mutation burden in FGFR3, KDM6A, and TBC1D12, C2 exhibited high tumor mutation burden in TP53 and FLG (Figures 2D–G). The C4 has the highest total tumor mutation load (TMB) (Figure 2H). There are previous studies that reported FGFR3 mutations are strongly associated with earlier stage and longer survival time (van Rhijn et al., 2020). While TP53 mutations have opposite outcomes (Sjödahl et al., 2020). These researches coincided with our findings.

After that, we explored TME of C1-C4. C2 had a highest estimate score, immune score and stromal score while C3 was lowest both in TCGA and GSE32894 cohorts (Figures 2I–K–KFigures 2I–K, 3D). Depending on seven immune infiltration methods of analysis, C2 had more immune infiltration which followed by C1 subtype. Conversely, C3 has least immune infiltration (Figure 4A). Immunomodulators are important for assessing the immune microenvironment of tumors. The gene expression of immunomodulators was elevated mainly in C1 and C2, further validating the high immune activation status of C1 and C2, especially in C2 (Figure 4B). In addition, C2 had stronger cancer immunity cycle activity and there was more activation of anti-tumor immune cells (Figure 5A). The C2 immune checkpoint expression was significantly higher (Figure 5B). These results suggested that the C2 may have a poor clinical outcome due to significant tumor immune escape though with a stronger immune response. The C2 phenotype was also closely associated with high M2 macrophage infiltration (Figure 5C). There was substantial evidence that M2 macrophage polarization suppresses the immune microenvironment and causes a poorer outcome (Komohara et al., 2016; Wei et al., 2022). This further substantiated our perspective that C2 phenotype is associated with immune evasion. Further validation was also obtained in the validation queue (Supplementary Figures S2A, B, S3A, B).

Next, we compared these four subtypes with pan-cancer immune subtypes and BLCA molecular subtypes. We obtained the pan-cancer subtypes of TCGA BLCA samples from a previous study (Thorsson et al., 2018). The results showed that C1 and C2 were significantly associated with pan-C2 (IFN-γ) subtype, exhibiting the greatest amount of immune infiltration and poor prognosis, which is consistent with our study. The C3 was clearly associated with pan-C4 (Lymphocyte Depleted) subtype and shows lymphocyte depletion (Figure 5E). In comparison with the currently accepted molecular typing of bladder cancer, the C2 was more of the basal type, the C3 and C4 was more luminal subtype, and the C1 might be intermediate (Figure 5F). These results were similar for the validation cohort (Supplementary Figure S2C).

Cisplatin, Gemcitabine, Paclitaxel and Doxorubicin are the most common chemotherapy drugs for BLCA. All four of these drugs showed the best efficacy in C1 and C2 (Supplementary Figures S4A–D). Therefore, the C1 and C2 subtype may have better benefit with chemotherapy drugs. Higher TIDE scores indicate that patients will benefit more from immunotherapy. Thus, we speculation C4 may more suitable for immunotherapy (Figure 5D).

To better assess the prognosis of BLCA patients for clinical benefit, we constructed a ERS related prognostic signature in TCGA training database. The 258 ERS-related genes were analyzed by Univariate Cox analysis and 75 prognosis-related genes were obtained (Supplementary Figure S5). After LASSO analysis, 21 genes were included in the Multivariable Cox analysis (Figures 6A, B). Ultimately, we developed the 4-gene prognostic signature. E R S   s c o r e = − A T P 2 A 3 * 0.29859 − S T I M 2 * 0.58523 + V W F * 0.31260 + P 4 H B * 0.35490 . We divided the sample into high ERS score and low ERS score groups according to the median scores (Figure 6E). Prognosis was predicted by Kaplan-Meier curves for OS in high and low ERS score groups (Figure 6C). The results showed that patients in the high ERS score group had a worse prognosis (p < 0.001). For the OS survival prediction of ERS signature, the AUC of the ROC curves for 1, 3 and 5, years were 0.698, 0.687, and 0.702 (Figure 6D). It demonstrated that ERS signature we constructed has strong predictive value for OS prognosis in the training set.

To further demonstrate that ERS scores could steadily predict prognosis, two validation sets GSE32894 and GSE13507 were integrated in our study. Similarly, we also divided the sample into high and low ERS score groups according to the median of the ERS scores (Figures 6J, K). The OS Kaplan-Meier curves for both validation sets indicate that high ERS score group has the worse OS prognosis (Figures 6F, G). The AUC of the ROC curves for 1, 3, and 5 years in the GSE32894 were 0.762, 0.744, and 0.745 (Figure 6H). The AUC of the ROC curves for 1, 3, and 5 years in the GSE13507 were 0.646, 0.712, and 0.690 (Figure 6I). The results were evidence of the powerful prognostic capability of the ERS signature.

Further, we explored the correlation between ERS score and clinical characteristics. Univariate and Multivariable Cox analysis showed that ERS score, stage and age were independent risk factors for the prognosis of BLCA (Supplementary Figures S6A, B). In addition, ERS scores were higher in patients older than 70 years, with higher stage, poorer OS prognosis, higher T-stage, and higher N-stage. However, there were no statistically significant differences in ERS scores by gender and M-stage (Supplementary Figures S6C–I). A clinically stratified analysis of gender, age, stage, T-stage, N-stage, and M-stage was performed in the TCGA database to analyze ERS signature prognostic performance. The results showed that patients in the low ERS scores cohort had improved survival outcomes compared with those in the high ERS score cohort in gender, age, stage, T-stage, and N0 (p < 0.05, Supplementary Figures S7A–I). In contrast, no prognostic differences were observed between the low and high ERS score cohort in the clinical stratification of N1-3 and M stage (Supplementary Figures S7J–L). To enhance clinical applicability, we constructed nomogram using stage, age, and scores (Figure 7A). Calibration plots and decision curves for 1-, 3-, and 5-year survival prediction indicated that nomogram has good predictive accuracy and benefit (Figures 7B–E).

To explore the characteristics between the high and low ERS score groups, we explored the differences in immune microenvironment and immunotherapy effects between the two groups. M0 macrophages, M2 macrophages, mast cells, and neutrophil infiltration were evident in the high ERS score group. In contrast, immune cells with tumor-killing properties such as CD8T cells was enriched in low-ERS score group (Figure 8A). The estimate score, immune score and stromal score were also higher in high ERS score group than in low ERS score group (Figures 8B–D). To evaluate the effect of immunotherapy, we introduced the IMvigor210 cohort and the IPS score for the study. The results showed that the prognosis and complete response rate of patients in low ERS score group were better than those in high ERS score group after PD-L1 blockage therapy (Figures 8E, F). Furthermore, IPS scores, IPS-CTLA4 blocker scores, IPS-PD1 blocker scores, and IPS-CTLA4 and PD1 blocker scores were higher in samples with low ERS scores (p < 0.05; Figures 8G–J), indicating that BLCA samples with low ERS scores may be suitable for anti PD-1 and CTLA-4 immunotherapy. Inspiringly, we were surprised to find that the ability of the ERS signature to predict immunotherapy efficacy was robust in our TMU-BLCA cohort. Patients with low ERS scores had good immunotherapy outcomes. These results further confirm the clinical importance of the ERS signature (Supplementary Figure S10F). Finally, we found that most patients in the high ERS score cohort were associated with C2 and advanced stage, whereas patients with low ERS score were C3 and earlier stage (Figure 8K).

Through TISCH website, we accessed to three single-cell databases GSE130001, GSE145281 and GSE149652 (Supplementary Figures S8A–C). P4HB expressed in a variety of cells, mainly expressed in stromal cells, including endothelial, fibroblasts, myofibroblasts and epithelial (Supplementary Figures S8D, S9A). VWF mainly expressed in endothelial (Supplementary Figures S8E, S9B). ATP2A3 mainly expressed in immune cells such as B-cell, while STIM2 expressed in a few amounts in a variety of cells (Supplementary Figures S8F–G, S9C–D).

We verified the expression patterns of ERS model genes (ATP2A3, STIM2, VWF, and P4HB) in BLCA cell lines and normal uroepithelial cell lines by RT-qPCR (Supplementary Figure S10A). Compared with the expression levels in SV-HUC-1 uroepithelial cells, ATP2A3 and VWF were meaningfully low expressed in 4 bladder cancer cell lines (T24, 253J-BV, EJ, UMUC3), while P4HB were high. STIM2 were low expressed in T24, 253J-BV and UMUC3 cell lines but not significant in EJ cells. In addition, we found downregulated expression of ATP2A3, VWF in tumor tissues (Supplementary Figures S10B–E). Conversely, P4HB were upregulated. The above results implied that ERS model genes may play an important potential role in bladder cancer progression.