Single-cell atlases link macrophages and CD8+ T-cell subpopulations to disease progression and immunotherapy response in urothelial carcinoma

Rationale: Immune checkpoint inhibitors (ICIs) have revolutionized the management of locally advanced or metastatic urothelial carcinoma. Strikingly, compared to urothelial carcinoma of the bladder (UCB), upper tract urothelial carcinoma (UTUC) has a higher response rate to ICIs. The stratification of patients most likely to benefit from ICI therapy remains a major clinical challenge. Methods: In this study, we performed the first single-cell RNA sequencing (scRNA-seq) study of 13 surgical tissue specimens from 12 patients with UTUC. The key results were validated by the analysis of two independent cohorts with bulk RNA-seq data for UCB (n = 404) and UTUC (n = 158) and one cohort of patients with metastatic urothelial carcinoma (mUC) who were treated with atezolizumab (n = 348). Results: Using scRNA-seq, we observed a higher proportion of tumor-infiltrating immune cells in locally advanced UTUC. Similar prognostically relevant intrinsic basal and luminal-like epithelial subtypes were found in both UTUC and UCB, although UTUC is predominantly of the luminal subtype. We also discovered that immunosuppressive macrophages and exhausted T-cell subpopulations were enriched in the basal subtype and showed enhanced interactions. Furthermore, we developed a gene expression signature (Macro-C3 score) capturing the immunosuppressive macrophages that better predicts outcomes than the currently established subtypes. We also developed a computational method to model immune evasion, and the Macro-C3 score predicted therapeutic response of mUC treated with first-line anti-PD-L1 inhibitors in patients with lower basal scores. Conclusions: Overall, the distinct microenvironment and Macro-C3 score provide an explanation for ICI efficacy in urothelial carcinoma and reveal new candidate regulators of immune evasion, suggesting potential therapeutic targets for improving antitumor immunity in the basal subtype.


Introduction
Urothelial carcinoma (UC) is the fourth most common tumor type, and upper tract urothelial carcinoma (UTUC) is a distinctly aggressive genitourinary entity of UC. UTUC arises from the urothelial lining of the upper urinary tract, which includes the renal pelvis and ureter [1]. Despite the common histologic appearance of urothelial carcinoma of the bladder (UCB), approximately 60% of UTUCs and 15-20% of UCBs are diagnosed as muscle-invasive (MI), which correlates with a worse prognosis than non-muscle-invasive (NMI) cases with an excellent overall prognosis [2]. Genomic and epigenetic subtyping of UCB and UTUC, including that completed by our group, has revolutionized the current understanding of UC pathogenesis [3][4][5][6][7]. However, even the most recent consensus molecular classification study of UTUC [8] does not provide compelling evidence for molecular pathogenesis and thus lacks the ability to predict the therapeutic response.
Tumor phenotype is represented by cell state measured as RNA expression. Therefore, cancer subtyping based on gene expression profiles can provide valuable insights into cancer biology and important clues for treatment planning. Notably, a previous study identified two intrinsic subtypes of muscle-invasive UCB, termed "luminal" and "basallike", which have the characteristics of different stages of urothelial differentiation. Moreover, a 47-gene predictor, bladder cancer analysis of subtypes by gene expression (BASE47), was developed to subtype UCB into prognostic subtypes [9]. However, it remains unclear whether these bulk RNA-seq data mask heterogeneity at the single-cell level derived from the presence of both malignant and nonmalignant cells. Moreover, it remains unclear whether similar prognostically relevant intrinsic subtypes exist in UTUC, and thus, our understanding of their relationship to therapeutic response is limited.
Cell state is an integration of inputs from cell-intrinsic (e.g., mutational background, epigenetic state) and cell-extrinsic (e.g., cell-to-cell interactions, tumor microenvironment) sources. Dynamic interactions between the tumor and tumor microenvironment (TME) are essential for cancer cell heterogeneity and therapeutic response [10]. In particular, distinct populations of T cells influence tumorigenesis by performing pro-and anti-tumorigenic functions. T cell exhaustion is the main reason why they fail to eliminate tumor cells and various forms of cancer immunotherapy are now effective in re-establishing and promoting T cell anti-tumor immunity [11,12]. By suppressing T cell function, other immune cell types maintain immune equilibrium, such as M2 polarized macrophages and tolerogenic dendritic cells [13]. Tumor-associated macrophages (TAMs) account for a significant fraction in human solid malignancies and influence tumor progression by providing nutritional support to malignant cells, now constituting promising targets for novel anticancer agents [14]. Advancements in single-cell RNA sequencing (scRNA-seq) have led to the ability to simultaneously evaluate the cell states of both malignant and nonmalignant cells in TME and thus provide unique opportunities to explore how tumor cell subtypes are shaped by cell-intrinsic and cell-extrinsic factors. In this context, we hypothesized that comprehensive profiling of malignant and nonmalignant cells in UTUC at the single-cell level would enable the deconvolution of current molecular subtypes into their constituent parts, which would consequently enable the development of more effective prognostic and predictive tools.
In this study, we performed single-cell RNA sequencing of 13 surgical tissue specimens obtained from 12 patients with UTUC. We observed that MI UTUC had a higher abundance of tumor-infiltrating immune cells than NMI UTUC. We also found that similar prognostically relevant intrinsic basal and luminal-like subtypes exist in UTUC. Notably, we identified distinct subtype-specific immunological phenotypes in both UTUC and UCB and determined their regulation by the tumor microenvironment and their response to immune checkpoint therapy.

