Myeloid Cell–associated Resistance to PD-1/PD-L1 Blockade in Urothelial Cancer Revealed Through Bulk and Single-cell RNA Sequencing

Standard treatment for metastatic urothelial cancer of the bladder has historically been limited to platinum-based chemotherapy. However, the treatment landscape has experienced a major shift with the introduction of several PD-1/PD-L1 immune checkpoint inhibitors (CPI) into the armamentarium (1–5). These therapies are characterized by durable responses, often measured in years, but achieved in only a subset of approximately 15%–25% of patients. This therapeutic profile has led to intensive investigation into dominant mechanisms of intrinsic resistance in pursuit of biomarkers and combination strategies to extend the benefits of CPIs to additional patients.

Responses to CPIs are thought to be predicated on a preexisting antitumor T-cell response restrained because of adaptive immune resistance (6). Measures of T-cell infiltration, IFNγ-related gene signatures, and PD-L1 expression, colloquially referred to as reflecting “hot” or “inflamed” tumors, have all been correlated with response to CPIs in patients with urothelial and other cancers (1, 3). While classifying tumors as “inflamed” versus “noninflamed” provides a convenient framework for conceptualizing the immunobiology underlying sensitivity and resistance to CPIs, this nomenclature fails to distinguish antitumor versus protumorigenic inflammation. Protumorigenic inflammation, a “hallmark of cancer” pathogenesis (7), involves a tumor microenvironment (TME) shaped by activated fibroblasts, endothelial cells, and innate immune cells, particularly myeloid phagocytic cells, which promote cancer growth and progression at least, in part, by impairing antitumor immunity (7–9). Antitumor immunity and protumorigenic inflammation coexist in delicate spatiotemporal balance in individual TMEs, complicating identification of tumors in the clinic which are resistant to CPIs due to the latter and obfuscating identification of cellular populations or signaling interactions for prioritization as therapeutic targets to overcome such resistance.

To identify dominant molecular and cellular interactions in the urothelial cancer TME associated with CPI resistance, we pursued an unbiased approach (Fig. 1). We first used pre-CPI treatment bulk RNA sequencing (RNA-seq) data from a large clinical trial cohort and identified two gene signatures independently associated with CPI outcomes beyond tumor mutational burden (TMB) alone: (i) a signature enriched in adaptive immune response genes and associated with better CPI outcomes, dubbed the adaptive immune response signature and (ii) a signature enriched in inflammation and innate immune genes and associated with worse CPI outcomes, dubbed the protumorigenic inflammation signature. Consistent with the notion that antitumor immunity and protumorigenic inflammation coexist within individual TMEs, the adaptive immune response:protumorigenic inflammation gene expression ratio, coined the 2IR score (adaptive immune:protumorigenic inflammation ratio), had the largest effect on clinical outcomes, and was validated in an independent urothelial cancer clinical trial cohort. We then generated single-cell RNA sequencing (scRNA-seq) data from urothelial cancer bladder specimens to uncover the cellular composition underlying the gene signatures, revealing diverse cellular populations contributing to both the adaptive immune response and protumorigenic inflammation gene signatures. Although the protumorigenic inflammation signature was expressed prominently by myeloid phagocytic cells as a whole, diverse expression of the adaptive immune response and protumorigenic signature genes was observed across individual macrophages/monocytes and neutrophils leading to application of the 2IR score to each individual cell (myeloid single cell immune:protumorigenic inflammation ratio or M_sc2IR score). Myeloid phagocytic cells with low M_sc2IR scores demonstrated upregulation of proinflammatory cytokines (e.g., IL1B) and chemokines (e.g., CCL20) and downregulation of antigen presentation genes, could not be discerned on the basis of M1 versus M2 classification, and were enriched in the pretreatment blood samples of patients with CPI-resistant metastatic urothelial cancer. Thus, the balance of adaptive immunity and protumorigenic inflammation is associated with CPI outcomes in urothelial cancer and resistance associated with the latter may be driven by interactions among diverse cell types in the TME and linked to a proinflammatory cellular state of myeloid phagocytic cells detectable in both the TME and peripheral blood.

Three datasets including bulk RNA-seq data from patients with urothelial cancer were analyzed in this study (Fig. 1; Supplementary Table S1): IMvigor 210, The Cancer Genome Atlas (TCGA) urothelial bladder cancer dataset, and the CheckMate 275 study.

