Abstract Despite remarkable advancements in therapeutic strategies, a considerable proportion of patients with bladder cancer (BC) still experience disease progression and unfavorable prognosis. The heterogeneity and biological functions of tumor endothelial cells (ECs) during BC progression remain poorly understood. We collected scRNA-seq data from BC samples and identified two EC subpopulations through hierarchical clustering analysis. The activity of signaling pathways in distinct EC subpopulations was assessed utilizing AUCell analysis. Gene regulatory networks (GRN) were constructed and analyzed for different EC subpopulations using the pySCENIC algorithm. Additionally, we investigated the association between the abundance of EC subpopulations and both clinical prognosis and immune cell infiltration. The biological effects of ESM1 protein on BC cells were further validated through EdU and Transwell assays. We analyzed 7,519 CD45-negative single cells from BC tissues and discerned two distinct EC subpopulations. The two subpopulations were characterized by high expression of ESM1 (S1 ECs) and CXCL2 (S2 ECs), respectively. In S1 ECs, we observed significant activation of signaling pathways involved in tumor promotion, including angiogenesis and cell proliferation. Additionally, our GRN analysis uncovered notable differences in transcription factor activity between S1 and S2 ECs. Moreover, ESM1 protein promoted proliferation and migration of BC cells. Patients with higher abundance of the S1 EC subpopulation exhibited more unfavorable clinical outcomes and increased infiltration of inhibitory immune cells. Our findings elucidate the transcriptional profiles and biological roles of the high ESM1-expression endothelial cell subpopulation in BC. This subpopulation is associated with poor prognosis and immunosuppressive tumor microenvironment. Accordingly, targeting endothelial cells with high ESM1 expression may offer a novel therapeutic strategy for patients with BC. Supplementary Information The online version contains supplementary material available at 10.1038/s41598-025-95731-2. Keywords: Single-cell RNA sequencing, Bladder cancer, Prognosis, Endothelial cell subpopulation, Immune microenvironment Subject terms: Cancer, Urological cancer Background In 2022, approximately 610,000 new cases of bladder cancer (BC) were diagnosed globally, with over 220,000 deaths attributed to the disease^[34]1. Tumors confined to the mucosa or submucosa are classified as non-muscle invasive bladder cancer (NMIBC), whereas those extending into the muscle layer are categorized as muscle invasive bladder cancer (MIBC)^[35]2,[36]3. Despite undergoing transurethral resection and subsequent intravesical treatment with chemotherapy or immunotherapy, many patients with NMIBC eventually experience recurrence and progression to MIBC, resulting in unfavorable prognosis and poor quality of life^[37]4. For patients with MIBC, cisplatin-based chemotherapy remains the mainstay treatment; however, its efficacy is relatively limited^[38]5. Despite evidence highlighting the efficacy of immune checkpoint inhibitors in treating advanced BC, a substantial proportion of patients still exhibit an unsatisfactory response to PD1/PD-L1 inhibition immunotherapy^[39]6,[40]7. Accordingly, elucidating the mechanisms underlying the recurrence and progression of BC is critical for improving the clinical outcomes of BC patients and exploring novel therapeutic targets. A growing body of evidence indicates that non-cancerous cells and their secreted factors in the tumor microenvironment (TME) play indispensable roles in cancer progression, metastasis, and the development of treatment resistance^[41]8. Endothelial cells (ECs) represent a critical stromal cell type in the TME and play a crucial role in angiogenesis and tumor metastasis^[42]9. Tumor ECs exhibit distinct morphological and molecular features compared with normal ECs^[43]10. In addition, tumor ECs in the TME display significant heterogeneity, and distinct EC subpopulations exert markedly different effects on tumor progression^[44]11. The interplay between cancer cells and ECs is critical for cancer metastasis and progression. A previous study has demonstrated that lung cancer cells facilitate tumor angiogenesis and endothelial cell migration via paracrine signaling mediated by extracellular vesicles^[45]12. Additionally, cancer cells have been reported to promote angiogenesis by activating the PI3K/Akt or Notch signaling pathways in ECs^[46]13,[47]14. On the other hand, tumor ECs can directly induce epithelial- mesenchymal transition in cancer cells and maintain tumor stemness through the production and release of EGF^[48]15. Biglycan, a secretory proteoglycan derived from tumor ECs, contributes to metastasis through the activation of the Toll-like receptor and NF-κB/ERK signaling pathways^[49]16. Furthermore, specific tumor EC subpopulations modulate the immune response and contribute to the formation of an immunosuppressive TME^[50]17. However, the heterogeneity of tumor-associated endothelial cells and their interactions with cancer cells in BC remain elusive. To elucidate the heterogeneity of tumor-associated endothelial cells (ECs) in BC, we conducted a comprehensive analysis of single-cell RNA-sequencing (scRNA-seq) data. Our analysis uncovered a previously unidentified EC subpopulation that is characterized by elevated expression of ESM1. We utilized bulk RNA-seq data to examine the clinical significance of high ESM1-expression EC subpopulation. In addition, recombinant ESM1 protein was employed to assess the impact of ESM1 on the proliferation and migration of BC cells in vitro. Furthermore, additional scRNA-seq and single-nucleus RNA-sequencing (snRNA-seq) datasets were employed to validate our findings. Our data indicate that targeting tumor-associated ECs may represent a promising therapeutic strategy for patients with BC. Methods Collection and processing of data We obtained scRNA-seq data from CD45-negative stromal cells, which were sorted by flow cytometry, from four BC tumor samples ([51]GSE211388). To further validate our results, we integrated three additional, independent scRNA-seq datasets ([52]GSE129845, [53]GSE135337, [54]GSE222315) that were not subjected to flow cytometry sorting. These datasets included 16 bladder cancer and 8 normal bladder samples. Single-nucleus RNA-seq data were obtained from [55]GSE169379, including 25 MIBC and four normal bladder samples. Bulk RNA-seq data and clinical features from The Cancer Genome Atlas (TCGA) dataset were accessed via the UCSC Xena platform ([56]https://xenabrowser.net/). FPKM values were converted into TPM values for further analysis. Transcriptome data and clinical information for patients with metastatic BC in the IMvigor210 immunotherapy cohort were obtained using the “IMvigor210CoreBiologies” R package. Microarray expression data and corresponding clinical information were retrieved from the datasets [57]GSE31684, [58]GSE13507, and [59]GSE32548. The main information for these datasets is presented in Supplementary Tables 1 and 2. scRNA-seq and snRNA-seq data processing The scRNA-seq data were analyzed using the “Seurat” package (version 4.1.1). Cells with a mitochondrial gene proportion exceeding 10% or with fewer than 200 or more than 5,000 detected genes were classified as low-quality and excluded from further analysis. To ensure analytical accuracy, the DecontX algorithm was employed to estimate and remove contamination from freely floating RNA molecules in individual cells^[60]18. Subsequently, we normalized and scaled the expression data. Batch effects were removed using the harmony algorithm^[61]19. The top 2000 highly variable genes were selected for analysis. For dimensionality reduction, principal component analysis (PCA) was conducted and the top 30 principal components were utilized for further analysis. Clustering analysis was conducted using a resolution parameter of 0.5. We manually annotated the main cell types according to the following canonical cell markers: epithelial cells (EPCAM, KRT8, KRT18), fibroblasts (DCN, COL1A1), myofibroblasts (ACTA2, RGS5), and endothelial cells (PECAM1,VWF, CLDN5)^[62]20–[63]22. The aforementioned methods were also employed to process and analyze snRNA-seq data in the present study. Differential expression analysis We identified marker genes that are specifically expressed in two distinct EC subpopulations utilizing the FindAllMarkers function. Genes with log2(fold-change) ≥ 1 and adjusted p-value < 0.05 were considered significantly upregulated. When comparing the gene expression profiles of the same cell subtype between two distinct groups, single-cell differential expression (DE) analysis methods typically yield more false discoveries than pseudobulk approaches^[64]23. To minimize false discoveries, we employed the pseudobulk method to conduct DE analysis between MIBC ECs and normal ECs in snRNA-seq data using the “muscat” package^[65]24. Genes with |log2(fold-change)|≥2 and adjusted p-value < 0.05 were identified as differentially expressed genes (DEGs) associated with MIBC ECs. Pathway enrichment analysis We employed the “ClusterProfiler” package to discover significantly enriched Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways in two distinct tumor EC subpopulations^[66]25–[67]27. Moreover, AUCell analysis was conducted to assess the activity scores of various functional pathways in each cell. The gene sets were downloaded from the Molecular Signatures Database ([68]https://www.gsea-msigdb.org/gsea/msigdb). Pseudotime trajectory analysis The “Monocle” package was employed to elucidate the developmental trajectory of ECs in BC^[69]28. DEGs were identified using the differentialGeneTest function with a q-value threshold of 0.05. Dimensionality reduction was conducted using DDRTree method. To elucidate differences in gene expression profiles between two distinct branches, we applied Branch Expression Analysis Modeling (BEAM). Gene regulatory network (GRN) analysis We employed pySCENIC (version 0.12.0) to analyze transcription factor (TF) regulatory networks in distinct tumor EC subpopulations. The GRNBoost method was utilized to construct a gene co-expression network. Regulons, comprising TFs and their target genes, were identified using the RcisTarget method. To evaluate the transcriptional regulatory activity of distinct TFs in ECs, the AUCell package was applied to calculate regulon activity scores (RAS). Definition and evaluation of EC subpopulation signature score The S1 and S2 EC signatures comprised the top 20 marker genes specific to their respective subpopulations. The single-sample gene set enrichment analysis (ssGSEA) algorithm was employed to calculate EC signature scores in both snRNA-seq data and bulk RNA-seq datasets. Analysis of immune cell infiltration in TCGA We analyzed the proportions of 22 immune cell types in the TME using the CIBERSORT algorithm, with gene sets obtained from the CIBERSORT website ([70]https://cibersortx.stanford.edu/). In addition, we evaluated the abundance of 28 tumor-infiltrating immune cell types using the ssGSEA algorithm^[71]29. Furthermore, the xCell algorithm was employed to assess the abundance of immune and stromal cells in the TME^[72]30. Establishment and validation of a MIBC-EC-related prognostic model First, we performed univariate Cox regression analysis on DEGs between MIBC ECs and normal ECs to select genes significantly associated with overall survival in bladder cancer patients. Subsequently, LASSO regression analysis was employed to diminish overfitting and screen prognostic genes using the glmnet package. Furthermore, we utilized stepwise regression analysis to establish a prognostic model according to the Akaike information criterion. This analysis was performed using the “MASS” package. The risk score was calculated using the following formula: mRNA1×COEF1 + mRNA2×COEF2 + mRNA3×COEF3+…+mRNAn×COEFn, in which “COEF” and “mRNA” represent the coefficient and expression value of each prognostic gene in the risk model, respectively. Patients were classified into high-risk and low-risk groups based on the median value of their risk scores. The receiver operating characteristic (ROC) curve was employed to evaluate the predictive capacity of the prognostic model using the “timeROC” package. To examine differences in overall survival between high-risk and low-risk groups, we conducted Kaplan–Meier analysis and performed the log-rank test using the “survival” package. To validate the prognostic model, we calculated risk scores for the additional GEO cohorts. In the validation cohorts, patients were categorized into high-risk and low-risk groups based on the optimal cut-off value, which was determined using the “survminer” package. Construction and validation of a nomogram Univariate and multivariate Cox regression analyses were performed on the risk scores and other potential clinical risk factors. Subsequently, we constructed a nomogram based on three key predictors: risk score, pathological T stage, and age. We calculated the nomogram score for each patient using the “nomogramFormula” package. The predictive capacity of the nomogram was assessed using ROC curves and calibration plots. Validation of ESM1 protein expression To investigate the distribution and expression levels of ESM1 protein in BC and normal bladder tissues, we obtained immunohistochemistry (IHC) staining images and relevant clinicopathological data from the Human Protein Atlas database ([73]https://www.proteinatlas.org/). Cell lines and reagents Human BC cell lines UMUC-3 and 5637 were purchased from Pricella Biotechnology Co., Ltd (Wuhan, China). These cell lines were cultured in DMEM (Gibco) supplemented with 10% fetal bovine serum (FBS). Recombinant human ESM-1 protein was purchased from MedChemExpress. EdU incorporation assay We performed the EdU incorporation assay to evaluate the proliferative capacity of BC cells using the Cell-Light EdU Apollo567 In Vitro Kit (RiboBio). Images were captured using a fluorescence microscope (ZEISS) at 200x magnification. The percentage of EdU-positive cells (with red fluorescence) relative to all Hoechst-labeled cells (with blue fluorescence) was calculated in each well. Transwell migration assay We conducted the Transwell migration assay to assess the effect of ESM1 protein on the migratory capacity of BC cells. BC cells (3 × 10^4) were resuspended in 300 µL of serum-free medium and plated into the upper compartment of 24-well Transwell inserts with an 8-µm pore size (Corning). The lower chamber was filled with 700 µL of medium containing 10% FBS and supplemented with recombinant human ESM1 protein at a final concentration of 100 ng/mL or with PBS as a control^[74]31,[75]32. After incubation for 24 h, BC cells were fixed with 4% paraformaldehyde solution and subsequently stained with 0.1% crystal violet solution for 15 min. Cells remaining in the upper compartment were removed using cotton swabs. Images of the migrated cells were captured at 200x magnification using a microscope (ZEISS). The number of migrated BC cells was quantified using ImageJ software. Statistical analysis In this study, the Student’s t-test and the Wilcoxon rank-sum test were employed to assess the differences in continuous variables between two distinct groups. All statistical analyses were conducted utilizing R software (version 4.2.0) and GraphPad Prism (version 9.5.1). The workflow of our research is illustrated in Fig. [76]1. Fig. 1. [77]Fig. 1 [78]Open in a new tab The workflow of this study. Results Single-cell analysis uncovers two distinct EC subpopulations in BC The xCell algorithm was applied to assess the abundance of ECs in BC patients from the TCGA cohort. Patients with higher EC levels exhibited poorer overall survival compared to those with lower EC levels (Supplementary Fig. 1A). To decipher the transcriptional and functional heterogeneity of distinct EC subpopulations in BC, we analyzed scRNA-seq data from CD45-negative stromal cells. After strict filtering, 7,519 high-quality CD45-negative single cells were included for further analysis. We identified twelve distinct clusters using an unsupervised clustering method and manually annotated four major cell types based on their marker genes (Fig. [79]2A-B; Supplementary Fig. 1B-F). A total of 705 endothelial cells were identified and categorized into four distinct clusters (Fig. [80]2C). Subsequently, hierarchical clustering analysis was performed, and the endothelial cells were categorized into two distinct subpopulations: S1 ECs and S2 ECs (Fig. [81]2D-E). S1 ECs exhibited elevated expression of angiogenesis-related genes, including ESM1, INSR, FLT1, and CXCL12 (Fig. [82]2F; Supplementary Table 3). S2 ECs displayed upregulated expression of immune-related genes, including CCL2, CXCL2, HLA-DRB1, and CD74 (Fig. [83]2F; Supplementary Table 4). Therefore, based on the expression of marker genes, S1 ECs were identified as the high ESM1-expression subpopulation, while S2 ECs were defined as the high CXCL2-expression subpopulation. To explore the biological functions of the two distinct EC subpopulations, KEGG and GO enrichment analyses were conducted on the marker genes of S1 and S2 ECs. KEGG analysis revealed that upregulated genes in S1 ECs were closely associated with PI3K-Akt signaling, focal adhesion, and extracellular matrix receptor interaction (Fig. [84]2G). S2 ECs exhibited significant enrichment in pathways related to immune-related diseases and pathogen infections (Fig. [85]2G). GO analysis showed that upregulated genes in S1 ECs were enriched in pathways associated with the proliferation and migration of epithelial and endothelial cells (Fig. [86]2H). Conversely, pathways related to antigen processing and presentation were markedly activated in S2 ECs (Fig. [87]2H). Fig. 2. [88]Fig. 2 [89]Open in a new tab Single-cell analysis identifying two distinct EC subpopulations in BC. (A) t-SNE plot of 7,519 CD45-negative single cells, color-coded by main cell types. (B) Dot plot displaying the average expression levels of marker genes in distinct cell types. (C) t-SNE plot of 701 endothelial cells, color-coded by clusters. (D) Dendrogram depicting hierarchical clustering of endothelial cells. (E) t-SNE plot of endothelial cells colored according to subpopulations. (F) Heatmap displaying the normalized average expression of the top 20 marker genes in two distinct EC subpopulations at the single-cell level. (G) Bar plot showing significantly enriched KEGG pathways in two distinct EC subpopulations. (H) Bar plot showing significantly enriched GO functional pathways in two distinct EC subpopulations. To further validate our findings, we integrated three additional independent scRNA-seq datasets that were not subjected to flow cytometry sorting. These datasets included a total of 16 bladder cancer and 8 normal bladder samples. Following quality control and filtering procedures, a total of 138,593 single cells were analyzed and classified into eight major cell types (Supplementary Fig. 2A). Among these, 6,146 endothelial cells were identified through the expression of specific marker genes, including PECAM1, VWF, and CLDN5 (Supplementary Fig. 2B-D). The endothelial cells, derived from both tumor and normal tissues, were further stratified into six distinct clusters (Supplementary Fig. 2E-F). The composition of endothelial cell clusters exhibited significant differences between tumor and normal tissues (Supplementary Fig. 2G). Specifically, cluster 2 constituted the highest proportion of endothelial cells in normal tissues, while cluster 0 was the most abundant in tumor tissues (Supplementary Fig. 2G; Supplementary Table 5). In contrast, the proportion of cluster 1 was similar between tumor and normal tissues (Supplementary Fig. 2G; Supplementary Table 5). The signature scores for the S1 and S2 EC subpopulations were computed using the ssGSEA algorithm, based on the top 20 marker genes for each subpopulation. The S1 EC signature score was significantly higher in cluster 0 compared to other clusters (Supplementary Figs. 2 H-I). Conversely, the S2 EC signature score was significantly higher in cluster 2 (Supplementary Figs. 2 K-L). Additionally, the S1 EC signature score was significantly elevated in tumor-associated endothelial cells (Supplementary Fig. 2J), while the S2 EC signature score was significantly elevated in normal endothelial cells (Supplementary Fig. 2M). Our analysis revealed that the six distinct endothelial cell clusters displayed significantly different transcriptomic profiles. The cluster 0 exhibited significant upregulation of S1 EC marker genes, including ESM1, IGFBP3, INSR, FLT1, MIR4435-2FG, and CA2 (Supplementary Fig. 2N). In contrast, S2 marker genes were significantly upregulated in clusters 1 and 2 (Supplementary Fig. 2N). Functional enrichment analysis demonstrated that highly expressed genes in cluster 0 were associated with angiogenesis and epithelial cell proliferation (Supplementary Fig. 2O). Conversely, highly expressed genes in clusters 1 and 2 were linked to immune response and antigen presentation processes (Supplementary Fig. 2O). The transcriptional and functional features of the high ESM1-expression EC subpopulation To examine the activity of various signaling pathways, we calculated AUCell scores for each cell. The activity scores of angiogenesis, NOTCH signaling, VEGF/VEGFR interaction, and Wnt signaling pathways were significantly increased in S1 ECs (Fig. [90]3A-B). In contrast, S2 ECs exhibited markedly elevated activity scores in immune response-related pathways, including T cell receptor signaling and activation, NK cell-mediated cytotoxicity, and B cell receptor signaling (Fig. [91]3A-B). Fig. 3. [92]Fig. 3 [93]Open in a new tab Comprehensive functional characterization of S1 and S2 EC subpopulations. (A) Lollipop chart showing differences in signaling pathway activity between S1 and S2 EC subpopulations. The Limma method was employed to compare activity scores between two distinct subpopulations and to compute the corresponding t-values. (B) Boxplots comparing AUCell scores between two distinct EC subpopulations, with statistical significance determined using a two-tailed Wilcoxon test. (*p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001.) (C–F) Pseudotime trajectory of endothelial cells, colored by pseudotime value (C), state (D), cluster (E) and subpopulation (F). (G) Heatmap illustrating dynamic changes in expression levels of DEGs between two EC subpopulations along the pseudotime trajectory. (H) Scatter plots displaying expression levels of DEGs between two EC subpopulations along the pseudotime trajectory. (I) Heatmap showing average regulon activity scores for the top five significantly activated regulons in S1 and S2 EC subpopulations, as estimated using the pySCENIC algorithm and AUCell analysis. (J) Violin plots comparing regulon activity scores between two distinct EC subpopulations, with statistical significance determined using a two-tailed Wilcoxon test. (*p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001.) To further explore the dynamic alterations in transcriptional profiles across distinct EC subpopulations, the pseudotime trajectory of ECs was constructed using the Monocle algorithm. This analysis showed that the S2 EC subpopulation was predominantly located at the initial stage of the pseudo-time trajectory, whereas the S1 EC subpopulation was mainly positioned at the middle and terminal stages (Fig. [94]3C-F). In addition, we observed that the expression levels of ESM1, IGFBP3, INSR, and SAT1 were significantly upregulated, whereas those of CD74, CXCL2, SELE, NFKBIA, ACKR1, and CCL2 were markedly downregulated as pseudo-time progressed (Fig. [95]3G-H). As shown in Fig. [96]3D, the pseudotime trajectory bifurcated into two main branches at the branch point 1, corresponding to states 3 and 4. To elucidate the differences in gene expression profiles between these two branches, we employed Branch Expression Analysis Modeling (BEAM) and categorized endothelial cells into three groups: pre-branch, state 3, and state 4 (Supplementary Fig. 3). The pre-branch group comprised cells extending from the starting point to branch point 1. Our analysis revealed that S1 EC marker genes, such as ESM1, IGFBP3, and FLT1, were significantly upregulated in state 4 compared to state 3 and were expressed at very low levels in the pre-branch group (Supplementary Fig. 3). State 3 exhibited significant upregulation of MALAT1, MCAM, HSPG2, and MT-ND5 compared to state 4 (Supplementary Fig. 3). In the pre-branch group, the expression levels of S2 ECs marker genes were significantly elevated, including SELE, CXCL2, and CD74 (Supplementary Fig. 3). Functional enrichment analysis demonstrated that genes upregulated in state 4 were associated with angiogenesis, cell proliferation and migration, and protein secretion (Supplementary Fig. 3). Genes highly expressed in state 3 were linked to cell mitosis and angiogenesis, whereas those upregulated in the pre-branch group were associated with immune response (Supplementary Fig. 3). To further validate our findings, we conducted pseudotime analysis on additional integrated scRNA-seq datasets, which included 16 bladder cancer and 8 normal bladder samples. As depicted in Supplementary Fig. 4A, endothelial cells with the smallest pseudotime values were defined as the origin of the pseudotime trajectory. Our analysis demonstrated that normal endothelial cells were predominantly located at the initial stage of the pseudotime trajectory, whereas tumor endothelial cells were primarily positioned at the middle and terminal stages of the pseudotime trajectory (Supplementary Fig. 4B). The pseudotime trajectory bifurcated into two primary branches at branch point 1, corresponding to states 2 and 3 (Supplementary Fig. 4C). Notably, S1 EC signature scores were significantly increased in state 2, whereas S2 EC signature scores were elevated in states 1 and 3 (Supplementary Figs. 4D-E). To further explore the differences between the two branches (states 2 and 3), BEAM analysis was conducted, and endothelial cells were categorized into three distinct groups: pre-branch, state 2, and state 3 (Supplementary Fig. 4F). Our analysis showed that S1 EC marker genes, such as ESM1, IGFBP3, and CA2, were significantly upregulated in state 2 (Supplementary Fig. 4F). State 3 exhibited substantial upregulation of immune-related genes, such as HLA-DRB1, HLA-DRA, and CD74 (Supplementary Fig. 4F). In the pre-branch group, the expression levels of S2 EC marker genes, including SELE, CCL2, and VCAM1, were significantly elevated (Supplementary Fig. 4F). Functional enrichment analysis revealed that genes upregulated in state 2 were associated with ECM remodeling and angiogenesis (Supplementary Fig. 4F). Genes highly expressed in state 3 were linked to the activation of immune responses, whereas those upregulated in the pre-branch group were associated with responses to bacteria and tumor necrosis factor (Supplementary Fig. 4F). To decipher the transcriptional regulatory patterns in the two distinct EC subpopulations, pySCENIC and AUCell analyses were employed to construct regulons and compute the activity scores of these regulons. The results demonstrated that the activity scores of TFs associated with tumor progression were remarkably increased in S1 ECs, including ETV5, TEAD4, HOXD10, ETS1, and ZNF384 (Fig. [97]3I-J). Conversely, inflammation-related TFs, such as ATF3, STAT1, IRF1, FOSB, and SPI1, were highly activated in S2 ECs (Fig. [98]3I-J). The abundance of the high ESM1-expression EC subpopulation is associated with prognosis and immune cell infiltration in BC To evaluate the expression levels of the S1 EC signature in tumor and normal ECs, we analyzed snRNA-seq data from MIBC and normal bladder samples. We annotated the primary cell types and identified endothelial cells using specific marker genes, including PECAM1 and VWF (Supplementary Fig. 5A-D). The S1 and S2 EC signature scores were calculated for each cell using the ssGSEA algorithm, based on the top 20 marker genes specific to each EC subpopulation (Supplementary Table 6). We observed that the S1 EC signature scores were markedly upregulated in MIBC ECs compared with normal ECs (Fig. [99]4A-C). To investigate the impact of distinct EC subpopulations on the prognosis and immune infiltration in BC, we computed the S1 and S2 EC signature scores for the TCGA and IMvigor210 cohorts. Kaplan–Meier analysis revealed that patients with higher S1 EC abundance exhibited markedly unfavorable overall survival compared with those with lower S1 EC abundance in both the TCGA and IMvigor210 cohorts (Fig. [100]4D-E). However, There is no significant difference in the objective response rate to immunotherapy between patients with higher and lower S1 EC abundance (Chi-square test p-value = 0.2844) (Fig. [101]4F). We observed no significant difference in overall survival between patients with higher and lower S2 EC abundance (Supplementary Fig. 6). The infiltration levels of various immune cells in the TME were assessed using the CIBERSORT, ssGSEA, and Xcell algorithms. The results of ssGSEA demonstrated that patients with higher S1 EC signature scores exhibited markedly elevated infiltration of macrophages, myeloid-derived suppressive cells (MDSCs), and regulatory T cells compared with patients with lower S1 EC signature scores (Fig. [102]4G). Moreover, CIBERSORT analysis showed that the proportion of M2-macrophages and neutrophils was significantly increased in patients with higher S1 EC abundance (Fig. [103]4H). Furthermore, patients with higher S1 EC abundance exhibited a significantly elevated abundance of M2 macrophages and regulatory T cells (Fig. [104]4I-J). These findings indicate that the high ESM1-expression EC subpopulation may facilitate recruitment of immune cells and contribute to the development of an immunosuppressive TME in BC. Fig. 4. [105]Fig. 4 [106]Open in a new tab S1 EC signature scores are associated with prognosis and immune cell infiltration in BC. (A) t-SNE plot of endothelial cells derived from MIBC and normal bladder tissues. (B) t-SNE plot of endothelial cells, with cells colored according to S1 EC signature scores. (C) Boxplot comparing S1 EC signature scores between tumor ECs and normal ECs, with statistical significance determined using a two-tailed Wilcoxon test. (*p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001.) (D) Kaplan–Meier curves showing that patients with higher S1 ECs signature scores exhibit shorter overall survival in the TCGA cohort. (E) Kaplan–Meier curves showing that patients with higher S1 EC signature scores exhibit shorter overall survival in the IMvigor210 cohort. (F) Bar plot displaying the objective response rate to anti-PD-L1 immunotherapy in patients with high versus low S1 EC signature scores. Statistical significance was assessed using a Chi-square test. (G) Comparison of ssGSEA scores for 28 immune cell types between patients with high and low S1 EC signature scores. (H) Comparison of the proportions of 22 tumor-infiltrating immune cell types, as estimated by CIBERSORT, between patients with high and low S1 EC signature scores. (I,J) Boxplots comparing abundance scores of M2 macrophages (I) and regulatory T cells (J), as assessed by the xCell algorithm. ESM1 is predominantly expressed in the endothelial cells of BC The snRNA-seq data showed that the expression of ESM1 was remarkably upregulated in MIBC ECs compared with normal ECs (Fig. [107]5A, C). Furthermore, ESM1 was identified as a common gene at the intersection of MIBC-EC-related DEGs, S1 EC marker genes, and ETV5 target genes (Fig. [108]5B). The expression of ESM1 was significantly higher in S1 ECs (Fig. [109]5D). To validate the localization and expression levels of the ESM1 protein, IHC staining results for ESM1 expression in BC and normal bladder tissues were obtained from the HPA database. IHC analysis demonstrated that ESM1 was predominantly expressed in endothelial cells rather than in tumor cells in BC tissues (Fig. [110]5E). Fig. 5. [111]Fig. 5 [112]Open in a new tab Validation of expression levels and biological function of ESM1 in BC. (A) Heatmap displaying expression levels of DEGs between tumor ECs and normal ECs from snRNA-seq data. (B) Venn diagram depicting the overlap of common genes among upregulated genes in MIBC ECs, S1 EC marker genes, and ETV5 target genes. (C) Violin plot and t-SNE plot showing expression levels of ESM1 in tumor ECs and normal ECs. (D) Violin plot and t-SNE plot showing expression levels of ESM1 in S1 and S2 EC subpopulation. (E) Immunohistochemical staining for ESM1 in normal bladder (n = 1) and BC tissues (n = 10). Arrows highlight abundant ESM1 expression in the vasculature. Scale bars: 200 μm (upper) and 50 μm (lower). (F) EdU incorporation assays in 5637 cells (upper panel) and UMUC-3 cells (lower panel) treated with PBS or recombinant human ESM1 protein (100 µg/ml) for 24 h. Scale bars: 100 μm. Statistical significance was determined using a two-tailed Student’s test. (*p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001.) (G) Transwell migration assays in UMUC-3 cells (upper panel) and 5637 cells (lower panel) treated with PBS or recombinant human ESM1 protein (100 µg/ml) for 24 h. Scale bar: 100 μm. Statistical significance was determined using a two-tailed Student’s test. (*p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001.) Recombinant human ESM1 protein promotes the proliferation and migration of BC cells To investigate the direct effects of the secretory proteoglycan ESM1 on the malignant phenotype of BC cells, we treated these cells with recombinant human ESM1. The concentration of ESM1 protein utilized in the present study was determined according to previous research^[113]31,[114]32. EdU incorporation assays showed that treatment with rhESM1 (100ng/ml) remarkably facilitated the proliferation of 5637 and UMUC-3 cells (Fig. [115]5F). Moreover, Transwell assays revealed that rhESM1 (100ng/ml) markedly promoted the migration of 5637 and UMUC-3 cells (Fig. [116]5G). Establishment and validation of the MIBC-EC-related prognostic model in BC Differential expression analysis was conducted on snRNA-seq data to compare gene expression profiles between tumor ECs and normal ECs. This analysis identified 244 DEGs, of which 117 genes exhibited significant upregulation and 126 showed significant downregulation in tumor ECs compared with normal ECs (Fig. [117]6A; Supplementary Table 7). Subsequently, we identified 55 genes that were closely correlated with overall survival in BC through univariate Cox regression analysis (Supplementary Table 8). To mitigate overfitting, we conducted LASSO regression analysis and identified 10 prognostic genes (Fig. [118]6B-C; Supplementary Table 9). Finally, we established a risk model using stepwise Cox regression analysis based on six prognostic genes, including RBP7, EVC2, SLC7A11, ELN, SCD, and TBX1 (Fig. [119]6D; Supplementary Table 10). We calculated the risk score using the following formula: Risk score = (RBP7 × 0.095) + (EVC2 × 0.190) + (SLC7A11 × 0.084) + (ELN × 0.120) + (SCD × 0.215) + (TBX1 × − 0.090). In the TCGA cohort, patients were categorized into high-risk and low-risk groups according to the median risk score. Kaplan–Meier analysis revealed that patients in the high-risk group exhibited remarkably shorter overall survival compared with those in the low-risk group (Fig. [120]6E). In addition, we generated ROC curves to assess the predictive accuracy of MIBC-EC-related risk score for overall survival (Fig. [121]6H). Furthermore, we validated the MIBC-EC-related prognostic model in two additional cohorts: [122]GSE13507 and [123]GSE31684. In these two cohorts, patients in high-risk group had markedly unfavorable overall survival compared with those in low-risk group (Fig. [124]6F-G). The predictive performance of the prognostic model was further assessed using ROC curves in the two validation cohorts (Fig. [125]6I-J). Fig. 6. [126]Fig. 6 [127]Open in a new tab Establishment and validation of a MIBC-EC-related prognostic model. (A) Volcano plot displaying DEGs between tumor ECs and normal ECs, with 117 genes upregulated and 126 genes downregulated in tumor ECs. (B) Visualization of LASSO regression analysis. The optimal parameter Lambda.min (left dotted line) was selected for the analysis. (C) Visualization of LASSO regression coefficient profiles. (D) Forest plot presenting six prognostic genes in the MIBC-EC-related risk model, identified through stepwise Cox regression analysis. (E–G) Kaplan–Meier curves showing that patients in the high-risk group exhibit shorter overall survival compared to those in the low-risk group in the TCGA (E), [128]GSE13507 (F), and [129]GSE31684 (G) cohorts. (H–J) Time-dependent ROC curves displaying the predictive performance of the MIBC-EC-related prognostic model for 1-, 3-, and 5-year overall survival in the TCGA (H), [130]GSE13507 (I), and [131]GSE31684 (J) cohorts. Construction and validation of a novel nomogram Univariate and multivariate Cox regression analyses demonstrated that the MIBC-EC-related risk score, pathological T stage, and age were independent predictors of overall survival in BC (Fig. [132]7A-B). Therefore, we constructed a novel nomogram incorporating these three independent prognostic risk factors (Fig. [133]7C). The predictive accuracy of the nomogram was evaluated using ROC curves in both the training and validation cohorts (Fig. [134]7D-E). Furthermore, the calibration curves were generated to assess the predictive capability of the nomogram in both the training and validation cohorts (Fig. [135]7F-G). Fig. 7. [136]Fig. 7 [137]Open in a new tab Construction and validation of a nomogram. (A,B) Forest plots showing the results of univariate (A) and multivariate (B) Cox regression analyses for MIBC-EC-related risk scores and other clinicopathological variables. (C) Construction of a nomogram incorporating three independent prognostic factors: MIBC-EC-related risk scores, pathological T stage, and age. (D,E) Time-dependent ROC curves displaying the predictive accuracy of the nomogram for 1-, 3-, and 5-year overall survival in the training (TCGA) and validation ([138]GSE31684) cohorts. (F,G) Calibration curves displaying the predictive accuracy of the nomogram for 1-, 3-, and 5-year overall survival in the training (F) and validation (G) cohorts. Discussion Recent advancements in scRNA-seq have shed light on the functional heterogeneity of tumor-associated endothelial cells and their contributions to tumor progression in lung and prostate cancers^[139]11,[140]33. A recent single-cell sequencing study has elucidated the immunosuppressive characteristics of the TME in urothelial carcinoma^[141]34. However, the phenotypic characteristics and biological functions of distinct EC subpopulations in BC have yet to be fully elucidated. In the present study, we identified a previously unrecognized EC subpopulation characterized by elevated ESM1 expression and activation of signaling pathways involved in angiogenesis and cancer cell proliferation. In addition, our analysis revealed that the signature score of the high ESM1-expression EC subpopulation was markedly increased in MIBC ECs compared with normal ECs. In line with our findings, a recent study in prostate cancer identified an activated EC subtype marked by elevated expression of cancer-associated fibroblast (CAF) markers and activation of signaling pathways implicated in ECM–receptor interaction and focal adhesion^[142]35. The tip EC subtype, characterized by upregulated expression of collagen remodeling-related genes, was closely associated with clinical prognosis and sensitivity to VEGF inhibitor therapy in lung cancer^[143]11. Moreover, a recent study identified an arterial EC subtype, characterized by elevated CXCL12 expression, as a potential therapeutic target for inhibiting angiogenesis in prostate cancer^[144]33. The interactions between ECs and tumor cells are critical for initiating tumor metastasis and angiogenesis. ECs have been shown to induce autophagy in tumor cells and facilitate metastasis via the activation of the CCL5/CCR5 paracrine signaling pathway in prostate cancer^[145]36. Increased expression and S-nitrosylation of VCAM-1 in ECs are critical for tumor cell adhesion to ECs and subsequent extravasation, thereby promoting breast cancer metastasis^[146]37. Moreover, ECs have been shown to enhance the proliferation and invasion of bladder cancer cells by secreting elevated levels of epidermal growth factor (EGF)^[147]38. A recent study has revealed that activation of Notch1 signaling in ECs enhances tumor metastasis by facilitating tumor cell transmigration across the vascular wall^[148]39. Consistently, our results indicate that ECs may promote the progression of bladder cancer by increasing the secretion of ESM1. Notably, our results demonstrate that patients with higher abundance of the high ESM1-expression EC subpopulation exhibit more unfavorable clinical prognosis compared to those with lower abundance. In addition, the signature scores of the high ESM1-expression EC subpopulation are positively correlated with the infiltration levels of inhibitory immune cells, including regulatory T cells and M2 macrophages. Mounting evidence suggests that ECs play paramount roles in modulating immune responses^[149]17. It has been reported that patients with elevated expression of DGKG in ECs exhibit unfavorable prognosis in liver cancer^[150]40. Elevated expression of DGKG in tumor-associated ECs facilitates angiogenesis and regulatory T cell differentiation, thereby diminishing the efficacy of combined immunotherapy and anti-angiogenic treatments^[151]40. Additionally, a recent study has identified a venous EC subtype with immunomodulatory functions, which influences the infiltration of immunosuppressive tumor-associated macrophages in breast cancer^[152]41. Moreover, tumor-associated endothelial cells can induce alternative macrophage polarization by secreting elevated levels of interleukin-6 (IL-6) and colony-stimulating factor-1 (CSF-1), thereby facilitating immune evasion and tumor progression in glioblastoma^[153]42. GRN analysis revealed that the activity of the transcription factor ETV5 was remarkably increased in the high ESM1-expression EC subpopulation. Additionally, ESM1 is a target gene of ETV5. These findings suggest that ETV5 may play a critical role in regulating ESM1 transcription and promoting disease progression in bladder cancer. ETV5 can promote angiogenesis by enhancing the expression and secretion of PDGFR-BB in colorectal cancer^[154]43. Moreover, elevated expression levels of ETV5 are closely associated with resistance to anti-angiogenic therapy with bevacizumab in colorectal cancer^[155]44. Our findings demonstrate that ESM1 is predominantly expressed in tumor-associated ECs rather than in tumor cells. Additionally, treatment with rhESM1 significantly enhances the proliferation and migration of BC cells in vitro. Consistent with our findings, it has been demonstrated that elevated expression of ESM1 in tumor-associated ECs is associated with reduced recurrence-free survival in patients with invasive bladder cancer^[156]45. Silencing ESM1 in tumor-associated ECs substantially inhibited angiogenesis driven by VEGF-A^[157]45. Notably, serum ESM1 concentrations are significantly elevated in patients with BC and are positively associated with more advanced pathological stages^[158]46. Downregulation of ESM1 suppresses tumor growth and metastasis in bladder cancer^[159]47. Furthermore, urine samples from patients with BC exhibited markedly elevated levels of ESM1 compared to those from healthy individuals, suggesting that ESM1 may serve as a promising biomarker for the noninvasive diagnosis of bladder cancer^[160]48. Recently, overexpression of ESM1 has been shown to facilitate the progression and metastasis of several malignancies, including cervical and gastric cancers^[161]49,[162]50. Recent studies have identified several endothelial cell-related prognostic signatures in BC^[163]51–[164]53. However, these studies have not elucidated the differences in gene expression between tumor-associated ECs and normal ECs. In this study, we conducted a comprehensive analysis of gene expression profiles in ECs derived from MIBC and normal tissues. Furthermore, we screened six prognostic genes and established a MIBC-EC-related prognostic model. Our findings suggest that this prognostic model may serve as a novel biomarker for predicting clinical outcomes in BC. The current study is subject to several limitations that remain to be fully addressed. First, we identified and characterized two distinct endothelial cell subpopulations in bladder cancer using scRNA-seq data from CD45-negative stromal cells. However, our results may be influenced by potential technical and biological biases, including sequencing depth, sampling variability, and batch effects. To mitigate these biases, we have employed rigorous quality control, data normalization, and batch-effect correction techniques. Additionally, we validated the findings using several independent scRNA-seq and single-nucleus RNA-sequencing datasets. Future research that integrates single-cell data across various disease stages and treatments, combined with spatial transcriptomics, could provide a more accurate exploration of the phenotypic and functional heterogeneity of endothelial cells in bladder cancer. Second, the clinical significance of the high ESM1-expression endothelial cell subpopulation in bladder cancer progression and immunotherapy response has been validated only in public datasets, including TCGA and IMvigor210 cohorts. To further confirm the generalizability of these findings, additional large-scale and prospective clinical cohorts are necessary. Third, we conducted only in vitro experiments to validate the influence of the ESM1 protein on bladder cancer cells. Future research should employ genetically engineered models, such as endothelial cell-specific ESM1 conditional knockout mice, to validate the impact of high ESM1-expression endothelial cells on bladder cancer progression. In conclusion, our investigation revealed that the high ESM1-expression EC subpopulation is intricately associated with clinical outcomes and immune characteristics in bladder cancer. These findings suggest that targeting ESM1 in tumor-associated endothelial cells may serve as a promising strategy for treating bladder cancer and enhancing anti-tumor immunity. In the future, well-designed animal experiments and clinical trials are essential to assess its potential as a therapeutic intervention. In addition, our results demonstrated that the ESM1 protein could enhance the proliferation and migration of bladder cancer cells in vitro. Accordingly, further cellular and molecular experiments are necessary to elucidate the mechanisms through which ESM1 drives bladder cancer progression. Furthermore, we have constructed and validated a prognostic signature associated with MIBC endothelial cells. Additional cohorts are warranted to further assess its potential as a prognostic biomarker in bladder cancer. Electronic supplementary material Below is the link to the electronic supplementary material. [165]Supplementary Material 1^ (1.8MB, pdf) Author contributions QH, KS, and BHL contributed to conceptualization and design. YFQ contributed to the acquisition of publicly available data, bioinformatics analysis and interpretation. YFQ and YHW performed cell experiments and statistical analyses. YFQ, YHW, and JHL wrote the manuscript. YFQ, QH reviewed and edited the manuscript. The manuscript was reviewed and approved by all authors. Funding This work was supported by HaiYa Young Scientist Foundation of Shenzhen University General Hospital; the National Natural Science Foundation of China (82002716); Shenzhen Scientific Research Foundation for Excellent Returned Scholars (000493); Natural Science Foundation of Shenzhen University General Hospital (SUGH2020QD005); Disciple gathering teaching project of Shenzhen University; Shenzhen Key Laboratory Foundation (ZDSYS20200811143757022); Teaching Reform Research Project of Shenzhen University (YXBJG202339); Shenzhen International Cooperation Research Project (GJHZ20220913143004008); Educational Reform Research Project of Clinical Teaching Base of Guangdong Province (2021JD207); Medical and Health Science and Technology Project (support) of Longgang District, Shenzhen (LGKCYLWS2021000032); Key Project of Medical and health technology of Longgang District, Shenzhen (LGKCYLWS2022034); Shenzhen Municipal Commission of Science and Technology Innovation (JCYJ20220818100009020, JCYJ20220818100016035). Data availability The scRNA-seq data ([166]GSE211388, [167]GSE129845, [168]GSE135337, [169]GSE222315) and snRNA-seq data ([170]GSE169379) are publicly available from the GEO database ([171]https://www.ncbi.nlm.nih.gov/geo). The bulk RNA-seq data used in this study include the TCGA-BLCA datasets ([172]http://xena.ucsc.edu/), the IMvigor210 datasets ([173]http://research-pub.gene.com/IMvigor210CoreBiologies), and the Affymetrix microarray datasets ([174]GSE31684, [175]GSE13507, [176]GSE32548). Declarations Competing interests The authors declare no competing interests. Footnotes Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Contributor Information Baohua Liu, Email: ppliew@szu.edu.cn. Kai Sun, Email: Henrysk@163.com. Qi Hou, Email: qi_hou@foxmail.com. References