Single-cell transcriptomic landscape in UTUC
To characterize intratumoral heterogeneity, we collected 13 surgical tissue samples, including 12 primary tumors and a paired tumor thrombus, from 12 UTUC patients for scRNA-seq (Supplementary Table S1). We also included two cohorts with bulk RNA-seq for UCB and UTUC, and one cohort of atezolizumab-treated patients with metastatic urothelial carcinoma (IMvigor210) [15,16] for further validation ( Figure 1A). Quality control of the scRNA-seq data resulted in the inclusion of 67,392 cells (Supplementary Figure S1A), which were then integrated into an unbatched dataset - (Figure 1B Figure S1E). Furthermore, we examined a larger cohort of UTUC patients [8] with available bulk transcriptome data and observed significantly increased proportions of immune cell subsets in the MI tumors compared to the NMI tumors ( Figure 1G & Supplementary Figure S1F). These findings suggest that crosstalk between tumor cells and the immune environment may play a role in shaping muscle invasiveness.

Intrinsic subtypes of MI UTUC reflect the luminal and basal subtypes of MI UCB
To better understand the cellular programs active in cancer cells that may function together with immune cells, we next sought to identify the intrinsic developmental subtypes of the malignant cells. We first confirmed 31,152 malignant cells based on large-scale chromosomal copy number variations (CNVs) in each cell (Supplementary Figure S2A-C) inferred by the inferCNV tool [17]. Considering the rarity of UTUC, we used the BASE47 panel, a 47-gene predictor defining malignant cells of MI UCB into prognostic subtypes [9], and identified two distinct expression profiles, termed luminal and basal, in the MI UTUC samples (Figure 2A & Supplementary  Table S2). However, with the limited number of cases, all three NMI UTUCs were classified as the luminal subtype, which is associated with a better clinical outcome (Supplementary Figure S2D). Notably, among the 47 markers, the serine peptidase inhibitor Kazal type I (SPINK1) showed enhanced expression in luminal UTUC but not in basal tumors (Supplementary Figure S2E). Immunohistochemistry showed higher expression of SPINK1, a potential luminal subtype-specific biomarker for UTUC ( Figure  2B). Using the 47 genes, we defined two computational scores, the basal score and the luminal score, to quantify the two subtypes. In the Japanese UTUC cohort, the basal score was significantly associated with decreased disease-specific survival (DSS), whereas the luminal score was significantly associated with increased DSS in the univariate analysis ( Figure 2C). Furthermore, the basal score and luminal score can generate "poor" or "good" prognosis calls for MI UCB in The Cancer Genome Atlas (TCGA) (Supplementary Figure S2F).
To further understand the gene expression patterns that differentiate the intrinsic subtypes of MI UTUC, we identified 521 differentially expressed genes (DEGs) in the tumor cells of these two subtypes ( Figure 2D, 2E & Supplementary Table S3). Gene Ontology (GO) analysis showed that the upregulated genes in the luminal subtype were enriched for pathways related to the epithelium development. The upregulated genes in the basal subtype were enriched in immune regulation pathways, such as the response to interferon beta and the response to cytokine ( Figure  2F). Interestingly, we observed enrichment of genes involved in regulating myeloid cell differentiation with upregulated expression of the myeloid chemotaxis signature, and immune evasion genes in the basal subtype compared to the luminal subtype ( Figure 2F & 2G). We found that the basal subtype had a stronger immune evasion phenotype than the luminal subtype.

