Longitudinal Immune Profiling Reveals Unique Myeloid and T-cell Phenotypes Associated with Spontaneous Immunoediting in a Prostate Tumor Model

The theory of cancer immunoediting describes how evolutionary pressures exerted by host immune responses shape the immunogenicity of emergent tumors (1–3). In turn, prosurvival adaptations by tumors shape the composition and character of host immune cells in a complex back-and-forth that evolves in the host over space and time. The outcome of this interplay may ultimately determine whether an evolving neoplasm becomes clinically apparent.

Three distinct phases of cancer immunoediting have been described: elimination, equilibrium, and escape (3). Central to the elimination phase are CD103⁺ dendritic cells (DC) driven by type I interferons to cross-present tumor antigens (4–6), polyfunctional tumor-specific CD8⁺ T cells elaborating IFNγ; refs. 7, 8), T_H1-polarized CD4⁺ T cells (9, 10), and innate inflammatory processes (11). These mechanisms have in part been defined using models driven to express model antigens for immunogenicity, or in response to immunomodulators that invigorate antitumor immunity (1). By contrast, the escape phase is characterized by T-cell dysfunction manifested in numeric deficiency and/or functional exhaustion (8, 12), insufficient antigen presentation due to DC skewing or exclusion, increased regulatory T-cell (Treg) function (13–15), an increased frequency of suppressive macrophages and myeloid-derived suppressor cells (MDSC; ref. 16), metabolic perturbations in nutrient deprivation and/or hypoxia (17, 18), and increased local expression of coinhibitory checkpoint ligands (19). The equilibrium phase, however, has proven much more difficult to study. The most informative model to date has been the 3′methylcholanthrene carcinogenesis model, which generates palpable stable masses that take upward of 200 days to form progressive tumors (20). This model shows that equilibrium is based on adaptive immunity to restrain tumor growth. However, significant gaps in the knowledge of what occurs at this phase remain, and models to study the equilibrium phase are therefore needed.

To address this knowledge gap, we studied a transplantable murine model of prostate cancer that we found to spontaneously and consistently recapitulate the three phases of immunoediting. This model, NPK-C1, was originally derived from a lung metastasis of NPK mice (Nkx3.1^CreERT2/+; Pten^flox/flox; Kras^LSL-G12D/+; R26R-^{LSL-\underline {\curr YFP}/+}) that generate bone-metastatic, castrate-resistant prostate tumors (21, 22). To define immune subsets with critical roles throughout the entire course of immunoediting, we utilized longitudinal 30-parameter spectral flow cytometry for phenotyping and performed in vivo depletion studies at early and late phases for functional validation. We defined subsets enriched in tumors that transition from the equilibrium to the escape phase versus those that do not. Finally, we explored unique cellular phenotypes (clusters) consistently associated with functional immunity at both the early and late phases of immunoediting. Together, these studies introduce a unique tumor model and provide a detailed view into the complex dynamics of T-cell and myeloid subpopulations over the course of immune editing in vivo.

Male C57BL/6J, B6 Albino, and J/NU mice (5–6 weeks old) were purchased from The Jackson Laboratory. Mice were 6 to 8 weeks old at time of use. All animals were housed in strict accordance with NIH and American Association of Laboratory Animal Care regulations. All experiments and procedures for this study were approved by the Columbia University Medical Center Institutional Animal Care and Use Committee (IACUC).

