Abstract
The tumor immune microenvironment (TIME) is commonly infiltrated by diverse collections of myeloid cells. Yet, the complexity of myeloid-cell identity and plasticity has challenged efforts to define bona fide populations and determine their connections to T-cell function and their relationship to patient outcome. Here, we have leveraged single-cell RNA-sequencing analysis of several mouse and human tumors and found that monocyte–macrophage diversity is characterized by a combination of conserved lineage states as well as transcriptional programs accessed along the differentiation trajectory. We also found in mouse models that tumor monocyte-to-macrophage progression was profoundly tied to regulatory T cell (Treg) abundance. In human kidney cancer, heterogeneity in macrophage accumulation and myeloid composition corresponded to variance in, not only Treg density, but also the quality of infiltrating CD8+ T cells. In this way, holistic analysis of monocyte-to-macrophage differentiation creates a framework for critically different immune states.
Introduction
A key component of most immune responses, including those to cancers, are mononuclear phagocyte cell populations, which share common features of phagocytosis, tissue repair, and immunoregulation but diverge in functional specialization. Conventional dendritic cells (cDC) are positioned in tissues to initiate and sustain adaptive T-cell responses (1), whereas macrophages engage in high rates of phagocytosis and tissue remodeling (2). Self-renewing tissue-resident macrophages are seeded during embryonic development (3), whereas inflammatory stimuli prompt infiltration of adult hematopoietic stem cell–derived monocytes that give rise to tumor macrophages (4–7). These monocyte-derived macrophages preferentially accumulate as tumors progress (8) and may predominate in regulating the ongoing antitumor T-cell response (9).
Macrophages consist of numerous subset populations that have been identified across tissues (10–13). Therapeutic blockade of key epigenetic and signaling pathways has demonstrated their amenability to transcriptional reprogramming (14), but how phenotypic diversity arises remains poorly understood. Recruited bloodborne monocytes exhibit plasticity in differentiation potential and can acquire features of macrophages and/or DCs depending on the inflammatory setting (5, 6, 10, 15–17). In addition, early studies demonstrated that macrophage exposure to type 1– or type 2–associated cytokines induces “M1” or “M2” cellular programs, respectively, and a model was put forth in which myeloid cells are polarized to be proinflammatory (“M1”) or anti- (“M2”) inflammatory (18–20). Although this nomenclature was thereafter understood to require nuance to account for additional plasticity (21), it remains undetermined whether these binary programs are applicable to describe tumor macrophage differentiation in vivo.
Myeloid phenotypic diversity has challenged efforts to utilize myeloid populations as biomarkers for patient treatment options and outcome. cDCs are critical for coordinating antitumor T-cell immunity (22–25) and higher cDC abundance is broadly associated with improved cancer patient survival, although additional tumor immune microenvironment (TIME) features may inform functionality (23, 24, 26). In contrast, macrophages have largely been considered to be protumorigenic (2, 14) and monocytes have often been described as myeloid-derived suppressor cells (MDSC; ref. 27). Yet, several studies have indicated that macrophages are not consistently a negative predictor of patient prognosis (28–31), and increased levels of circulating monocytes were unexpectedly linked to patient responsiveness to immune checkpoint blockade (ICB; ref. 32). These contrary findings speak to the need for improved resolution of myeloid-cell categorization and phenotype to dissect heterogenous responses among patients with cancer.
Using single-cell RNA sequencing (scRNA-seq), we uncovered transcriptional heterogeneity among tumor-infiltrating myeloid cells and distinguished monocyte and macrophage lineage- and activation-induced programs shared between multiple mouse tumor models and human kidney cancer samples. Monocyte differentiation is dynamically regulated, and we found that regulatory T cell (Treg) density was one immunoregulatory axis capable of modulating macrophage density. Further comprehensive analysis of key myeloid populations revealed distinct network connections between different myeloid-cell types and T-cell subsets, including Tregs and effector T cells. This is consistent with an archetypal organization of immune systems in tumors—collections of cell types that move together as modules (33)—and improved classification of patients such that we could identify those with effective antitumor T-cell responses.
Materials and Methods
Mice
The following mice were housed and/or bred under specific pathogen-free conditions at the University of California, San Francisco Animal Barrier Facility: C57BL/6J (The Jackson Laboratory), MMTV-PyMT-mCherry-OVA transgenic (34), and Foxp3-DTR (The Jackson Laboratory). All mice were handled in accordance with NIH and American Association of Laboratory Animal Care standards, and experiments were approved by the Institutional Animal Care and Use Committee of the University of California, San Francisco (San Francisco, CA).
Human tumor samples
Renal cell carcinoma (RCC), melanoma, and head and neck tumor samples were transported from various cancer operating rooms or outpatient clinics. All patients provided informed written consent to the UCSF IPI clinical coordinator group for tissue collection under a UCSF Institutional Review Board (IRB)-approved protocol (UCSF IRB# 20-31740) in accordance with the Declaration of Helsinki guidelines. Patients were selected without regard to prior treatment and 21 RCC (all primary tumors), 22 melanoma (20 primary and 2 metastasis tumors), and four head and neck tumors (all primary tumors) were collected. All samples were defined as primary tumor or metastasis by pathology assistants. Samples were obtained after surgical excision with biopsies taken to confirm the presence of tumor cells. Freshly resected samples were placed in ice-cold Dulbecco's Phosphate Buffered Saline (DPBS; Thermo Fisher Scientific, catalog no. 14190144) or Leibovitz's L-15 medium (Thermo Fisher Scientific, catalog no. 11415064) in a 50 mL conical tube and immediately transported to the laboratory for sample labeling and processing. As described below (Human Tissue Processing and Flow Cytometry Analysis), the whole tissue underwent digestion and processing to generate a single-cell suspension. In the event that part of the tissue was sliced and preserved for imaging analysis, the remaining portion of the tissue sample was used for flow cytometry analysis.
Tumor cell lines
B16-F10 cells (ATCC, CRL-6475) were purchased in 2015 and stock vials were generated from an initial thaw. Cells in use were taken from early passages and maintained consistent and homogenous morphologic characteristics, during which time they were tested for Mycoplasma. B16-F10-ZsGreen was previously generated in our laboratory as described previously (35) and tested for maintained expression of ZsGreen by flow cytometry. After thawing, tumor cells were cultured at 37°C in 5% CO2 in DMEM (Thermo Fisher Scientific, catalog no. A4192101), 10% FCS (Benchmark, catalog no. 100-106), penicillin, streptomycin, and l-Glutamine (Thermo Fisher Scientific, catalog no. 10378016). Cells were generally cultured for 2 to 5 days (0–1 passages) before use for subcutaneous injection.
Mouse tumor cell injections and growth
Prior to injection, adherent B16-F10 or B16-ZsGreen tumor cells were dissociated with 0.05% Trypsin-EDTA (Thermo Fisher Scientific, catalog no. 25300120) and washed 2–3× with DPBS (Thermo Fisher Scientific, catalog no. 14190144). A total of 1.0 × 105 to 2.5 × 105 cells were resuspended in DPBS and mixed 1:1 with Matrigel GFR (Corning, catalog no. 356231). Mice were injected subcutaneously with a volume of 50 μL either unilaterally or bilaterally depending on the experimental setup. Tumor tissue was harvested approximately 12 to 16 days later.
MMTV-PyMT-mCherry-OVA transgenic mice were bred and genotyped for the transgene by PCR. Spontaneous tumor growth was monitored by measuring with electronic calipers. Tumors were harvested when the mice were approximately 20 to 30 weeks of age such that palpable tumors had developed but had not exceeded a size of 100 mm2.
In vivo mouse treatments
To deplete Tregs, Foxp3-DTR and control mice were injected intraperitoneally with 500 ng of unnicked diptheria toxin (List Biologics, catalog no. 150). Mice were typically injected 9, 10, and 12 days following initial inoculation with tumor cells (see Mouse Tumor Cell Injections and Growth).
For specified experiments, wild-type mice were injected intraperitoneally 7, 9, 10, 11, and 13 days following tumor injection with 250 μg of anti-mouse CTLA-4 IgG2c (modified clone 9D9, Bristol-Myers-Squibb), mouse IgG2C isotype, anti-mouse CTLA-4 IgG1 (modified clone 9D9, Bristol Myers Squibb), or mouse IgG1 isotype.
Mouse tissue processing and flow cytometry analysis
Mouse tumor tissue was harvested and enzymatically digested with 0.2 mg/mL DNase I (Sigma-Aldrich, catalog no. D5025), 100U/mL Collagenase I (Worthington Biochemical, catalog no. LS004197), and 500 U/mL Collagenase Type IV (Worthington Biochemical, catalog no. LS004189) for 30 minutes at 37°C while under constant agitation. Blood was collected via cardiac puncture from mice that were euthanized by overdose with 2.5% Avertin. Blood samples were treated with 175 mmol/L NH4Cl for 5 minutes at room temperature to lyse red blood cells.
Samples were filtered, washed with stain media (DPBS, 2% FCS), and resuspended in stain media. Cells from this single-cell suspension were washed with DPBS and stained with Zombie NIR fixable viability dye (BioLegend, catalog no. 423106) for 30 minutes at 4°C. Cells were washed and resuspended in stain media containing anti-CD16/32 (BioXCell, catalog no. BE0307), 2% rat serum (Thermo Fisher Scientific, catalog no. 10710C), 2% Armenian hamster serum (Innovative Research, catalog no. IGHMA-SER), and antibodies against surface proteins of interest. Cells were stained for 30 minutes at 4°C. In some experiments, cells were then washed and stained for intracellular protein levels, for which they were fixed, permeabilized, and stained according to BD Cytofix/Cytoperm Kit (BD Biosciences, catalog no. 554655) or the FoxP3/Transcription Factor Staining Buffer Set (Thermo Fisher Scientific, catalog no. 00-5523-00).
The following antibodies were from BioLegend: anti-mouse CD45 (clone A20, catalog no. 110727), anti-mouse Ly-6C (clone HK1.4, catalog no. 128037), anti-mouse CD11b (clone M1/70, catalog no. 101257), anti-mouse CD11c (clone N418, catalog no. 117339), anti-mouse MHC-II (clone M5/114.15.2, catalog no. 107622), anti-mouse F4/80 (clone BM8, catalog no. 123135, 123107, or 123131), anti-mouse CD24 (clone M1/69, catalog no. 101822), anti-mouse Ly-6G (clone IA8, catalog no. 127645), anti-mouse NK1.1 (clone PK136, catalog no. 108749), anti-mouse CD90.2 (clone 30-H12, catalog no. 105331), anti-mouse/human CD45R/B220 (clone RA3-6B2, catalog no. 103246), anti-mouse CD301b (clone URA-1, catalog no. 146814 or 146803), anti-mouse CD64 (clone X54-5/7.1, catalog no. 139306), anti-mouse CD127 (clone A7R34, 135031). The following antibodies were from BD Biosciences: anti-mouse Siglec-F (clone E50-2440, catalog no. 740956), anti-mouse CD106 (clone 429, catalog no. 745672). The following antibodies were from R&D: anti-mouse/human ARG1 (polyclonal, catalog no. IC5868F), normal sheep IgG (polyclonal, catalog IC016F). The following antibodies were from Thermo Fisher Scientific: anti-mouse FoxP3 (clone FJK-16s, catalog no. 48-5773-82).
Following staining, cells were washed, resuspended in stain media, and analyzed on a BD Biosciences Fortessa or sorted with a BD Biosciences FACSAria Fusion. Flow cytometry data were analyzed using FlowJo software (BD Biosciences, version 9 or 10).
Human tissue processing and flow cytometry analysis
Human tumor or metastatic tissue was thoroughly chopped with surgical scissors and transferred to gentleMACS C Tubes (Miltenyi Biotec) containing 20 μL/mL Liberase TL (5 mg/mL, Roche, catalog no. 5401020001) and 50 U/mL DNAse I (Roche, catalog no. 10104159001) in RPMI1640 (Thermo Fisher Scientific, catalog no. 11875093) per 0.3 g tissue. gentleMACS C Tubes were installed onto the gentleMACS Octo Dissociator (Miltenyi Biotec) and incubated for 45 minutes according to the manufacturer's instructions. Samples were then quenched with 15 mL of sort buffer (DPBS, 2% FCS, 2 mmol/L EDTA), filtered through 100 μm filters, and spun down. Red blood cell lysis was performed with 175 mmol/L NH4Cl if needed. Cells were incubated with Human FcX (BioLegend, catalog no. 422301) to prevent nonspecific antibody binding. Cells were then washed in DPBS and incubated with Zombie Aqua Fixable Viability Dye (Thermo Fisher Scientific, catalog no. L34957). Following viability dye, cells were washed with sort buffer and incubated with cell surface antibodies that had been diluted in the BV stain buffer (BD Biosciences, catalog no. 563794) for 30 minutes on ice. Cells were subsequently fixed in either Fixation Buffer (BD Biosciences, catalog no. 554655) or in Foxp3/Transcription Factor Staining Buffer Set (Thermo Fisher Scientific, catalog no. 00-5523-00) if intracellular staining was required.
The following antibodies were from BD Biosciences: anti-human HLA-DR (clone G46-6, catalog no. 564040), anti-human CD56 (clone NCAM16.2, catalog no. 564448), anti-human CD127 (clone HIL-7R-M21, catalog no. 563225), anti-human CD25 (clone 2A3, catalog no. 340939), anti-human CD45RO (clone UCHL1, catalog no. 561889), anti-human PD-1 (clone EH12, catalog no.563789), anti-human CTLA-4 (clone BNI3, catalog no. 565931), and anti-human CD64 (clone 10.1, catalog no. 564425). The following antibodies were from Thermo Fisher Scientific: anti-human CD45 (clone HI30, catalog no. 47-0459-42), anti-human CD3ε (clone OKT3, catalog no. 46-0037-42), anti-human FoxP3 (clone 236A/E7, catalog no. 25-4777-41), anti-human Ki-67 (SolA15, catalog no. 11-5698-82), anti-human CD19 (clone H1B19, catalog no. 45-0199-42), anti-human CD20 (clone 2H7, catalog no. 45-0209-42), anti-human CD56 (clone CMSSB, catalog no. 46-0567-42), and anti-human CD11c (clone 3.9, catalog no. 56-0116-42). The following antibodies were from BioLegend: anti-human CD4 (clone S3.5, catalog no. 100455), anti-human CD8α (clone RPA-T8, catalog no. 301039), anti-human CD38 (clone HIT2, catalog no. 303523), anti-human CD16 (clone 3G8, catalog no. 302039), CD1C/BDCA-1 (clone L161, catalog no. 331515), anti-human CD14 (clone M5E2, catalog no. 301837), anti-human CD304 (clone 12C2, catalog no. 354503), and streptavidin. Anti-human BDCA-3 (clone AD5-14H12, catalog no. 130-098-843) was purchased from Miltenyi Biotec.
Stained cells were washed and analyzed on a BD Biosciences Fortessa or sorted with a BD Biosciences FACSAria Fusion. Flow cytometry data were analyzed using FlowJo software (BD Biosciences, version 10.6). Clustering and heatmap analyses were performed using Morpheus (Broad Institute).
scRNA-seq data generation
Mouse B16 samples were generated over two independent experiments. Samples were pooled from at least 5 mice per experiment to ensure representation across a cohort of tumor-bearing mice. In the first experiment, CD45+CD90−B220−NK1.1−Ly6G− cells that were Ly6C−MHC-II+ or Ly6C+CD11b+ were sorted as one bulk myeloid sample. Individual monocyte (CD45+CD90−B220−NK1.1−Ly6G−Ly6C+CD11b+) and macrophage (CD45+CD90−B220−NK1.1−Ly6G−Ly6C−MHC-II+F4/80+CD24loCD11clo/hi) populations from these tumors were also sorted. Each of these samples was processed separately, but in parallel, for scRNA-seq analysis. In the second experiment, B16 tumor myeloid cells (CD45+CD90−B220−NK1.1−Ly6G− cells that were CD11b+ and/or CD11c+) were sorted from control or Foxp3-DTR mice. Blood myeloid cells (CD45+CD90B220−NK1.1−Ly6G− cells that were CD11b+ and/or CD11c+) were also sorted from these B16 tumor-bearing wild-type mice. Each of these three samples was processed separately, but in parallel, for scRNA-seq analysis. In a subsequent experiment, myeloid cells (CD45+CD90−B220−NK1.1−Ly6G−CD11b+ and/or CD11c+) from a mouse PyMT tumor were sorted. In addition, myeloid cells (CD45+CD3ε−CD19−CD20−CD56−HLA-DRdim/hi) from a total of one RCC, four head and neck, and six melanoma biopsies were sorted and processed individually as they became available.
Once sorted, cells were resuspended at a concentration of 1 × 103 cells/μL in media (DPBS, 0.04% BSA) and loaded onto the Chromium Controller (10× Genomics). Samples underwent single-cell encapsulation and cDNA library preparation using the Chromium Single Cell 3′ v1 or v2 Reagent Kits (10× Genomics, catalog no. 120237). The cDNA library was sequenced on an Ilumina HiSeq 4000 or Illumina Novaseq.
scRNA-seq data processing
Sequencing data were processed using the 10× Genomics Cell Ranger V1.2 pipeline. Fastq files were generated from Ilumina bcl files with the Cell Ranger subroutine mkfastq. Fastq files were then processed with Cell Ranger's count to align RNA reads against UCSC mm10 or GRCh38 genomics for mouse and human cells, respectively, using the aligner STAR (36). Redundant unique molecular identifiers (UMI) were filtered, and a gene–cell barcode matrix was generated with count. Mkfastq and count were run with default parameters.
For mouse B16 tumor samples, the gene–cell barcode matrix was passed to the R software package Seurat (v2.3.0; ref. 37) for all downstream analyses. Genes that were expressed in at least three cells were included. Cells that did not express at least 200 genes were excluded, as were those that contained >5% reads associated with cell-cycle genes (38, 39). For mouse PyMT and human tumor samples, raw feature–barcode matrices were loaded into Seurat (v3.1.5; ref. 40) and genes with fewer than three UMIs were dropped from the analyses. Matrices were further filtered to remove events with greater than 20% percent mitochondrial content, events with greater than 50% ribosomal content, or events with fewer than 200 total genes. The cell-cycle state of each cell was assessed using a published set of genes associated with various stages of human mitosis (41).
Using Seurat's ScaleData function, read counts were log2 transformed and scaled using each cell's proportion of cell-cycle genes as a nuisance factor. A set of highly variable genes was generated by Seurat's FindVariableGenes function, which were used for principal component analysis (PCA). Genes associated with principal components (selected following visualization with scree plots) were used for graph-based cluster identification and dimensionality reduction using t-distributed stochastic neighbor embedding (t-SNE) or Uniform Manifold Approximation and Projection (UMAP) analysis. Seurat's FindAllMarkers function was used for subsequent cluster-based analyses, including cluster marker identification and diffe-rentially expressed (DE) gene analyses.
scRNA-seq signature generation
To generate mouse monocyte- and macrophage-specific gene signatures, sorted monocyte, CD11clo macrophage, and CD11chi macrophage samples were aggregated, log2 transformed, and scaled using Seurat. DE gene analysis was performed using FindMarkers with the parameters log N fold change > 1.5 and a min.pct of 0.25. Genes were selected by ranked fold change and the criteria that min.pct1/min.pct2 > 1.5. Genes used for cell-cycle regression analysis were excluded. The top 10 genes (or fewer if less remained) were median normalized and aggregated, scaled 0–1, and plotted on specific t-SNE plots.
Gene signatures for cellular programs such as metabolism (42), “M1” and “M2” polarization (43), and MHC-II–associated genes (GSEA, REACTOME_MHC_CLASS_II_ANTIGEN_PRESENTATION), previously published cell types (44, 45), and cell populations identified here were generated by taking the mean of log-normalized expression for a particular set of genes related to the specific pathway or phenotype. To visualize the distribution of these scores across cells, we binarized the distribution of the score at the 70th percentile unless specified otherwise and overlaid on the calculated t-SNE coordinates.
For correlation analysis of “M1” and “M2” genes, the expression of each gene in the signatures was calculated for each B16 tumor Csf1r+Mafb+ cluster cell and binarized at the median. A cross-correlation gene–gene matrix was obtained using the R cor function with method = “pearson.”
scRNA-seq sample aggregation
Pairwise comparison analyses were performed between B16 tumor myeloid-cell clusters from wild-type and Treg-depleted mice. For this, one sample from wild-type mice and one sample from Foxp3-DTR mice were used, along with an additional wild-type sample that had been generated in a previous independent experiment. The 3 objects were first transformed from Seurat v2 to Seurat v3. The raw UMI counts were renormalized using person residuals from “regularized negative binomial regression,” with sequencing depth a covariate in a generalized linear model via the R sctransform package (46). Pairwise “anchor” cells were identified between the 3 objects using the original wild-type mouse sample as a reference via the Seurat FindIntegrationAnchors function. Briefly, each pair of samples was reduced to a lower dimensional space using diagonalized canonical correlation analysis using the top 3,000 genes with the highest dispersions. The canonical correlation vectors were normalized using L2-normalization. Multiple nearest neighbors across datasets were identified for each cell in each dataset and cell–cell similarities are used as anchors to integrate the datasets together using the Seurat IntegrateData function.
For the integration of all human samples, the individually processed samples were normalized, and variance stabilized using negative binomial regression via the scTransform method offered by Seurat (46). Counts matrices were merged into a single Seurat object and the batch (or library) of origin was stored in the metadata of the object. The log-normalized counts were reduced to a lower dimension using PCA and the individual libraries were aligned in the shared PCA space in a batch-aware manner (each individual library was considered a batch) using the Harmony algorithm (47). The resulting Harmony components were used to generate a batch-corrected UMAP, and to identify clusters of transcriptionally similar cells across each of the 11 samples.
scRNA-seq pseudotime analysis
Raw UMI counts from the cleaned and processed Seurat objects for the control and Treg-depleted mouse experiment were extracted and coerced into Monocle2 (48, 49). CellDataSet objects were generated, normalizing the data using a negative binomial distribution with fixed variance (negbinom.size). Each object was independently processed to identify a pseudotime trajectory. Briefly, each object was processed to estimate per-cell coverage and sequencing depth (estimateSizeFactors) and gene dispersions (estimateDispersions). Cells were further filtered to retain high-quality cells with ≥500 genes and genes were filtered to retain only those in at least 10 cells. The dataset was reduced to two dimensions using the DDRTree algorithm and the marker genes that differentiated the Ly6c2+Hp+ monocytes and C1qa+ macrophage clusters from other cell types were used to guide the trajectory inference. Relative pseudotime was obtained through a linear transformation relative to the cells with the lowest and highest pseudotimes (1-min_pseudotime)/max_pseudotime. The “wave” plots were constructed using the Seurat LogNormalized counts and the relative pseudotime described above for the top five genes per cluster as identified by scRNA-seq.
Human samples were analyzed with Monocle3 (48, 49), and the cell_data_set object along with the batch-corrected PCA and UMAP embeddings were imported directly from the Seurat object. Each cell-specific trajectory was inferred by reverse embedding the UMAP coordinates using the DDRTree method. Relative pseudotime was obtained as described for the mouse tumor samples.
The Cancer Genome Atlas survival analyses
Tumor RNA-seq counts and transcripts-per-million (TPM) values for 985 kidney renal clear cell carcinoma (KIRC) samples from the Toil recompute data in The Cancer Genome Atlas (TCGA) Pan-Cancer cohort (50) were downloaded from the UCSC Xena browser (51). A gene signature score for each patient was calculated using the gene signature score method below. The feature gene signature scores were calculated using an m × n matrix where m represented the TPM normalized log2 counts per million expression of the feature signature genes and n represented the selected sample set (52). The expression of each gene was converted to percentile ranks across the samples using the SciPy Python module (53). The top and bottom 30 percentile were then used to define low and high patients for each respective signature unless otherwise noted.
Statistical analysis and data visualization
Comparisons between groups were analyzed using GraphPad Prism software. Experimental group allocation was determined by genotype or by random designation when all wild-type mice were used. Error bars represent mean ± SEM calculated with Prism unless otherwise noted. Comparisons between two groups were analyzed with Student t test. For statistical measures between more than two groups, one-way ANOVA was performed unless otherwise noted. Nonsignificant comparisons are not shown. Investigators were not blinded to experiment group assignment during experimental data generation or analyses. The R packages Seurat and ggplot2 were used to generate figures.
Data availability
The expression matrices for the scRNA-seq samples reported in this article can be found with the following Gene Expression Omnibus accession numbers: GSE188548, GSE184096, GSE159913, and GSE184398. Code used to generate the scRNA-seq analyses is included on Github (https://github.com/UCSF-DSCOLAB/mujal_et_al_MonoMac_2021).
Results
Mouse B16 tumors harbor a diversity of myeloid states
Subcutaneous implantation of B16 melanoma cells is a well-established mouse tumor model with abundant infiltration of monocytes, macrophages, and cDCs (23). To study these cells along their differentiation trajectories, we used conventional markers to sort bulk myeloid populations, along with reference populations of Ly6C+ monocytes and two tumor-associated macrophage (TAM) populations, distinguished on the basis of level of expression of CD11c and MHC-II (Fig. 1A; Supplementary Fig. S1A; ref. 23). Each of these samples was then subjected to scRNA-seq analysis.
Within the bulk myeloid population, t-SNE clustering yielded eight transcriptionally distinct cell populations (Fig. 1B; Supplementary Fig. S1B and S1C), including three Flt3+Kit+ cDC populations (clusters 4, 6, 7), which were marked by signatures specific to cDC1s, cDC2s, and conserved cDC activation programs (Supplementary Fig. S1D; Supplementary Table S1; refs. 23, 24, 26). The remaining myeloid cells (clusters 0, 1, 2, 3, 5) broadly expressed Csf1r and Mafb (Fig. 1C), indicative of monocytes and macrophages. In addition, these cells broadly expressed Ccr2 and modest but appreciable levels of Cx3cr1 (Supplementary Fig. S1E). Having focused on the stimulatory capacity of cDCs in previous work (23, 24, 45), here we focused on the diversity of monocytes and macrophages as it related to the TIME.
To align transcriptional cell-type categorization with flow cytometry analysis, we generated cell type–specific gene signatures from the scRNA-seq analysis of the FACS-sorted monocytes and TAMs (Fig. 1A; Supplementary Fig. S1F). When applied (Fig. 1D), these indicated that four Csf1r+Mafb+ populations (clusters 0, 1, 2, 5) expressed monocyte-specific genes. The four monocyte populations expressed Ly6c2, but varied in levels of other monocyte-associated genes (e.g., Hp, Chil3) and, as found in cluster 0, also expressed TAM-associated genes (e.g., H2-Ab1, C1qa, Ms4a7; Fig. 1E and F). Monocyte-like clusters were differentiated from one another by cellular activation programs. For example, cluster 1 (“IFN-responsive”) was specifically enriched for IFN-inducible genes such as Cxcl10, Gbp2, and IFIT-family members. Cluster 2 (“stress-responsive”) cells expressed Arg1 and were enriched for cellular stress processes, including oxidative stress–responsive genes and HSP genes such as Hmox1, Hspa1a, Hilpda, Bnip3, Ero1L, and Ndrg1 (Fig. 1F; Supplementary Fig. S1G). In contrast to the heterogeneity observed amongst monocytes, signatures for both populations of TAMs localized within cluster 3 (Fig. 1D).
We applied pseudotime analysis (48) to generate a model of tumor monocyte-to-macrophage differentiation (Fig. 1G and H; Supplementary Fig. S1H). This model placed cluster 5 Ly6c2+Hp+ monocytes and cluster 3 C1qa+ TAMs at opposite ends of a linear trajectory, consistent with our expectations. Cluster 0 monocytes occupied the continuum between them and expressed a combination of both monocyte- and TAM-associated gene signatures such that we designated these cells “Intermediate monocytes” (“Mono-Int”). Kinetic analysis of cluster-enriched genes confirmed gradual downregulation of Ly6c2+Hp+ monocyte-associated genes and upregulation of “Mono-Int”- and TAM-associated genes along the pseudotime trajectory (Fig. 1I). This transcriptional model thus supported a framework of progressive monocyte-to-TAM differentiation, in which Ly6C downregulation is paired with up-regulation of CD64, MHC-II, and F4/80 (Supplementary Fig. S1I and S1J; ref. 15). In contrast, IFN- and stress-responsive cells occupied distinct branches that diverged from the dominant differentiation trajectory at intermediate timepoints (Fig. 1H and I; Supplementary Fig S1H).
Heterogeneous acquisition of “stress-” and “IFN-responsive” cellular programs
To gain higher resolution on the differentiation trajectories within the monocyte/macrophage lineage, we performed cluster analysis on the sorted monocyte and TAM samples. Sorted monocytes expressed Ly6c2 and contained clusters similar to those identified within the bulk myeloid-cell sample (Fig. 2A; Supplementary Fig S2A), indicating that these cells may not consist purely of early-stage monocytes but also include some cells that have acquired macrophage attributes (Supplementary Fig. S1J). Cluster analysis of CD11clo and CD11chi TAMs, however, resolved diversity beyond the C1qa+ TAM signature (Fig. 2B; Supplementary Fig. S2B and S2C) including identifying clusters enriched for cell cycle–related genes, and an Mgl2+ TAM subset that expressed immune modulators such as Ccl6, Il1b, and Retnla as compared with the C1qa+ cluster, which more highly expressed genes such as Ms4a7. Although these cells had not formed a distinct population in our original analysis of bulk myeloid cells (Fig. 1), we did retrospectively detect Mgl2+ cells in that scRNA-seq dataset, as well as by flow cytometry (Supplementary Fig. S2D). TAM-subset clusters were also accompanied by an Arg1+ stress-responsive cluster akin to that found in the sorted monocytes (Fig. 2B; Supplementary Fig. S2B and S2C). Indeed, reclustering of the entire stress-responsive cluster from the bulk tumor myeloid sample revealed that this program was acquired by monocytes, “Mono-Int” and TAMs (Fig. 2C; Supplementary Fig. S2E).
Segregated expression of stress-responsive genes and canonical TAM-associated genes suggested divergent transcriptional programs and we sought to determine whether these populations could also be distinguished by flow cytometry. Differential gene expression analysis of the stress-responsive and C1qa+ TAM clusters from our bulk myeloid-cell sample revealed cluster-specific expression of cell surface genes Il7r and Vcam1, respectively (Fig. 2D). Using the same gating as in Supplementary Fig. S1A, we confirmed this split in both “Mono-Int” and TAMs (Fig. 2E) and we found enriched arginase 1 (ARG1) expression in both IL-7Rα+ populations (Fig. 2E and F; Supplementary Fig. S2F). As expected from the single-cell transcriptional analysis, VCAM1+ cells were more abundantly found within TAMs (Fig. 2E and F; Supplementary Fig. S2F).
Together, this dissection of sorted cell populations lent support to a model in which monocytes and TAMs exist in a differentiation trajectory, along which cells can adopt specialized cellular programs (Fig. 2G and H). Some programs, such as those associated with Mgl2+ or Vcam1+ TAMs, selectively emerged later, in mature TAMs. Others, such as IFN-induced signaling or stress-responsiveness may be more universally accessible across differentiation stages. In addition, we detected populations of IFN-responsive monocytes in the peripheral blood of B16 tumor-bearing mice (Supplementary Fig. S2G and SSH), perhaps suggesting that systemic IFN signaling, or other induction of this program, may define monocytes prior to tumor entry. In contrast, stress-responsive populations were not detected in the blood, suggesting that microenvironmental cues in the TIME likely induce this activation program locally. Further studies are warranted to explore whether these programs directly influence monocyte differentiation processes or act as “layers” that accessorize a canonical differentiation trajectory.
Mouse tumor macrophage subset heterogeneity does not reflect “M1/M2” polarization
Macrophage exposure to type-1 or type-2 cytokines in vitro results in “M1” and “M2” transcriptional signatures that are often used to describe “proinflammatory,” or “anti-inflammatory” and wound-healing processes, respectively (18–20). To address whether “M1/M2” polarization was a useful construct to define tumor macrophage diversity in vivo, we tested how “M1” and “M2” gene signatures (43) corresponded to the tumor myeloid-cell subsets profiled here. Using correlation and clustering analyses (Fig. 3A; Supplementary Fig. S3A), we found that, contrary to in vitro findings, tumor myeloid cells were marked by broad expression of both “M1”- and “M2”-associated genes, and we did not observe substantial correlation of gene expression within “M1” or “M2” gene groups across single cells. These data suggest that although tumor myeloid cells can express individual “M1” and “M2” genes, they rarely do so in any distinguishably consistent way during unperturbed tumor growth. We next processed tumor myeloid cells from MMTV-PyMT spontaneous mammary carcinomas for scRNA-seq analysis, sorting on Lin−CD11c+ and/or CD11b+ cells to capture the full cadre of myeloid populations including MHCII+/− cells (Supplementary Fig. S1A). We found that MMTV-PyMT tumors shared populations with the identical signatures as those defined for B16 tumors in Fig. 1, albeit in different proportions, and also showed a lack of coassociation between “M1” and “M2” signatures amongst the clusters (Fig 3B and C; Supplementary Fig. S3B–S3E).
While myeloid-cell populations appeared to be largely defined by differentiation stage and activation programs, we considered whether other core cellular features could help to further distinguish subsets across diverse microenvironments. It is now increasingly appreciated that metabolic reprogramming accompanies differentiation of immune cells, including macrophages (54). Indeed, assessment of metabolism-related genes (42) demonstrated that glycolysis-associated genes were specifically enriched in the stress-responsive cell cluster, whereas genes pertaining to oxidative phosphorylation were specifically enriched in C1qa+ TAMs in two distinct mouse models (Fig. 3D and E; Supplementary Fig. S3F). This suggests that these populations have additional important biological features in common—namely those coupled to distinct bioenergetic processes and demands.
Together, our data provide compelling evidence that “M1” and “M2” pathways have limited use in defining in vivo tumor myeloid-cell differentiation and subset plasticity during normal tumor development. Rather, common microenvironmentally induced programs and associated metabolic programs may yield greater insight in efforts to transcriptionally define and selectively target monocyte/TAM subsets.
Human RCC-infiltrating monocytes and macrophages mirror murine populations
We next assessed how the mouse monocyte/macrophage transcriptional programs we identified might compare with those from human cancers. We performed scRNA-seq analysis on HLA-DRdim/+Lin− myeloid cells sorted from an RCC sample, which are described to have substantial myeloid-cell infiltration (55), as well as six melanoma and four head and neck cancer samples (Fig. 4A and B; Supplementary Fig. S4A). Signatures derived from previously described blood myeloid-cell populations (11, 44) guided cluster identification and exclusion of cDCs (Supplementary Fig. S4B). Analysis of the CSF1R+MAFB+ clusters revealed a heterogenous collection of monocytes and macrophages with varying levels of CD14 and CD16 (Fig 4C–E).
As in mouse models, we detected early-stage CD14+S100A8+ classical monocytes along with terminally differentiated C1QC+ TAMs (Fig. 4E and F; Supplementary Table S2). Another population was CD14+ and differentially expressed LYPD3 and MHC-II genes, consistent with intermediate differentiation of monocytes towards TAM (“Mono-Int”; Fig. 4E and F). A population of FCGR3A+ non-classical monocytes also expressed IFN-stimulated genes and thus appears to functionally represent “IFN-responsive” cells (Fig. 4E and F; Supplementary Fig. S4C).
Finally, we found that there were a mix of cells on the monocyte–macrophage trajectory that expressed the stress-responsive program identified in mice, including the gene SPP1 (Fig. 4F; Supplementary Fig. S4C). When compared further with C1Q+ TAMs, this SPP1+ cluster was less mature based on higher expression of monocytic markers (i.e., S100A genes) and lower expression of MHC-II–related genes (Supplementary Fig. S4D). This analysis revealed also a population marked by expression of the antioxidant gene SEPP1 (Fig. 4F). These cells largely resembled C1Q+ TAMs but were enriched for FOLR2 (Fig. 4F; Supplementary Fig. S4E). FOLR2 is a marker previously associated with tissue-resident macrophages in breast cancer (56). Pseudotime analysis recapitulated a monocyte-to-macrophage differentiation trajectory (Fig. 4G and H) but did not connect the SEPP1+ cluster to the other monocyte–macrophage subsets, potentially due to the distinct ontogeny of tissue-resident macrophages, and thus this cluster was not considered for further trajectory analysis (Fig. 4G and H).
As in the mouse samples, the stress- and IFN-responsive programs aligned over the monocyte-to-macrophage trajectory, although in these samples, IFN-responsive monocytes appeared more advanced in differentiation stage. Again, there was broad coexpression of “M1”- and “M2”-associated genes across the populations (Supplementary Fig. S4G). Also, as in mice, there was a striking enrichment of a glycolytic signature (42) within the stress-responsive (SPP1+) cluster as compared with the C1Q+ TAMs, supporting the notion that these cells were functional orthologs in the two species (Fig. 4I). Altogether, these data confirm the limitation of “M1” and “M2” applicability in human tumors and illustrate the ability of other pathways to better define myeloid-cell subsets in vivo.
Treg depletion elicits reprogramming of the tumor myeloid-cell compartment
Myeloid-cell density can vary across patients (55), but how myeloid-cell infiltration and differentiation is collectively regulated in human cancer is still not well understood. When we quantified myeloid-cell populations in 20 RCC patient biopsies using flow cytometry, we found that the proportion of myeloid cells among live immune cells was increased in tumors of greater size and later stages (Supplementary Fig. S5A). Closer examination revealed that the ratio of macrophages-to-monocytes was also specifically increased in more advanced tumors (Fig. 5A, top). This suggested that the balance between monocytes and macrophages is dynamically regulated and that tumor growth is tied to higher macrophage density.
We thus sought other immunosuppressive parameters that might work in concert with increased macrophage accumulation. Tregs are a potent immunosuppressive force in the TIME and ablation can result in tumor clearance (24, 57). We found that Tregs accumulated as kidney tumor size increased (Fig. 5A, bottom), and that Treg abundance correlated well with macrophage-to-monocyte ratios in RCC as well as melanoma (Fig. 5B). The positive correlation between Treg and macrophage density spurred us to ask whether one caused the other. Using Foxp3-DTR mice, we found that depletion of Tregs dramatically reduced the macrophage-to-monocyte ratio in mouse B16 tumors (Fig. 5C, top) as assessed by flow cytometry. These results were phenocopied by treatment of mice with an anti–CTLA-4 that specifically depletes tumor Tregs (Fig. 5C, bottom; Supplementary Fig. S5B–S5C; refs. 24, 58). The shift in macrophage-to-monocyte ratio observed following both methods of Treg loss preceded subsequent tumor growth control (24, 58).
To further examine how Tregs may be influencing monocyte and macrophage proportions, we performed scRNA-seq analysis on mouse tumor Lin−CD11b+ and/or CD11c+ myeloid cells from B16 tumor-bearing control and FoxP3-DTR mice. Csf1r+Mafb+ clusters from this experiment were aggregated with those from the original wild-type B16 tumor sample in Fig. 1 and we observed similar cell populations across both experiments and treatment conditions (Fig. 5D; Supplementary Fig. S5D and S5E). Cluster proportions were modestly shifted with Treg loss (Fig. 5D), but cells from control and Treg-depleted tumors shared similar differentiation trajectories (Fig. 5E). However, Monocle analysis revealed differences in the accumulation of cells along the trajectory. Namely, whereas tumor monocytes, “Mono-Int,” and TAMs from the control sample acquired progressively increased pseudotime scores, “Mono-Int,” and TAM populations in the Foxp3-DTR sample did not exhibit sequential increases in pseudotime scores (Fig. 5E). In effect, TAM progression appeared stunted following depletion of Tregs.
In addition to increased expression of inflammatory and immunomodulatory genes (e.g., Ccl24, Arg1, Retnla, Mmp12, Mmp13, Nos2), expression of monocyte-associated genes was sustained in TAMs from Treg-depleted tumors (Fig. 5F and G; Supplementary Table S3). Moreover, expression of genes tied to macrophage differentiation (e.g., C1qa, H2-Ab1, Apoe, Ms4a7) was decreased across stages of differentiation (Fig. 5H; Supplementary Fig. S5F), further indicating these TAMs were more immature. Our analysis suggests that Treg depletion may impair implementation of TAM transcriptional programs, a remodeling detected early during tumor monocyte differentiation. Altogether these findings support a model in which Treg abundance promotes an accumulation of terminally differentiated TAMs in both mouse and human tumors.
Multiparametric immune-cell analysis improves classification of patients with kidney cancer
Given this association between T-cell subset density and TAM maturation, we sought to further explore how features of tumor macrophage infiltration could be harnessed to reliably inform features of patient outcome, such as survival. Analysis of TCGA KIRC samples using a myeloid gene signature from CiberSort (59) demonstrated that patients with varying levels of overall myeloid-cell density did not significantly differ in their survival (Fig. 6A, left). We next stratified TCGA patient data based on levels of the monocyte–macrophage lineage genes CSF1R and MAFB, finding that patients with higher levels of these had modest improvements in outcome (Fig. 6A, middle). As these genes are not strictly macrophage specific, we leveraged our scRNA-seq analyses of human RCC samples to generate signature scores based on the ratio between macrophage and monocyte (Fig. 4E). However, no significant differences in survival were revealed using this metric (Fig. 6A, right).
As TAM density did not appear to robustly inform patient outcome, we sought to test how TAM abundance corresponds with other immune parameters and may stratify patients with kidney cancer, using flow cytometry analysis of their biopsies. We thus performed unbiased clustering analysis using measurements of myeloid-cell, Treg, and conventional T cell (Tconv) frequencies. This revealed three groups of patients that were characterized by nearly binary (all high or all low) levels of CD8+ T cells, CD4+ Tconv, cDC2s, and cDC1s (Fig. 6B). These groups could be further parsed by varied macrophage–monocyte density: low (pink), moderate (yellow), and high (red; Fig. 6B, right). Consistent with our previous finding (Fig. 5), the group with the highest degree of macrophage differentiation was associated with the consistent presence of Treg infiltration (Treg-Mp). Of two groups with lowest macrophage differentiation, one (CD8-Mo-cDC1; pink) was distinguished by notable infiltration of cDC1s, which are critical for CD8+ T-cell responses (22, 23), and that group presented with uniformly high CD8+ T-cell infiltration. A second group had higher frequencies of both CD4+ T cells and cDC2s along with a variable amount of Tregs, consistent with the demonstrated role of cDC2s in supporting both CD4+ Tconv and Treg responses (24). Furthermore, the same patients clustered together again based on phenotypic analysis of the CD8+ T-cell compartment (Fig. 6C). Notably, CD8+ T cells from the CD8-Mo-cDC1 group (pink) expressed low levels of the exhaustion markers PD-1 and CD38 (Fig. 6C and D) and were also distinguished by higher expression of the checkpoint regulator CTLA-4, which may indicate ongoing activation (Fig. 6C and D; ref. 60). In contrast, the Treg-Mp (red) group showed the highest levels of both exhaustion markers and proliferative capacity (i.e., Ki-67).
In testament to the heightened antitumor CD8+ T-cell profile associated with low macrophage and Treg abundance but high cDC1 density, the subset of patients with these attributes (pink) had dramatically improved survival, showing no mortality for over 3 years (Fig. 6E). This multiparametric clustering parsed patients with the highest survival rates more profoundly from our dataset than the sole metric of cDC1 infiltration (Supplementary Fig. S6A). Similarly, although sole measurement of a cDC1 signature alone corresponded to higher survival rates among patients with TCGA (Supplementary Fig. S6B) in support of previous studies (23, 45), a combined measurement of the ratio of cDC1s to macrophages through combined gene signatures using the C1Q+ macrophage gene signature (Supplementary Fig. S6C and S6D) allowed for identification of patients with kidney cancer with better survival. Thus, fine-tuned stratification of the kidney cancer TIME provided the resolution critical for identifying three biologically distinct patient classes including this CD8-Mo-cDC1 group, which defines patients with the best CD8+ T-cell infiltration and outcome.
Discussion
Here we undertook scRNA-seq analysis of tumor monocytes and macrophages to determine the key hallmarks of their transcriptional diversity. We found two types of differentiation manifest during tumor development. On the one hand, we found a classical lineage differentiation trajectory that progresses from monocytes-to-macrophages in a way that has been long appreciated (61) with a discernable “intermediate” monocyte (“Mono-Int”) cell population. A “Mono-Int” population is, for reference, well described in other settings. For example, Randolph and colleagues detect “intermediate” monocytes in lymphoid and nonlymphoid tissue in steady-state conditions (62), and fluorescent real-time lineage tracing identifies cells undergoing that transition during allergic challenge (63).
On the other hand, we found two differentiation layers—“stress-responsive” and “IFN-responsive”—that coexist along that trajectory and that were shared across multiple mouse models as well as a profiled human RCC biopsy (Figs. 1 and 4). These programs were also present in other recently published studies (11, 13, 64, 65). For example, in a pan-cancer study, Cheng and colleagues discern myeloid populations whose primary distinction is their expression of IFN-induced genes (e.g., ISG15+ TAMs; ref. 11). Similarly, we noted that the stress-responsive population shares characteristics with cells historically contained within MDSCs (i.e., Arg1 expression and glycolytic programming; ref. 27). A notable difference in our interpretation compared with these previous reports lies in our incorporation of these layers within the monocyte–macrophage differentiation axis, rather than proposing them as a unique trajectory. Through independent profiling of purified monocytes and macrophages in our study and pseudotime analysis (Figs. 1 and 2), we find the stress-responsive signatures evident in both cell populations and indeed across them. In additional support of such a view, we found that an IFN-responsive signature was enriched among monocytes in one mouse model and macrophages in another (Figs. 1 and 3). We believe that this indicates that macrophages can differentiate in two dimensions—progression through the classical lineage as well as acquisition of specialized states characterized by examples of IFN or stress exposure. For these reasons, we prefer employing a nomenclature system that integrates the degree of monocyte-to-macrophage differentiation first, followed by additional transcriptional and functional qualities. Intuitively, this is similar to CD4+ T cells that can differentiate along a naïve–effector–memory axis while also being able to layer on Th1/Th2/Th17 programs.
Despite the comparison with CD4+ T cells, we do not find any populations, nor indeed any cells, that have an exclusively “M1” or “M2” signature (Fig. 3). Individual genes such as Arg1 are associated with certain clusters, as some have observed previously (11), but both correlation and signature analyses fail to identify any of the described “M1” or “M2” genes as either being selectively linked with one another in single cells, or as key classifiers of cell clusters. To this extent, the “M1/M2” nomenclature has provided direction in the fruitful study of myeloid-cell signaling and differentiation but does not appear to be accurately categorize distinct differentiation states, at least for tumors in vivo. We note the absence of data to the contrary of this conclusion in other recent reports (43, 55, 65), although of course individual nomenclature (e.g., “M2-like”) is clearly a matter of choice and needs only discussion as to which part of the in vitro signature might be biologically relevant.
One important aspect of myeloid-cell biology that requires further elaboration is how to identify IFN- and stress-responsive phenotypes. For example, Gubin and colleagues use iNOS as a marker by flow cytometry to define the IFN-stimulated population induced by ICB (12), whereas Cheng and colleagues utilize ISG15 (11). Particularly in the former study, which studied macrophage identity following ICB therapies, varied levels of type I and II IFNs may also modulate properties of this differentiation layer. In the case of “stress-responsive” populations, our data also point to IL7Rα expression, which may indicate involvement of TSLP signaling through heterodimeric pairing with TSLPR (66). An important set of conserved genes for “stress-responsive” macrophages, taken from our article, is their consistent and significant enrichment for glycolytic genes, particularly in comparison with conventional C1q “mature” TAMs. Given that HIF-1 is known to induce glycolytic genes under inflammatory and/or hypoxic conditions (54), this finding raises the question of whether these cells are selected for hypoxic environments where oxidative phosphorylation may not proceed, as well as their specific function. Going forward, the use of multiplexed imaging technologies such as ion beam imaging (MIBI) and single-cell spatial transcriptomics will enable this question to be addressed.
Our investigation of monocyte–macrophage differentiation led us to explore how its regulation could inform our understanding of antitumor immunity. Analysis of RCC and melanoma patient cohorts revealed an increase in macrophage-to-monocyte ratios with tumor grade, a rise that coincided with Treg density and was Treg dependent. Tregs exert potent immunosuppression and are thought to restrain T-cell activity and antitumor responses through modulation of DC stimulatory capacity, production of immunosuppressive cytokines and substrates, and competitive usage of growth factors and metabolic byproducts (24, 67, 68). It is becoming clear now that tumor Tregs also strongly influence the monocyte–macrophage lineage, likely through multiple mechanisms. In a recent study, tumor Tregs promoted tumor macrophage numbers by supporting their mitochondrial capacity and viability (69). Here, our scRNA-seq data demonstrates that early-stage monocytes and “Mono-Int” cells are already unable to properly implement TAM-associated transcriptional programs in the absence of Tregs, indicating that Tregs also fuel macrophage differentiation processes. This liaison between Tregs and macrophages mirrors one identified in the adipose fat of lean mice, where Tregs are thought to actively maintain homeostasis and hold inflammatory macrophages at bay (70, 71). Similarly, during the resolution of injury and inflammation in skeletal muscle and heart tissue, a transition from proinflammatory to anti-inflammatory macrophages occurs in a manner that appears to rely on Treg accumulation (72, 73). That Tregs may act on tumor macrophages in a similar fashion offers another example of how the TIME can exploit immune programs of “accommodation” that are otherwise in place to achieve tissue homeostasis in the face of perturbations (74).
Accumulation of a broad swath of macrophages in the TIME has previously been implicated with poor outcome (75). Consistent with this but at higher resolution, we detected a group of patients with kidney cancer for whom high macrophage-to-monocyte abundance was associated with diminished T-cell infiltration and exhaustion of those cells detected, concurring with other reports (55, 64). Our article thus points to an emerging trio of Tregs, macrophages, and exhausted T cells, whereby effector T cells may be corrupted through direct cellular interactions with TAMs, as has been suggested by observations of TAM–CD8+ T-cell colocalization in clear-cell RCC (ccRCC; ref. 64), or indirectly through macrophage-induced Treg expansion and activity (8, 76) or DC suppression (24, 77).
Yet, high myeloid-cell infiltration or skewed macrophage-to-monocyte ratios alone were not prognostic for KIRC patient survival. Indeed, although macrophages have often been found to be negatively associated with patient outcome, macrophage abundance as a sole biomarker has not been universally useful with prior studies similarly reporting instances in which macrophage abundance is not informative for patient cohorts with specific cancer subtypes, treatment regimens, or tumor stage (78–81). Clustering analysis of kidney TIME composition using comprehensive immune parameters, however, uncovered an archetype characterized by low macrophage-to-monocyte differentiation in conjunction with high cDC1 infiltration. These patients (CD8-Mo-cDC1) had elevated infiltration of CD8+ T cells with low surface expression of proteins associated with exhaustion and highly enhanced survival rates (Fig. 6, pink). Notably, recent work focused on ascertaining the different immune archetypes across solid tumors suggests that these patient groups, though most frequent in kidney cancer, span cancer types including frequent representation in colorectal and bladder tumors (33).
Identification of a CD8–Mo–cDC1archetype emphasizes the value of integrating multiparametric biomarkers as a means to better parse patient outcome and to establish principles of TIME organization. Given that T-cell activity appears to be collectively influenced by multiple immune cell populations with distinct partnering patterns, our analysis suggests that dual targeting of TIME axes may elicit the best CD8+ T-cell responses. For example, reprogramming and/or depletion of macrophages may relieve active suppression (2, 14) and strategies that boost cDC1 recruitment and survival (1) may further benefit even those with favorable macrophage-to-monocyte density. It is also notable that this protective archetype is specifically enriched for monocytes. Indeed, monocyte differentiation into macrophages may not be inevitable and accumulation of “Mono-Int” cells has been detected in multiple forms of inflammation (10, 16, 82, 83). In addition, the potential importance of monocytes is indicated by their increased numbers in the blood of ICB responsive as compared with patients with nonresponsive melanoma (32). In patients with ccRCC, IFN-responsive TAMs exhibited lower levels of HLA-DR, reminiscent of the “Mono-Int” cells described here, and higher levels of these ISGhi TAMs were predictive of survival after tyrosine kinase inhibitor treatment (79). Such a relationship opens questions across cancer type; namely, whether “Mono-Int” are distinct in their antitumor function, and how might monocytes be additive or synergistic with cDC1s to drive antitumor CD8+ T cells?
Altogether these findings contribute to the endeavor of clarifying useful distinctions in myeloid-cell gene expression and highlight settings in which multiparametric analysis of tumor myeloid-cell composition can inform patient immune archetype and guide development of relevant therapies.
Authors' Disclosures
A.M. Mujal reports grants from NIH and Cancer Research Institute during the conduct of the study. A.J. Combes reports grants from Pfizer, Abbvie, Amgen, Genentech, and Bristol Myers Squibb during the conduct of the study. J. Tsui reports other support from Pfizer, Amgen, Abbvie, and nonfinancial support and other support from Bristol Myers Squibb during the conduct of the study. S.P. Porten reports other support from Photocure and Coloplast outside the submitted work. M.F. Krummel reports grants from NIH, Abbvie, Amgen, Bristol Myers Squibb, and Pfizer during the conduct of the study; other support from Pionyr Immunotherapeutics and Founder Innovations|Studios outside the submitted work. No disclosures were reported by the other authors.
Authors' Contributions
A.M. Mujal: Conceptualization, formal analysis, investigation, visualization, writing–original draft, writing–review and editing. A.J. Combes: Conceptualization, formal analysis, investigation, visualization, writing–original draft, writing–review and editing. A.A. Rao: Data curation, software, formal analysis, visualization. M. Binnewies: Conceptualization, investigation, writing–review and editing. B. Samad: Formal analysis, investigation. J. Tsui: Investigation. A. Boissonnas: Investigation, writing–review and editing. J.L. Pollack: Data curation, software, formal analysis, visualization. R.J. Argüello: Investigation, writing–review and editing. M.V. Meng: Resources. S.P. Porten: Resources. M.K. Ruhland: Investigation. K.C. Barry: Investigation. V. Chan: Supervision, project administration. M.F. Krummel: Conceptualization, resources, supervision, funding acquisition, writing–original draft, project administration, writing–review and editing.
Acknowledgments
This work was supported by the U.S. NIH (R01CA197363 and U01CA217864, to M.F. Krummel). A.M. Mujal is supported by the Cancer Research Institute as a Cancer Research Institute/Amgen Fellow. Acquisition and analysis of certain human samples was partially funded by contributions from Abbvie, Amgen, Bristol Myers Squibb, and Pfizer as part of the UCSF Immunoprofiler Initiative.
We would like to thank E. Wan and the Institute for Human Genetics at UCSF for helping prepare samples for scRNA-seq as well as the UCSF Flow Cytometry Core for maintenance of flow cytometers and sorters. We would like to also thank J.J. Engelhardt and Bristol Myers Squibb for Fc-modified anti-CTLA-4 antibodies.
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.