An immunosuppressive macrophage population is enriched in the basal subtype and correlates with poor clinical outcome
We next characterized myeloid cells in UTUC since the tumor cells of the basal and luminal subtypes showed a discrepancy in the regulation of  Table S4). Macrophages were the most dominant myeloid cell type in UTUC and expressed TAM markers, whereas subcluster-specific expression of conventional proinflammatory (M1) and anti-inflammatory (M2) macrophage markers were not observed ( Figure 3C). Interestingly, we identified a subcluster of macrophages (Macro-C3) that was significantly enriched in the basal subtype ( Figure 3D). This subcluster highly expressed CD8 + T-cell recruitment-related genes, such as CXCL10 and CXCL11, as well as immune checkpoint and evasion genes ( Figure 3B & Supplementary Figure S3A). Compared with other macrophages, Macro-C3 showed the upregulation of genes that could positively affect the induction of immune tolerance, indicating a central role in the ecosystem of UTUC ( Figure 3E). This subcluster resembled recently identified immunosuppressive macrophages that consisted predominantly of genes associated with immunosuppression [18].
After identifying a subset of subtype-relevant macrophages, we explored their clinical implications. For this, a Macro-C3 signature was developed and consisted predominantly of genes associated with immunosuppression, such as IDO1 and CD274 ( Figure 3F, Supplementary Table S2 & S4). Consistent with the observations in our scRNA-seq dataset, evaluation of the signature in the Japanese UTUC cohort revealed that the basal subtypes had significantly higher Macro-C3 signature scores than the NMI and luminal subtypes ( Figure 3G). Additionally, higher Macro-C3 signature scores in the Japanese UTUC cohort were significantly correlated with reduced DSS and progression-free survival (PFS) ( Figure 3H). Collectively, these findings suggest a shift toward enhanced immunosuppression of macrophages in basal UTUC.

