Abstract Bladder carcinoma (BLCA) is characterized by a high rate of post-surgery recurrence and multifocality. Multifocal tumors have a higher risk of recurrence compared to single tumors, significantly impacting bladder cancer-specific mortality. However, the interregional or intraregional heterogeneity within both primary and recurrent tumors remains poorly understood. Here, we employed single-cell RNA sequencing to analyze tumor lesions from five multifocal bladder cancer patients comprising three primary tumors and two recurrent tumors. Our findings revealed that malignant cells derived from recurrent multifocal bladder cancer exhibited higher interregional transcriptional similarity and consistent cellular communication. Furthermore, our analysis uncovered that malignant cells from recurrent tumors may evade immune destruction by suppressing cytokine responses and natural killer cell activity. Notably, we identified a preference for the expression of the tryptophan metabolic enzyme IL4I1 on SPP1^+ macrophages in recurrent tumors. Functional analyses have revealed that IL4I1 may promotes tumor progression in recurrent tumors by activating the aryl hydrocarbon receptor (AHR) and recruiting regulatory T cells to suppress adaptive immunity. Taken together, our study provides a comprehensive understanding of primary and recurrent multifocal bladder tumors, offering valuable resources for analyzing the multifocality and recurrence of bladder cancer. Subject terms: Cancer microenvironment, Data mining, Bladder cancer __________________________________________________________________ By shedding light on the complexities of primary and recurrent multifocal bladder tumors, this research provides valuable insights and resources for further exploration of the multifocality and recurrence of bladder cancer. Introduction Bladder cancer (BLCA) is a significant health concern worldwide^[36]1. Non-muscle-invasive bladder cancer (NMIBC) accounts for 70%–75% of these cases and is typically treated through surgical resection^[37]2,[38]3. Despite efforts to reduce recurrence and progression using intravesical chemotherapy or bacillus Calmette-Guérin infusion after tumor removal, a substantial number of patients still experience disease recurrence and progression to muscle-invasive bladder cancer (MIBC), especially those with multifocal tumors^[39]4,[40]5. Multifocal tumors carry a 40% higher risk of recurrence compared to single tumors, and distinguishing them from recurrence based on clinical symptoms alone can be challenging^[41]6. Therefore, there is a critical need for more precise diagnostic and treatment approaches. Previous multi-level analysis of human cancer tissue unveils a vast degree of molecular heterogeneity among cells and molecules within each tumor lesion^[42]7,[43]8. Investigating the tumor microenvironments (TME) of multifocal and recurrent tumors may contribute to unveiling the variation between multifocal tumors within primary or recurrent tumors, and will aid in the development of more precise diagnostic and treatment approaches. Single-cell RNA sequencing (scRNA-seq) has emerged as a powerful technology for investigating the intratumoral heterogeneity of cell composition and communication within the tumor microenvironment (TME)^[44]9–[45]11. Recent studies utilizing scRNA-seq to explore bladder cancer (BLCA) have uncovered the diverse nature of the TME and identified specific cell types that play crucial roles in either promoting or inhibiting tumor progression, as well as influencing immunotherapy response^[46]12. For example, the combination of scRNA-seq and T-cell receptor (TCR) sequencing has shed light on the existence of distinct tumor-specific and cytotoxic CD4^+ T-cell states capable of selectively eradicating autologous tumors through MHC class II-dependent mechanisms^[47]13. This finding provides predictive value for anti-PD-L1 therapy in bladder cancer^[48]14,[49]15. Furthermore, stromal cells, such as SLC14A1^+ cancer-associated fibroblasts (CAFs), have been implicated in promoting BLCA stemness through the paracrine pathway involving WNT5A^[50]16. While previous studies have reported various cell types involved in bladder cancer progression and immunotherapy response, they have often overlooked the heterogeneity of malignant cells between different lesions in primary or recurrent tumors^[51]17,[52]18. The evolution and heterogeneity of malignant cells can contribute to intratumor heterogeneity by exploiting different cell types within the tumor tissue to establish a supportive microenvironment for tumor growth or by competing for nutrients and inducing metabolic reprogramming in non-cancer cells, thereby facilitating tumor recurrence^[53]19–[54]21. Therefore, it is crucial to unravel the intratumor and intertumoral heterogeneity underlying primary and recurrent tumors to gain a comprehensive understanding of bladder cancer. In this study, we analyzed distinct tumor lesions from five bladder patients, encompassing three primary and two recurrent cases. We employed scRNA-seq to examine cellular heterogeneity and intercellular crosstalk within the primary and recurrent tumor microenvironments. Our findings indicate that the expression patterns and communication networks of malignant cells derived from distinct lesions are more similar in recurrent tumors compared to primary tumors. Additionally, through a comparison of primary and recurrent-derived malignant cells, we observed distinct transcriptomic and metabolic features in recurrent malignant cells when compared to their primary counterparts. Finally, we discussed the characteristic differences in immune cells between recurrent and primary tumors. Overall, our study significantly advances our understanding of the variations within the tumor microenvironment of primary and recurrent bladder tumors, offering valuable insights into the identification of recurrent biomarkers and potential therapeutic targets. Results Malignant cells in recurrent patients exhibited higher interregional similarity To investigate the heterogeneity of tumor microenvironment in multifocal bladder cancer between primary and recurrent tumors, we collected two separate surgical tumor specimens from five treatment-naive patients with multifocal bladder cancer. This included three primary and two recurrent patients, and the specimens were subjected to scRNA-seq analysis (Fig. [55]1A). Detailed clinical and pathological information, such as tumor grade, TNM stage, and tumor site, can be found in Table [56]S1. In total, we analyzed 10 samples, and the distinct lesion-derived tumor tissues from one patient are referred to as R1 and R2. After implementing quality control measures and removing doublets, we identified 121,554 high-quality cells and annotated them into 10 major cell types (Fig. [57]1B, Supplementary Fig. [58]1A–C). Notably, we observed that malignant cells derived from recurrent patients formed patient-specific clusters (Supplementary Fig. [59]1D). In contrast, malignant cells from primary tumors exhibited greater diversity and were clustered based on both regions and patients, indicating variations in interregional heterogeneity between recurrent and primary tumors, as well as the potential presence of batch effects. To reduce technical noise, we corrected for these batch effects and found that the cells were grouped according to cell type (Fig. [60]1C and Supplementary Fig. [61]1E). Additionally, the proportion of malignant cells varied among different patients but showed larger similarity within regions derived from the same patient, particularly in the case of recurrent patients (Fig. [62]1D). Furthermore, hierarchical clustering analysis of the regions revealed that malignant cells from the same patients tended to cluster together, indicating a smaller interregional heterogeneity compared to intertumoral heterogeneity (Fig. [63]1E). Fig. 1. Multilesional single-cell transcriptome profiling of primary and recurrent bladder cancer. [64]Fig. 1 [65]Open in a new tab A Workflow of multilesional tumor tissue collection, processing, scRNA-seq, and data analysis, Image created with BioRender.com, with permission. B Heatmap shows the marker genes expression of all cell types. C Uniform manifold approximation and projection (UMAP) plots of all malignant and epithelial cells, colored by the cell type annotation. D Histogram indicating the proportion of cells in each sample. E Hierarchical clustering of malignant cells from each tumor lesion across all samples. F The distribution of pair-wise correlations of malignant cells within each tumor lesion (intraregion), across tumor lesions within each patient (interregion) and across patients (intertumor). Pearson’s correlation coefficient was applied. Solid and dashed gray lines indicate the mean and standard deviation of all intraregional correlation values. G Boxplot showing the correlation of malignant cells across tumor lesions within patients. H UMAP of non-malignant cells, colored by cell type. I Boxplot showing the correlation of each cell type cells across tumor lesions within patients. J Histogram indicating the proportion of each cell type in each sample. To quantitively assess the heterogeneity of malignant cells, we calculated the intraregional similarity (IAS) score by examining the pairwise correlation distribution of malignant cells within each specific region^[66]7. Additionally, we determined the interregional similarity (IRS) score by measuring the correlation of malignant cells between two regions within the same patient. For comparison, we also calculated the intertumoral similarity (ITS) score among five patients. Remarkably, the tumor regions within recurrent patients exhibited higher IAS scores compared to primary tumors, while regional differences within each tumor still existed (Fig. [67]1F). Consistent with high IAS scores, malignant cells from recurrent patients demonstrated significantly higher IRS scores than those from primary tumors (Fig. [68]1G). The ITS scores were considerably lower than the IAS and IRS scores (Fig. [69]1F). In conclusion, recurrent tumors exhibited higher intraregional and interregional similarity in terms of tumor biopsies compared to primary tumors, suggesting the presence of potentially shared tumor clones that may contribute to recurrence. Variations of the non-malignant cells in multifocal primary and recurrent bladder cancer To investigate the variation among non-malignant cells, we visualized all cells by uniform manifold approximation and projection (UMAP). Interestingly, while malignant and epithelial cells were primarily clustered by sample or patient, immune and stromal cells exhibited distinct groupings based on their cell types. We further annotated the non-malignant cells into 8 major cell types, including T-cells, B-cells, plasma cells, myeloid cells, mast cells, endothelial cells, fibroblast, and myofibroblast cells using specific cell type marker genes (Fig. [70]1B, H, and Supplementary Fig. [71]1F). Next, we evaluated the interregional similarity of non-malignant cells separately in primary and recurrent patients. Notably, mast cells displayed a remarkably high IRS score across different tumor regions within a patient, suggesting a consistent and stable transcriptomic status among diverse tumor regions (Fig. [72]1I). In contrast, plasma cells exhibited the lowest IRS score compared to all other cell types, indicating substantial variation among different sampling regions within a single patient (Fig. [73]1I). These findings align with previous studies that reported higher variations among B lymphoid cells in different cases^[74]7. Interestingly, the observed variations were primarily attributed to plasma cells rather than B-cells, which displayed high IRS scores in both primary and recurrent tumors. Furthermore, we compared the cellular composition within each tumor region among different patients. Again, distinct tumor regions within each patient exhibited similar cellular patterns of non-malignant cells, except for plasma cells (Fig. [75]1J). Taken together, our findings indicated that plasma cells exhibited higher proportional and functional variation among tumor regions, while primary and recurrent samples shared similar major immune cell types. Distinct transcriptomic and metabolic features of malignant cells between primary and recurrent tumors Next, we sought to investigate the transcriptomic variations of malignant cells between primary and recurrent tumors. We characterized these cells based on their origins and conducted a differential expressed analysis across patients. Our functional enrichment analysis revealed that malignant cells from primary tumors exhibited upregulated genes associated with the P53 pathway, leukocyte trans-endothelial migration, and regulation of lymphocyte activation (Fig. [76]2A). Strikingly, malignant cells in recurrent tumors not only displayed heightened expression of genes linked to hypoxia, apoptosis, and the cell cycle, but also exhibited elevated levels of antigen processing and presentation genes, and interferon-gamma response pathway. This finding is consistent with previous studies that reported a subset of tumor cells exhibiting an epithelial-immune dual feature in nasopharyngeal carcinoma^[77]22,[78]23. Fig. 2. Deciphering expression programs revealed the epithelial-immune dual feature of malignant cells in recurrent tumors. [79]Fig. 2 [80]Open in a new tab A Two-sided bar graph showing the enriched pathways in primary and recurrent tumor-derived malignant cells. B Hierarchical clustering of intra-tumor expression programs, defined by Non-negative factorization (NMF), to identify the meta-programs based on the Jaccard index. C Heatmap showing the enriched pathways for each meta-program. D Boxplot showing the signature scores of meta-programs among malignant cells from primary and recurrent tumors. E Boxplot showing the expression of GATA3 in malignant cells from primary and recurrent tumors. F Venn diagram illustrating the shared genes between IFN response-associated genes and the upregulated gene of recurrent tumor-derived malignant cells. G Immunofluorescence images showing the localization of CD74 and the malignant cells marker gene CK-P in recurrent bladder tumors. Scale bar, 0.05 mm. H Boxplot showing the expression of CD74 across epithelial and malignant cells. I The rank of driver transcription regulators in the CD74^+ and CD74^- malignant cells. The regulators are ranked by the regulatory importance from the SCENIC result. J Boxplot showing the expression of STAT1 in epithelial and malignant cells. Furthermore, our analysis uncovered significant enrichment of metabolic-associated pathways in malignant cells derived from recurrent tumors, such as glycolysis and gluconeogenesis, indicating metabolic alterations between primary and recurrent tumors. To validate these metabolic changes, we estimated the metabolic activity of epithelial and malignant cells using gene set variation analysis (GSVA)^[81]24. Consistent with prior research, we observed that all malignant cells exhibited specific enrichments in energy production-associated metabolic pathways, particularly oxidative phosphorylation^[82]25,[83]26. This observation indicates an increased energy demand to support proliferation and disease progression (Supplementary Fig. [84]2A). Despite the overall upregulation of energy production-associated pathways in malignant cells compared to normal epithelial cells, we also detected variations in energy supply within different tumor regions of individual patients. Specifically, our differential metabolic activity analysis revealed that malignant cells from recurrent patients upregulated nitrogen metabolism-associated genes, such as hypoxia-sensor genes CA9 and CA12, which aligns with the significant enrichment of the hypoxia pathway in these cells^[85]27,[86]28 (Supplementary Fig. [87]2B–D). In summary, our findings highlight distinct transcriptomic and metabolic features of malignant cells in recurrent tumors. These cells demonstrate a stronger association with hypoxia and exhibit a notable epithelial-immune dual feature. Recurrent malignant cells exhibit an elevated IFN response and a reduced cytokine/NK response To further explore the transcriptomic diversity of malignant cells, we conducted non-negative matrix factorization (NMF) analyses on 10 different tumors separately^[88]29,[89]30. This analytical approach allowed us to decipher intrinsic biological programs and mitigate batch effects among malignant cells. Previous studies have indicated that intratumor programs could be shared among different tumors^[90]29. Consequently, we calculated the Jaccard index for each pair of intratumor programs, which led us to identify eight prominent meta-programs that were consistently present across distinct tumors (Fig. [91]2B). These meta-programs were characterized by the functional enrichment results of selected commonly shared genes and were labeled based on associated functions, including cell cycle, stress response, cytokine/NK response, hypoxia, estrogen response, IFN response, and EMT (Fig. [92]2C). Encouragingly, we observed that malignant cells derived from recurrent tumors exhibited an up-regulation of hypoxia-related genes, aligning with the findings from the differential expression analysis (Fig. [93]2D). Recurrent tumors also displayed higher levels of EMT, stress response, and IFN response programs, whereas primary tumors showed significantly higher cytokine/NK response and estrogen response programs (Fig. [94]2D, Supplementary Fig. [95]2E, F). Interestingly, the increased IFN response and decreased cytokine/NK response programs were observed exclusively in malignant cells from recurrent tumors, not in normal epithelial cells (Supplementary Fig. [96]2F). Given hamper T and NK cell function in recurrent tumors, we further examined whether the recurrent tumors-derived malignant cells could modulate cytokine/NK response by upregulated immune inhibitors, such as the transcription factors and cytokines. As expected, the transcription factor GATA3, inducing T cell dysfunction independent of activation, and immunosuppressive cytokine TGFB1 significantly up-regulated in malignant cells from recurrent tumors compared to primary tumors^[97]30–[98]32 (Fig. [99]2D). Altogether, these results collectively suggest that recurrent tumors may evade immune surveillance by elevating the levels of immune inhibitors, thereby reducing cytokine/NK response. It has been reported that excessive and prolonged levels of IFN could orchestrate tumor progression and immune escape in solid tumors^[100]32. Therefore, we sought to investigate the impact of elevated IFN response on recurrent tumor-derived malignant cells. First, we performed cross-referencing and observed that malignant cells from recurrent tumors exhibited intrinsic upregulation of IFN-stimulated genes, including CD74, IFI27, IFITM3, HLA-DRA, and HLA-DRB1 (Fig. [101]2F). For the top differentially expressed gene CD74, we further confirmed its co-expression with malignant cell markers CK-P and CD74 in individual malignant cells using multiplex immunohistochemistry (mIHC) (Fig. [102]2G, H). This finding was consistent with previous studies indicating that IFNG enhances CD74 expression in melanoma. Moreover, we identified a positive correlation between IFNG and CD74 expression both in single-cell and bulk RNA-seq datasets (Supplementary Fig. [103]2H)^[104]33. CD74 plays a crucial role not only in antigen presentation but also as a receptor for macrophage migration inhibitory factor (MIF), which activates the PI3K/AKT pathway. PI3K/AKT further degrades p53 and activates mTOR to promote tumorigenesis and metastasis^[105]34,[106]35. Given the oncogenic significance of CD74, we confirmed elevated levels of PI3K/AKT in malignant cells derived from recurrent tumors compared to primary tumors (Supplementary Fig. [107]2I). To further investigate the transcription factors responsible for mediating CD74 expression in the presence of excessive IFN levels, we performed SCENIC analyses comparing CD74^+ and CD74^− BLCA malignant cells. Our analysis revealed that STAT1 had the highest regulatory importance, consistent with previous studies demonstrating the involvement of IFNG/STAT1 signaling in CD74 regulation (Fig. [108]2I). Additionally, we observed significantly higher expression of STAT1 in recurrent tumors (Fig. [109]2J). Interestingly, CD74^+ malignant cells displayed enhanced energy production-associated pathways, such as pyruvate metabolism, fatty acid degradation, and glycolysis. This suggests that these tumor cells are more metabolically active and may exert a stronger influence on tumor growth (Supplementary Fig. [110]2J). In conclusion, our findings suggest that malignant cells derived from recurrent tumors could reduce cytokine/NK response and activate the PI3K/AKT pathway through CD74-MIF interaction. This mechanism may contribute to tumor escape and progression. Enhanced immune suppression and cellular plasticity of the recurrent TME To gain a deeper understanding of the diverse infiltrating immune cells, we conducted a comprehensive analysis of two major cell lineages: lymphocytes and myeloid cells. In the case of lymphocytes, we meticulously addressed any batch effects, followed by re-clustering to identify and annotate 13 distinct lymphocyte subtypes using canonical marker genes (Fig. [111]3A and Supplementary Fig. [112]3A, B). To characterize the functional state of each cell type, we performed enrichment analyses using cytotoxic and regulatory signatures^[113]36. Lymphocytes expressing high levels of GZMB and GZMA, along with increased cytotoxic scores, were annotated as CD8Teff_GAMB and CD8Teff_GZMA, respectively. Cell types exhibiting elevated regulatory scores were identified as Treg (Fig. [114]3B). Naïve T-cells, such as Tnaive_TCF7, were defined by their lowest cytotoxic and regulatory scores. Interestingly, we observed a decrease in cytotoxic scores in T-cells derived from recurrent tumors, accompanied by an increase in regulatory and exhausted scores, indicating a more pronounced immunosuppressive state of T-cells within the recurrent tumor microenvironment (TME) (Fig. [115]3B). Proportional analyses of different cell types revealed a consistent cellular composition within patients across various tumor regions (Fig. [116]3C and Supplementary Fig. [117]3C). Notably, recurrent tumor samples exhibited significant enrichment of Treg_TNFRSF4 and T_CXCL13, both demonstrating high exhausted scores (Supplementary Fig. [118]3D). Taken together, these findings strongly suggest that recurrent tumors exhibit a more immunosuppressive TME, characterized by reduced cytotoxic T-cell activity and an increased presence of regulatory T-cells. Fig. 3. Immune cell heterogeneity of primary and recurrent bladder tumors. [119]Fig. 3 [120]Open in a new tab A UMAP visualization of the distribution of lymphocytes, colored by cell types. B Scatter plot showing the cytotoxic, exhausted, and regulatory signature score for each lymphocyte. C Histogram indicating the proportion of each cell type in each sample. D UMAP visualization of the distribution of myeloid cells, colored by cell types. E Boxplot showing the infiltration levels of each cell type in primary and recurrent tumors. F Scatter plot showing the M1/2 polarization of all macrophages and monocytes. G Bubble heatmap showing the interaction strength of gene pairs between macrophage and other cell types. These scores are normalized expression levels, and the sizes of the bubbles indicate the significance of the interactions, calculated by CellPhoneDB. H Kaplan–Meier curves showing the clinical effect of SPP1^+ macrophages in public NMIBC dataset. For myeloid cells, we identified 14 subtypes based on Louvain clustering and expressed markers, encompassing four clusters for macrophages, two for monocytes, four for dendritic cells (DCs), two for mast cells, one for neutrophils, and one for cycling cells (Fig. [121]3D and Supplementary Fig. [122]3E). Distribution analysis revealed that cDC_LAMP3 and cDC1_CLEC9A were highly enriched in recurrent tumor samples, while CAP3^+ and CTSG^+ mast cells were more prevalent in primary tumor samples (Fig. [123]3E and Supplementary Fig. [124]3F). To describe the polarization of macrophages and monocytes, we utilized M1 and M2 signature scores^[125]37. Intriguingly, our observations indicated that macrophages and monocytes derived from recurrent tumors exhibited simultaneous M1 and M2 polarization, exemplified by the presence of Macro_SPP1 and Mono_FCN1 (Fig. [126]3F). Furthermore, we employed phagocytosis, angiogenesis, and myeloid-derived suppressor cell (MDSC) signature scores to characterize the heterogeneity of monocytes and macrophages^[127]38. Notably, Mono_FCN1 displayed an angiogenesis and MDSC phenotype, while Macro_SPP1 exhibited the highest M2 and angiogenesis signature scores, suggesting their potential roles in regulating immunity as tumor-associated macrophages (Supplementary Fig. [128]3F, TAMs). Interestingly, we observed that SPP1^+ macrophages showed a higher propensity for interactions with lymphocytes through SPP1, specifically involving SPP1-CD44 and SPP1-PTGER4. Previous studies have demonstrated that the interaction of SPP1-CD44 and SPP1-PTGER4 has the potential to inhibit T-cell proliferation^[129]7,[130]39 (Fig. [131]3G and Supplementary Fig. [132]3G). Consistently, the elevated expression of the signature gene of Macro_SPP1 was associated with poorer clinical outcomes in public NMIBC cohort (Fig. [133]3H). Overall, our analyses suggest the presence of immunosuppressive tumor-associated macrophages (TAMs) in both primary and recurrent tumors, with TAMs associated with recurrence showing increased cellular plasticity. IL4I1^+SPP1^+ macrophage promotes AHR-driven cancer cell motility and suppresses adaptive immunity through recruiting Tregs Metabolic reprogramming plays a crucial role in cell adaptation and influences the cellular state within the TME. Given the remarkable cellular plasticity of myeloid cells, we investigated the metabolic activity of different myeloid subtypes using GSVA. Our findings revealed that Macro_SPP1 displayed the highest activity in tryptophan (Trp) metabolism compared to other myeloid cells (Fig. [134]4A, Supplementary Fig. [135]4A). Previous studies have implicated that Trp metabolism can induce PD-L1 expression, facilitating immune evasion in human pancreatic beta cells, indicating a potential connection between upregulated Trp metabolism and tumor occurrence^[136]40,[137]41. Further gene set analyses indicated a significant upregulation of the Trp metabolic-associated enzyme IL4I1 in Macro_SPP1 and a subset of mature cDC_LAMP3 (Fig. [138]4B, C and Supplementary Fig. [139]4B, C)^[140]42. To validate these findings, we performed mIHC staining on BLCA samples and confirmed the colocalization of IL4I1 and SPP1 in CD14^+ macrophages (Fig. [141]4D and Supplementary Fig. [142]4D). Interestingly, the expression of IL41 is enhanced in the Macro_SPP1 cells derived from recurrent tumors (Fig. [143]4E). Moreover, the co-expression of IL4I1 and SPP1 in Macro_SPP1 cells was also observed in public kidney renal clear cell carcinoma (KIRC) and liver hepatocellular carcinoma (LIHC) datasets, suggesting a potential pan-cancer regulatory mechanism involving IL4I1 and SPP1 (Supplementary Fig. [144]4E). While a subset of cDC_LAMP3 also expressed IL4I1, the Trp metabolic activity in these cells was mainly mediated by IDO1 and TDO2^[145]43 (Supplementary Fig. [146]4B, C). Notably, IL4I1 has been reported to have a stronger association with Aryl hydrocarbon receptor (AHR) activity compared to IDO1 and TDO2^[147]43,[148]44. Therefore, IL4I1-mediated Trp catabolism could drive AHR activation, promoting cancer cell motility, impairing T-cell proliferation, and recruiting regulatory T-cells^[149]44. Ultimately, this contributes to enhanced tumor malignancy and suppression of anti-tumor immunity. Fig. 4. IL4I1 activates the aryl hydrocarbon receptor (AHR) and suppresses adaptive immunity. [150]Fig. 4 [151]Open in a new tab A Boxplot showing the activity of tryptophan metabolism in SPP1^+ macrophages and other macrophages. B Venn diagram illustrating the shared genes between Tryptophan metabolism-associated genes and the SPP1^+ macrophage upregulated genes. C Contour line delineating expression of IL41 on myeloid cells. D Immunofluorescence images showing the colocalization of CD14, SPP1, and IL4I1. Scale bar, 0.05 mm. E Dot plot depicting the expression of tryptophan metabolism-associated enzymes across all cell types. F Density plot illustrating the AHR signature score between malignant cells from primary and recurrent tumors. G Dot plot showing the correlation between the proportion of SPP1^+ macrophage and TNFRSF4^+ regulatory T cells. H Dot plot showing the correlation between the signature score of SPP1^+ macrophages and infiltrated regulatory T cells in the TCGA BLCA cohort. I Top ligands inferred to regulate the SPP1^+ macrophage according to NicheNet (middle). Heatmap showing the expression of ligands mentioned in the middle panel across cell types (left). Ligands ranked by the Pearson correlation (right). J Boxplot showing the CytoSig predicted TNFA and IFNG activities between primary and recurrent tumors (ns, not significant; *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001). To gain further insights into the role of IL4I1, we calculated the AHR signature score for each malignant cell. The results revealed that malignant cells from recurrent samples had significantly higher AHR scores compared to their primary counterparts (Fig. [152]4F). This observation aligns with the elevated expression of IL4I1 in Macro_SPP1 cells from recurrent tumors, suggesting that IL4I1 derived from TAMs is associated with AHR-driven tumor progression, especially in recurrent cases. Besides malignant cells, we found a significant positive correlation between the proportion of Macro_SPP1 and activated Treg_TNFRSF4 (Fig. [153]4G). To validate this, we analyzed a large TCGA BLCA dataset and confirmed the association between Macro_SPP1 and regulatory T-cell infiltration, indicating the role of IL4I1^+SPP1^+ macrophages in modulating tumor immunity by recruiting regulatory T-cells (Fig. [154]4H). Given the known immunosuppressive properties of Macro_SPP1, we investigated whether the enrichment of regulatory T-cells in recurrent tumor samples was primarily influenced by IL4I1 or SPP1. Our analysis revealed that IL4I1 expression exhibited a significantly stronger positive correlation with the level of regulatory T-cell infiltration across different cancer types when compared to SPP1 (Supplementary Fig. [155]4F). Furthermore, IL4I1 overexpression was frequently observed across various cancer types, and heightened levels of IL4I1 were associated with worse overall survival in public NMIBC dataset (Supplementary Fig. [156]4G, H). These findings collectively suggest that IL4I1 may promote tumor progression at a pan-cancer level. To identify potential upstream regulators of IL4I1 expression in Macro_SPP1, we conducted NicheNet analysis and found that IFNG and TNF, predominantly released by GZMA^+ T-cells, MS4A1^+ B-cells, and CLEC9A^+ dendritic cells, were highly expressed in recurrent BLCA tumors and held the highest potential to influence the expression of IL4I1 and SPP1 in Macro_SPP1^[157]45 (Fig. [158]4I, J). This indicates that cellular cross-talk plays a role in shaping the function of TAMs within the tumor microenvironment. In summary, our analyses unveiled that IL4I1 derived from SPP1^+ macrophages drives AHR-mediated cancer cell motility while suppressing adaptive immunity, particularly in recurrent tumors. ECM-associated RGS5^+ myofibroblasts are enriched in recurrent tumors Stromal cells, which play a central role in the TME, exhibit high heterogeneity and are involved in various functions, including extracellular matrix (ECM) remodeling, angiogenesis, immune regulation, and antigen presentation^[159]46–[160]48. In our study, we focused on characterizing stromal cells and identified 11 distinct subtypes based on clustering and marker gene expression analysis (Fig. [161]5A). Notably, we observed that almost all marker genes for both pan-cancer and MIBC-specific stromal cells were expressed in NMIBC-derived stromal cells^[162]16,[163]47, indicating that while stromal cells in NMIBC may not be in direct contact with malignant cells, they are still influenced by the products of the tumor microenvironment (Supplementary Fig. [164]5A, B). To explore the immune regulation of stromal cells, we employed CellphoneDB to infer cell communication and TIMER to estimate the regulation of stromal cells on immune infiltration^[165]49,[166]50. Our findings revealed that Fibro_CTHRC1, Fibro_CXCL1, and Fibro_POSTN were more likely to interact with SPP1^+ macrophages through MDK-LRP1 and WNT5A-ANTXR1. These interactions are known to promote immunosuppressive macrophage differentiation and M2 polarization, indicating that stromal cells can modulate the functional states of macrophages^[167]51 (Fig. [168]5B). Furthermore, we observed a significant positive correlation between the signature score of Fibro_CCL5 and the estimated proportion of CD8^+ T cells, CD4^+ T cells, and neutrophils in TCGA-BLCA datasets (Fig. [169]5C). Pathway enrichment analysis and functional signature scores revealed that CCL5^+ fibroblasts were not only involved in antigen processing and presentation in cells from recurrent tumors but also contributed to leukocyte transendothelial migration. These findings suggested their potential role in recruiting immune cells and consistent with previous studies reported that antigen-presenting associated fibroblast can facilitate immune cell infiltration (Fig. [170]5D, E)^[171]52. Importantly, several stromal cell types, including Fibro_CTHRC1, Fibro_POSTN, Myfibro_RGS5, and Myfibro_ACTA2 were found to be highly enriched in functions related to EMT and ECM remodeling. This indicates that besides their immune regulation roles, these stromal cells actively participate in remodeling the ECM surrounding tumors. Next, we examined the distribution of stromal cell subtypes and observed the proportions of these cell types varied across different samples, but showed more consistency within regions of the same patient (Fig. [172]5F). Specifically, ECM-associated RGS5^+ myofibroblasts were found to be significantly enriched in recurrent tumors, while ECM-associated fibroblasts like Fibro_POSTN and Fibro_CTHRC1 tended to be slightly more enriched in primary tumor samples (Supplementary Fig. [173]5C). In conclusion, our study highlights the presence of diverse ECM-associated fibroblast subtypes in recurrent tumors and primary tumors and underscores the critical role of stromal cells in modulating tumor immunity by influencing immune cell infiltration and regulating macrophage polarization. Fig. 5. Stromal cell heterogeneity of primary and recurrent bladder tumors. [174]Fig. 5 [175]Open in a new tab A UMAP visualization of the distribution of stromal cells, colored by cell types. B Bubble heatmap showing the interaction strength of gene pairs between myeloid cells and fibroblast cells. C Dot plot showing the correlation between the estimated immune cell infiltration and the signature score of CCL5^+ fibroblast cells. D Heatmap showing the enriched pathways for each stromal cell type. E Heatmap showing different expression patterns of function-associated signature genes among all stromal cell types. F Histogram indicating the proportion of each cell type in each sample. Higher consistency of interregional communication network in recurrent tumors In light of the marked transcriptomic heterogeneity between primary and recurrent tumor lesions described earlier, we put forth the hypothesis that each tumor lesion might harbor distinct molecular communication networks between malignant and non-malignant cells. To comprehensively investigate these intratumoral communication networks, we inferred ligand-receptor interactions among all cell types within each sample. Subsequently, we integrated the interaction results of all samples to compare the differences between tumor lesions in primary and recurrent tumors. Notably, our analysis unveiled that recurrent tumors exhibited fewer but more consistent intercellular interactions compared to primary tumors. This implies that the tumor lesion-specific communication networks are more conserved in recurrent tumor patients (Fig. [176]6A, B). These findings align with the higher intra-regional similarity observed in recurrent tumors, providing support for the notion that malignant cells play a pivotal role in shaping the unique communication network^[177]7. Furthermore, we endeavored to identify ligand-receptor pairs specific to their cellular sources and observed enrichment in the EGF, CSPG4, CD40, and TENASCIN signaling pathways in recurrent tumor samples (Fig. [178]6C, D, Supplementary Fig. [179]6A). Notably, we discovered that the interaction between HBEGF and EGFR from the EGF pathway, which has been reported to drive carcinogenesis and tumor growth in lung cancer and gliomas, was exclusively enriched in malignant cells derived from recurrent patients^[180]53–[181]55. Hence, it suggests the involvement of various immune and stromal cell types in stimulating growth factors within recurrent tumors to support the growth of malignant cells (Fig. [182]6E and Supplementary Fig. [183]6B). Additionally, myofibroblast cells, particularly those derived from recurrent tumors, released CSPG4, which could promote malignant cell growth and stimulate endothelial cell migration by modulating integrin function in the TME (Fig. [184]6E and Supplementary Fig. [185]6C)^[186]54. Taken together, our results indicate that the recurrent tumor may possess relatively stable communication networks among distinct regions compared to primary tumors. Moreover, the TME in recurrent tumors appears to rely more on growth factors to drive tumor progression. Fig. 6. Communication of malignant cells and non-malignant cells. [187]Fig. 6 [188]Open in a new tab A Bar plot showing the number of inferred interactions in primary and recurrent tumors. B Similarity of ligand-receptor interactions among tumor lesions of different patients. Zero indicates no overlap of ligand-receptor interactions while 1 means a full overlap of ligand–receptor interactions between samples. C Stacked histograms of the significant pathways were ranked based on differences in the overall information flow within the inferred networks between primary and recurrent tumor samples. D Heatmap showing the outgoing signaling patterns of source-specific interaction pathways among all cell types. E Circos plots showing signaling interaction between malignant and non-malignant cells in recurrent tumors. Discussion The presence of molecular changes in multi-regional tumors, both within individual patients and across different patients, contributes to a considerable level of heterogeneity compared to single-regional tumors. Furthermore, although advances in surveillance and treatment strategies have improved overall survival rates, the clinical outcome for multifocal BLCA patients remains low, particularly due to its early recurrence. Recent studies have focused on analyzing tumor heterogeneity, the interaction between malignant cells and fibroblasts, and the response to immunotherapy using techniques such as scRNA-seq and spatial transcriptome datasets in primary and recurrent BLCA tumors. However, a comprehensive understanding of the heterogeneity among tumor lesions within a single patient and the differences in the TME between primary and recurrent multifocal tumors is still lacking. In this study, we present a comprehensive single-cell level atlas to characterize the TME in primary and recurrent multifocal BLCA. Our analyses reveal significant differences in malignant cells, immune cells, and the stromal context between recurrent and primary tumors. This thorough understanding of the interregional or intraregional dynamics of tumor heterogeneity not only provides valuable insights but also offers potential novel targeting strategies for the treatment of recurrent tumors. Our analyses revealed that recurrent and primary tumors exhibited a smaller degree of interregional heterogeneity compared to intertumoral heterogeneity. Additionally, malignant cells derived from recurrent tumors showed even higher interregional similarities, suggesting that multifocal tumors in recurrent cases may originate from a common set of “seed” cells. Furthermore, we observed that recurrent tumor-derived cells demonstrated more consistent communication networks within each region. Upon functional analysis of malignant cells from recurrent tumors, we observed a suppressed cytokine and NK cell response, along with an elevated expression of IFN-stimulated genes, including CD74, HLA-DRA, and HLA-DRB1, which are involved in antigen presentation. Of particular interest, the expression of MHC II-associated genes, typically associated with canonical antigen-presenting cells, was detected in malignant cells, indicating their dual epithelial-immune nature. Previous studies have reported the expression of MHC II-associated genes in epithelial cells of the respiratory tract or skin, as well as in malignant cells from nasopharyngeal carcinoma, melanoma, and colon cancer. However, the impact of MHC II expression on patient survival and response to immunotherapy remains controversial across different cancer type^[189]23. For instance, tumor-specific MHC II expression has been linked to favorable outcomes in melanoma, colon cancer, and breast cancer, but correlated with shorter overall survival in nasopharyngeal cancer patients. Our findings suggest that the increased expression of CD74 in malignant cells may activate the oncogenic PI3K/AKT pathways, contributing to tumor progression in recurrent bladder tumors. Furthermore, our investigation into the regulatory mechanisms of CD74 revealed its upstream regulator STAT1, a key mediator of IFN signaling, in malignant cells from recurrent tumors. Considering the well-established role of the IFNG signaling pathway as a critical immune signaling pathway in bladder cancer, mediating response to immunotherapy such as Bacillus Calmette Guérin (BCG) and immune checkpoint inhibitors (ICI)^[190]4,[191]5, we hypothesize that CD74 could serve as a promising therapeutic target for combined therapy in the treatment of recurrent bladder tumors. Macrophages can adopt distinct polarization states dictated by the TME, either promoting or suppressing tumor growth. In our study, we observed that SPP1^+ macrophages derived from recurrent tumors exhibited increased activity in tryptophan metabolism, which was attributed to the upregulation of IL4I1 compared to their counterparts in primary tumors. IL4I1, a secreted enzyme, plays a role in catalyzing the conversion of tryptophan^[192]42. This metabolic pathway activation can activate the AHR, leading to tumor progression and the establishment of immune-suppressive microenvironments. Notably, we also observed a significant elevation in IL4I1 expression in malignant cells derived from recurrent tumors, resulting in a pronounced increase in the AHR signature score. This, in turn, promotes tumor progression by enhancing cancer cell motility. Furthermore, our findings revealed a positive correlation between SPP1^+ macrophage proportion and regulatory T-cells, indicating their potential to stimulate angiogenesis, foster tumor progression, and suppress adaptive immunity through their associated metabolic activities. Considering the critical role of tryptophan metabolism in cancer therapy, our discovery provides valuable insights into the potential development of novel pharmacological targets for the treatment of BLCA. Different from analyses of other cancer types, NMIBC includes stromal cells that may not directly connect with adjacent tumor cells or other TME cells, such as those in the Ta stage. According to the American Joint Committee on Cancer (AJCC) staging system, NMIBC is classified into Ta and T1 stages. In T1 stage bladder cancer, the tumor invades the lamina propria, a region primarily consisting of loose connective tissue with numerous elastic fibers, stromal cells, and delicate smooth muscle bundles. While the majority of our collected samples were from patients at the T1 stage, we also included patients at the Ta stage. Importantly, almost all stromal cell marker genes associated with pan-cancer and bladder cancer were highly expressed in our stromal cells at both T1 and Ta stages. These findings suggest that stromal cells, either directly or indirectly, interact with malignant cells and other TME components through cytokine signaling, indicating a complex interplay within the tumor microenvironment. In conclusion, our study presents a valuable contribution by providing comprehensive insights into the heterogeneity and variation of the TME in multi-tumor BLCA lesions from both primary and recurrent tumors. The distinctive features observed in recurrent tumors not only offer unique therapeutic targets but also provide indications of the co-evolution between malignant cells and the TME. Despite the relatively small sample size resulting from sampling limitations, our findings offer an initial glimpse into the remarkable differences between primary and recurrent tumors, with potential implications for studying other cancer types in the future. Methods Sample collection and processing A total of five urothelial cancer patients from Shanghai Tenth People’s Hospital were enrolled in this study. All patients received either CT or MRI scanning preoperatively and were clinically diagnosed as multi-region non-muscle invasive urothelial tract cancer. Among them, three patients were primary and two were recurrent urothelial cancer. Transurethral resection was performed as the golden standard for bladder cancer patients. For one patient (P03) with concurrent bladder and ureter cancer, a radical nephroureterectomy was also performed to remove the ureter cancer. Samples from two distinct parts of the urothelial tract were harvested for each patient, for example, the right and left bladder wall, ureter, and bladder. Each sample was measured about 5 mm diameter in size before single-cell library preparation. All patients were pathologically confirmed as non-muscle-invasive urothelial carcinoma. Sample collection was performed with written informed consent from patients. All ethical regulations relevant to human research participants were followed. Tissue dissociation and preparation The fresh tissues were stored in the sCelLiveTM Tissue Preservation Solution (Singleron Bio Com, Nanjing, China) on ice after the surgery within 30 min. The specimens were washed with Hanks Balanced Salt Solution (HBSS) 3 times and then digested with 2 ml sCelLiveTM Tissue Dissociation Solution (Singleron) by Singleron PythoN™ Automated Tissue Dissociation System (Singleron) at 37 °C for 15 min. Afterward, the GEXSCOPE red blood cell lysis buffer (Singleron, 2 ml) was added, and cells were incubated at 25 °C for another 10 min to remove red blood cells. The solution was then centrifuged at 500 g for 5 min and suspended softly with PBS. Finally, the samples were stained with trypan blue (Sigma, United States) and the cellular viability was evaluated microscopically. Library preparation and scRNA-seq Single-cell suspensions (1 × 10^5 cells/ml) with PBS (HyClone) were loaded into microfluidic devices using the Singleron Matrix Single Cell Processing System (Singleron). Subsequently, the scRNA-seq libraries were constructed according to the protocol of the GEXSCOPE Single Cell RNA Library Kits (Singleron^[193]55). Individual libraries were diluted to 4 nM and pooled for sequencing. At last, pools were sequenced on Illumina HiSeq X with 150 bp paired-end reads. Immunofluorescence staining Protein expression was evaluated by immunofluorescence staining (IF) of bladder cancer tissue harvested from patients. Briefly, 5-μm paraffin-embedded cross-sections of tissues were fixed with 4% paraformaldehyde, permeabilized with 0.1% Triton-100, and blocked with 5% bovine serum albumin (BSA). Subsequently, tissues were incubated with anti-CD74 (1: 500, Proteintech) and anti-pan-keratin (1:200, CST) or anti-SPP1(1:400, Proteintech), anti-IL4I1(1:400, Proteintech) and anti-CD14 (1:1000, Proteintech) overnight at 4 °C. Then, the cells were treated with the corresponding secondary antibody after washing three times. DAPI (Beyotime) staining was used for nuclear localization. Images were captured with a confocal microscope (Leica Microsystems, Mannheim, Germany). Single-cell RNA-Seq data processing The CeleScope pipeline from Singleron was employed to process the fastq file. The data process encompassed filtering out low-quality reads, aligning reads to the GRCh38 reference genome, assigning cell barcodes, and extracting the unique molecular identifier (UMI) barcode. After the initial processing, we employed Scrublet to remove cells likely to be doublets^[194]56. Scrublet was applied to each sequencing library, targeting potential doublets with an expected doublet rate of 6%. Then, the output expression matrix of each sample was combined by MAESTRO^[195]57. Next, we imported the merged count matrix into the Seurat toolkit to filter genes detected in fewer than 10 cells, as well as cells with fewer than 500 genes or exceeding 15% mitochondrial gene counts^[196]58. The raw UMI counts were normalized by the total counts per cell (library size) and followed by scaling to a factor of 10^6 and logarithmic transformation. Dimension reduction and cell clustering for scRNA-seq data To perform dimension reduction and unsupervised clustering on scRNA-seq data, we followed the Seurat standard workflow. The top 3000 highly variable genes (HVGs) were selected based on their average expression and dispersion level using the FindVariableGenes function. To reduce noise, we employed the normalized and scaled HVGs’ expression profiles to do principal component analysis (PCA). From a total of 75 components, we retained the top 50 components for further use for non-linear dimension reduction, resulting in the generation of the Uniform Manifold Approximation and Projection (UMAP) plot for visualization. Then, the Seurat function “FindClusters” was used to cluster cells with a resolution parameter set to 1.0. Additionally, we calculated differentially expressed genes for each cluster using the FindAllMarkers function. Next, we annotated each cluster based on the expression of canonical marker genes, such as EPCAM and KRT14 for epithelial cells, PECAM1 and VWF for endothelial cells, and MS4A1 and CD79A for B-cells. To identify malignant cells, we leveraged the CopyKat R package to estimate the degree of clonal large-scale chromosomal copy number variations (CNV) in each cell within each sample^[197]59. According to the CNV clusters result, we divided all epithelial cells into malignant cells and epithelial cells. From the distribution of cell type in the UMAP, we observed the cells were clustering according to the cell types, except for epithelial cells and malignant cells, in which cells were clustered according to the samples instead of cell types (Supplementary Fig. [198]1D). To enhance the visualization of cell type distributions, we partitioned the cells into two groups: one comprising epithelial and malignant cells and the other encompassing all other cell types. Subsequently, we repeated the dimensionality reduction, graph clustering, and UMAP visualization steps separately for each of these two cell groups. Deciphering intratumor expression programs and meta-programs of BLCA To investigate the intratumor heterogeneity and capture the biology variation to the full extent, we employed a non-negative factorization (NMF) algorithm to identify underlying intratumor expression programs^[199]60. For each sample, the normalized expression matrix of malignant cells was subjected to NMF analysis, resulting in the extraction of twelve putative programs. Then, for each program, we ordered the features by weight and selected the top 30 features to represent the programs. To evaluate program similarity, we calculated the average similarity (Jaccard index) between each program pair and leveraged hierarchical clustering based on similarity to identify the recurring programs, which were subsequently termed “meta-programs”. Finally, we named the meta-program based on their functional enrichment result. Batch effect correction At the major cell type levels, we observed cells, except for malignant and epithelial cells, were clustering according to their cell types. For the malignant and epithelial cells, the clustering pattern was influenced by a combination of batch effects and the inherent heterogeneity of malignant cell populations. Consequently, we employed canonical correlation analysis (CCA) to correct batch effects and assess the proportion at the major cell type level. To keep the biology variation while mitigating the batch effect, we used NMF to elucidate the expression program within samples. Next, to investigate the function and state of tumor microenvironment cells, we categorized cells into three major-cell lineages: lymphoid cells (including T-cells, B-cells, and plasma cells), myeloid cells (encompassing monocytes, macrophages, and, mast cells), and stromal cells (comprising endothelial cells, fibroblasts, and myofibroblasts). For each major cell lineage, we repeated the dimensionality reduction, graph clustering, and UMAP visualization. Notably, cells from the same samples and patients tend to cluster together, indicating the clustering results were influenced by the batch effects (Supplementary Fig. [200]2A). To rectify this issue, for each major lineage, we employed CCA from the Seurat package to eliminate the batch effect. Subsequently, all downstream analyses were conducted based on the batch-corrected results, ensuring not confounded by technical noise. Metabolic activity From the Kyoto Encyclopedia of Genes and Genomes (KEGG) database, we collected the metabolic pathways, including 85 pathways and 1566 genes grouped into 11 categories based on KEGG classifications^[201]61. Due to the scRNA-seq data with the high drop-out rate, we first employ the Markov Affinity-based Graph Imputation of Cells (MAGIC) method to imputation expression profiles, aiming to recover the missed expression values for enzymes associated with metabolism^[202]62. To ensure the accurate representation of metabolic enzyme expression, we accounted for interdependencies among enzymes involved in the same or separate reactions. Then, we used the Gene Set Variation Analysis (GSVA) to quantify the metabolic activities of each metabolic pathway. Differential pathway activation within each subtype was identified through Wilcoxon rank sum tests, comparing the specific subtype against other cell types within the same lineage. Function analysis Pathway enrichment analysis: To depict the functional role of gene sets, such as the gene list derived from the NMF programs and the cell type-specific up-regulated genes, we collected the cancer hallmark gene sets and KEGG pathways from the Molecular Signatures Database (MSigDB v6.1)^[203]48. Then, the hypergeometric test from the clusterProfiler package was used to estimate the overrepresented of interested gene sets in known pathways^[204]63. Significantly enriched pathways were defined as those exhibiting a log-fold-change > 0.05 and an adjusted P-value < 0.01 based on the query gene set. Scoring cell types by using function-associated signature genes: To describe the functional diversity of cell types, we collected the signature gene lists from the previously published studies (Table [205]S2). For T and NK cells, we used the AddModuleScore function within Seurat to compute the cytotoxic, exhausted, and regulatory scores for each cell, thus characterizing the functional states across different regions and validating subtype annotations. Additionally, we assessed the functional attributes of monocyte and macrophage cells by evaluating M1/M2 polarization, angiogenesis, phagocytosis, and pro-/anti-inflammatory capabilities. Cell–cell interaction analysis To analyze cell-cell communication, we employed CellPhoneDB to infer interactions based on the expression of known ligand-receptor pairs^[206]50. CellPhoneDB was used to predict the potential interaction strength between cell types based on the expression of ligand-receptor pairs. The significance of ligand-receptor interactions between two cell types was calculated through permutation tests. Ligand-receptor pairs exhibiting statistically significant interactions (P-value < 0.05) were extracted. To delineate interactions specific to primary or recurrent tumors, we followed the standard CellChat workflow with default parameters. Then, we merged the interaction results from primary and recurrent tumors using the mergeCellChat function, enabling us to visualize the number and strength of inferred interactions for each source (Fig. [207]6A). Utilizing the merged Cellchat object, we further examined and visualized shared and source-specific signaling patterns between any two sources using the netAnalysis_signalingChanges_scatter function (Fig. [208]6C). Survival analysis To evaluate the clinical effect of interest cell types or NMF-derived programs, we collected and downloaded the bladder cancer expression and clinical cohort from public NMIBC dataset. For each NMF-derived intra-tumor meta-program, we selected the genes that appear in more than half the number of programs and performed univariable Cox proportional hazards regression to associate the meta-program with overall survival. A hazard ratio higher than one means an increased mortality risk, while a hazard ratio lower than one suggests a decreased mortality risk. For the tumor microenvironment cells, we directly used the top 50 upregulated genes to represent the cell type and calculated the cell type signature score by gene sat variation analysis (GSVA) in NMIBC patients. Patients were then stratified into high and low groups based on the median values of these signature scores. Kaplan–Meier survival curves were plotted to show the differences in survival time, with the statistical significance determined using log-rank p-values reported by Cox regression models implemented in the R package ‘survival’. Gene regulatory analysis for scRNA-seq data To identify potential transcription factors (TFs) inducing the CD74 expression on malignant cells, we employed LISA, which constructs an epigenetic model and utilizes histone mark ChIP-seq data and chromatin accessibility profiles^[209]62. The up-regulated genes of CD74 positive and negative cells were imported into LISA with default parameters to identify the top TF candidates (Fig. [210]2L). NicheNet, a powerful tool for predicting ligands that drive transcriptomic change in target cell types, was used to identify potential ligands responsible for inducing IL4I1 expression in SPP1^+ macrophages. Highly expressed genes in SPP1^+ macrophage with log2FC > 0.25 and adjusted P-value < 0.05 were used as the input, with all macrophage expressed genes serving as background genes. The ligand-receptor interactions were narrowed down by the expression profile of target cells and sender cells. Then, the function predict_ligand_activities from the R package “nichenetr” was used to predict and rank the potential ligands, and the top potential ligands were visualized. Statistics and reproducibility Statistical analysis and visualization were conducted using R. The statistical methods employed for each analysis are explained in the text and figure legends. Statistical significance is denoted by specific symbols: p-value  <  0.1, *p-value  <  0.05, **p-value  <  0.01, and ***p-value  <  0.001. Graphs were created using various R packages, including ggplot2 (version 3.4.3), ggpubr (version 0.6.0), ggsci (version 3.0.0), ggrepel (version 0.9.3), Seurat (version 4.3.0.1), pheatmap (version 1.0.12), VennDiagram (version 1.7.3), nichenetr (version 1.0.0), and CellChat (version 1.1.3). Each boxplot displays the 25% (lower hinge), 50% (median), and 75% (upper hinge) quantiles, with the lower and upper whiskers representing the smallest and largest observations that are greater than or equal to the lower hinge minus 1.5 times the interquartile range (IQR) and less than or equal to the upper hinge plus 1.5 times the IQR, respectively, as defined by the default settings in the geom_boxplot function of ggplot2. Reporting summary Further information on research design is available in the [211]Nature Portfolio Reporting Summary linked to this article. Supplementary information [212]Transparent Peer Review file^ (908.1KB, pdf) [213]Supplementary Information^ (37.8MB, pdf) [214]42003_2024_7343_MOESM3_ESM.pdf^ (29KB, pdf) Description of Additional Supplementary Files [215]Supplementary data 1^ (16.5MB, xlsx) [216]Supplementary data 2^ (14MB, xlsx) [217]Reporting Summary^ (1.9MB, pdf) Acknowledgements