IMvigor 210 was a single-arm phase II study investigating PD-L1 inhibition with atezolizumab (1,200 mg, intravenously, every 3 weeks) in patients with metastatic urothelial cancer (NCT01208652 and NCT02951767). The primary endpoint of the trial was the objective response rate according to RECIST v1.1. Patients with metastatic urothelial cancer progressing despite prior platinum-based chemotherapy, or chemotherapy-naïve patients who were not eligible for cisplatin-based chemotherapy, were eligible. The results of IMvigor 210 have been reported previously (1, 10). Patients enrolled on IMvigor 210 were required to have archival tumor tissue obtained within 2 years of study entry submitted for analysis, which included bulk RNA-seq, as well as targeted next-generation sequencing–based genomic profiling for 395 cancer-related genes (FoundationOne, Foundation Medicine). For this study, RNA-seq data, TMB (“FMOne mutation burden per MB”), objective response rate, and survival data for 348 unique patients were extracted using the R package IMvigor210CoreBiologies (http://research-pub.gene.com/IMvigor210CoreBiologies/).

TCGA bladder cancer dataset includes patients with clinically localized muscle-invasive urothelial cancer of the bladder who underwent radical cystectomy. This cohort has been described in detail previously (11) and RNA-seq data (“Level_3_RSEM_genes_normalized”) for 408 unique patients were downloaded from Firehose (2016_01_28) at the Broad Institute (https://gdac.broadinstitute.org/). The updated clinical data were downloaded from an integrated TCGA pan-cancer clinical data resource (12).

CheckMate 275 was a single-arm phase II study investigating PD-1 inhibition with nivolumab (3 mg/kg, intravenously, every 3 weeks) in patients with metastatic urothelial cancer (NCT02387996). The primary endpoint of the trial was the objective response rate according to RECIST v1.1. Patients with metastatic urothelial cancer progressing despite prior platinum-based chemotherapy were eligible. The results of CheckMate 275 have been reported previously (3). Patients enrolled on CheckMate 275 were required to have archival tumor tissue submitted for analysis, which included bulk RNA-seq and whole-exome sequencing. Patients who were consented for genomic studies and had tumor material that passed quality control (QC) were included in this analysis. RNA-seq and TMB data were provided by Bristol Myers-Squibb and the latter was calculated as the missense mutation count. TMB (n = 139) and/or RNA-seq (n = 72) data, objective response rate, and survival data were available, with both TMB and RNA-seq data available for 54 patients.

For the IMvigor 210 dataset, only genes with a read count >1 in more than 10% of the samples were considered. The raw read count data from the IMvigor and CheckMate 275 datasets were first transformed to RPKM, and then scaled patient-wise such that the 75% quantile of each sample was equal to 1,000 (similar to the RSEM normalization; ref. 13). To facilitate analysis across datasets, only 16,339 genes common among the three datasets were analyzed in this study. Batch effects were removed across the three datasets using R package ComBat (14).

With the goal of identifying consistently coexpressed gene modules (CCGM) associated with survival, we first identified genes nominally associated with better overall survival (OS) outcomes in the IMvigor 210 dataset. A bivariable Cox regression model was used to estimate the association between the expression of each gene, |$Gene_i$|⁠, with the OS conditional on TMB: |$Surv (Event, Time) \sim Gene_i \ \plus \ {\rm {log}}\it (TMB)$|⁠. We identified 1,193 genes for which higher expression was associated with better survival outcomes (nominal P value of two-sided Wald test <0.05). We employed a lenient P value cutoff to be as inclusive as possible at this initial gene selection step, and then identified CCGMs to enrich for true signals and filter out possible noise. A CCGM is defined as a list of genes that are coregulated in multiple datasets. Using weighted correlation network analysis (15) among these 1,193 genes, we identified one coexpression module in the IMvigor 210 dataset (735 genes) and one coexpression module in TCGA dataset (600 genes). Significant overlap (575 genes) existed between the modules identified in the IMvigor 210 and TCGA datasets (P < 1e-16 by two-sided Fisher exact test) and these 575 overlapping genes were considered a CCGM, referred to as the “adaptive immune response” signature.

Next, we identified genes associated with worse survival outcomes conditioned on both TMB and the adaptive immune response signature genes (⁠|$M_{adaptive\ immune})$|⁠. Specifically, we assessed the association of each gene (⁠|$Gene_i$|⁠) with OS using a multivariable Cox regression model |$Surv(Event, Time) \sim Gene_i \ \plus\ M_{adaptive \ immune} \ \plus \ {\rm {log}}(TMB)$|, where |$M_{adaptive \ immune}$| was calculated for each sample by averaging expression of the adaptive immune response signature genes. A total of 1,498 genes were associated with worse survival outcomes (nominal P value of two-sided Wald test <0.05). The weighted correlation network analysis was conducted for the 1,498 genes in both the IMvigor 210 and TCGA datasets, followed by the overlapping analysis analogous to the methodology described for derivation of the adaptive immune response signature which resulted in two CCGMs, that is, the “protumorigenic inflammation” signature (437 genes) and stromal signature (287 genes). A third CCGM (50 genes) enriched with HALLMARK_MYC_targets was not further pursued in this study given its small size. We further updated the “adaptive immune response” signature by excluding genes associated with worse survival in this three-variate Cox regression model (Z > 1.5), resulting in 483 genes in the module.

To investigate the pathways enriched in each CCGM, we compared the signature genes with the HALLMARK and canonical gene sets in the Molecular Signatures Database (software.broadinstitute.org/gsea/msigdb; ref. 16) using Fisher exact test (nominal P value of two-sided test <1e-5).

Calculation of the 2IR score, identification of top-ranked genes, comparison with other biomarkers, and multiplex IHC are detailed in the Supplementary Materials and Methods.

Cox proportional hazard regression models [coxph() function] were performed using the R package survival to evaluate the association between the gene signatures and TMB with OS. When signature expression and TMB were treated as continuous variables, they were standardized to N(0,1) before entering the Cox regression model to estimate HR and confidence interval, and the significance testing was performed by Wald test. When the 2IR score was discretized into tertiles, the R package survminor was used to plot the Kaplan–Meier curve, and the significance testing for differences in OS was performed using the log-rank test. Logistic regression models were performed to evaluate the association between the gene modules and TMB with objective response. In the logistic regression, a complete response (CR) or partial response (PR) was treated as 1, and stable disease (SD) or progressive disease (PD) was treated as 0. The signature expression and TMB were similarly standardized before entering the logistic regression model to estimate the coefficient, and the significance testing was performed by Wald test. All statistical analyses and figures were generated in R version 3.6.3.

Primary urothelial bladder cancer tumor tissue was obtained after obtaining informed consent in the context of an Institutional Review Board (IRB)-approved genitourinary cancer clinical database and specimen collection protocol (IRB #10-1180) at the Tisch Cancer Institute, Icahn School of Medicine at Mount Sinai (New York, NY). Patients undergoing transurethral resection of bladder tumor had a portion of their tumor placed immediately into RPMI medium in the operating room. The specimen was then transferred to the laboratory for further processing. Patients undergoing radical cystectomy and lymph node dissection had their bladder and lymph nodes sent directly to the pathology suite upon completion of the lymph node dissection. The bladder was bivalved and a portion of visible tumor was then placed in media as above. Adjacent normal tissue was identified in a subset of specimens on the basis of visual inspection. The specimen was then transferred to the laboratory for further processing.

Tissue specimens were processed immediately upon receipt and dissociated into single-cell suspensions using the GentleMACS Octodissociator with kit matched to the tissue type (Miltenyi Biotec) following the manufacturer's instructions. scRNA-seq was performed on these samples using the Chromium Platform (10x Genomics) with the 3′ gene expression V3 kit, using an input of approximately 10,000 cells. Briefly, Gel-bead in emulsions (GEM) were generated on the sample chip in the Chromium controller. Barcoded cDNA was extracted from the GEMs by post-GEM RT-cleanup and amplified for 12 cycles. Amplified cDNA was fragmented and subjected to end-repair, poly A-tailing, adapter ligation, and 10 ×–specific sample indexing following the manufacturer's protocol. Libraries were quantified using Bioanalyzer (Agilent Technologies) and QuBit (Thermo Fisher Scientific) analysis. Libraries were sequenced in paired-end mode on a NovaSeq Instrument (Illumina) targeting a depth of 5 × 10⁴–1 × 10⁵ reads per cell. Sequencing data were aligned and quantified using the Cell Ranger Single-Cell Software Suite (version 3.0, 10x Genomics) against the provided GRCh38 human reference genome.

Seurat (version 3.0; ref. 17) was used to process the scRNA-seq data. After filtering cells with a high percentage (>20%) of mitochondrial reads and cells with <200 or >6,000 genes detected, as well as potential doublets uncovered during subsequent analysis steps, 19,708 cells from 10 samples and 22,175 genes with nonzero read counts in >5 cells were included for further analysis.

After the read count data were log normalized, the most variable 2,000 genes were selected. Then the effect of the unique molecular identifier count and percentage of mitochondria per cell were regressed out, followed by dimensionality reduction using principal component analysis. Finally, the cells were clustered using the K-nearest neighbors graph-based methods as implemented in Seurat (with the top 20 principal component and resolution = 0.5). Cells were grouped into nine major cell clusters on the basis of the canonical cell type–specific markers: T/NK (“CD3E”), B/plasma (“MS4A1,” “MZB1,” and “CD79A”), DC (dendritic cell; “HLA-DQA1” and “HLA-DQB1”), mast (“MS4A2”), macrophage/monocyte (“C1QA” and “LYZ”), endothelial (“PLVAP”), fibroblast related (“DCN” and “ACTA2”), epithelial (“KRT19”), and neuronal cells (NNAT). The identification of minor cell clusters, M_sc2IR score calculation, and NicheNet analysis are detailed in the Supplementary Materials and Methods.

For each of the nine major cell clusters, we identified cell type–overexpressed genes using FindAllMarker() function in Seurat package (with default parameters). The overlap between each of the adaptive immune response, protumorigenic inflammation, and stromal signature genes and overexpressed genes among the major cell clusters were assessed using OR and P value (two-sided Fisher exact test).

scRNA-seq data for 10 frozen peripheral blood mononuclear cell (PBMC) samples derived from pretreatment peripheral blood of 5 patients with metastatic urothelial cancer who achieved an objective response to treatment with atezolizumab and 5 patients with metastatic urothelial cancer who did not achieve an objective response to treatment with atezolizumab in the setting of the IMvigor 210 study were downloaded from Gene Expression Omnibus (GSE145281). The peripheral blood scRNA-seq cohort and analysis are detailed in the Supplementary Materials and Methods.

To identify molecular features associated with survival in CPI-treated patients with metastatic urothelial cancer, we utilized bulk RNA-seq and TMB data from the IMvigor 210 study, a large single-arm phase II trial testing the PD-L1 inhibitor, atezolizumab (Fig. 1A; refs. 1, 10, 18). This cohort has been described previously and additional details are provided in Supplementary Table S1; RNA-seq and TMB data were available for 348 and 272 patients, respectively (18). We pursued stepwise identification of CCGMs, which focused on identifying gene modules associated with OS and utilizing gene modularity to enrich for true signals (Fig. 1B; Supplementary Fig. S1; see Materials and Methods). Given the correlation between TMB and response to CPI in urothelial cancer (18, 19), we explored genes associated with OS conditioning on TMB (see Materials and Methods) and identified a signature consisting of 1,193 genes associated with longer OS. To refine this signature, we performed meta-analysis of coexpression patterns (20, 21) using both the IMvigor 210 and TCGA urothelial cancer datasets and identified a CCGM comprising of 483 genes (see Materials and Methods; Supplementary Data S1). Gene set enrichment analysis revealed that this module was highly enriched in adaptive immune response–related genes (Fig. 1C; Supplementary Fig. S2) and was, therefore, labeled the adaptive immune response signature.

In the second step, we further analyzed the IMvigor 210 dataset to identify genes associated with survival conditioning on both TMB and the adaptive immune response signature (Fig. 1B). We identified 1,498 genes associated with shorter OS. We again applied meta-analysis of coexpression patterns (22, 23) in the IMvigor 210 and TCGA urothelial cancer datasets and identified two CCGMs for further analysis. The first module consisted of 437 genes, was enriched in inflammation and innate immune genes (Fig. 1C; Supplementary Fig. S2; Supplementary Data S1), and associated with shorter OS and was, therefore, labeled the protumorigenic inflammation signature. The second module associated with shorter OS, consisted of 287 genes, was enriched in epithelial-to-mesenchymal transition (EMT)- and extracellular matrix–related genes (Fig. 1C; Supplementary Data S1), and consistent with our prior work (24) was named the stromal signature. Importantly, expression of the adaptive immune response and stromal signatures was both positively correlated with the adaptive immune response module (Supplementary Fig. S3) such that their disparate impact on OS was only revealed using our stepwise approach (Supplementary Fig. S4) and suggested that the balance of these features in individual tumors may impact CPI outcomes.

We next sought to define the independent contribution of the three gene signatures to outcomes with CPIs in the IMvigor 210 cohort. Models combining both the adaptive immune response and protumorigenic inflammation signatures (see Materials and Methods), particularly the 2IR score, demonstrated the largest effect size on OS and similar findings were observed with objective response rate (Fig. 2A–C; Supplementary Table S2). Importantly, when both the protumorigenic inflammation and stromal signatures were entered into a multivariable model along with the adaptive immune response signature and TMB, the stromal module was no longer independently associated with OS (Fig. 2A; Supplementary Table S2). These findings indicated that (i) the balance of cellular and molecular events underlying the adaptive immune response versus protumorigenic inflammation signatures within an individual urothelial cancer TME may dictate outcomes with CPIs and (ii) the negative impact of the stromal signature on outcomes may be more indirect and mediated via the events underlying the protumorigenic inflammation signature (Fig. 2A).

Given the potential practical advantages of smaller sets of genes for validation and clinical biomarker development, we identified the top-ranked genes within the three signatures (see Materials and Methods; Supplementary Table S3; Supplementary Fig. S5). Signature scores derived from these smaller gene sets demonstrated similar associations with OS compared with scores derived from the full gene sets (Supplementary Table S4).

For validation, we utilized TMB and RNA-seq data from the CheckMate 275 study, a single-arm phase II trial evaluating the PD-1 inhibitor, nivolumab, in patients with metastatic urothelial cancer (Fig. 1D; ref. 3). This cohort has been described previously, with further detail provided in Supplementary Table S1; RNA-seq and TMB data were available for 72 and 139 patients, respectively (3). The adaptive immune response, protumorigenic inflammation, and stromal gene signatures demonstrated similar associations with OS, progression-free survival, and response rate in the CheckMate 275 cohort (Supplementary Table S5). As observed in the IMvigor 210 cohort, the 2IR score in the CheckMate 275 cohort demonstrated the largest effect on CPI outcomes (Fig. 2E and F). Furthermore, the 2IR score remained significantly associated with OS after several clinical prognostic factors were taken into consideration (Supplementary Results; Supplementary Table S6).

Other cancer cell intrinsic and TME-related features have been associated with CPI outcomes in urothelial cancer and other tumors (18, 25–27). The 2IR score demonstrated favorable performance characteristics relative to such features, including PD-L1 protein expression, the tumor immune dysfunction and exclusion and CD8 effector T-cell gene signatures, ARID1A mutation status, and CXCL13, TGFB1, or CXCL8 (IL8) gene expression (Fig. 2D and G; see Supplementary Fig. S7 for correlation between these features and the 2IR score), in both the IMvigor 210 and CheckMate 275 cohorts. Less dramatic findings were observed applying the 2IR score to TCGA urothelial cancer dataset, a cohort of patients with muscle-invasive urothelial cancer of the bladder treated with curative intent cystectomy (Supplementary Fig. S6), suggesting that the 2IR score may impart predictive rather than solely prognostic information. Taken together, the 2IR score, representing the balance of expression of the adaptive immune response and protumorigenic inflammation gene signatures within individual TMEs, was associated with objective response and OS in CPI-treated patients with metastatic urothelial cancer in two clinical trial cohorts and conveys information beyond that achieved with previously identified features.

We next sought to characterize the relationship between the 2IR score and the cellular organization of the urothelial cancer TME, particularly the spatial localization of T cells, based on prior work from our group and others linking T-cell spatial localization with CPI resistance in urothelial cancer (18, 24). We employed a tissue profiling approach known as multiplexed IHC consecutive staining on a single slide (MICSSS; refs. 28, 29) on a subset of 19 specimens from the CheckMate 275 cohort with matched RNA-seq data. Notably, MICSSS revealed that specimens with higher 2IR scores exhibited occasional tertiary lymphoid-like structures (Fig. 3A and B; Supplementary Fig. S8) consistent with prior findings linking such structures with improved CPI outcomes (30, 31). To quantify the localization of T cells according to gene signature expression, we defined cancer cell and stromal zones on the basis of pan-cytokeratin staining using a machine learning segmentation tool and examined CD8⁺ expression in 76 regions of interest across the 19 specimens (see Materials and Methods; Fig. 3E). Lower 2IR scores correlated with decreased CD8⁺ T-cell enumeration in cancer cell nests, with T cells more prominently localized to the stromal regions, suggestive of a T-cell “excluded” phenotype (Fig. 3C–G; refs. 18, 24). These findings suggested that CPI resistance associated with lower 2IR scores may be related to impairment of T-cell trafficking and/or function, prompting us to further probe the cellular populations and interactions underlying the gene signatures.

Our gene signatures were derived from bulk RNA-seq data from archival urothelial cancer specimens obtained pretreatment with CPIs, the vast majority of which represented invasive primary tumors (Supplementary Table S1). Therefore, to map the cellular origins of the three gene signatures (Fig. 4A), we performed droplet-based encapsulation scRNA-seq on an analogous set of eight freshly resected invasive urothelial cancer bladder specimens, as well as two specimens derived from adjacent grossly normal urothelium using the 10x Genomics Chromium System (see Materials and Methods; cohort characteristics in Supplementary Table S7). The characteristics of the cohort are detailed in Supplementary Table S7. After excluding cells not passing QC (see Supplementary Fig. S9 for QC plots), 19,708 cells from the 10 samples were analyzed. A median of 1,456 genes were detected per cell. We performed graph-based clustering as implemented in the Seurat package (17). A two-stage clustering approach was employed in which cells were first grouped into major populations and subsequently further partitioned into minor populations (see Materials and Methods).

Canonical marker genes revealed nine major cell populations identified by scRNA-seq, including T and natural killer (NK) cells, B cells, myeloid lineage cells, nonhematopoietic stromal cells, and epithelial cells (Fig. 4B and C; Supplementary Data S2). To determine the origins of the adaptive immune response, protumorigenic inflammation, and stromal signatures, the expression pattern of the signature genes was assessed among these major cell populations (Fig. 4D). While adaptive immune response gene expression was prominent among T and B cells (OR = 7.65 and 11.19, respectively) and protumorigenic inflammation gene expression was prominent among monocytes/macrophages (OR = 12.49), diversity in expression of signature genes among cell populations was observed (Fig. 4D). For example, there were protumorigenic inflammation signature genes overexpressed in T cells (OR = 5.01) and adaptive immune response signature genes overexpressed in monocytes/macrophages (OR = 2.55). Notably, expression of stromal signature genes demonstrated much less cellular diversity and was detected predominantly from cancer-associated fibroblast (CAF; OR = 22.03) and endothelial (OR = 2.78) cell populations.

To determine whether the adaptive immune response and protumorigenic inflammation signature genes expressed by a given major cell type were arising from discrete cellular subpopulations (e.g., adaptive immune response genes from one subset of macrophages and protumorigenic inflammation genes from another subset of macrophages), we subjected each major cluster to a second round of partitioning revealing a total of 50 minor cell clusters (described in detail in Supplementary Results; Supplementary Fig. S10–S16). Unexpectedly, we observed expression of both adaptive immune response and protumorigenic signature genes within most minor cell populations (Fig. 4E). Hence, the adaptive immune response and protumorigenic inflammation gene signatures are contributed by diverse cell types within the TME and may be linked to underlying cellular states rather than discrete cellular subpopulations.

Myeloid phagocytic cells demonstrated the most prominent expression of the protumorigenic inflammation signature genes that were associated with CPI resistance in our clinical trial cohorts, yet also expressed some adaptive immune response signature genes (Fig. 4E). Consequently, the expression level of protumorigenic inflammation signature per cell [as assessed by the AddModuleScore() function in the Seurat package] was highest in myeloid phagocytic cells (Fig. 4F). In addition, compared with other types of cells in the scRNA-seq dataset, myeloid phagocytic cells had much higher variance in the expression of the protumorigenic inflammation signature genes (also known as heterogeneity of molecular state; Fig. 4F), yet comparable variance in adaptive immune signature genes (Supplementary Fig. S4F). Thus, we turned further attention to the myeloid phagocytic cells. We identified seven minor monocyte/macrophage populations and one neutrophil population by scRNA-seq analysis (Fig. 5A and B; Supplementary Fig. S17A). While some of these minor cell populations demonstrated higher expression of M1 versus M2 signature genes (Fig. 5C), or vice versa, heterogeneity of monocyte/macrophage minor populations was observed beyond classical M1 and M2 polarization as has been documented in prior analyses (32, 33). The macrophage populations resembled previously described “tumor-associated macrophage (TAM)-like macrophages” with increased expression of APOE, C1QA, C1QB, SLC40A1, and TREM2 (32, 34). Two populations demonstrated higher expression of S100A family genes, but lower M1 and M2 signature gene expression, and were annotated as monocyte-Jun and monocyte-LYZ. These clusters resembled previously described “myeloid-derived suppressor cell (MDSC)-like macrophages” with overexpression of THBS1, S100A8, FCN1, and VCAN (32). A population annotated as MM-CCL2 shared marker genes with both “TAM-like” and “MDSC-like” macrophages.

We next sought to better characterize the myeloid cells that might be linked to CPI resistance. Myeloid phagocytic cells are highly plastic, educated by cellular and signaling interactions in the TME, and play diverse roles in promoting and restraining anticancer immunity (35). Our single-cell characterization of urothelial cancer specimens revealed diversity in expression of the protumorigenic inflammation and adaptive immune response signature genes across individual macrophages/monocytes and neutrophils (Fig. 5D). Intrigued by the observation that the balance of protumorigenic inflammation and adaptive immunity plays a key role in bulk urothelial cancer specimens, and reasoning that such balance might be relevant at the level of individual cells, we extended this concept to the monocyte/macrophage and neutrophil population and calculated a 2IR score for each individual cell, termed the myeloid single-cell 2IR score (M_sc2IR score; see Fig. 5D and E and Materials and Methods). While low M_sc2IR score cells were observed across all minor myeloid phagocytic cell populations and were not correlated with M1 or M2 signatures (Supplementary Fig. S17B), these cells were highly enriched in the MM-CCL20 minor population (OR = 11.0; P < 1e-16 by Fisher exact test) and underrepresented in the macrophage-C1QA minor population (OR = 0.14; P < −16 by Fisher exact test; Fig. 5F).

Differential gene expression and gene set enrichment analysis of myeloid phagocytic cells with low versus high M_sc2IR scores revealed upregulation of proinflammatory pathways and top-ranking genes, such as IL1B, CXCL8 (IL8), SPP1, and CCL20, in the former, while the latter demonstrated upregulation of genes and pathways related to antigen presentation and the T-cell recruiting chemokines CXCL9 and CXCL10 (Fig. 5E and G; Supplementary Fig. S18).

Protumorigenic monocytes derived from patients with renal cancer have been shown previously to express proinflammatory cytokines and chemokines, including IL1B, CCL20, and CXCL8 (IL8), via an IL1β-dependent mechanism (36). We sought to define putative therapeutic targets implicated in polarizing myeloid phagocytic cells with low M_sc2IR scores and not restrict our analysis to genes overexpressed in low M_sc2IR score cells, but rather seek upstream ligands implicated in driving the expression of such genes. We, therefore, used NicheNet (37), an approach that predicts ligands that modulate target gene expression by leveraging prior knowledge of signaling pathways and transcriptional regulatory networks (see Materials and Methods). Indeed, this analysis revealed that IL1α and IL1β were the top-ranked ligands inferred to regulate genes overexpressed in low M_sc2IR score cells (Fig. 5H). Both IL1α and IL1β were also predominantly expressed by myeloid phagocytic cells in our single-cell cohort (Supplementary Fig. S19).

Thus, the M_sc2IR score, reflecting the balance of adaptive immune response and protumorigenic inflammation gene expression in individual myeloid phagocytic cells, may reflect the plasticity of these cells in the TME (Fig. 5E). Low M_sc2IR score monocytes/macrophages and neutrophils, with upregulation of proinflammatory genes and downregulation of antigen presentation genes and not delineated by classical M1 versus M2 polarization or graph-based unsupervised cell clustering, may define a cellular state of myeloid phagocytic cells contributing to CPI resistance.

We next asked whether similar heterogeneity in M_sc2IR scores was present in monocytes in the peripheral blood of patients with metastatic urothelial cancer and whether these populations were associated with CPI resistance. scRNA-seq data from PBMCs collected prior to the initiation of treatment with anti-PD-L1 CPI from 5 patients who achieved an objective response, and 5 patients who did not achieve an objective response, were utilized (see Materials and Methods; Supplementary Fig. S20). We calculated M_sc2IR scores in individual monocytes identifying low, intermediate, and high M_sc2IR score populations. Monocytes with low M_sc2IR scores were significantly enriched in the peripheral blood of patients with CPI-resistant versus CPI-responsive metastatic urothelial cancer (Fig. 6A; P = 0.0048 by two-sided t test). Alternatively, the 5 patients who responded to CPI could not be readily distinguished from the 5 patients with CPI-resistant metastatic urothelial cancer using monocyte minor populations identified by graph-based unsupervised cell clustering (Fig. 6B), individual genes such as CXCL8 (IL8; Fig. 6C), or M1 and M2 signatures (Supplementary Fig. S21). Similar to our findings in the urothelial cancer TME, low M_sc2IR score monocytes in the pretreatment peripheral blood of patients with metastatic urothelial cancer demonstrated upregulation of proinflammatory genes and downregulation of antigen presentation genes (Fig. 6D) and IL1α and IL1β were the top-ranked ligands inferred to regulate this gene expression program. Therefore, low M_sc2IR score myeloid cells are present in both the TME and peripheral blood of patients with urothelial cancer and are associated with CPI resistance.

Protumorigenic inflammation is recognized as a “hallmark of cancer” pathogenesis (7, 38, 39). However, antitumor immunity and tumor-promoting inflammation coexist in delicate balance complicating dissecting the role of the latter in mediating CPI resistance in studies using human specimens. Using an unbiased approach, we identified an adaptive immune response gene signature associated with better CPI outcomes and protumorigenic inflammation and stromal gene signatures associated with worse CPI outcomes in patients with metastatic urothelial cancer. We further demonstrated that: (i) expression of the three gene signatures was positively correlated with one another, consistent with the coexistence and balance between antitumor immunity and tumor-promoting inflammation, (ii) the stromal gene signature, linked to activated CAFs, did not convey independent information related to CPI outcomes beyond the protumorigenic inflammation signature, suggesting a more indirect role (e.g., recruitment and education of myeloid cells), (iii) the 2IR score, reflecting the balance of antitumor immunity and tumor-promoting inflammation in individual TMEs, best correlated with CPI outcomes and was reflective of diverse cell types in the TME, and (iv) low M_sc2IR score myeloid phagocytic cells were characterized by increased expression of proinflammatory genes and decreased expression of antigen presentation genes, could not be discerned on the basis of classical M1 versus M2 polarization, and were enriched in the pretreatment blood samples of patients with metastatic urothelial cancer resistant to CPI. Together, our findings define a myeloid phagocytic cell state associated with CPI resistance, highlight potential approaches to identify patients potentially best suited for therapies seeking to overcome protumorigenic inflammation–related CPI resistance, and delineate putative therapeutic targets for further study.

Our overarching goal was to identify dominant clinically relevant features correlated with CPI resistance that might be linked to underlying immunobiology and associated therapeutic targets for prioritization for further preclinical and clinical testing as CPI-based combination strategies. In addition, with further refinement and validation, the identified tissue- and blood-based features could prove valuable in establishing proof of concept in early-phase clinical development of combination regimens targeting myeloid-related CPI resistance, through associations with clinical outcomes and/or pharmacodynamics monitoring (e.g., serial changes in low Msc2IR score monocytes in peripheral blood).

Other studies have linked aspects of the TME with response/resistance to CPI in urothelial cancer and other malignancies. Features reflecting an activated tumor stroma, including EMT- and TGFβ-related gene signatures, have been correlated with poor outcomes with CPI treatment (18, 24). In this study, our stromal gene signature was no longer independently associated with CPI outcomes when the protumorigenic inflammation signature was included in multivariable models, suggesting the former may play a more indirect role. Multiple studies have correlated T-cell gene signatures, or related measures of adaptive immune resistance, with sensitivity to CPIs (1, 3, 40). However, analyses utilizing human tumor specimens demonstrating an independent association between gene signatures, reflecting protumorigenic inflammation, or related myeloid cells, and CPI resistance, have been much more limited. This disconnect is likely, in part, related to the positive correlation between gene signatures, reflecting the presence of T cells, and other stromal and immune cells, despite a disparate impact on outcomes complicating dissecting the role of the latter and due to the plasticity of immune cells in the TME. Recently, elevated levels of serum/plasma IL8, traced primarily to the myeloid cell compartment, were associated with decreased efficacy of CPIs across several tumor types (26, 41). These studies likely relate to similar biology exposed in our analysis, although we extend these findings by: (i) contextualizing the importance of protumorigenic inflammation in contributing to CPI resistance in urothelial cancer relative to other features defined through an unbiased approach, (ii) providing a refined understanding of the myeloid cellular state associated with protumorigenic inflammation characterized by expression of a number of inflammatory cytokines and chemokines, beyond IL8 alone, suggesting that targeting upstream regulators may be required for optimal therapeutic modulation, and (iii) defining tumor tissue- and blood-based measures using bulk or scRNA-seq to identify tumors for which protumorigenic inflammation may be contributing to CPI resistance in urothelial cancer.

Myeloid phagocytic cells have been linked to suppression of antitumor immunity across a range of malignancies via a variety of mechanisms, although clinically tractable approaches to target myeloid cell–related CPI resistance have remained elusive (42–44). IL1 was among the top-ranked ligands inferred to regulate the low M_sc2IR score myeloid phagocytic cell gene program, in line with prior experimental data demonstrating that inflammatory cytokine and chemokine production from protumorigenic monocytes in patients with renal carcinoma was IL1β dependent (36). IL1β has been considered a “master regulator” of inflammation involved in the tumor-promoting and immunosuppressive function of myeloid cells, anti-IL1β combined with anti-PD-1 therapy abrogated tumor growth in model systems, and anti-IL1β has been associated with lower cancer mortality in human studies (36, 45–47). Building on this collective work, our findings raise the hypothesis that targeting IL1 may reverse the inflammatory phenotype of low M_sc2IR score myeloid phagocytic cells and may represent a rational combination strategy to overcome CPI resistance in a defined subset of patients with urothelial cancer. Additional studies are required to refine the role of IL1α versus IL1β in this context, although IL1β is not present in cells from healthy individuals, and is a product of limited cell types, such as myeloid phagocytic cells, whereas IL1α is more ubiquitously expressed. Clinical trials combining CPIs and anti-IL1 therapies have already been initiated (NCT03631199 and NCT03742349). Furthermore, the gene expression program of low M_sc2IR score myeloid phagocytic cells in the TME reveals several additional putative targets previously linked to inflammatory disorders worthy of further investigation (48, 49).

There are potential limitations to our study. While our study is among the first to characterize the urothelial cancer TME at the single-cell RNA level, the quantity of single cells from each specimen was variable, with two specimens contributing a large proportion of cells; a larger cohort is required to establish a definitive cellular atlas of urothelial cancer specimens. Still, the main goal of our scRNA-seq cohort in this study was to uncover the cellular origins of our gene signatures derived from bulk RNA-seq data. Features of urothelial cancer cells are likely associated with sensitivity and resistance to CPIs. However, beyond TMB, on which our three gene signatures were conditioned, this analysis was focused on the TME given the expression of the module genes when projected onto our scRNA-seq data. Cancer cell intrinsic features that contribute to immune escape and ultimately shape the protumorigenic inflamed TME require further study. Although we linked low Msc2IR score monocytes in peripheral blood with resistance to CPIs in patients with metastatic urothelial cancer, we did not have paired scRNA-seq data from the matched primary tumors to directly explore the association between the TME and circulating immune cells. Together, these considerations underscore the need for additional studies of urothelial cancer specimens profiled at single-cell resolution and linked to CPI treatment outcomes.

Our study identified and validated key gene signatures associated with sensitivity or resistance to CPIs in patients with metastatic urothelial cancer related to adaptive immunity and protumorigenic inflammation, defined the 2IR score as reflecting such balance in individual urothelial cancer TMEs, established the M_sc2IR score as reflecting the cellular state of myeloid phagocytic cells linked to CPI resistance, and identified putative therapeutic targets to overcome resistance. Future work exploring the 2IR and M_sc2IR scores in clinical trials seeking to overcome myeloid-related CPI resistance, further defining and credentialing “master regulators” of low M_sc2IR score myeloid cell polarization as putative therapeutic targets, and dissecting the dominant mechanisms of immune suppression related to these cells, may help facilitate extension of the benefits of CPIs to additional patients with urothelial cancer.

L. Wang reports a patent for tech no. (P-200303-US1) pending. R. Sebra reports other from Sema4 outside the submitted work. A.M. Farkas reports grants from NIH during the conduct of the study. S. Gnjatic reports personal fees from Merck and OncoMed, and grants from Regeneron, Immune Design, Agenus, Genentech, Takeda, Pfizer, and Janssen R&D outside the submitted work. W.K. Oh reports other from Sema4 and personal fees from Astellas, Sanofi, AstraZeneca, Janssen, Merck, and Bayer outside the submitted work. P. Szabo reports other from BMS outside the submitted work. M. Wind-Rotolo reports other from Bristol Myers Squibb during the conduct of the study and outside the submitted work. E. Schadt reports grants from NIH during the conduct of the study and other from Sema4 outside the submitted work. P. Sharma reports personal fees from Achelois, BioAtla, Jounce, Neon, Glympse, Infinity Pharma, Lava, Lytix, Earli, Polaris, Dragonfly, Oncolytics, Hummingbird, and Marker Therapeutics outside the submitted work. N. Bhardwaj reports grants from NIH during the conduct of the study, as well as grants from Regeneron, Harbour Biomedical, and Dragonfly Therapeutics, grants and personal fees from Parker Institute for Cancer Immunotherapy, and personal fees from Neon Therapeutics, Novartis, Avidea, Boehringer Ingelheim, Rome Therapeutics, Roswell Park Comprehensive Cancer Center, BreakBio, Carisma Therapeutics, CureVac, Genotwin, BioNTech, Gilead, Tempest Therapeutics, and Cancer Research Institute outside the submitted work. J. Zhu reports grants from NCI during the conduct of the study, other from Sema4 outside the submitted work, and plans for filing a patent pending. M.D. Galsky reports grants from Bristol Myers Squibb during the conduct of the study, as well as grants from Novartis and Dendreon, grants and personal fees from Bristol Myers Squibb, Merck, Genentech, and AstraZeneca, and personal fees from Pfizer, EMD Serono, Seattle Genetics, Janssen, Numab, Dragonfly, GlaxoSmithKline, Basilea, UroGen, and Rappta Therapeutics outside the submitted work, and has a patent for immune checkpoint inhibitor resistance in urothelial cancer pending. No disclosures were reported by the other authors.

L. Wang: Conceptualization, data curation, software, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. J.P. Sfakianos: Resources, data curation, investigation, writing–review and editing. K.G. Beaumont: Resources, data curation, investigation, writing–review and editing. G. Akturk: Conceptualization, data curation, software, formal analysis, investigation, visualization, methodology, writing–original draft, writing–review and editing. A. Horowitz: Data curation, writing–review and editing. R.P. Sebra: Resources, data curation, writing–review and editing. A.M. Farkas: Data curation, investigation, writing–review and editing. S. Gnjatic: Resources, data curation, investigation, writing–review and editing. A. Hake: Data curation, writing–review and editing.S. Izadmehr: Data curation, investigation, writing–review and editing. P. Wiklund: Investigation, writing–review and editing. W.K. Oh: Investigation, writing–review and editing. P.M. Szabo: Data curation, writing–review and editing. M. Wind-Rotolo: Data curation, writing–review and editing. K. Unsal-Kacmaz: Data curation, writing–review and editing. X. Yao: Investigation, writing–review and editing. E. Schadt: Investigation, writing–review and editing. P. Sharma: Investigation, writing–review and editing. N. Bhardwaj: Conceptualization, resources, data curation, supervision, funding acquisition, validation, investigation, project administration, writing–review and editing. J. Zhu: Conceptualization, resources, supervision, funding acquisition, investigation, visualization, methodology, writing–original draft, project administration, writing–review and editing. M.D. Galsky: Conceptualization, resources, data curation, supervision, funding acquisition, investigation, visualization, methodology, writing–original draft, project administration, writing–review and editing.

The results shown here are, in part, based upon data generated by TCGA research network: http://cancergenome.nih.gov/. We thank Jill Gregory at the Icahn School of Medicine Instructional Technology Group for assistance with graphic design. This work was supported by CA196521 (to M.D. Galsky, N. Bhardwaj, and W.K. Oh) and NIH Loan Repayment Program (to S. Izadmehr).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.