Expansion of exhausted T cells in the basal subtype
Next, we performed unsupervised clustering of T and NK cells and obtained eight subclusters, including a population of naïve T cells, two clusters of CD4 + T cells, three clusters of CD8 + T cells, and two of NK cells ( Figure Figure  S4B). However, the CD4 + T cells were mainly regulatory (CD4-C1) or exhausted (CD4-C2) ( Figure  4B). To ascertain the developmental trajectory of the T cells, we performed a pseudotime analysis. The CD4 + T-cell branches originated from naïve T cells, which differentiated into regulatory T cells (CD4-C1) and exhausted CD4 + cells (CD4-C2) (Supplementary Figure S4C). The CD8 + T cells, however, showed a different developmental trajectory, starting from the naïve state, passing through the cytotoxic phase (CD8-C2), and ending in exhaustion (CD8-C1) (Supplementary Figure S4C & S4D). Next, we sought to characterize the T lymphocytes in distinct UTUC subtypes. We found that the relative proportion of naïve T cells was significantly lower in the basal subtype than in the other two subtypes (Figure 4D & Supplementary  Figure S4E). Moreover, the exhausted CD8 + T-cell population (CD8-C1) showed significant enrichment in the basal subtype, which indicated increased immune exhaustion and CD8 + T-cell dysfunction in the basal tumor microenvironment ( Figure 4E). A high proportion of exhausted CD8 + T cells was significantly associated with the poor survival of patients with UCB and UTUC according to single-cell deconvolution in the TCGA-BLCA and Japanese UTUC cohorts, respectively ( Figure 4F). Therefore, we concluded that CD8 + T cells in basal UTUC show the highest rate of exhaustion, suggesting weakened immune surveillance.  (E) Relative percentage of CD8-C1 cells (exhausted CD8 + T cells) among the NMI, luminal, and basal subtypes. P10 was excluded due to the presence of <100 T/NK cells. The relative percentage of CD8-C1 cells is defined as the proportion of CD8-C1 cells in the T/NK populations. Error bars represent mean ± standard error of the mean. P values were calculated with a t test. (F) Kaplan-Meier curves of OS and progression free interval (PFI) with the decomposed CD8-C1 (exhausted CD8 + T cells) proportions in the TCGA-BLCA cohort, and DSS in the Japanese UTUC cohort. Samples were categorized as high (red, top 50%) or low (blue, bottom 50%) CD8-C1 abundance. P values were calculated using the log-rank test.

The interaction between myeloid and lymphoid compartments shapes an immunosuppressive microenvironment in basal UTUC
Given the common trends of immunosuppressive phenotypes in the malignant cells, macrophages, and T cells of the basal subtype compared to the luminal subtype, we hypothesized that different cell populations participate in complex crosstalk. Therefore, we firstly investigated possible intercellular interactions influencing TAMs. We found that macrophages had the most frequent communication with epithelial cells over other myeloid cells, especially in the basal subtype ( Figure 5A). Between the two cell types, intense signaling mediated by MIF and CD74 was identified, which has been reported to activate the recruitment of macrophages [20] and promote the secretion of growth factors by TAMs [21] ( Figure 5B). Moreover, basal UTUC exhibited higher CD47-SIRPA and C3-C3AR1 interactions, which are known to inhibit phagocytosis [21] and associated with TAM infiltration [22], respectively ( Figure 5B). The active crosstalk between basal malignant cells and macrophages explained increased immunosuppression in this subtype, yet the process ruling T cell recruitment remains elusive. Therefore, we screened receptor-ligand pairs and found a significant positive correlation between CD8 + T-cell abundance, exhaustion and macrophage-derived CXCL expression ( Figure 5C & Supplementary Figure S5A-C). Moreover, in the basal subtype we identified stronger macrophage-CD8 + T-cell interactions mediated by CXCL and their receptors (Figure 5B), which were essential for the recruitment of T cells into tumors [23]. This association was further validated in the TCGA-BLCA and Japanese UTUC cohorts, where the basal subtype showed significantly enhanced CXCL10 expression and a greater number of tumor-infiltrating CD8 + T cells and macrophages (Supplementary Figure S5D). Mechanistically, multiplex immunofluorescence staining revealed the colocalization of CD68-expressing macrophages (CXCL10 + ) and CD8-expressing CD8 + T cells (CXCR3 + ) in the basal subtype ( Figure 5D) with a closer average distance compared to luminal UTUC (Supplementary Figure  S5E), which further indicated that macrophages may play an important role in the recruitment of CD8 + T cells to UTUC.
Next, we explored the mechanism underlying the exhaustion of recruited CD8 + T cells. In the basal subtype, both tumor cells and macrophages showed enhanced LGALS9-HAVCR2, SPP1-CD44 interactions with CD8 + T cells (Figure 5B), which are known to suppress T-cell functions [24,25]. Furthermore, macrophages of the basal subtype expressed PDCD1LG2, CD274, CD80, and CD86, which target the immune checkpoints CTLA-4 and PDCD1 to inhibit CD8 + T-cell activation ( Figure 5B). Interestingly, the exhausted CD8 + T cell population, CD8-C1, showed intensive signaling in the basal subtype with macrophages ( Figure 5E), particularly Macro-C3 ( Figure 5F). Consistently, the proportions of Macro-C3 and CD8-C1 were positively correlated (Figure 5G), implying a pivotal role of immunosuppressive macrophages in driving the exhaustion of CD8 + T cells. This association was further validated by analyzing the Japanese UTUC cohort, in which exhaustion scores and Macro-C3 signature scores also showed a strong correlation ( Figure 5H). Collectively, this evidence reveals intense interactions between the myeloid and lymphoid compartments in UTUC and suggests that epithelial cells, macrophages, and T lymphocytes cooperate to shape the immunosuppressive microenvironment in the basal subtype with high infiltration of exhausted T lymphocytes ( Figure  5I).

