Forty paraffin histopathological sections were obtained from the pathology department of the First Affiliated Hospital of Guangxi Medical University. All participants had not received any tumor therapy before enrolment, and the diagnosis of MIBC was confirmed by two experienced pathologists. The N, L and T region were delineated by the pathologist. Written informed consent was provided by all participants. This study was approved by the ethics committee of the First Affiliated Hospital of Guangxi Medical University (2024-E443-01). The following cohort were used for pan-cancers analysis: TCGA-BLCA, TCGA- MESO, TCGA- GBM, TCGA-LGG, TCGA-STAD, TCGA-UVM, TCGA-KIRC, TCGA-KIRP, TCGA-LUSC, GSE32894, GSE188715. The bulk RNA-seq data and clinical information of TCGA datasets can be downloaded from https://portal.gdc.cancer.gov/. As for GEO datasets, they can be found at https://www.ncbi.nlm.nih.gov/geo/.

IMC was performed as previously described (17). Briefly, after dewaxing the tissue sections (thickness of 5µm), antigen retrieval was performed using Tris-EDTA antigen retrieval solution (FC16FA0005, Sangon Biotech (Shanghai) Co.,Ltd.,China). The sections were then blocked with a 3% BSA solution at room temperature and incubated overnight at 4°C with antibody cocktail, as shown in Supplementary Table S1 . The next day, the sections were stained with Cell-IDTM Intercalator-Ir (Fluidigm) for DNA and subjected to analysis.

