Dynamic single-cell RNA sequencing identifies immunotherapy persister cells following PD-1 blockade

Drug-tolerant persister cells (DTPCs) mediate resistance to oncogene-targeted cancer therapies by adopting a quiescent and apoptosis-resistant state (1–3). Analogous to the concept of “antibiotic-tolerant” bacterial persister cells, these preexisting subpopulations of cancer cells evade drug killing by slowing cell growth and then recovering upon drug withdrawal (4). Furthermore, DTPCs commonly avoid drug-induced cell death by adopting features of the epithelial-mesenchymal transition (EMT) (5) and activating prosurvival cytokine signaling pathways such as IL-6 (6, 7).

Discovery of similar subpopulations that evade immune checkpoint blockade has been hampered by a lack of functional models that recapitulate the tumor immune microenvironment (TIME) because characterization of persister cells in vivo is challenging and requires accurate repeat biopsies (8, 9). However, emerging evidence suggests that EMT may be linked to an immunosuppressive microenvironment that confers resistance to immune checkpoint blockade. For example, breast carcinoma models with high expression of the transcription factor Snail and activation of a mesenchymal transcriptional program promote an immune suppressive microenvironment and resistance to CTLA4 blockade (10). Similarly, activation of Snail in murine KRAS-driven lung cancer cells fosters infiltration of neutrophils, which can promote resistance to PD-1 blockade (11, 12). Yet, it remains unknown whether preexisting cancer cell populations avoid cancer immunotherapy response by adopting this mesenchymal cell state.

Features of EMT have also been linked with a tumor-initiating, stem cell–like state in cancer (13, 14). Increasing evidence suggests that activation of EMT-promoting transcription factors can maintain cells in a hybrid state that is intermediate between a full epithelial or mesenchymal program, promoting cancer stem cell–like behavior (15), which could also be involved in immune evasion. Indeed, tissue stem cells from multiple origins can adopt a quiescent cell state that avoids immune recognition and is likely an evolutionary adaptation to protect these critical cells (16).

Emerging technologies such as single-cell RNA-Seq (scRNA-Seq) and enhanced primary tumor culture methods provide a new opportunity to probe cancer biology in a high-resolution and dynamic manner (17, 18). A major challenge for technologies such as organoid culture has been to maintain original features of the TIME during the time needed to expand cancer cells. Recently, we demonstrated the feasibility of assessing ex vivo response to immune checkpoint blockade using short-term microfluidic culture of patient- or murine-derived organotypic tumor spheroids (PDOTS/MDOTS) embedded in a 3D collagen matrix (19). Detailed characterization of these spheroids from patients with multiple tumor types and mice with a variety of syngeneic tumor backgrounds revealed that they retained many of the key components of the TIME, including PD-1–expressing T cells, multiple myeloid cell populations and dendritic cells, as well as fibroblasts and endothelial cells. Although we observed specific CD8⁺ T cell–dependent killing of MDOTS derived from well-established syngeneic tumor models such as MC38, phenotypic analyses were limited to live/dead cell staining and cytokine profiling of conditioned media (CM). Here, we explored the potential of bulk RNA-Seq and scRNA-Seq of the dynamic response to PD-1 blockade in this system to uncover potentially more fundamental biology of immune escape after CD8⁺ T cell activation.