Macro-C3 score predicts the immunotherapy response in mUC
To further evaluate the determinants of response and resistance to immunotherapy, we examined whether subtype-specific immunological phenotypes can predict the treatment response. Pretreatment tumor samples from a large phase II trial (IMvigor210) investigating the clinical activity of PD-L1 blockade with atezolizumab in metastatic urothelial cancer (mUC) were used for this integrated evaluation. We categorized patient tumors into the NMI, luminal, and basal subtypes. In the IMvigor210 dataset, a larger proportion of basal tumors (38%) than luminal (20%) or NMI (11%) tumors exhibited an inflamed phenotype ( Figure 6A). Immunohistochemical staining further confirmed the inflamed phenotype in the basal subtype and the excluded phenotype in the luminal subtype ( Figure 6B & Supplementary Figure S6A). PD-L1 expression on immune cells (more than 5% of immune cells in mUC were detected by the PD-L1 antibody SP142) and the Macro-C3 score were associated with the response in the NMI and luminal subtypes but not in the basal subtype ( Figure 6C). A nonlinear relationship between the log(basal/ luminal) score and fraction of responders to ICI treatment was identified (Supplementary Figure  S6B). Thus, we hypothesized that immunosuppressive interactions between macrophages and T cells may contribute to immune evasion in the basal subtype.
Next, we examined the most significant interactions between macrophages and CD8 + T cells, which were positively and negatively associated with the abundance of exhausted CD8 + T cells. We found multiple known interactions, such as CD274-PDCD1, CD86-CTLA4, and LGALS9-HAVCR2, among the positively associated interactions that promoted T-cell exhaustion [26] (Figure 6D & Supplementary Figure  S6C). Consistently, we identified multiple known interactions, such as IFNG-IFNGR1 and IL-7-IL7R [27], among the negatively associated interactions that reversed T-cell exhaustion (Figure 6D). Moreover, we found that the expression of immunosuppressive Macro-C3 signature genes and T-cell cytotoxicity genes increased with an enhanced basal signature (Figure 6E & Supplementary Figure  S6D). Therefore, we developed a computational method to model the extent of immune evasion using the Macro-C3 score and cytotoxic T cells using the T-cell cytotoxicity score (Figure 6F). We identified an intersection point between the z score of T-cell cytotoxicity and the Macro-C3 score ( Figure 6F); based on the log(basal/luminal) value of this intersection, we divided the IMvigor210 cohort into two groups. Below this threshold, the Macro-C3 score and the fraction of responders to PD-L1 blockade treatment were significantly correlated; thus, we denoted the Macro-C3 signature as the ICI response score (Figure 6G). Above this threshold, the Macro-C3 score represented immune evasion capacity, and thus, we assigned the minimal Macro-C3 score of the cohort as their ICI response score. To evaluate the predictive performance regarding the ICI response, we compared the ICI response score with the Macro-C3 score and a published signature for predicting the clinical response to immune checkpoint blockade, termed the tumor immune dysfunction and exclusion (TIDE) score [28]. The receiver operating characteristic (ROC) curves showed that the ICI response score had a higher AUC (0.664) than the TIDE score (0.552) ( Figure 6H). Overall, the Macro-C3 score predicted the ICI response, especially in the NMI and luminal subtypes, whereas the Macro-C3 score predicted ICI resistance in the basal subtype, likely based on the enhanced T-cell dysfunction induced by immunosuppressive macrophages.