The MCD viewer software was used to open the files generated by the IMC system, and 16-bit OME-TIFF files were exported. The 16-bit OME-TIFF files were imported into the CellProfiler software, where channel information was set and signal intensities of each channel were extracted. The DNA (Ir191 and Ir193) and cell surface markers were used to identify the cell nuclei and cell boundaries, respectively. Single cells were then segmented and TIFF files were exported. The TIFF files exported from CellProfiler were imported into the HistoCAT software. The phenograph algorithm was applied for clustering, while the tSNE algorithm was used for high-dimensional data reduction. Neighbor analysis was employed to examine the interactions between different cell clusters. The immune cells and stromal cells were gated in cytobank (https://community.cytobank.org) by specific markers: stromal cell (Collagen I, vimentin), immune cell (CD45), B cell (CD20), macrophage (CD68), CD8+ T cell (CD3, CD8) and CD4+ T cell (CD3, CD4).

During cell division, the X/Y coordinates of each cell type were determined. Similarly to previously reported methods (18), we utilize the euclidean distance to calculate the closest distance between each cell and different cell types. On a per-ROI basis, we then computed the average closest distance between cell types within each ROI. Cellular neighborhood identification: For the 1,303,975 cells in our experiment, we set the window size to 26, consisting of one central cell and its 25 nearest neighbors, determined by euclidean distance in X/Y coordinates. We chose a 25-cell radius as a rough approximation, where the central cell distance in each direction is visually determined to be a good indicator of local functional activity.

The scRNA-seq data from eight bladder cancer patients were downloaded from the cited literature (19). As a quality control measure, we excluded sequenced cells that had less than 200 detected genes and more than 10% mitochondrial genes. Additionally, genes expressed in fewer than three cells were also excluded. We used the FindClusters function in Seurat (version: 4.0.1) to identify 35 clusters and visualize them using a t-distributed stochastic neighbor embedding (t-SNE) plot. Cluster-specific genes were used to annotate cell types with classic markers described in previous studies (19): tumor epithelial cells (EPCAM+/KRT17+); endothelial cells (PECAM1+/ENG+); fibroblasts (COL1A1+/RGS5+); myeloid cells (LYZ+/CD68+); T cells (CD3D+/CD2+); B cells (CD79A+/MZB1+); and mast cells (TPSAB1+/CPA3+).The “ComplexHeatmap” R package was used to visualize the average expression values of marker genes for each cell cluster.

For the analysis of interactions between cell types, we extracted the expression profiles of each cell type from eight bladder cancer tissues. The expression matrix, along with cell type metadata, was used as input for running cellphonedb with the statistical analysis method in Python. This study detected statistically significant interactions (p < 0.05) between receptors and ligands by comparing them against a curated L-R database. The “ktplots” R package was used to visualize the significant L-R pairs identified by cellphonedb. Transcription factor prediction of S3 up-regulated gene was performed by BART, an online analytical tool (http://bartweb.org/).

Bladder cancer ST data can be obtained from published literature (13). The “SpaCET” R package was utilized to compute the proportion of tumor cells and other cell types for each spot. Spots with more than 60% tumor cells, FB, or endothelial cells were designated as tumor region, FB region, or endothelial region, respectively (20). We estimated S3 enrichment in spots by scoring each spot of the ST data using S3 up-regulated genes. Positive L-R spots need to meet the following criteria: 1. They should be located at the boundary between the FB and endothelial region. 2. The ligand and receptor should be queryable in the cellponedb database. 3. The spot should express both the corresponding ligand and receptor genes simultaneously.

The gene sets used for ST scoring can be obtained from the GSEA official website (https://www.gsea-msigdb.org/gsea/downloads.jsp). The GSVA package (version: 1.38.2) was employed to calculate the enrichment score for each spot.

Based on the expression of specific markers (Collagen1, CD90, Vimentin, and α-SMA), FB can be classified into FB S3, FB S4, and FB SX. “FindMarkers” was utilized to perform differential gene expression analysis for FB subtypes, with logFC. threshold >= 0.25 and p_val_adj < 0.05 considered as significantly differentially expressed genes. Further filtering of the characteristic genes specific to FB S3 was performed using xgboost.

Statistical analysis was conducted using two-sided Student’s t-test and Kruskal-Wallis rank sum test. A significance level of p<0.05 was considered statistically significant. Kaplan-Meier analysis was used to estimate the overall survival curves of patients with MIBC, and the log-rank test was employed to compare the curves. GraphPad Prism 8 software was utilized for statistical analysis.

To explore the complex spatial characteristics of the MIBC microenvironment and enhance our understanding of its role in promoting metastasis. We performed IMC with a panel of 33 TME-associated markers to detect 100 tumor core regions (referred to as T), 62 leading-edge regions (referred to as L, the boundary of tumor cell invaded into muscular layer), and 23 nontumor regions (referred to as N) from paraffin-embedded tissue sections of 40 patients with MIBC, yielding 6,105 high-quality images ( Figure 1 and Supplementary Figure S1 ). Through cellular identification and segmentation of these images, we obtained single-cell expression profiles of 33 proteins ( Figures 2A, B ). Segmentation obtained data from 185 ROIs, ranging from 1562 to more than 10,000 cells per ROI (median 7317 cells/ROI; Supplementary Figure S2 ).

We found that the spatial structure of N, L, and T regions showed significant heterogeneity ( Figures 2C, D ). By cell clustering and spatial reconstruction, we found that cells with similar phenotypes cluster nearby within the T and N regions. On the other hand, cells in the L region exhibit a more disorganized distribution, suggesting a more complicated microenvironment in the L region compared with that in the T and N regions ( Figure 2E ). The deep analysis using Cytobank (an online analytic platform) revealed that the frequencies of stromal cells decreased from N to L to T regions progressively. The frequencies of immune cells including B cells, macrophages, CD4+ T cells, and CD8+ T cells increased in the L region than in the N and T regions ( Figure 2F ). Altogether, these results suggest that the microenvironments of the N, L, and T regions are spatially heterogeneous in MIBC.

Subsequent analyses were performed to determine the spatial heterogeneity of the microenvironment in N, L, and T regions. In the N region, numbers of B cells were recruited, whereas T cells and macrophages were less prevalent. The L regions displayed a higher density of both B cells and T cells, which were positioned near each other, with macrophages more dispersed near the tumor cells. The T regions showed a decreased number of immune cells, which were primarily located in the stromal areas rather than infiltrating the cancer nests ( Figure 3A ). The collagen structures in these regions also exhibited regional heterogeneity. The L regions were characterized by disordered collagen configurations, and some collagen appeared as dots ( Figure 3B ). Collagen was observed as curved structures in the N regions, loosely encasing small segments of muscle tissue ( Figure 3C ). Within the T regions, the following four unique types of collagen structures were identified: curved collagen encapsulating the tumor nests, parallel collagen aligning along the cancer nest margins, directionally arranged collagen, and disorderly arranged collagen with unclear boundaries relative to the tumor ( Figure 3D ). Additionally, tumor budding cells were found in the invasive region, which were dispersed within the tumor stroma as single cells or small clusters ( Figure 3E ).

To explore the important components of the TME involved in the recurrence and metastasis of MIBC, we clustered cells from 100 T regions, 62 L regions, and 23 N regions into 100 cell clusters with distinct phenotypes ( Figure 4A ). Then, we consolidated these 100 clusters into 17 groups based on the expression of specific markers for tumor cells, immune cells, stromal cells, vascular endothelial cells, and muscle cells. Among these were the following four unique FB clusters: FB S1 (α-SMA− Collagen I+ Vimentin- CD90+), FB S2 (α-SMA+ Collagen I+ Vimentin- CD90+), FB S3 (α-SMA+ Collagen I+ Vimentin+ CD90+), and FB S4 (α-SMA- Collagen I+ Vimentin+ CD90+) ( Figure 4B ). These FB clusters exhibited unique spatial distributions. The L regions predominantly contained two types of FB, with nearly 75% being FB S3, referred to as S3 hereafter. Furthermore, all four FB clusters were present in the T regions, whereas the N regions exclusively contained FB S4 ( Figure 4C ).

The L regions, a critical boundary for tumor invasion and metastasis, emerged as the most dynamic area, with S3 as the predominant FB type. This suggests that S3 exhibits high activity and functional characteristics, facilitating tumor progression. Survival analysis showed that MIBC patients with S3-positive tumors had shorter overall survival than those with S3-negative tumors ( Figures 4D, E ), indicating a crucial role of S3 in MIBC progression. Investigating the TME heterogeneity between S3-positive and S3-negative tumors, we found that S3-positive tumors exhibited a higher density of blood vessels, which included more extensive and severely distorted blood vessels with abnormal dilation. Furthermore, these tumors showed high expression of VEGFR on vessel walls and outward expansion ( Figures 4F, G ), suggesting that S3 may contribute to abnormal angiogenesis. Subsequent investigation using IMC spatial analysis determined the interactions between S3 and endothelial cells. By evaluating spatial coordinates, we determined the distances from each FB cluster to the nearest endothelial cells and found that the proximity of S3 to endothelial cells was shorter compared with other FB types ( Figures 4H–J ). Neighbor component analysis among the four FB clusters showed that endothelial cells contained 2.38%, 1.43%, 8.39%, and 2.79% of the neighbors for FB S1, S2, S3, and S4, respectively ( Figure 4K ). These findings reveal an important close interaction between S3 and endothelial cells, highlighting the significant effect of S3 on MIBC angiogenesis.

To elucidate the mechanisms underlying S3, we analyzed publicly available scRNA-seq data from eight bladder cancer samples. Through cell clustering and annotation, we found 35 unique cell clusters characterized by their specific markers, which included B cells, endothelial cells, tumor epithelial cells, fibroblasts, mast cells, myeloid cells, and T cells ( Supplementary Figures S3A, B ). Subsequent refinement of our analysis by reclustering the fibroblasts yielded 14 fibroblast clusters ( Supplementary Figure S3C ). Using the expression of stromal cell-specific markers, including Collagen 1, α-SMA, Vimentin, and CD90, we identified the previously characterized FB clusters S3 and S4 from IMC analysis and labeled the remaining FB clusters as FB SX ( Figures 5A, B ). Differential expression analysis showed unique molecular profiles among these FB types ( Supplementary Figure S3D ). The gene expression upregulated in S3 was enriched in pathways related to angiogenesis, cell migration, epithelial cell proliferation and differentiation, endothelial cell proliferation, and fibroblast proliferation ( Figure 5C ). Analysis of cellular communication among clusters via single-cell L–R interaction showed that S3 engaged in more L–R pair interactions with other cells compared to FB S4 and SX. Specifically, S3 exhibited a significantly higher number of L–R pairs with endothelial cells, showing strong communicative activity ( Figures 5D, E ). These findings indicate that S3 maintains active communication with different cell types and also engages in close interactions with endothelial cells.

To elucidate the mechanisms between S3 and endothelial cells, we analyzed L–R pairs unique to S3 and endothelial cells ( Figure 5F ). While scRNA-seq offers insights into the cell–cell interactions based on L–R gene expression, it does not account for cell localization. Spatial transcriptomics (ST) analysis, which profiles spatial gene expression within tissues, allows the detection of cell–cell interactions across different regions. By analyzing publicly available ST data from bladder cancer, we identified specific regions corresponding to tumor cells, FB, and endothelial cells. Furthermore, the L–R pair NOTCH1–JAG2 was co-located in the endothelial spots adjacent to the FB region, indicating that S3 may predominantly interact with endothelial cells through the NOTCH1–JAG2 pathway ( Figure 5G ).

Further investigation into the regulatory mechanisms of S3, we employed the BART transcription factor prediction tool to analyze transcription factors associated with genes upregulated in S3 relative to other FB clusters. YAP1 emerged as one of the top predicted transcription factors ( Figure 6A ). IMC imaging confirmed high YAP1 expression in the microenvironment where S3 is located, contrasting with lower expressions in microenvironments dominated by FB S1, S2, and S4 ( Figures 6B, C ), suggesting a potential regulatory link between S3 and YAP1. To further determine the spatial characteristics and mechanisms associated with S3, we analyzed ST data from bladder cancer. ssGSEA results showed that S3-enriched samples showed higher angiogenesis signature scores, consistent with increased YAP1 expression compared with low-S3-enriched samples ( Figures 6D–G ). Altogether, these findings show that S3 may promote angiogenesis and could be regulated by YAP1.

To investigate the relationship between the presence of S3 and poor prognosis in MIBC, we utilized the XGBoost algorithm to analysis the gene expression profiles of S3 in scRNA-seq and construct an S3-related signature ( Supplementary Table S2 ). Using this signature, we calculated the Fibroblast Cluster S3 Score (FCS3S) in the TCGA-BLCA cohort with the gene set variation analysis (GSVA) algorithm. Our results revealed that patients with high FCS3S had significantly poorer overall survival ( Figure 7A ). This finding was supported by two independent BLCA cohorts (GSE188715 and GSE32894) ( Figures 7B, C ). Subsequent analysis showed that patients with high FCS3S exhibited increased activation of pathways associated with malignant tumor progression, including angiogenesis, tumor stemness, cellular exosomes, immune suppression, extracellular matrix (ECM) remodeling, neuroendocrine differentiation, and epithelial–mesenchymal transition (EMT) in the TCGA-BLCA cohort, indicating a higher degree of malignancy in these patients ( Figure 7D ). ESTIMATE analysis showed that patients with high FCS3S had increased immune and stromal scores, as well as overall ESTIMATE scores, compared to those with low FCS3S ( Figure 7E ). Notably, the immune suppression score was significantly elevated in high FCS3S patients, suggesting greater immune suppression in this group ( Figure 7D ). Furthermore, we evaluated the prognostic value of FCS3S across various cancers and found that it can effectively identify patients with poor outcomes and the following eight cancer types: glioblastoma multiforme (GBM), kidney renal clear cell carcinoma (KIRC), kidney renal papillary cell carcinoma (KIRP), lung squamous cell carcinoma (LUSC), uveal melanoma (UVM), lower-grade glioma (LGG), stomach adenocarcinoma (STAD), and mesothelioma (MESO) ( Figure 7F ). To summarize, our findings verified that S3 is complicated and related to the progression of MIBC, and FCS3S serves as a robust tool for predicting patient outcomes across different cancer types.