The NPK-C1 cell line (originally LM7304) was provided by Dr. Cory Abate-Shen at Columbia University. NPK-C1 cells were originally derived in 2011 to 2013 from a lung metastasis of NPK mice (for Nkx3.1^CreERT2/+; Pten^flox/flox; Kras^LSL-G12D/+; R26R-^{LSL-\underline {\curr YFP}/+}) that utilize an inducible Cre under control of the Nkx3.1 homeobox gene promoter for temporal and spatial regulation of gene recombination in luminal prostate epithelial cells, leading to highly bone-metastatic, castrate-resistant prostate tumors (21). Prior to in vivo use, the cell line was validated for contamination using the Charles River CLEAR Panel through the Columbia University Institute for Comparative Medicine. NPK-C1 cells were maintained in complete RPMI medium (Corning); supplemented with 10% fetal bovine serum (HyClone), penicillin (100 U/mL), and streptomycin sulfate (100 mg/mL; Gibco). NPK-C1 cells at passage number 10 were thawed for all in vivo experiments, with implantations occurring within two to three passages after thawing. NPK-C1 cells were engineered to express luciferase (termed NPK-C1-Luc cells), by lentiviral transduction of a modified pHAGE PGK-GFP-IRES-LUC-W plasmid (Addgene # 46793), in which blue fluorescent-protein (BFP) was inserted in place of GFP by Gibson cloning (New England Biolabs, E2611) in order to allow for flow sorting and selection of BFP-IRES-Luciferase cells. Lentiviral particles were generated in HEK-293 cells (ATCC), using second-generation packaging vectors (psPAX2 and pMD2.G; Addgene).

NPK-C1 cells at 70% to 90% confluence were harvested with 0.05% trypsin (Gibco), washed with PBS, counted, and resuspended at 10 × 10⁶ cells/mL in ice-cold PBS. On day 0, 6- to 8-week-old mice were implanted on the right flank with 1 × 10⁶ NPK-C1 cells. Beginning on day 6 after implantation, tumor measurements were recorded every 2 to 3 days by digital caliper. For in vivo depletion experiments, mice received 250 μg/mL anti-CD4 (GK1.5 IgG2b) or anti-CD8 (2.43 IgG2b), alternating 50 μg/mL anti-Ly6G (1A8 IgG2a) and anti-rat IgGκ (MAR18.5 IgG2a), or isotype antibody (LTF-2 IgG2b) by intraperitoneal injection beginning either day 5 or 17 after implantation. All antibodies were purchased from Bio X Cell. Please see respective figures for details on timing of implantation and harvest. For reimplantation studies, NPK-C1 or -C2 tumors were harvested and disaggregated as described for flow cytometry preparation. Tumor cell suspensions were cultured for 2 weeks in standard RPMI medium as above, and then were cryopreserved. In vivo–passaged cell lines were thawed and implanted in mice as described above. For cumulative growth curves, the average tumor volume ± SD for each group was plotted until the first mouse in that treatment group was euthanized. Mice were sacrificed upon reaching either of the following criteria: (i) tumor volume equal to or greater than 1,000 mm³ or (ii) development of tumor ulceration >5 mm in diameter.

For orthotopic implantation of NPK-C1 and NPK-C1-Luc cells, J/NU (nude), C57BL/6, or B6-Albino mice were anesthetized by intraperitoneal injection of ketamine (100 mg/kg) and xylazine (16 mg/kg; Patterson Vet Supply), and an incision was made in the lower abdomen under IACUC-approved protocols. Tumor cells were resuspended in PBS and 10⁵ cells were injected in 30 μL into the left anterior prostate using a 30-guage ½ inch needle on a Hamilton syringe. The incision site in the muscle layer was closed with surgical suture (Ethicon J385H Vicryl absorbable suture, size 5–0) and the skin incision by wound clips. The orthotopic growth of NPK-C1-Luc was monitored twice weekly by bioluminescence imaging using an IVIS Spectrum Optical Imaging System (PerkinElmer). Ten minutes prior to imaging, mice were injected intraperitoneally with 150 mg/kg D-luciferin (PerkinElmer). Images were generated and quantified using Living Image Software (PerkinElmer).

Following euthanasia, NPK-C1 tumors were measured with a digital caliper, harvested, weighed, and diced into 70-μm filters (Thermo Fisher) in 6-well plates containing 5 mL of X-VIVO 15 media (Lonza) supplemented with 1 mg/mL collagenase H (Sigma) and 160 μg/mL DNAse (Roche). Tumors were digested in a 37^o incubator for 40 minutes. Digested tumors were physically passed through the 70-μm filter using syringe plungers and washed with complete RPMI media to obtain single-cell suspensions. Cells were counted, washed with PBS and then stained with the Zombie NIR Fixable Viability Kit at 1:1,000 in PBS (BioLegend) for 15 minutes. Samples were washed with FACS buffer (consisting of PBS, 3% FBS, 1% pen/strep, 1 mmol/L EDTA, and 10 mmol/L HEPES; Gibco), blocked with 1:50 TruStain Mouse FcX Reagent (BioLegend) in FACS buffer on ice for 15 minutes, and then stained with surface antibodies in FACS buffer (150 μL total volume) supplemented with 50 μL Brilliant Stain Buffer (BD) for 45 minutes on ice in the dark. After washing with FACS buffer, cells were fixed using the FoxP3 Fix/Perm Kit (eBioscience) for 30 to 60 minutes on ice in the dark. Following fixation, samples were stained with intracellular antibodies in 1× FoxP3 kit perm wash for 45 minutes on ice in the dark. Following washing with 1× perm wash and FACS buffer, samples were filtered through 70-μm blue cap tubes (Corning) and acquired on a 3-laser (RBV) Cytek Aurora spectral flow cytometer (Cytek). All antibodies were purchased from BioLegend, BD, or Thermo Fisher as described in Supplementary Table S1. Ultracomp eBeads and ArC Amine Reactive beads (Thermo Fisher) were used for single reference controls. All downstream data analysis was performed using FlowJo v10.6 (Treestar). Sample gating trees with cell population–defining markers used for conventional flow cytometry data analysis are described in Supplementary Data.

NPK-C1 tumors were harvested and dissociated to single-cell suspensions as described above. 5–10 × 10⁶ total cells were resuspended in 50 μL X-VIVO15 medium (Lonza) and plated in a 96-well round bottom plate with 50 μL of X-VIVO15 medium containing phorbol myristate acetate (PMA) and ionomycin (InvivoGen) to final concentrations of 50 ng/mL PMA and 500 ng/mL ionomycin. After 30 minutes, 100 μL of X-VIVO15 medium containing PMA/ionomycin and Brefeldin A at a final concentration of 5 μg/mL (BioLegend) were added. After 3.5 hours (4 hours total stimulation), cells were washed and fixed with the FoxP3 Fix/Perm Kit (eBioscience) for 30 to 60 minutes on ice in the dark. Cells were then stained for surface and intracellular markers as described and acquired on a Cytek Aurora spectral flow cytometer (Cytek). All antibodies used are listed in Supplementary Table S1.

For longitudinal immune profiling, two independent replicate experiments were performed; each included n = 5–15 mice per group; day 10 = 15 mice, day 15 = 15 mice, day 20 = 15 mice, and day 24 = 25 mice. For in vivo depletion studies, two independent replicate experiments were performed consisting of 5 to 10 mice per group. Student t test or two-way ANOVA with Tukey correction was used for individual frequency comparisons, Holm–Sidak correction was used for multiple comparisons of tumor growth curves, and standard linear regression was used for all correlations. To control for multiple linear regression comparisons, Bonferroni and Benjamini–Hochberg corrections were used as shown in each respective figure. These statistics were run using GraphPad Prism version 9.

For high-dimensional visualization and clustering of flow cytometry data, t stochastic neighbor embedding (tSNE) was performed within FlowJo, with 3,000 iterations and k = 30 for all projections. Clustering was performed on all CD45⁺ cells using the FlowSOM implementation within FlowJo, and selected metaclusters were isolated for successive tSNE projections and FlowSOM clustering as described in respective figures. For longitudinal cluster frequency comparisons, statistical significance was calculated using two-way ANOVA with Tukey correction for multiple comparisons. To generate phenotyping heat maps, each cluster’s median fluorescence intensity (MFI) for all phenotyping markers was exported from FlowJo. Expression normalization was used to scale each MFI to the minimum and maximum cluster MFI for each marker in a given heat map. The following equation was used: Normalized MFI = [(cluster MFI) − MIN(all cluster MFIs)]/[MAX(all cluster MFIs) − MIN(all cluster MFIs)]. For markers with universally high or low MFIs within a given compartment (for CD8: TCR-β, CD8, CD4, KLRG1, CX3CR1, CD86; for CD4: TCR-β, CD8, CD4, CD62L), scaling was performed against the minimum and maximum of all immune clusters. For bootstrapping, 10 mice in a group were sampled at random, and the average frequency of a given cluster was compared with that from 10 random mice chosen from the comparison group (e.g., day 10 group vs. day 15 group). This procedure was repeated 100 times without replacement to generate average fold-change values with an approximated deviation. Bootstrapping and differential fluorescence analysis was performed in R version 3.6.1.

The NPK-C1 cell line was derived from a lung metastasis of an NPK mouse that utilizes an inducible Cre expressed under control of the Nkx3.1 homeobox gene promoter for temporal and spatial regulation of gene recombination in luminal prostate epithelial cells (21). As a genetically engineered mouse model, the NPK model generates castrate-resistant prostate tumors that recapitulate the clinical tendency of prostate cancer to metastasize to bone (22). To study the growth dynamics of NPK-C1 tumors in immune-competent mice, we implanted NPK-C1 cells subcutaneously into C57BL/6J mice, as well as into T cell–deficient athymic J/NU (nude) mice. We observed a clear difference in tumor growth dynamics between strains. Whereas NPK-C1 grew progressively in immune-deficient J/NU mice, tumors spontaneously regressed around day 10 in C57BL/6 mice, consistent with immune-mediated regression (23). However, these tumors were not immediately eliminated; they remained at a size of ∼100 to 200 mm³ for at least 10 days, at which point roughly 60% of tumors progressed rapidly, around 30% of tumors remained in prolonged equilibrium but ultimately escaped, and roughly 10% of tumors ultimately underwent rejection (Fig. 1A; Supplementary Fig. S1A and S1B). These growth kinetics were also observed when NPK-C1 was orthotopically implanted into the anterior prostate (Supplementary Fig. S1C and S1D). Reimplantation of escaping NPK-C1 tumors into naïve mice yields progressively growing tumors, further suggesting NPK-C1 tumors undergo immune editing in vivo (Supplementary Fig. S1E). Taken together, these data support the notion that the in vivo growth kinetics of NPK-C1 tumors in immunocompetent mice are consistent with the canonical “three E’s” of cancer immunoediting: elimination, equilibrium, and escape (Fig. 1A; ref. 3).

To quantify the immune composition of NPK-C1 tumors across the immunoediting process, we harvested subcutaneously implanted NPK-C1 tumors for flow cytometry at days 10, 15, 20, and 24, representing tumors at elimination, equilibrium, early and late escape phases, respectively (Fig. 1B). To facilitate these studies, we developed a 28-color antibody panel designed for tumor–immune phenotyping using the Cytek Aurora spectral flow cytometer (Supplementary Table S1). With this panel, we defined key canonical immune populations of interest in NPK-C1 tumors using standard biaxial gating (Supplementary Fig. S2A). On day 10 (elimination phase), NPK-C1 tumors contained roughly equal frequencies of CD8⁺ T cells, CD4⁺ conventional T (Tconv) cells, and natural killer (NK) cells (∼10% of CD45 each; Fig. 1C). The myeloid compartment was dominated by granulocytes (∼40% of CD45), followed by monocytes (∼15% of CD45) and dendritic cells (DC; ∼8% of CD45). As NPK-C1 tumors transitioned to equilibrium phase at day 15, we observed an expansion of the granulocyte compartment (to ∼55% of CD45; P < 0.0001), with corresponding contractions in the CD4⁺ Tconv compartment (to ∼2% of CD45; P = 0.005) and monocyte compartment (to 10% of CD45; P = 0.0073). CD8⁺ T-cell frequencies did not significantly change at this stage. On days 20 and 24, as NPK-C1 tumors progressed to the escape phase, the granulocyte compartment retracted to 20% of CD45 (P < 0.0001). At this later point, a significant CD8⁺ T-cell expansion (up to ∼25% of CD45, P < 0.0001), together with moderate increases in NK-cell frequencies (to ∼15% of CD45; P < 0.0001) and a trend toward an increase in CD4⁺ Tconv cells, was observed. However, no significant changes in overall frequency of CD45⁺ cells between each phase was seen, although we did observe a greater deviation from the mean at later time points (Supplementary Fig. S2B). Thus, these conventional immunophenotyping studies showed that CD4⁺ Tconv cells were present at high frequencies during the elimination phase but were lost during transition to equilibrium. This coincided with a temporary expansion of granulocytes, which were numerically replaced in part by CD8⁺ T cells during the escape phase.

Based on the above observations, we tested whether CD8⁺ T cells, CD4⁺ T cells, and/or granulocytes had functional significance in restraining tumor progression during immunoediting of NPK-C1. For these studies, we performed antibody-mediated in vivo depletion experiments at both early and late stages (Fig. 2A and C). First, we administered depleting anti-CD4 or anti-CD8 on days 5 and 8 after implantation to determine each cell type’s contribution to elimination. We verified >95% depletion of both populations in spleen and tumor by flow cytometry (Supplementary Fig. S3). We found that depletion of CD8⁺ T cells or CD4⁺ T cells derepressed NPK-C1 growth kinetics, indicating both CD8⁺ and CD4⁺ T cells were required to suppress NPK-C1 progression to the escape phase (Fig. 2B). By contrast, depletion of granulocytes using alternating daily injections of anti-Ly6G (clone 1A8) and a secondary depleting antibody (24) showed no significant impact on tumor growth, suggesting that granulocytes were not necessary for elimination of NPK-C1 (Supplementary Fig. S4).

To test the requirement for these cells during the equilibrium and early escape phases, we administered depleting antibodies on days 17 and 20. Again, we found that CD8⁺ or CD4⁺ T-cell depletion promoted NPK-C1 growth, indicating that both populations are required to restrain tumor progression at the late equilibrium/early elimination phase (Fig. 2D). By contrast, granulocyte depletion beginning on day 17 yielded no clear phenotype, suggesting granulocytes were not involved NPK-C1 control or escape, although we acknowledge the relative difficulty of fully depleting granulocytes with 1A8 (Supplementary Fig. S3 and S4). Regarding the T-cell compartment, these data support the well-established role of CD8⁺ T cells as critical effectors of tumor control throughout immunoediting and are consistent with a growing body of literature implicating CD4⁺ T cells as critical effectors of the antitumor immune response (25).

As shown in Fig. 1, at approximately day 20, tumor growth in the NPK-C1 model bifurcated between progressing tumors and those that remained in equilibrium. This delayed emergence of progression, which was not observed in immunocompromised mice, afforded a unique opportunity to compare immunologic escape to sustained equilibrium in synchronously implanted, untreated tumors. For such analyses, we focused on day 24 samples, where we observed a roughly even split between progressing tumors and those which were stable or regressing, labeled together as “regressors” (Fig. 3A and B). Our initial analyses focused on immune cell populations defined by canonical gating strategies (Supplementary Fig. S2A); in later analyses (Figs. 4–6), we performed unsupervised clustering to define novel populations. These experiments showed that regressing tumors on average contained higher frequencies of CD45⁺ cells compared with progressors (P = 0.039; Fig. 3C). We also found a statistically significant direct correlation between tumor size and CD45⁺ frequency (P = 0.0129; Fig. 3D). However, this correlation was only significant on days 24 and 20, and not at earlier (day 10 or 15) time points (Supplementary Fig. S5). Using conventionally defined (canonical) populations, we found that the only immune cell population with significant enrichment in regressors relative to progressors was Tregs (Fig. 3E; P < 0.003). Tregs were the only population whose frequency correlated significantly with decreased tumor size on day 24 (Fig. 3F and G; P = 0.012). Of note, Tregs here were defined broadly, thus it is possible that subpopulations within the FoxP3⁺ T-cell subset could show the expected direct correlation with tumor mass. Granulocytes showed a trend toward correlation with larger tumors; however, this did not reach statistical significance after correcting for multiple comparisons (Fig. 3F and G; P = 0.089). Functionally, CD8⁺ and CD4⁺ T-cell cytokine expression potential was unchanged on average throughout immunoediting phases, as measured by intracellular cytokine staining upon ex vivo restimulation. However, on day 24, we observed significant negative correlations between expression of IFNγ, granzyme B, and/or IL2 and tumor mass (Supplementary Fig. S6). These findings suggest bulk population analyses were insufficient to identify functionally significant immune subsets involved in the equilibrium to escape transition in NPK-C1 tumors.

To leverage the high dimensionality of this longitudinal flow data set, we used dimensionality reduction and minimally supervised clustering algorithms. Specifically, we concatenated all live CD45⁺ cells from each tumor and then performed tSNE to generate a two-dimensional projection of all CD45⁺ cells arranged according to their relative similarity in 28-parameter space. For those clustering studies, we utilized the FlowSOM algorithm (26). After identifying common bulk immune populations as metaclusters in an unbiased fashion (Fig. 4A), we next sought to identify phenotypically unique subpopulations within key metaclusters by performing an additional round of unbiased, algorithm-based clustering on each population of interest. Thus, first we isolated the CD8⁺ T-cell metacluster, then performed tSNE and FlowSOM using markers with variable expression in the CD8⁺ T-cell compartment. This analysis yielded a CD8⁺ T-cell tSNE map with 7 subclusters designated Tc1-Tc7 (Fig. 4B). We next used the Mann–Whitney U test to verify that each cluster possessed at least one phenotypic marker that was statistically differentially expressed relative to all other cells (Supplementary Table S2). We implemented this into our analysis to ensure clusters were meaningful and to exclude any potential clustering based primarily on noise in the high-dimensional data. For each cluster, we next enumerated the number of cells within the cluster on a per-tumor basis to calculate the frequency of each cluster as a percentage of CD45⁺ cells over time in the longitudinal data set (Fig. 4C). Finally, we calculated the normalized MFI for each marker of interest in each cluster to facilitate interpretation of cell phenotype and inferred functional state (Fig. 4C). We profiled both subcutaneous and orthotopically implanted NPK-C1 tumors by this approach, and, finding them largely concordant, focused our analyses on our larger subcutaneous tumor data set (Supplementary Fig. S7). These studies showed that CD8⁺ T-cell clusters Tc1–Tc4 were relatively numerous in NPK-C1 tumors, whereas Tc5–Tc7 were rare cells. In general, Tc1–Tc4 increased significantly on days 20 and/or 24 and could be differentiated, in part, by expression of PD-1, Tim-3, and Ly6C. Correlating each cluster frequency with tumor mass showed that only two of the CD8⁺ T-cell clusters correlated significantly (Fig. 4D and E). One cluster, Tc7, was a granzyme B⁻ CD8⁺ population with moderate PD-1 and Tim-3 expression and low expression of CD44 (Fig. 4C). Although undetectable in tumors <25 mg, Tc7 represented ≥1% of CD45 in tumors with a mass >50 mg (Fig. 4E) and may represent an “exhausted” population. By contrast, Tc5 CD8⁺ T cells were more common, and their frequency correlated significantly with reduced tumor size. Tc5 expressed CD62L, Ly6C, and granzyme B, while lacking PD-1 (Fig. 4C) and appeared to represent an activated central memory-like CD8⁺ cluster. Together, these data highlight the significant heterogeneity within CD8⁺ T cells that infiltrate NPK-C1 tumors and suggest that central memory-like CD8⁺ T cells are associated with productive antitumor immunity that suppresses tumor escape, whereas progressive tumors are enriched with granzyme B^neg potentially exhausted CD8⁺ T cells.

Given our previous observations that CD4⁺ T cells were critical throughout the immunoediting process (Fig. 3E and F), we applied the analysis above to the CD4⁺ T-cell population. As shown in Fig. 4A, FlowSOM clustering of the entire CD45⁺ compartment showed two CD4⁺ T-cell metaclusters, grossly corresponding to CD4⁺ Tconv cells and Tregs. To further resolve these compartments, we performed tSNE and FlowSOM analysis of each metacluster independently, using the analytical techniques described for Fig. 4. These analyses revealed seven unique subclusters within the CD4⁺ Tconv compartment, labeled Tn1–Tn7. Expression data (Fig. 5B) suggested that Tn1 and Tn2 were FoxP3⁺Helios⁺ cells that clustered with CD4⁺ Tconv cells, although we cannot rule out activation-induced expression of these markers (Fig. 5B). Each CD4⁺ Tconv cluster was present at a relatively high frequency at day 10, at reduced frequencies at days 15 and 20, with partial recovery by day 24 (Fig. 5B). Clusters Tn3–Tn7 were predominantly differentiated by expression of PD-1, Tim-3, and Ly6C, and Tn4 was the sole CD4⁺ Tconv subpopulation that correlated significantly with reduced tumor mass on day 24. Although Tn2 also correlated with smaller tumors, it phenotypically resembled a Treg phenotype (discussed below). Among CD4⁺ clusters, Tn4 cells showed the highest expression of PD-1, Tim-3, and CD44 (Fig. 5C), suggesting that this population may include the activated CD4⁺ T cells presumed to be involved in preventing escape (Fig. 2), although their lack of granzyme B expression is likely inconsistent with direct cytotoxic activity.

Similar analyses of the FoxP3⁺ (presumptive Treg) population showed phenotypic diversity, and we identified 10 discrete clusters, Tr1–8 (and Tn1–2) that varied significantly in their expression of Helios, CD103, KLRG1, Tim-3, and Ly6C (Fig. 5D). As above, we validated the discreteness of these clusters by Mann–Whitney U Test to exclude noise-based or artificial overclustering. The majority of these clusters were exceedingly rare (<0.25% of CD45) at days 10 and 15 in NPK-C1 tumors; however, approximately half of these expanded through days 20 and 24 to reach 0.5% to 1% of CD45 (Fig. 5D). Two putative Treg clusters, Tn2 and Tr4, were significantly correlated with reduced tumor mass on day 24. Both Tn2 and Tr4 expressed Helios, CD44, and Tim-3, had relatively elevated Ki67 expression, but differed in their expression of KLRG1 and CD103 (Fig. 5C and D). The divergence between the putative Treg populations was notable, where Tr8 correlated with larger tumor mass (escape) and Tr4 and Tn2 correlated with smaller tumors. Tr8 expressed KLRG1 and CD103, markers which we (27) and others previously associated with suppressive capacity, whereas Tn2 expressed granzyme B, possibly consistent with cytotoxic capacity. Taken together, these data showed heterogeneity in tumor-infiltrating CD4⁺ Tconv cells and Treg populations, a finding not well appreciated when immune cell populations are analyzed according to prespecified phenotypes.

Applying the FlowSOM algorithm to the non–T cell CD45⁺ compartment yielded four clusters that comprised a single myeloid metacluster in tSNE space (Figs. 4A and 6A). We thus concatenated these clusters together for downstream analyses, using an additional round of tSNE dimensionality reduction and FlowSOM clustering as described above. Again, clusters were statistically validated by Mann–Whitney U test to avoid noise-based overclustering. These high-dimensional analyses yielded nine unique myeloid subclusters M1–M9 (Fig. 6B). M1–M4 represented monocyte populations, which ranged in frequency from 2% to 10% of CD45 on average (Fig. 6C). M6–M9 resembled DC-like cells having moderate to high CD11c and low CD11b expression and were generally rare (<2% of CD45) with the exception of M8. M5 expressed both CD11b and CD11c and thus was difficult to define (Fig. 6C). Most of these did not significantly vary over time, although M4 showed an atypical up-down–up-down behavior and M2 was significantly reduced after day 10 (Fig. 6C). Focusing first on clusters significantly correlated with smaller NPK-C1 tumors, we found M1 (a proliferating monocyte population), M6 (an iNOS⁺ DC-like population), and M8 (a classic CD11c⁺MHC-II^hi DC population). The enrichment of proliferating monocytes in regressing tumors (M1) corresponded with an enrichment of nonproliferative monocytes in progressing NPK-C1 tumors (M2; Fig. 6D and E). Thus, in the context of NPK-C1, local monocyte proliferation versus quiescence on day 24 separated regressing from progressing tumors. Focusing further on clusters enriched in progressing NPK-C1 tumors, we found two MHC-II^lo DC-like populations in M7 and M5, and CD103⁺ DCs in M9. Given the established importance of Batf3-lineage CD103⁺ DCs in driving antitumor immunity, this association was unexpected. These data suggest that myeloid phenotypes are complex and that the particular populations associated with disease progression or regression may have context-dependent roles over the course of immunoediting.

Immunoediting involves two main inflection points: the elimination–equilibrium transition and the equilibrium–escape transition. To determine whether compositional changes occurring at the primary elimination–equilibrium transition were similar or different from those occurring at the secondary equilibrium–escape transition, we used a bootstrapping procedure to generate average fold-change values for all computationally defined subclusters between elimination (day 10) tumors and equilibrium (day 15) tumors (Fig. 7A), as well as between day 24 regressor and progressor tumors (Fig. 7B). We observed a number of clusters to be enriched in opposite directions in the early comparison versus the late comparison. For example, Tc7 (granzyme B⁻ putative exhausted CD8) existed at higher frequency in “eliminating” day 10 tumors relative to “equilibrating” day 15 tumors, yet Tc7 was reduced in “controlled” day 24 tumors relative to “escaping” day 24 tumors (Fig. 7B). These data implied that either (i) the role of Tc7 was dichotomous in these two conditions or (ii) changes that occurred as a tumor transitioned from elimination to equilibrium were unique from what happened during the switch from equilibrium to escape.

To define clusters exhibiting statistically significant enrichment in opposite directions between these two comparisons, we dichotomized each fold-change value depending on whether it was greater than or less than 1, then used Fisher exact test to identify significant differences (Fig. 7C). Overall, we found CD4⁺ Tconv-cell enrichment patterns to be largely consistent in early versus late comparisons; in other words, CD4⁺ Tconv clusters tended to be enriched in the tumor microenvironment (TME) where a productive response was occurring, regardless of the inflection point (Fig. 7C). By contrast, we observed broad divergence of myeloid clusters. The only myeloid clusters enriched in inflamed tumors at both inflection points were M6 (CD11c⁺iNOS⁺) and M8 (CD11c⁺MHC-II^hiCD103⁻) DC-like cells. Proliferating monocytes, although enriched in regressing day 24 tumors, were contrarily enriched in the elimination–equilibrium transition. The opposite was true for nonproliferative monocytes, suggesting that the association between monocyte proliferation and NPK response was only valid in later stages of tumor outgrowth. Finally, we also found the CD103⁺ DC cluster M9 to be similarly divergent in early versus late comparisons (Fig. 7A–C). Together, these data highlight the dynamic nature of the TME, in that immune cell populations associated with a favorable phenotype at one time point may show an opposite association at another. These data also suggest canonically defined phenotypes are likely insufficient to capture the immune cell diversity in the TME.

In this study, we performed high parameter longitudinal immune profiling and cellular depletion assays to define the cellular determinants of immunoediting in NPK-C1, a transplantable murine model of prostate cancer that spontaneously models the distinct stages of immunoediting: elimination, equilibrium, and escape. This is, to our knowledge, one of the few transplantable models that consistently exhibit such behavior in C57BL/6 mice, with a tumor volume and time-scale that allows for facile longitudinal immunophenotyping at each stage. Although the paucity of alternative models restricted our study to this single-cell line, our observation that NPK-C1 exhibited immune-mediated editing at both subcutaneous and orthotopic sites with relatively minor differences in immune contexture lends credence to the relevance of this model. Overall, these studies revealed four key insights: (i) CD4⁺ T cells are required for productive antitumor immunity during elimination phase and play active roles in restraining tumor growth during equilibrium phase; (ii) Tregs exhibit diverse phenotypes in the TME and can be associated with favorable or unfavorable response depending on phenotype; (iii) myeloid phenotypes are complex and context dependent, as CD103⁺ DCs appear to be associated with progression to escape phase in this model; and (iv) to understand complex immune biology within the TME, data sets with increasing dimensionality in both marker number and time points may be critical, along with unbiased high-dimensional analysis techniques.

Our data substantiate a growing body of literature emphasizing the role of CD4⁺ T cells in antitumor immunity (25). Our study is reminiscent in approach to that used by Wei et al, who used CyTOF profiling to find an ICOS⁺ T_H1 effector CD4⁺ T-cell subset enriched in mice with MC38 and B16 tumors treated with CTLA-4 blockade (28). Alspach and colleagues identified class II–restricted neoantigens in 3-methylcholanthrene (MCA) sarcomas and showed neoantigen-reactive CD4⁺ T cells are required at the tumor site for spontaneous and therapeutically augmented antitumor immunity (9). A report by Oh and colleagues used single-cell sequencing from a cohort of human bladder cancer specimens to identify cytotoxic CD4⁺ T cells within the TME and found this cytotoxic CD4 gene signature predicts response to anti–PD-L1 in a validation cohort (10). In relation to these studies, our data set is limited in its ability to phenotype CD4⁺ T cells in that we could not specifically delineate cells recognizing tumor antigen. However, our model uniquely allowed us to demonstrate that CD4⁺ Tconv activation was associated with productive antitumor immunity throughout each stage of immunoediting. Although we did not observe substantial granzyme B positivity in the CD4⁺ Tconv compartment in NPK-C1, future studies are required to detail how cytokine secretion and/or cytotoxicity contribute to the antitumor activities of CD4⁺ Tconv cells.

Treg frequency as a prognostic biomarker in cancer is controversial. Although Treg infiltration is associated with poor outcomes in a majority of cancers, in colorectal and some other cancers, Tregs associate with favorable response (29). Our data mirror these seemingly contradictory findings. We also found that Tregs were phenotypically diverse, as we identified 10 distinct Treg phenotypes in NPK-C1 tumors using FoxP3, Helios, KLRG1, CD103, Tim-3, and Ly6C markers alone. Although we have previously shown Helios and CD103 expression associate with suppressive tumor Tregs (27), these data suggest heterogeneity among Tregs bearing these markers; Tr4 expressed Helios and CD103 and correlated with tumor regression, whereas Tr8 expressed CD103 and correlated with tumor progression in this model. Little is known about what drives these specific phenotypes and to what extent there are functional differences in Tregs at this level of granularity. Although these data clearly require further investigation, an important implication of this result is that high-dimensional approaches will be of increasing importance in defining Treg phenotypes with consistent prognostic or predictive ability in cancer, particularly in human studies.

Our data corroborate the widely acknowledged understanding that myeloid cell behavior is complex and context dependent in the TME. Myeloid cells can play active roles in immunoediting, both in antagonizing or promoting tumor escape (11, 16). Regarding tumor antagonistic function, O’Sullivan and colleagues showed IFNγ-stimulated M1 macrophages expressing iNOS and MHC-II mediate immunoediting in RAG2^−/−γc^−/− mice (11). Using single-cell RNA sequencing and CyTOF of checkpoint blockade treated MCA sarcomas, Gubin and colleagues elegantly established a direct association between tumor regression and iNOS⁺ M1-like macrophage polarization, whereas tumor progression was associated with CD206⁺ M2-like skewing (30). Although our data included few macrophages, we observed the iNOS⁺ phenotype in DC-like cells enriched in regressor tumors, confirming the use of iNOS as a marker defining myeloid populations associated with productive antitumor immunity in mice. We observed in progressor tumors a higher frequency of CD103⁺ DCs, which are generally thought to be required for productive antitumor immunity (4, 31–33). This finding was unexpected, given the role of CD103⁺ DCs in tumor antigen cross-presentation to prime CD8⁺ T-cell responses (31, 33). We note, though, that CD103⁺ DCs were enriched in elimination phase tumors relative to equilibrium phase tumors, as expected. These data may possibly be explained in two ways: (i) a high frequency of CD103⁺ DCs in the TME following initial elimination and equilibrium phases indicates failed DC migration to draining lymph nodes, which results in insufficient tumor-specific CD8 priming (32, 34); or (ii) CD103⁺ DCs can undergo immunosuppressive polarization in the TME, analogous to the mature regulatory DC phenotype described by Maier and colleagues (35). Further investigation into the transcriptional profile of these cells, as well as their cross-presenting functional capacity, is required to understand the curious enrichment of CD103⁺ DCs in the TME upon progression to escape in our model.

In summary, our study introduces the utility of the prostate tumor model NPK-C1 for studying mechanisms of immunoediting in C57BL/6 mice. Using high-dimensional longitudinal immune profiling and in vivo depletion experiments, we found that central memory–like CD8⁺ T cells, activated CD4⁺ effectors, proliferative Helios⁺ Tregs, and iNOS⁺ myeloid cells associated with immunity involved in tumor escape, whereas exhausted CD8⁺ T cells, Helios⁻CD103⁺KLRG1⁺ Tregs, and CD103⁺ DCs were enriched in progressing tumors. We also found compositional changes that were unique between each immunoediting inflection point, highlighting the constantly evolving nature of tumor–immune interactions over time. These observations illustrate the importance of relevant models and high-dimensional phenotyping in understanding the complex immune biology underlying cancer immunoediting.

A. Califano reports other support from DarwinHealth, Inc. outside the submitted work; patent 10,790,040 issued and with royalties paid; and is founder of, equity holder in, and consultant for DarwinHealth, Inc., a company that has licensed some of the algorithms used in this article from Columbia University. Columbia University is also an equity holder in DarwinHealth, Inc. C.G. Drake reports personal fees from AstraZeneca/MedImmune, from Bristol-Myers Squibb, Compugen, F-Star, Genocea, Janssen, Merck, Harpoon, Tizona, Werewolf, Genentech, and Rubius outside the submitted work, and is currently an employee of Janssen Pharmaceuticals. No disclosures were reported by the other authors.

C.R. Ager: Conceptualization, data curation, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. A.Z. Obradovic: Formal analysis, methodology, writing–review and editing. J.M. Arriaga: Resources, writing–review and editing. M.G. Chaimowitz: Project administration, writing–review and editing. A. Califano: Resources, supervision, writing–review and editing. C. Abate-Shen: Resources, supervision, writing–review and editing. C.G. Drake: Conceptualization, resources, supervision, funding acquisition, writing–review and editing.

This research was funded in part through the NIH/NCI Cancer Center Support Grant P30CA013696 (C.G. Drake), the Prostate Cancer Foundation (C.G. Drake), the Movember Challenge Award (C.G. Drake), the NIH R01 CA183929-05 and R01 CA173481 (C. Abate-Shen), the NIH R35 CA197745 (A. Califano), the NIH U01 CA217858 (A. Califano), and the NIH S10 OD012351 and S10 OD021764 (A. Califano). J.M. Arriaga is currently supported by the Dean’s Precision Medicine Research Fellowship from the Irving Institute for Clinical and Translational Research at Columbia University Irving Medical Center (UL1TR001873). The authors also acknowledge Dr. Siu-hong Ho and Dr. PJ Stewart-Hutchinson at Cytek Biosciences for assistance and training on spectral flow cytometry. They also thank the CCTI at CUIMC for departmental support.

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.