Discussion
Herein, to the best of our knowledge, we have performed the first comprehensive profiling of UTUC at the single-cell level, which allowed us to simultaneously elucidate the cell state of malignant and nonmalignant cells and to derive more therapeutically relevant molecular signatures with higher resolution. Although UTUC is more frequently diagnosed in the locally advanced stage than UCB, we surprisingly found that the intrinsic subtype of tumor cells and subtype-specific immune contexture were similar in UTUC and UCB. Thus, it is reasonable to hypothesize that there are similar epithelial and/or nonepithelial subpopulations with the ability to predict the therapeutic response in both UTUC and UCB. Consistently, in a single-arm phase II study of atezolizumab in mUC (IMvigor210) [15], we identified a similar prediction signature for anti-PD-L1 therapy in the NMI, basal, and luminal subtypes.
It has been shown that among UTUC and UCB patients who receive atezolizumab (anti-PD-L1 therapy), the response rate is higher in UTUC than in UCB (39% vs. 17%) [29,30]. However, they found no significant differences in baseline covariates, including anatomic sites of metastases, tumor mutation load, T-effector gene expression, TCGA subtype, and baseline tumor burden, between the UTUC and UCB [29]. The higher response rate in UTUC may be at least partially explained by its predominant luminal subtype with a lower inflamed immunosuppressive contexture than the basal subtype [31] (Figure 2G, 5E & Supplementary Figure  S7). Additionally, anatomic differences, but not intrinsic subtype or immune contexture, may account for much of this disparity in muscle-invasive UTUC and UCB at diagnosis because the thinner smooth muscle covering of the upper tract may allow for more rapid progression to non-organ-confined disease.
Although the application of bulk gene expression signatures has shown promise in identifying subgroups of patients responsive to immunotherapy [16], yet they provided a limited mechanistic understanding of the cell types responsible for mediating clinical benefit. In contrast, at the single-cell resolution, we found progressive T-cell dysfunction and/or exhaustion, likely through enhanced interactions between immunosuppressive macrophages and T-cell subpopulations. We identified a macrophage population, Macro-C3, that strongly expressed several immunosuppressive effectors, showed the strongest putative interaction with exhausted CD8 + T cells and defined poor patient prognosis in UTUC.
More importantly, we identified diverse cell-cell interactions between T cells and macrophages that promoted or reversed T-cell exhaustion ( Figure 6D). Thus, strategies to inhibit TAM functions, especially those of immunosuppressive Macro-C3, in combination with immunotherapy may have great potential for the treatment of UC patients. Indeed, both preclinical and clinical strategies targeting the tumor-promoting functions of TAMs in cancer are being developed [32]. It has been reported that M2-like macrophages impede ICIs that were based on the PD-1/PD-L1 axis [33]. Consistently, PD-L1 was highly expressed in Macro-C3 cells. Meanwhile, other genes, such as CXCL9, CXCL11, and IDO1, derived from the Macro-C3 signature may also serve as therapeutic targets to ameliorate T-cell exhaustion. Then, we developed a computational method based on the Macro-C3 score to model the two primary mechanisms of tumor immune evasion and T-cell cytotoxicity. Of note, compared to the widely used ICI response markers, PD-L1 levels and the TIDE signature, the Macro-C3 score achieved better performance in the mUC cohort. One explanation for the better performance of the Macro-C3 score is that this score utilized cancer type-specific T-cell dysfunction regulators. A central limitation of this study was the relatively small cohort of patients employed to generate these signatures, and larger cohorts with scRNA-seq data of UC patients would be required. Although the prognostic and predictive values of the identified signatures were validated in multiple external cohorts, further studies to deconvolute other less-characterized immune cell types are necessary for identifying novel therapeutic targets that address immune dysfunction in UC.