Bulk RNA-Seq analysis of MDOTS after PD-1 blockade uncovers a distinct transcriptional state. To examine dynamic cellular transcriptomic changes after ex vivo culture and immune stimulation of well-defined MC38 MDOTS, we isolated bulk RNA from collagen-embedded spheroids after 6 days (Figure 1A). MC38 MDOTS were treated during this time period with isotype IgG control, IFN-γ, anti–PD-1 (αPD-1), or αPD-1 + CD8 neutralizing antibody (αCD8). As expected, this resulted in specific αPD-1–induced cell death that was CD8⁺ T cell dependent (Figure 1, B and C). Bulk RNA-Seq analysis of IFN-γ–treated MDOTS revealed highly significant upregulation of genes associated with IFN response, such as Gbp2, Gpb3, Irf1, and Cxcl10, compared with IgG (Figure 1D). Indeed, gene set enrichment analysis (GSEA) identified both IFN-γ and IFN-α signatures as the top 2 Hallmark pathways induced by IFN-γ treatment (Figure 1D and Supplemental Figure 1A; supplemental material available online with this article; https://doi.org/10.1172/JCI135038DS1), capturing the expected transcriptomic response and confirming the utility of this approach.

Figure 1

Bulk RNA-Seq of MDOTS reveals a unique transcriptomic response to immune checkpoint blockade. (A) Schematic of the workflow to profile surviving cells after treatment. αPD-1, anti–PD-1 therapy. (B) Representative immunofluorescent microscopy images of live (AO)/dead (PI) staining of MC38 MDOTS in the microfluidic devices after 6 days of treatment (scale bar: 200 μm). αCD8 = CD8 neutralizing antibody (n = 3). (C) Quantitative measurement of tumor cell viability. Data are mean ± SEM and were analyzed by multiple t tests with Bonferroni’s correction (n = 3). (D) Volcano plot of genes differentially expressed by MC38 MDOTS after 6 days of IFN-γ treatment compared with IgG isotype control. Genes with a log₂(fold change) greater or less than 0.5 are shown in red. The top 5 most significant Hallmark gene sets are shown. (E) Volcano plot of genes differentially expressed by MC38 MDOTS after 6 days of αPD-1 treatment compared with IgG isotype control. Genes with a log₂(fold change) greater or less than 0.5 are shown in red. The top 5 most significant Hallmark gene sets are shown. (F) Two-dimensional t-distributed stochastic neighbor embedding (t-SNE) plots of single-cell RNA-Seq (scRNA-Seq) performed on CD45/CD90⁺ and CD45/CD90^– cells isolated from MC38 tumors at the end of in vivo treatment with either αPD-1 or IgG isotype control. t-SNE plots are colored based on treatment condition (left) and type of cells (middle). Right, t-SNE plot showing projection of bulk EMT signature (E-M) onto scRNA-Seq clusters. *P < 0.05, **P < 0.01, ***P < 0.001.

We next analyzed bulk RNA-Seq data generated after αPD-1 treatment. Unsupervised gene expression clustering revealed that the bulk transcriptome of cells that survived ex vivo PD-1 blockade was highly distinct compared with IFN-γ–treated cells (Supplemental Figure 1B). Furthermore, these transcriptomic changes were completely abolished by αPD-1 + αCD8 treatment, which clustered most closely with IgG control treatment (Supplemental Figure 1B). Specifically, we observed that αPD-1–treated cells uniquely upregulated multiple genes involved in IL-6 and mesenchymal stem-like pathways, such as Socs3, Lif, Mmp2, and Snai1 (Figure 1E and Supplemental Table 1). We confirmed these results using qPCR (Supplemental Figure 1C) and observed increased levels of IL-6 and LIF in CM in response to PD-1 blockade (Supplemental Figure 1, D–F). GSEA of Hallmark gene sets identified TNF-α/NF-κB and EMT signatures as among the top pathways expressed in these cells (Figure 1E, Supplemental Figure 1G, and Supplemental Tables 2 and 3). Evaluation of this ex vivo RNA-Seq pipeline using CT26, another colorectal cancer (CRC) model with relative αPD-1 resistance (19), again demonstrated EMT among the top pathways enriched after PD-1 blockade, as well as substantial overlap with multiple other MC38-enriched pathways (Supplemental Figure 2 and Supplemental Tables 4 and 5). These data suggested that tumors contain specific cellular subpopulations that persist despite effective PD-1 blockade by activating programs similar to those engaged by DTPCs. However, we recognized that this transcriptional program could be activated in residual tumor cells or in other cells within the TIME of MDOTS, and therefore pursued scRNA-Seq to assess this response at higher resolution.

scRNA-Seq identifies unique cancer cell clusters resistant to αPD-1 therapy. We first sought to validate these results in vivo by treating syngeneic mice with established MC38 tumors with αPD-1 therapy (10 mg/kg i.v.) or negative control (isotype IgG) and then performing scRNA-Seq. To focus on cancer cell subpopulations from residual MC38 tumors after treatment, we sorted CD45/CD90 double-negative cells (enriched for tumor cells) and compared them with CD45/CD90 double-positive immune cells (Figure 1F). We used the specific Hallmark pathway EMT genes that were upregulated by bulk RNA-Seq analysis in αPD-1–treated MC38 MDOTS to generate a more specific signature (E-M⁺) for this resistant cell state (Supplemental Table 6). Similar to our findings in MDOTS, we observed strong enrichment of this E-M⁺ signature uniquely in αPD-1–treated cancer cell subpopulations (Figure 1F). To isolate the nature of these cell subpopulations further, we took advantage of the fact that MC38 MDOTS contained a minority of CD90⁺ mesenchymal cells by day 6 and that digesting MDOTS from collagen substantially depleted CD45⁺ immune cells (Supplemental Figure 3, A and B). We therefore extracted collagen-embedded residual cells from MC38 MDOTS at the end of treatment with IgG or αPD-1 and performed scRNA-Seq from this ex vivo model (Supplemental Figure 3B). Unsupervised clustering of 2543 IgG- and 2626 αPD-1–treated cells confirmed robust and distinct αPD-1–specific transcriptional states, including 4 unique clusters, with the majority of αPD-1–treated cells appearing in clusters 1 and 2 (Figure 2A). We next applied GSEA of Hallmark pathways to these clusters, which revealed that the 2 dominant αPD-1–resistant clusters 1 and 2 both shared prominent downregulation of IFN-γ and IFN-α signatures (Figure 2, B and D, and Supplemental Figure 4A). In contrast, clusters 1 and 2 were distinguished by highly divergent changes in cell-cycle gene expression, with upregulation of E2F targets and G2M-associated genes in cluster 1, versus strong negative enrichment of these genes in cluster 2, suggestive of quiescence (Figure 2, C and D, and Supplemental Figure 4A). Clusters 1 and 2 shared activation of oxidative phosphorylation pathway genes and notably were also distinguished by upregulation of the EMT Hallmark signature in cluster 2 (Figure 2D). Taken together, these findings suggest that a specific subset of MC38 cells evaded αPD-1 therapy by adopting a quiescent, mesenchymal phenotype.

Figure 3

Snai1 and Ly6a (Sca-1) mark a subpopulation of αPD-1 immunotherapy persister cells. (A) Two-dimensional t-SNE plot showing projection of bulk EMT signature (E-M) onto scRNA-Seq clusters of MC38 MDOTS. (B) Volcano plot of E-M⁺ (z-based score > 0.25 for bulk EMT signature) versus E-M^– cells in scRNA-Seq clusters of MC38 MDOTS. (C) Left, representative mRNA in situ hybridization (ISH) images showing Snai1 expression in perinecrotic regions at scale of 50 μm in hematoxylin-stained tumor specimens from MC38 mice 1 week after treatment with either negative control (isotype control IgG + vehicle) or αPD-1 monotherapy (n = 3 each). Right, summary of quantification of H score of Snai1 in the nonnecrotic and perinecrotic regions of tumors. Data are box-and-whisker graphs with box representing IQR, solid line representing median, and all the points ranging from minimum to maximum, and were analyzed by 1-way ANOVA followed by Bonferroni’s correction. (D) t-SNE plot showing relative SNAI1 expression (log₂ expression scale) localized to carcinoma cells from single-cell transcriptome analysis of an intrinsically αPD-1–resistant microsatellite instability–high (MSI-H) colorectal cancer (CRC) patient tumor. Gray color represents cells in which Snai1 reads were not detected. (E) Normalized bulk RNA-Seq SNAI1 expression scores for this αPD-1–resistant MSI-H CRC patient tumor compared with those of all CRC patient samples (n = 557) or MSI-H CRC patient samples (n = 73) available in The Cancer Genome Atlas (TCGA). (F) Normalized bulk RNA-Seq SNAI1 expression scores on pretreatment melanoma tumor specimens from patients who achieved versus did not achieve durable clinical benefit (see Methods for definition) from treatment with nivolumab. Dotted red line represents top 10th percentile. Data are number of patients and were analyzed by Fisher exact test (n = 51, shown in Table 1).

Figure 2

Single-cell RNA-Seq identifies unique cell clusters resistant to α–PD-1 therapy. (A) Two-dimensional t-SNE plots of single-cell RNA-Seq (scRNA-Seq) performed on MC38 MDOTS after 6 days of treatment. t-SNE plots are colored based on treatment condition (left) and unsupervised clustering distribution (middle). Right, bar graphs showing the number of cells in each cluster from αPD-1–treated (gray) or isotype IgG (black) groups. (B and C) t-SNE plots showing relative expression of Hallmark IFN-γ response (B) and E2F target (C) signature on scRNA-Seq clusters of MC38 MDOTS treated with αPD-1 (clusters 1–4, 5) versus isotype IgG control (clusters 5, 6–10). (D) Cluster pathway analysis using the Hallmark gene sets of the genes differentially expressed for scRNA-Seq clusters 1 and 2 of MC38 MDOTS, with their normalized enrichment scores.

Snai1 and Ly6a (Sca-1) mark a subpopulation of αPD-1 immunotherapy persister cells. We next compared expression differences in E-M⁺ signature high versus low cells to identify specific markers of cells expressing this feature. In consonance with the Hallmark pathway analysis, cells within the quiescence-associated cluster 2 most strongly expressed the top genes from this E-M⁺ signature (Figure 3A). Additionally, as expected, this analysis identified increased expression of genes included in the signature itself, such as Sdc1, Mmp2, and Dcn1 (Figure 3B and Supplemental Table 7). This analysis also uncovered Snai1 and Ly6a (encoding stem cell antigen 1, Sca-1) as among the most significantly upregulated transcripts in these E-M⁺ cells, using a stringent statistical threshold of log₂ fold change of more than 0.5 and –log₁₀-adjusted P value greater than 50. These genes were selected for further evaluation because of their biological plausibility as immunotherapy persister cell (IPC) markers. As discussed above, Snai1 has been linked to a similar resistant hybrid E-M state that is a feature of tumor-initiating cells (20). Similarly, Sca-1 is also a well-established marker of hematopoietic, tumor initiating, and tissue stem cells, which can evade immune surveillance (21, 22). Together, these findings further suggested the potential existence of a specific Snai1⁺ Sca-1⁺ anti–PD-1 IPC population that escapes CD8 T cell–mediated killing after effective PD-1 blockade.

To further validate enrichment of this IPC signature in vivo after PD-1 blockade by an orthogonal method, we evaluated Snai1 expression in MC38 tumors using mRNA in situ hybridization 1 week after αPD-1 therapy or negative control (isotype IgG + vehicle) to match ex vivo profiling studies. Consistent with our scRNA-Seq profiling results, clusters of tumor cells expressing high levels of Snai1 mRNA signal were observed and significantly enriched in the perinecrotic areas of all 3 αPD-1–treated MC38 tumors compared with isotype IgG-treated control tumors (Figure 3C). Evaluation of nonnecrotic areas did not show increased Snai1 mRNA signal.

Similar to syngeneic murine models such as MC38, microsatellite instability–high (MSI-H) status in patient tumors is a robust predictive biomarker for response to anti–PD-1 therapy in the clinic (23). Recently, genomic and transcriptional profiling of an MSI-H CRC patient tumor that was resistant to αPD-1 therapy identified biallelic loss of β₂-microglobulin as a potential mechanism for resistance to immune activation (24). In consonance with our findings, scRNA-Seq analysis also identified the presence of SNAI1⁺ cancer cell subpopulations preexistent within this αPD-1–resistant patient tumor (Figure 3D). Moreover, standardized SNAI1 expression from bulk RNA-Seq of this patient’s pretreatment tumor compared with available samples of CRC patients in The Cancer Genome Atlas (TCGA) was in the top 4.9 percentile of all samples and 8.2 percentile of MSI-H CRC tumor samples (Figure 3E), consistent with the intrinsic resistance to PD-1 blockade compared with typical MSI-H patients and our findings from the MC38 model. We also examined SNAI1 expression in pretreatment melanoma tumors and found that high (top 10th percentile) SNAI1 expression significantly enriched for lack of durable clinical benefit from treatment with nivolumab (Figure 3F and Table 1). In contrast, pretreatment expression of genes related to the TNF-α/NF-κB signature such as baculoviral IAP repeat-containing protein 2 (Birc2) and Birc3 were not predictive (Supplemental Figure 4B). However, the paucity of publicly available, large-scale, paired, pre- and post-progression RNA-Seq data sets limits robust analysis of on-treatment enrichment for each of these markers.

Table 1

Association of clinical benefit from nivolumab with Snai1 expression in pretreatment melanoma tumor samples in the Riaz cohort

IPCs that resist CD8⁺ T cell killing preexist in syngeneic cancer cell models. In contrast to Snai1, Ly6a (Sca-1) is a mouse-specific gene, but its utility as a cell-surface marker for hematopoietic and tissue stem cells has been well documented (21, 22, 25). We therefore used this marker to investigate the preexistence of IPCs in syngeneic murine cancer cell lines and their resistance to T cell killing. First, we generated ovalbumin-expressing MC38 cells (MC38-ova), cocultured them with ova-specific cytotoxic OT-I CD8⁺ T lymphocytes, and measured whether Sca-1⁺ cells emerge under immune selective pressure in this model system (26). As expected, increasing effector/target (E/T) ratios of OT-I T cells to MC38-ova cancer cells led to decreased viability of target cancer cells (Figure 4A). Notably, we observed that MC38-ova cancer cells indeed contained a small proportion of Sca-1⁺ cells that were highly enriched among the surviving cells at higher E/T ratios, becoming the dominant population after immune selection pressure. These data are consistent with the persistence of this MC38 tumor cell subpopulation after ex vivo PD-1 blockade and confirms their relative resistance to cytotoxic T cell killing.

Figure 4

IPCs that resist CD8⁺ T cell killing preexist in murine syngeneic cancer cell models. (A) Representative graph of 2 independent experiments (each done in triplicates) of coculture assay of target MC38 cells expressing ovalbumin antigen (MC38-ova) (T) and effector OT-I CD8⁺ T cells (E), showing viability of MC38-ova and percentage of Sca-1⁺ cells at increasing E/T ratios. (B) Representative flow cytometry plots and summary of experiments evaluating percentage of preexisting Sca-1⁺ CD44⁺ cells (hereafter Sca-1⁺ cells) in MC38 and CT26 cells (n = 5). (C) Cluster dendrogram evaluating correlation between whole-exome sequencing results for the 2 independently sorted MC38 Sca-1⁺ purified and Sca-1^– samples. (D) Schema of the following experiments conducted with Sca-1⁺ purified fraction of syngeneic cancer cells isolated by FACS. (E) Representative flow cytometry plots and summary of experiments evaluating for persistence of Sca-1⁺ cells in culture of Sca-1⁺ purified fractions of MC38 and CT26 at 24 and 94 hours (n = 3). (F) Heatmap showing log₂ fold change (log₂F) of differentially produced cytokines in conditioned media (CM) of Sca-1⁺ purified cell culture of MC38 versus CT26 cells at 24 and 96 hours (n = 2, each was run in duplicates). Triangle symbol represents above the level of detection. (A, B, and E) Data are mean ± SEM and were analyzed by multiple t tests with Bonferroni’s correction (A) or 2-tailed Student’s t test (B and E). *P < 0.05, **P < 0.01, ***P < 0.001.

Next, to characterize this preexistent Sca-1⁺ cell subpopulation further, we used flow cytometry to quantify their presence in CRC (MC38 and CT26) and other syngeneic cancer cell models. Notably, both CRC lines were predominantly CD44⁺, and contained a small subpopulation of Sca-1⁺ cells averaging approximately 1% in MC38 and approximately 3% in CT26 (Figure 4B and Supplemental Figure 5, A–C). Ras-mutant lung cancer models LLC and CMT167 had substantially higher proportions of Sca-1⁺ CD44⁺ cells (Supplemental Figure 5, B and C). Whole-exome sequencing of MC38 Sca-1⁺ and Sca-1^– populations, isolated by flow sorting, confirmed that they shared known MC38 cancer gene variants, but uncovered broader genomic heterogeneity across independently sorted samples, consistent with stochastic acquisition of this cell state (Figure 4C and Supplemental Figure 5, D and E). Of note, there was no significant difference in surface MHC class I H-2Kb expression between Sca-1⁺ and Sca-1^– cells in culture (Supplemental Figure 6A), indicating that resistance in the coculture assay to cytotoxic T cells was not simply driven by impaired antigen presentation.

We next characterized the behavior of these Sca-1⁺ IPCs over time in culture. In the absence of immune cell pressure, purified MC38 Sca-1⁺ cells quickly reverted to Sca-1^– cells, returning close to baseline proportion within 96 hours, whereas Sca-1⁺ cells from the more αPD-1–resistant CT26 and LLC models persisted over time (Figure 4, D and E, and Supplemental Figure 6B). Furthermore, Sca-1 depleted fractions also acquired Sca-1⁺ cells at a proportion commensurate with their baseline positivity in bulk cultures (Supplemental Figure 6C). These data further reinforce the stochastic nature and plasticity of this cell state.

We wondered whether differences in the proportion of MC38 and CT26 Sca-1⁺ cells could be related to certain growth factors or cytokines that support their expansion. Both Hgf and Fgf7 were upregulated in bulk RNA-Seq data after αPD-1 treatment and have been implicated in resistance to oncogene-targeted therapies (27, 28). However, culture of Sca-1⁺–sorted MC38 cells in the presence of recombinant HGF or FGF7 could not rescue their depletion over time (Supplemental Figure 6D). We therefore explored differences in autocrine cytokine/chemokines in CM from sorted MC38 versus CT26 Sca-1⁺ cells, as a more unbiased screen to identify factors that could explain the relatively greater persistence of Sca-1⁺ cells in CT26. Notably, among 32 potential candidates, IL-6 was consistently the most upregulated cytokine in CT26 Sca-1⁺ CM versus MC38 CM (Figure 4F and Supplemental Figure 7A). IL-6 was also found at higher levels in the CM of Sca-1⁺ cells versus Sca-1^– cells for both MC38 and CT26 (Supplemental Figure 7, B and C).

Sca-1⁺ IPCs expand in response to IL-6 and show differential thresholds to TNF cytotoxicity. We next performed a focused in vitro screen using sorted MC38 Sca-1⁺ cells to identify which among the top upregulated cytokines and growth factors could significantly expand these cells over time. In addition to IL-6, HGF, and FGF7, we included LIF, which was secreted by αPD-1–treated MC38 MDOTS and is known to support stem cell growth (29), as well as G-CSF, GM-CSF, M-CSF, and VEGF, which were also higher in CT26 Sca-1⁺ CM (Figure 4F). IFN-γ treatment was used as a positive control inducer of Sca-1 surface expression (21, 30). Notably, only IL-6 was capable of increasing the proportion of MC38 Sca-1⁺ cells, similar to positive control IFN-γ treatment (Figure 5A and Supplemental Figure 8A). However, in contrast to IFN-γ, IL-6 failed to upregulate MHC class I expression in Sca-1⁺ cells (Figure 5B). We further confirmed that IL-6 expanded Sca-1⁺ cells in MC38 MDOTS (Figure 5C). Similar to MC38, IL-6 boosted Sca-1⁺ cells in vitro in other models, including CT26, LLC, and CMT167 (Figure 5D). IL-6 also promoted expansion of Sca-1⁺ cells in Sca-1–depleted MC38 and CT26 cells, confirming a robust role for IL-6/STAT3 signaling in expanding this cell population (Supplemental Figure 8B). Thus, IL-6 signaling supports this cell state, while limiting IFN response and tumor antigenicity.

Figure 5

Sca-1⁺ IPCs expand in response to IL-6 and show differential thresholds to TNF cytotoxicity. (A) Summary of experiments showing proportion of Sca-1⁺ cells in 96-hour cultures of Sca-1⁺ purified MC38 cells supplemented with growth factors (100 ng/mL) versus media (negative control) and IFN-γ (positive control) (n = 3). (B) Left, representative histogram of MFI of H-2Kb of MC38 Sca-1⁺ cells in 96-hour culture of Sca-1⁺ purified fractions (media alone vs. IFN-γ vs. IL-6 stimulation [100 ng/mL]). Right, summary of experiments [n = 3]. (C) Representative immunofluorescence images and summary of experiments showing Sca-1 expression in MC38 MDOTS after treatment with IL-6 versus control (IgG) (scale bar: 200 μm). Bottom, data are median (solid line) with first and third quartiles (dashed lines) and were analyzed by Wilcoxon rank sum test. (D) Representative flow cytometry plots showing effect of stimulation by IL-6 on expansion of Sca-1⁺ cells in 96-hour cultures of Sca-1⁺ purified fractions. Right, summary of experiments (n = 4 for MC38/ CT26, n = 3 for LLC/CMT167). (E) Summary of experiments showing effect of stimulation by IL-6, TNF-α, or IL-6 + TNF-α (100 ng/mL) on fold change of proportion of Sca-1⁺ cells in 96-hour cultures of Sca-1⁺ purified fractions of MC38 and CT26 cells versus media alone (negative control) (n = 3). (F) Heatmap showing log₂ fold change (log₂F) of expression of genes involved in apoptosis pathway in MC38 Sca-1⁺ cells stimulated with TNF-α (100 ng/mL) for 2 and 6 hours versus media alone. Asterisks denote genes upregulated on RNA-Seq analysis of MC38 MDOTS after anti–PD-1 therapy. (A, B, D, and E) Data are mean ± SEM and were analyzed by multiple t tests with Bonferroni’s correction (A and E), 1-way ANOVA with Bonferroni’s correction (B), or 2-tailed Student’s t test (D). *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001.

We next compared IL-6 stimulation of MC38 Sca-1⁺ cells with TNF-α given that TNF-α/NF-κB signaling was also among the top αPD-1–enriched signature in MC38 MDOTS, but TNF-α was secreted at comparatively low levels in the CM (Figure 1E and Supplemental Figure 9A). In contrast to IL-6, exogenous TNF-α exposure potently inhibited MC38 Sca-1⁺ cell growth in culture and abrogated the ability of IL-6 or IFN-γ to expand these cells, consistent with its known dual pro-and antiapoptotic functions (Figure 5E, Supplemental Figure 9B, and Supplemental Figure 10A) (31, 32). Indeed, TNF-α treatment of Sca-1⁺ cells induced expression of multiple proapoptotic genes, including Fas and Traf1, countered by survival genes identified in the αPD-1–induced TNF-α/NF-κB signature, including Birc2 and Birc3 (Figure 5F, Supplemental Table 2, and Supplemental Figure 9C).

Evaluation of other syngeneic models in culture further revealed differential sensitivity of Sca-1⁺ cells to TNF-α. Although LLC Sca-1⁺ cells were especially sensitive, CT26 Sca-1⁺ cells were only modestly sensitive and CMT167 Sca-1⁺ cells were comparatively resistant (Figure 5E, Supplemental Figure 9D, and Supplemental Figure 10, A and B). Furthermore, in contrast to MC38, IL-6 supplementation was able to rescue Sca-1⁺ cells at least partially in these 3 models. Finally, consistent with these in vitro culture results, MC38 MDOTS showed significant reduction in live cells after culture for 6 days in the presence of TNF-α in a concentration-dependent manner, while less marked (and nonsignificant) effects were seen in CT26 MDOTS (Supplemental Figure 9E). These results support the presence of differential tumor TNF cytotoxicity thresholds within these stem-like populations.

Birc2/3 degradation by LCL161 further sensitizes IPCs to TNF-α and promotes durable αPD-1 response. Since Birc2 and Birc3 were induced in Sca-1⁺ IPCs, and have been previously implicated in CD8⁺ T cell TNF cytotoxicity thresholds and αPD-1 sensitivity (33, 34), we explored whether inhibition of Birc2/3 function could exploit this TNF-α–mediated vulnerability of Sca-1⁺ IPCs. Indeed, addition of the second mitochondria-derived activator of caspase (SMAC) mimetic, LCL161, to TNF-α treatment completely eradicated MC38 Sca-1⁺ cells, irrespective of supplementation with IL-6 (Figure 6A and Supplemental Figure 10A). Treatment with LCL161 also reduced enrichment of MC38 Sca-1⁺ cells under immune selection pressure with coculture of MC38-ova and OT-1 CD8⁺ T lymphocytes (Supplemental Figure 11A). Moreover, growth of comparatively TNF-α–resistant CT26 Sca-1⁺ cells was also potently inhibited by cotreatment with LCL161, though still partially rescued by IL-6 (Figure 6A and Supplemental Figure 10A). These differential sensitivities and combinatorial activity of TNF-α and LCL161 were further recapitulated in MC38 and CT26 MDOTS (Supplemental Figure 11, B–D). Further evaluation in lung syngeneic cancer models also confirmed the potent activity of cotreatment with TNF-α and LCL161 in eradicating Sca-1⁺ IPCs in culture, especially in CMT167, which was resistant to TNF-α alone (Supplemental Figure 10B).

Figure 6

Birc2/3 degradation by LCL161 further sensitizes IPCs to TNF-α and promotes durable αPD-1 response. (A) Representative flow cytometry plots showing effects on proportion of Sca-1⁺ cells in culture of Sca-1⁺ purified MC38 and CT26 cells at 96 hours with or without treatment with Birc2/3 antagonist LCL161 in the presence of TNF-α, IL-6, TNF-α + IL-6 (both at 100 ng/mL), or media alone (n = 3). (B) Representative immunofluorescence microscopic images and quantification of CD45 and DAPI-stained tumor specimens from MC38 tumor-containing mice 1 week after treatment with negative controls (isotype control IgG + vehicle), αPD-1 monotherapy, LCL161 monotherapy, or αPD-1 + LCL161 combination therapy (n = 6 mice per group, except 4 for IgG + vehicle) (scale bar: 100 μm). Data are box-and-whisker graphs, with box representing IQR, solid line representing median, and all points ranging from minimum to maximum, and were analyzed by Kruskal-Wallis test. (C and D) Summary of in vivo experiments with Kaplan-Meier curves showing tumor volumes (mm³) over time (days) in MC38 (C) or CT26 (D) tumor-containing mice treated weekly with negative controls (vehicle + IgG), αPD-1, LCL161, or αPD-1 + LCL161 combination therapy. Box with shaded region represents duration of treatment. Mice with complete response (CR)/total number of mice tested are presented for each group. (C) For MC38, n = 11 mice, except 10 in LCL161 monotherapy group. (D) For CT26, n = 8 mice, except 7 in LCL161 monotherapy group. Some mice were euthanized because of tumor necrosis before they reached threshold for progression (2 mice in vehicle + IgG, 5 in αPD-1, 6 in LCL161, and 2 in combination therapy groups); 1 mouse in combination therapy group that achieved CR had accidental death on day 15. *P < 0.05, **P < 0.01, ***P < 0.001.

To determine the potential therapeutic relevance of these findings, we next evaluated the efficacy of αPD-1 ± LCL161 treatment in vivo. Syngeneic mice with established MC38 and CT26 tumors were treated with αPD-1 therapy (10 mg/kg i.v. weekly), LCL161 (100 mg/kg weekly orally in 2 split doses), combination therapy with αPD-1 and LCL161, or negative controls (isotype IgG + vehicle). Histopathologic evaluation 1 week after the first treatment dose suggested enhancement of αPD-1 immune-mediated MC38 tumor cell killing by cotreatment with LCL161, with more complete eradication of tumor cells (Figure 6B). Indeed, over time, LCL161 and αPD-1 combination therapy yielded significant improvement in survival and complete responses compared with each monotherapy alone or isotype control in mice bearing MC38 tumors (Figure 6C and Supplemental Figure 12, A and B). These results were confirmed in 2 independent sets of experiments, with complete responses observed more frequently in the combination therapy group. Furthermore, consistent with the differential IPC TNF cytotoxicity thresholds observed in vitro and ex vivo, numerically higher but not statistically significant rates of complete response and durable survival were seen with combination therapy in mice bearing CT26 tumors compared with monotherapy with either αPD-1 or LCL161 (Figure 6D and Supplemental Figure 12, C and D). Taken together, these data reveal that Birc2/3 antagonism can co-opt this TNF-mediated vulnerability of IPCs and improve durable responsiveness to PD-1 blockade.

Here, we used dynamic ex vivo scRNA-Seq of organotypic tumor spheroids and identified a specific transcriptional program engaged by cells that persisted despite effective PD-1 blockade. We further confirmed that these IPCs represented preexisting subpopulations of cancer cells with a stem-like phenotype that resisted αPD-1–mediated CD8⁺ T cell reinvigoration. The specific features of these cells are highly analogous to the DTPCs that emerge after oncogene-directed targeted therapies, sharing properties of quiescence, apoptosis resistance, and EMT with tissue and cancer stem cells (13–16). In addition, IPCs downregulate IFN response, which has been associated with resistance to immune checkpoint blockade in multiple studies (33, 35–37). Despite this ability to escape T cell killing, these cells remain vulnerable given their reliance on activation of TNF-associated NF-κB survival signaling. By enhancing the TNF-α prodeath response that these cells exhibit and lowering their TNF cytotoxicity thresholds, for example by interfering with the key survival factors Birc2/3, it is possible to kill these cells and improve the ability of PD-1 blockade to induce durable responsiveness in vivo.

IPCs are characterized by expression of Snai1 and Ly6a (Sca-1). Snai1 is a master transcription factor that is known to drive a hybrid epithelial-mesenchymal state (20). EMT induction in the MCF10 epithelial mammary cell line by stably overexpressing SNAI1 resulted in increased activation of NF-κB/MAPK signaling and induction of IL-6/IL-8 on IL-1β stimulation and was associated with chemoresistance (38). Snai1 has also been linked with resistance to CTLA-4 blockade in murine breast carcinoma models (10). Knockdown of SNAI1 by shRNA in high-grade serous ovarian cancer models was associated with decreased expression of Nanog and Lin28, increased let-7 expression, and decreased self-renewal capacity reflected in tumor burden in orthotopic xenografts (39). Moreover, the role of metabolic pathways in promoting lung cancer metastasis has also been ultimately linked to enhanced stability of SNAI1 mRNA via depletion of UDP-glucose by UDP-glucose 6-dehydrogenase on EGFR activation (40). Sca-1, on the other hand, has been directly linked to hematopoietic and tissue stem cells (21, 22, 25, 41). For example, in the lung, Sca-1 expression marks bronchioalveolar stem cells that upon transformation give rise to adenocarcinoma (22). Besides normal stem cells, Sca-1 and other Ly6 family members have been linked with cancer stem cells (41). Of note, the tumor expression of human homologs of Sca-1, Ly6K, and Ly6E has been associated with poor overall survival in breast cancer with a mechanistic link to TGF-β signaling (42). These data lend support to our findings identifying stem cell–like subpopulations of cancer cells that persist despite immune activation and may ultimately drive tumor progression and/or relapses. Similar to the DTPC concept, these primitive populations provide a reservoir that can escape effective T cell killing and contribute to anti–PD-1 immune evasion.

The discovery that NF-κB survival signaling pathways and IL-6 in particular are integrally involved in expansion and survival of IPCs is further reminiscent of activation of these pathways in DTPCs (6, 7, 43). The identification of IL-6 as an IPC growth factor is also consistent with the role of IL-6 signaling in promoting EMT (44, 45) and supporting cancer stem cell expansion (46, 47). IL-6 has been shown to enrich for a stem cell phenotype after treatment with platinum-based chemotherapy in ovarian cancer (48) and non–small cell lung cancer (49). Specific induction of IL-6 in DTPCs has been observed in multiple oncogene-addicted lung cancer models and is sensitive to transcriptional inhibitors such as THZ1 (6). Activation of NF-κB signaling after EGFR tyrosine kinase inhibitor (TKI) treatment in lung adenocarcinoma induces formation of an EGFR-TRAF2-RIP1-IKK complex (43) and is associated with ER stress and STING activation (7), which promote DTPC survival. Thus, activation of NF-κB and induction of IL-6 may play a more general role in fostering expansion of these stem-like cancer cell subpopulations that resist not only chemo- and targeted therapies, but also anti–PD-1 immune checkpoint blockade.

Antiapoptotic adaptation by increased expression of specific regulators of mitochondrial membrane permeability are often employed by cancer cells in response to therapeutic stress (50). For example, inhibition of Bcl2 or Bcl-XL is capable of overcoming these prosurvival signals and triggering apoptosis in EGFR TKI DTPCs (7, 51). In the case of IPCs, we have identified Birc2 and Birc3 as alternate antiapoptotic nodes that are activated downstream of TNF-induced NF-κB signaling, providing a unique pharmacological target by which to trigger their apoptosis. This is in agreement with genetic identification of Birc2 in 2 independent genome-wide CRISPR/Cas9 loss-of-function screens as a target that augments activity of immunotherapy (33, 34). Furthermore, these findings suggest variation in antiapoptotic defense mechanisms adopted by cancer persister cells to resist targeted therapies versus anti–PD-1 therapies, potentially due to preferential induction of Birc2/3 dependency by TNF-α in the microenvironment in response to immune-activating therapies.

There remains an unmet clinical need to develop predictive biomarkers to guide precision medicine efforts involving PD-1 blockade and antiapoptotic drugs. There are at least 3 ongoing clinical trials that are evaluating combinations of SMAC mimetics/IAP antagonists with anti–PD-1 (NCT02890069 in advanced solid malignancies and NCT03111992 in multiple myeloma) or anti–PD-L1 (NCT03270176 in advanced solid malignancies) therapies. By uncovering IPCs, our study provides potential predictive biomarkers such as Snai1 expression that could enhance the precision of these trials. Furthermore, differential TNF cytotoxicity thresholds noted among the syngeneic CRC and lung cancer models in our study may potentially explain the variability in response anticipated in these clinical trials. IPCs with relatively lower thresholds (exemplified by MC38 and LLC models) were exquisitely sensitive to IAP antagonism in presence of TNF-α, whereas those with higher TNF cytotoxicity thresholds (exemplified by CT26 and CMT167 models) showed comparative resistance. Given the role of IL-6/STAT3 activation in expanding IPCs, evidence of IL-6 signaling could also represent a biomarker in pretreatment and/or on-treatment tumor biopsies, which would enhance robustness of predicting tumor response or resistance in patients. Finally, ex vivo modeling of therapeutic response and resistance in patient-derived organotypic tumor spheroids may provide a direct avenue to study these dynamic changes further without the burden of repeated in vivo biopsies.

Discovery of IPCs and their unique TNF cytotoxicity thresholds and antiapoptotic dependencies uncovers multiple additional potential therapeutic opportunities. Although our data demonstrated the feasibility of targeting one of these defense mechanisms to enhance durable anti–PD-1 responses, higher-order combinations with PD-1 blockade may be necessary to completely eradicate IPCs in vivo, especially in those with high TNF cytotoxicity thresholds. More generally, this study highlights the power of dynamic high-resolution scRNA-Seq using functional models of immunotherapy to understand more fundamental mechanisms of immune evasion, which could also be applied to dissect resistance to additional modes of treatment or combination immunotherapies.

Cell lines and primary cultures. MC38 murine colon adenocarcinoma cells were provided by Gordon Freeman (Dana-Farber Cancer Institute, DFCI) in 2015 under a material transfer agreement (MTA) from Jeffrey Schlom of the National Cancer Institute (NCI, Bethesda, Maryland, USA). CT26, LLC, and B16 cells were purchased from ATCC in 2015; 4T1 cells were purchased in 2016. ID8 and CMT167 cells were purchased from MilliporeSigma in 2016 and 2017, respectively. MC38 cells for in vivo experiments performed at Novartis Institute for Biomedical Research (Novartis) were received from the NCI (Rockville, Maryland, USA) under MTA 38699-15. All cells were routinely tested for mycoplasma and found to be free of contamination.

Murine in vivo experiments. Thawed syngeneic cancer cells were cultured for 3 passages in DMEM (MC38) or RPMI (CT26) supplemented with 10% FBS (Gemini Bioproducts). A cell aliquot was stained with acridine orange/propidium iodide (Nexcelom), and cell counts were performed prior to implantation using the Nexcelom Cellometer K2 image cytometer. For the MDOTS studies, MC38 or CT26 cells (5 × 10⁵ cells/mouse in 100 μL), resuspended in FBS-free media, were s.c. injected in the upper right dorsal flank of 8–12-week-old female C57BL/6 albino mice (The Jackson Laboratory, RRID:IMSR_JAX:000058) or BALB/c mice (The Jackson Laboratory, RRID:IMSR_JAX:000651), respectively. After the first week, postinjection tumor volumes were measured 3–4 times per week. Tumors were collected 2 to 3 weeks after implantation before they reached 600 mm³ (MC38) or 350 mm³ (CT26), or for humane reasons per IACUC regulations.

For in vivo therapeutic studies at Novartis, MC38 cells were thawed and cultured for 1 week prior to implant in DMEM supplemented with 10% heat-inactivated FBS; 10 mM HEPES; 1 mM sodium pyruvate; and 1× NEAA, during which they were split approximately 2–3 times. Next, 1 × 10⁶ MC38 cells/mouse were s.c. implanted on the upper-right dorsal flank of 5–7-week-old female C57BL/6 mice (Charles River Laboratories, RRID:MGI:5658459). When tumor volume reached approximately 100 mm³, mice were randomized by a computer using a Novartis program to control (isotype + vehicle) or treatment groups (monotherapy with anti–mouse PD-1, LCL161, or combination therapy with anti–PD-1 and LCL161); 17 mice were randomized per group. The vehicle control (30% 0.1 N HCL, 70% sodium acetate, pH 4.63) and LCL161 (50 mg/kg) were administered orally once a week in 2 split doses (7 hours apart). Isotype control (mIgG1, 10 mg/kg; (Novartis, clone MOPC-21), and anti–mouse PD-1 (10 mg/kg; (Novartis, clone 1D2) were i.v. administered once a week. The group with the combination regimen received anti–PD-1 therapy and LCL161 concurrently. The treatment was continued for 4 consecutive weeks except for 6 mice in the control, anti–PD-1, and combination groups (7 mice in LCL161 group), which were euthanized on days 7–8 to study treatment-related tumor changes. For the CT26 therapeutic study performed at DFCI, 0.25 × 10⁶ CT26 cells/mouse were s.c. implanted on the right flank of 8-week-old female BALB/cJ mice (The Jackson Laboratory). On day 9 after implantation, animals were randomized into the 4 treatment groups with n = 7–8 mice/group, based on tumor volume using Studylog software as vehicle control (30% 0.1 N HCl, 70% sodium acetate, pH 4.63 + rat IgG2a; Bio X Cell, clone 2A3), LCL161 (50 mg/kg orally once a week in 2 split doses; Novartis), anti–mouse PD-1 (10 mg/kg i.v. once a week; Bio X Cell, clone RMP1-14) or the combination of LCL161 and anti–PD-1 for 3 weekly treatments.

In therapeutic studies, the tumor volume was measured 2–3 times per week (length × width) with a digital caliper. The tumor volume was determined using the formula ([l × w²] × 3.14159)/6. Body weight was recorded at the same time. Survival endpoint was defined as either tumor volume exceeding 1000 mm³ or if the tumors became necrotic and mice were humanely euthanized. Mice with no measurable tumors, defined as complete responders, were followed for 55 days (MC38 model) or 71 days (CT26 model).

Spheroid preparation and microfluidic culture. MDOTS preparation and microfluidic culture (AIM Biotech) were performed as described previously (19, 52). Briefly, fresh tumor specimens were first minced using sterile forceps and a scalpel, followed by resuspension in DMEM (MC38) or RPMI (CT26) media containing 100 U/mL collagenase type IV (Life Technologies). Samples were resuspended in fresh media and strained sequentially through 100 μM and 40 μM filters to obtain S2 (40–100 μM) spheroid fraction. Finally, the mixture of S2 with type 1 rat tail collagen (Corning) was cultured in 3D microfluidic culture devices in the presence of the media. Treatments were performed by diluting drugs in full media and adding 300 μL to each well of the device. The treatments consisted of αPD-1 (10 μg/mL; Bio X Cell, clone RMP1-14), IgG2a (10 μg/mL; Bio X Cell, clone 2A3), αCD8a (10 μg/mL; Bio X Cell, clone 53-6.7), αIFN-γ (10 μg/mL; Bio X Cell, clone R4-6A2), IL-6 (100 ng/mL; R&D Systems, 406-ML), TNF-α (10-100 ng/mL; R&D Systems, 410-MT), and/or LCL161 (500 nM; Novartis).

Live/dead staining and quantification. Viability staining was performed by loading microfluidic devices with Hoechst 33342 and propidium iodide (Invitrogen) diluted in full media to the final concentrations of 10 μg/mL and 1 μg/mL, respectively. The dyes were incubated in the microfluidic device for 30 minutes in a 37°C incubator with 5% CO₂. Images were captured on a Nikon Eclipse 80i fluorescence microscope equipped with ProScan III-controlled Z-stack and linear stage (Prior), and Zyla 5.5 sCMOS camera (Andor). Image capture and analysis were performed using NIS-Elements AR software package. The whole well area was acquired using stitching of 4 fields of view captured with 4× objective. Three planes of focus 50 μm apart were acquired and an extended depth-of-focus image was generated. Live/dead measurements were done by quantifying total Hoechst-stained and dead propidium iodide–stained cell areas.

Immunofluorescence staining and cytokine analysis of microfluidic devices. For immunofluorescence studies, MC38 MDOTS in triplicate wells at day 6 of the culture with either IgG, αPD-1, or IL-6 were washed directly in the devices with PBS and blocked with 1% BSA in PBS for 15 minutes at room temperature. Directly conjugated antibodies were anti–mouse CD45-Alexa Fluor 647 (BioLegend, clone 30-F11; shown in green for better visualization), anti–mouse CD90.2-PE (BioLegend, clone 30-H12), and anti–mouse Ly6a/e (Sca-1)-PE (BioLegend, clone D7). Antibodies were diluted 1:50 in 10 μg/mL solution of Hoechst 33342 in PBS and loaded into microfluidic devices for a 45-minute incubation at room temperature in the dark. Spheroids were washed twice with PBS with 0.1% Tween 20 followed by PBS. Images were captured and whole-well image acquisition was performed as described in Live/dead staining and quantification. For Sca-1 fluorescence quantification of tumor cells, we manually placed regions of interest (ROI) frames in the areas containing only CD45^– cells. In these areas, mean Sca-1 fluorescence intensity was measured in all ROIs. The total number of ROIs was greater than 1000 for each treatment condition. Cytokine analysis of CM from MDOTS at the end of the experiment was performed using a Bio-plex pro mouse cytokine 23-plex assay (Bio-Rad), MILLIPLEX MAP mouse cytokine/chemokine magnetic bead panel premixed 32-plex (MilliporeSigma), or V-PLEX mouse TNF-α kit (MSD).

Bulk RNA-Seq from microfluidic devices. After treatment in the microfluidic devices, cells were lysed directly from the devices, and RNA was isolated using the Agencourt RNAdvance tissue isolation kit (with DNAse). RNA libraries were prepared from 250 ng total RNA using the Illumina exome capture kit per the manufacturer’s instructions. RNA-Seq was performed per the standard protocols at the DFCI Molecular Biology Core Facilities (Illumina NextSeq 500). A sample-to-sample correlation heatmap was generated using VIPER (visualization Pipeline for RNA-Seq analysis) and displayed hierarchical clustering of Spearman rank correlations across samples. Differential expression analysis of bulk RNA-Seq files was then performed using the R package DESeq2 (53). Volcano plots were generated using R with a log₂ fold-change cutoff of 1.5 and an adjusted P value cutoff of 1 × 10^–10. Pathway analysis of bulk RNA-Seq was performed using the DESeq2 results and run through the R package ClusterProfiler, using the Hallmark GeneSets from MSigDB (Broad Institute) (54). Plots of running enrichment score were generated with R.

scRNA-Seq of tumors treated in microfluidic devices and in vivo. MC38 MDOTS were isolated from the microfluidic device on day 6 of treatment after a 15-minute collagenase treatment. The spheroids were then treated with trypsin in a 37°C incubator for 5 minutes to obtain single-cell suspensions. The viable tumor cells were isolated via FACS, washed, and assessed for viability. For residual MC38 tumors after in vivo treatment, an additional step with FACS for CD45/CD90 double-positive and double-negative populations was performed. The cells from MDOTS or residual in vivo MC38 tumors were loaded onto a 10x chromium instrument (10x Genomics) per the manufacturer’s instructions. ScRNA libraries were generated using the single cell 3′ reagent kit (10x Genomics) per the user guide. Quality control of the completed libraries was performed using a bioanalyzer high sensitivity DNA kit (Agilent) and then sequenced using the Illumina NextSeq 500 platform.

The raw sequencing reads were processed using the 10x Genomics CellRanger bioinformatics pipeline v3.0.2. The assembled matrix was then fed into the standard workflow of the R package, Seurat v2.1.0. Genes that were expressed in at least 3 cells, and only cells that expressed at least 200 genes, were kept for downstream processing. Additionally, cells expressing more than 7000 genes and cells with more than 18% of UMIs mapping to mitochondrial genes were removed from the analysis. All the samples were prepared and sequenced together on the same platform.

The filtered matrix was log-normalized using global scaling in Seurat. UMI and mitochondrial transcript content were used as regression parameters. The normalized matrix was scaled and centered gene-wise, and then underwent dimensionality reduction using principal component analysis (PCA) on the highly varying genes. After visual inspection of the PCA elbow plot, the top 18 PCs were chosen for further analysis. Clustering was performed on the chosen PCs using the shared nearest neighbor algorithm in Seurat with default parameters. A t-distributed stochastic neighbor embedding (t-SNE) map was computed and plotted using the RunTSNE and TSNEPlot modules of Seurat. Cluster differential expression analysis was performed in Seurat using the FindMarkers command using the Wilcoxon rank sum test without thresholds. To further characterize these clusters, pathway analysis was carried out using the ClusterProfiler R package. FDR correction was done using Benjamini-Hochberg method. To estimate pathway scores for single cells, the z normalized expression values for pathway genes in each cell were summed. Scores using this method correlated highly with those generated using the AUCell R package (data not shown). Pathway scores were then added to the MetaData slot of Seurat objects (55).

E-M⁺ signature was generated using specific Hallmark pathway EMT genes that were upregulated by bulk RNA-Seq analysis in αPD-1–treated MC38 MDOTS. We used the statistical threshold of log₂ fold change of more than 0.5 or less than –0.5 and –log₁₀-adjusted P value greater than 50 for evaluation of differentially enriched genes in E-M⁺ cells. Candidate genes for further evaluation were selected based on biological plausibility and pretest probability.

Flow cytometry and FACS. Syngeneic cancer cells in culture were examined for surface expression of Sca-1–positive cells by flow cytometry with BD Biosciences LSRFortessa and the following antibodies: Ly6a/e (Sca-1)-FITC (3:100; Thermo Fisher Scientific, clone D7) and CD44-APC (3:100; BioLegend, clone IM7). MHC class I expression was measured by evaluating surface expression of H-2Kb molecule using the pacific blue–labeled antibody (2:100; BioLegend, clone AF6-88.5). CD45 expression was determined on pre-FACS and post-FACS tumor specimens by AF488 (1:100; BioLegend, clone 30-F11) and APC-Cy7 (1:100; BioLegend, clone 30-F11) antibodies, respectively. FACS was performed with BD FACSMelody device to isolate Sca-1⁺ purified and Sca-1 depleted (Sca-1^–) fractions by labeling MC38, CT26, LLC, and CMT167 cells with the above-mentioned Sca-1 antibody and CD44-PE-Cy7 (3:100; BioLegend, clone IM7). FlowJo v10 was used to perform analysis of flow cytometry raw data. Gating was based on forward- and side-scatter areas; positive and negative staining cell populations were determined using negative controls (rat IgG2a kappa-FITC, Thermo Fisher Scientific, clone eBR2a; rat IgG2b kappa-APC, BioLegend, clone RTK4530; rat IgG2b-PE-Cy7, BioLegend, clone RTK4530).

Coculture assay. The MC38 cell line was transduced to stably express OVA antigen from pLVX-lucOS-IRES-neo lentiviral vector (26). CD8⁺ T cells were isolated from 8- to 12-week-old C57BL/6-Tg (TcraTcrb)1100Mjb/J OT-I mice (The Jackson Laboratory, stock 003831) using magnetic separation and LS columns per the manufacturer’s protocol (Miltenyi Biotec, 130-049-401). OT-I T cells were activated with Dynabeads mouse T-activator CD3/CD28 beads (Life Technologies) for 24 hours before coculture. MC38-ova cells (target) were seeded in 24-well plates at a density of 50,000 cells per well in 1 mL of media and cocultured with OT-I CD8⁺ T cells (effector) at effector to target (E/T) ratios of 0, 0.5, 1, and 2 for 48 hours. In some experiments, LCL161 (500 nM) or DMSO control were added concurrently with OT-I T cells for 48 hours. Flow cytometry experiments were performed as described above. Viability was assessed with Live/Dead Fixable Zombie NIR (BioLegend). Dead cells were gated out before Sca-1 assessment to eliminate the possibility of nonspecific staining.

Bulk and single-cell sequencing for patient sample. The details of the MSI-H CRC patient’s case and RNA-Seq studies have been described previously (24). The comparison with other CRC patients (with both MSI-H and microsatellite-stable cancers) was performed using data from previously published data and TCGA program (56).

Human αPD-1 RNA-Seq data set. Bulk RNA-Seq data from pretreatment samples were obtained from the Riaz cohort (57). Patients were grouped into durable clinical benefit (DCB) and nondurable clinical benefit (NCB) response groups using a combination of response evaluation criteria in solid tumors (RECIST) and survival criteria. DCB patients included those who experienced complete response, partial response, or stable disease with an overall survival greater than 1 year. NCB patients included those who experienced progressive disease, stable disease with an overall survival less than 1 year, or patients with death prior to disease assessment. Gene expression data represent upper quartile normalized TPM counts. Fisher exact test was performed to evaluate enrichment in the top 10th percentile of expressers for analyzed genes of interest.

In vitro syngeneic cell cultures. Sca-1⁺ purified and Sca-1 depleted fractions were loaded at the same numbers and maintained in parallel in vitro cultures at 37°C in media in 6-well plates for 24 and 96 hours. Flow cytometry and cytokine/chemokine evaluation of CM were performed at the end of respective culture time points. Mouse cytokine/chemokine magnetic bead panel (MilliporeSigma, MCYTMAG-70K-PX32) was used to perform Luminex analysis of CM. The role of cytokine/chemokine supplementation on expansion of Sca-1⁺ population was evaluated by addition of the following cytokines at the final concentration of 100 ng/mL individually to parallel in vitro cultures on the day they were established using the same number of cells under identical conditions (except test cytokine): mouse IFN-γ (485-MI/CF), IL-6 (406-ML), LIF (8878-LF), VEGF (493-MV), G-CSF (414-CS), M-CSF (416-ML), GM-CSF (415-ML), FGF-7 (5028-KG), HGF (2207-HG) and TNF-α (410-MT); all R&D Systems). In additional experiments, LCL161 drug was added at the final concentration of 500 nM to MC38, CT26, LLC, and CMT167 cultures in the presence of TNF-α, IL-6, or both. Flow cytometry was performed after 96 hours of culture using Sca-1, CD44, and H-2Kb targeted fluorescently labeled antibodies as described above.

Whole-exome sequencing. DNA was isolated from MC38 Sca-1⁺ purified cells and Sca-1^– cells immediately after FACS, using DNeasy blood and tissue kit (QIAGEN) per the manufacturer’s instructions. Whole-exome sequencing was performed by GeneWiz (Illumina HiSeq 2X150PE). On average, 39 million 150 bp paired-end reads were generated for each sample, which were first checked for quality control, and then the adapters were trimmed. Reads were first mapped to the mouse genome (MGSCv37). The genome analysis toolkit (GATK) Best Practices workflow was used for variant calling followed by hard filtering on raw variants. Filtered SNP variants detected were annotated using SnpEff v.4.3i. Additional variant filtering was performed to remove known common variants. This was performed using the SnipSift v.4.3i package against dbSNP (v150), ExAC, and dbNSFP_3.5c databases. For 38,549 SNPs, the coverage of reads calculated for each SNP and those SNPs with coverage smaller than 10× were filtered out. For the remaining SNPs, additional filtering was performed for removing putative germline variants by comparing with C57BL_6NJ_R SNPs provided by Wellcome Sanger Institute (58). As a result, 2260 SNPs were selected for calculating the correlation between samples using BAMixChecker (59). Additionally, the previously characterized 7 exclusive somatic variations of the MC38 cell line were examined in all the samples (60).

Quantitative PCR with reverse transcription. MC38 Sca-1⁺ purified cells were cultured for 6 hours and subsequently stimulated with mouse TNF-α (100 ng/mL) at 37°C for 2 or 6 hours in parallel. The cells were washed with PBS and then extracted from 6-well plates using RLT Buffer supplemented with β-mercaptoethanol (QIAGEN RNeasy Mini Kit) and cell lifter (Corning). RNA isolation was performed per the manufacturer’s instructions. RT² Profiler PCR array mouse apoptosis kit (QIAGEN) was used to profile differential expression of 84 genes simultaneously in MC38 Sca-1⁺ cells on stimulation with TNF-α compared with no stimulation (control media). RT-qPCR was performed using Applied Biosystems StepOne Plus Real-Time PCR system, and data were analyzed in GeneGlobe Data Analysis Center. The results for Birc2 (forward primer, Fwd:5′-AGTAGATTTGCACATTCGTCACC-3′; reverse primer, Rev:5′-AGAATTAAGAGGGCTAGAGCACA-3′) and Birc3 (Fwd:5′-TGAAGAGTGCTGACACCTTTG-3′; Rev:5′-GGAAAAGCTGAATACGTGGACAA-3′) gene expression were validated subsequently by repeating the experiments twice. Additionally, the results of bulk RNA-Seq from MC38 MDOTS were validated by evaluating expression of Snai1 (Fwd:5′-CACACGCTGCCTTGTGTCT-3′; Rev:5′-GGTCAGCAAAAGCACGGTT-3′), Socs3 (Fwd:5′-ATGGTCACCCACAGCAAGTTT-3′; Rev:5′-TCCAGTAGAATCCGCTCTCCT-3′), Mmp2 (Fwd:5′-CAAGTTCCCCGGCGATGTC-3′; Rev:5′-TTCTGGTCAAGGTCACCTGTC-3′) and Lif (Fwd:5′-ATTGTGCCCTTACTGCTGCTG-3′; Rev:5′-GCCAGTTGATTCTTGATCTGGT-3′) genes among αPD-1– versus IgG-treated MC38 MDOTS by RT-qPCR analysis, as described before (19). The relative level of gene expression was normalized using the level of mouse 36B4 (Fwd:5′-AGATTCGGGATATGCTGTTG-3′; Rev:5′- CGGGTCCTAGACCAGTGTTC-3′).

Immunofluorescence staining of in vivo tumors. First, 5 μM FFPE MC38 tumor slides were deparaffinized and placed in a pressure cooker in antigen retrieval citrate (pH 6.0) buffer for 30 minutes at 115°C, followed by 15 minutes at 90°C. Slides were cooled for 30 minutes at room temperature and were then placed in Shandon Sequenza (Thermo Fisher Scientific) slide rack for immunostaining. Slides were first blocked with 1% BSA for 30 minutes. Anti–mouse Ly-6A/E (1:50; BioLegend, clone D7) was added to the cover plates and incubated overnight at 4°C. The slides were then washed 3 times with PBS with 0.1% Tween 20, incubated with Alexa Fluor 568 goat anti–rat IgG (1:100; Invitrogen, A-11077) for 30 minutes in the dark, and washed again 3 times as described above. Directly conjugated anti-mouse CD45-Alexa Fluor 647 (1:50; BioLegend, clone 30-F11) was subsequently added to the cover plates. After a 1-hour incubation in the dark, the slides were washed 4 times as described above, removed from the cover plates, rinsed in double-distilled H₂O, and dried before Vectashield mounting medium with DAPI (Vector Laboratories) was applied. Images were captured on a Nikon Eclipse 80i fluorescence microscope as described above; 5–6 areas (2560 × 2160 px) were captured per slide using a ×40 objective. Percentage of tumor and immune cells were evaluated based on CD45 staining and nuclear morphology by a pathologist who was blinded to the treatment groups.

mRNA in situ hybridization of in vivo tumors. mRNA in situ hybridization was performed on freshly cut 5 μm FFPE tumor sections using RNAscope 2.5 HD duplex assay (Advanced Cell Diagnostics) with Snai1 probe (red) and control Huwe1 probe (green). In situ hybridization scores were generated from perinecrotic and nonnecrotic regions by a board-certified pathologist at ×400 magnification (2–5 fields of view per tumor). RNAscope analysis was performed as a modified H score with the following signal intensity grades: 0 = no staining or less than 1 dot/10 cells; 1 = 1–3 dots per tumor cell; 2 = 4–9 dots per tumor cell, none or very few dot clusters; 3 = more than 10 dots per cell; and 4 = presence of dot clusters. Snail H score was calculated by multiplying the percentage of Snail-expressing tumor cells by signal intensity grade (61).

Data deposition statement. Bulk RNA-Seq and scRNA-Seq data sets have been deposited in NCBI’s Gene Expression Omnibus (GSE160228 and GSE160400).

Statistics. All graphs and scatter plots depict mean ± SEM unless otherwise indicated. Graphs were generated and statistical analysis was performed using GraphPad Prism (v8.0) or Microsoft Excel. Statistical tests have been mentioned in figure legends where applicable. P value was considered significant at less than 0.05.

Study approval. Informed written consent to participate in DFCI/Harvard Cancer Center IRB-approved research protocols was obtained from the human subjects (24). The study was conducted per World Medical Association Declaration of Helsinki and IRB-approved protocols. All animal experiments at DFCI were performed in compliance with established ethical regulations in Association for Assessment and Accreditation of Laboratory Animal Care International–accredited vivarium and were approved by the IACUC. In vivo therapeutic studies at Novartis were conducted in accordance with Novartis IACUC regulations and guidelines.