Conclusions
Taken together, our results show the common underlying transcriptional subtypes and immune microenvironments in UTUC and UCB, identifying potential opportunities for common management strategies, especially in terms of ICI for UC. Overall, our findings contribute to the understanding of the pathophysiology of UC and provide novel prognostic assessment strategies and individualized treatment recommendations for different UC subtypes.

Patients and samples
Tissue samples were obtained from the UTUC patients who received ureterectomy or radical nephroureterectomy at the Peking University First Hospital. Twelve primary UTUC tumor tissues (three low-stage and nine high-stage), along with one adjacent tumor thrombus, were included in this cohort (Supplementary Table S1). This study was approved by the Ethics Committee of Peking University First Hospital (Grant No. 2018[186]).

Single-cell suspension preparation and droplet-based single-cell RNA sequencing
Tissue samples were processed immediately after the surgical resection. Single-cell suspensions were obtained by mechanical and enzymatic dissociation. By following the manufacturer's protocol, we used the Single Cell 3' Library and Gel Bead Kit V3.1 (10x Genomics, 1000075) and Chromium Single Cell B Chip Kit (10x Genomics, 1000074) to prepare barcoded scRNA-seq. Paired-end 150 bp reads were then generated by sequencing the libraries on the Illumina NovaSeq6000 platform (performed by CapitalBio Technology, Beijing).

Quality control of the scRNA-seq data
CellRanger (version 4.0.0) was used to generate raw gene expression matrixes from the sequencing data of each sample based on the human reference version GRCh37. The R package Seurat (version 4.0.0) was used for downstream quality control [34]. Only genes expressed at >3 cells and cells with >200 genes detected were kept in the expression matrix. We then removed low-quality cells that met at least one of the following criteria: (i) >5000 expressed genes, (ii) <200 expressed genes, or (iii) >25% UMIs derived from the mitochondrial genome.

Integration of multiple scRNA-seq datasets
The R package Seurat (version 4.0.0) was used for the dataset integration. Briefly, the function NormalizeData was applied to each expression matrix for log-transformation and the function FindVariableFeatures was used to select the top 2000 genes with high cell-to-cell variation. Then, "anchors" between individual datasets were identified using the function FindIntegrationAnchors. Based on these anchors, an unbatched dataset was created using the function IntegrateData.

Unsupervised clustering and dimension reduction
We used the function ScaleData to scale and center the expression of 2000 variable genes, on which principal component analysis (PCA) was subsequently performed. Then, functions FindNeighbors and FindClusters were used for the first-round cluster on the first 20-50 PCs with a resolution of 0.1-0.6. Finally, the main cell clusters were identified with a resolution of 0.3, and nonlinear dimensional reduction was visualized using the two-dimensional uniform manifold approximation and projection (UMAP).

Cell type annotation and cluster marker identification
The major cell types were first annotated using the SingleR (version 1.4.1) R package based on the expression of canonical cell type markers [35]. Preferentially expressed genes in each cluster were identified using the FindAllMarkers function and differential expression testing was performed with the "MAST" method with sample identity as a latent variable [36]. Cells that expressed two canonical cell type markers were classified as doublet cells, which were excluded from further analyses.

CNV estimation of the malignant cells
CD45epithelial cells were analyzed using the inferCNV R package (version 1.11.1) to estimate copy number aberrations in scRNA-seq data [37]. The EPCAMimmune and stromal cells were used as references. A cut-off of 0.1 was used for the min average read counts per gene among reference cells.
The option "random_trees" was used to define the tumor sub-clusters.

Pathway analysis
Gene set variation analysis (GSVA) was performed using the GSVA R package (version 1.38.2) to estimate the enrichment scores of biological pathways [38]. The hallmark and gene ontology (GO) gene sets were obtained using the R package msigdbr (version 7.4.1) and from the MsigDB website (http://software.broadinstitute.org/gsea/msigdb) [39]. The limma R package was then used to evaluate the differential activities of pathways between different clusters or groups (P value < 0.05) [40].
DEGs between the two subtypes of cells were identified using the FindMarkers function of the Seurat R package with the same parameters used for FindAllMarkers. When analyzing the malignant epithelial cells of different subtypes, we overcame the intertumor heterogeneity by keeping only high-confidence DEGs for a subtype such that the expression in each sample had a >0.8 log2 fold change over the average expression of the other subtype. Overrepresentation analysis of gene ontology (GO) biological process (BP) terms was performed using the fora function of the fgsea R package (version 1.16.0) [41].

Diffusion map and pseudotime analysis
PCA embeddings were extracted from the integration produced by Seurat R package. The destiny R package (version 3.1.1) was then used to produce the diffusion trajectory from the predefined PCA embeddings for pseudotime analysis. Naïve T cells were chosen as the root, and the function DPT was used to estimate the diffusion pseudotime (DPT) for single cells.

Cell-cell communication analysis
We used CellphoneDB (version 2.0.0) [42] to explore the potential communication between different cell types based on the expression of ligand-receptor pairs. The receptor or ligand should be expressed by more than 10% of cells in a cluster to be included in the downstream analysis. To assess the significance of a ligand-receptor pair between two clusters, an empirical P value was determined by randomly assigning the cluster labels of each cell for 1,000 times. The R package CellChat (version 1.1.3) [43] was used to characterize the interaction number and strength with default parameters.

Prediction of ICI response
The putative ICI response of patients in the mUC cohort was predicted based on the TIDE as well as the in-house developed Macro-C3 gene signature. For the TIDE evaluation, gene expression in the TPM was supplied and force-normalized. The cancer type was set to "other".
To predict the response to ICI based on the Macro-C3 gene signature, enrichment scores were first estimated based on log-normalized TPM values using ssGSEA. For a cohort of size N comprising a group of R responders, the luminal, basal, cytotoxicity, and Macro-C3 gene signature scores were L = {l1, l2, …, lN}, B = {b1, b2, …, bN}, C = {c1, c2, …, cN}, and M = {m1, m2, …, mN}. The z-score of the cytotoxicity score is: We then fitted a regression line LC taking the log(basal/luminal) score as the independent variable and Z(C) as the response. The same procedure was repeated for Z(M), which produced another regression line, LM intersecting LC at (a, b). This allows us to define an optimum log(basal/luminal) threshold a such that patients below this threshold have a relatively higher Z(C) compared to Z(M) and vice versa. Subsequently, we denoted the Macro-C3 signature as the ICI response score for patients with a log(basal/luminal) score < a: To prove that the log(basal/luminal) threshold a had a clinical implication, we generated 256 evenly distributed cutoff points of the log(basal/luminal) score X = {x1, x2, …, x256}. For a given cut-off point xj, we wrote the fraction of responders as follows: The log(basal/luminal) threshold a was then visualized to compare with xj.
The R package pROC (version 1.17.0.1) was used to generate ROC curves visualizing the true-positive rates versus false-positive rates at various thresholds of the TIDE, Macro-C3 signature, and ICI response scores. The area under the ROC curve (AUC) was adopted as the quality metric of prediction.

Statistical analysis
Statistical analysis was performed using the R (version 4.0.3) and Python (version 3.8.8) and has been described in the Methods section and figure legends. Unless otherwise specified, the statistical analyses were performed in a two-sided manner. P < 0.05, P <