Multi-cellular phenotypic dynamics during the progression of an immunocompetent breast cancer model

Abstract

The breast tumor microenvironment (TME) has recently been profiled at high resolution by performing single-cell RNA sequencing (scRNAseq) on patient samples. However, from patients' samples, analyzing the temporal dynamics of the TME is ethically, practically and scientifically challenging. Revealing these dynamics could structure inter-tumor heterogeneity into a temporally ordered sequence of causes and consequences in cellular events. Here, we survey the dynamics of the TME by performing scRNAseq at different time points of the progression of a mouse breast tumor allograft model driven by the PyMT antigen. We find that multi-cellular phenotypic dynamics follow one of three possible temporal patterns: stable colonization, wave-like, or progressive increase. In particular, IFN-responsive cancer cells, GzmB+ cytotoxic T cells, as well as C1q macrophages, increase in parallel with tumor progression. These findings establish the single-cell types and phenotypes in a progressing breast tumor, and reveal when these cellular players enter and leave the TME.

Keywords: cancer; immunology; transcriptomics.

Graphical abstract

Figure 1

Performing scRNAseq on tumors across multiple growth time points reveals the temporal dynamics of the transcriptional heterogeneity of the TME (A) An allograft model of breast cancer was used for the experiment. 86R2.2 cells derived from an MMTV-PyMT genetically engineered mouse model in a C57BL/6 background were labeled with Luciferase and GFP as reporters. 500 cells were injected in the mammary fat pads of 30 immunocompetent C57BL/6 mice. Tumor growth was monitored non-invasively by bioluminescence imaging. Continuous lines represent longitudinal bioluminescence measurements from the same individual. Tumors were harvested at four time points and scRNAseq performed in duplicates. Y axis: tumor volume estimated from the bioluminescence signal (STAR Methods). (B) Hierarchical principal component analysis coupled with Louvain clustering (colors) identifies 9 cell types in this data. Cell type markers expressed in each cluster appear in each box. Axes represent the two first PCs. (C) The temporal dynamics of the cellular composition of the TME have a biphasic structure. Tumors are initially rich in cancer and stroma, and stroma is subsequently diminished as the immune system colonizes the TME. In a second phase, Th cells, Tregs, and group 1 ILCs recede while cancer cell abundance increases. All data points are represented (dots), and the data are also summarized as a mean (solid line).

Figure 2

The prevalence of interferon (IFN)-responsive cancer cells, cytotoxic Tc and C1q macrophages progressively increases along tumor growth (A–E and G) The transcriptional heterogeneity of cancer cells describes a linear continuum of proliferation-survival polarization and a gradient in interferon (IFN) signaling and response. Dots: cancer cells from all time points. Axes: first three PCs of cancer cell gene expression. Color bar: log transcription quotients (LTQs), a measure of gene expression. (A) At one end of the linear continuum, we find the high expression of Plk1, a cell division gene. (B) Mmp2, an ECM remodeling enzyme associated with tumor invasion, follows the inverse gradient. (C) Expression of Irf7 is orthogonal to the survival-proliferation continuum. (D) Temporal changes in the transcriptional heterogeneity of cancer cells follow the interferon (IFN)-response gradient. (E) Summary of the transcriptional heterogeneity of cancer cells. (F) The abundance of interferon (IFN)-responsive cancer cells, proliferative Tc, and C1q macrophages progressively increases over the four-week experiment. All data points are represented (dots), and the data are also summarized as mean (solid line) ± one standard deviation (error bars). (G) Interferon (IFN)-responsive cancer cells have the highest predicted fitness among cancer cells. Color bar: Cancer cell fitness score computed from the CRISPR experiment of Lawson et al. (2020). p-value: Wilcoxon’s rank-sum test on fitness scores of the 5% closest cells to archetype 3 versus other cancer cells (Figure S2B). (H–M) The transcriptional heterogeneity of macrophages and monocytes is a weighted average of four archetypal phenotypes. Dots: macrophages and monocytes from all time points. Axes: first three PCs of gene expression (H). Individual cells fall on a simplex in the first three PCs. p-value: t-ratio test (I) Single cell from macrophage/monocyte cluster can be described as a weighted average of the four endpoints (archetypes) of the 3D simplex: monocytes, M2-like/immune-suppressive, M1-like/immunostimulatory, C1q macrophages. (J–M) Cells close to the different archetypes express specific phenotypic markers. (J) Ly6c2, a monocyte marker, is specifically expressed in cells closest to archetype 1. (K) Mrc1, an M2 marker, is expressed in cells closest to archetype 2. (L) Cells closest to archetype 3 upregulate Ciita, part of the M1-like/immunostimulatory macrophage signature. (M) Cells closest to archetype 4 upregulate Trf, an iron uptake protein.

Figure 3

The transcriptional heterogeneity of cytotoxic T cells is characterized by the expression of TCR vs. NK lectin receptor genes and a cytotoxicity/proliferation gradient, while DCs adopt three established phenotypes (A–E) Tc cells and DCs form phenotypic clusters. Labels represent pathways significantly associated with PCs. Arrows represent gene expression gradients, with the length of the arrows representing the strength of the gradient. (A–D) Cytotoxic T cells can adopt two phenotypes characterized by the expression of TCR or NK lectin receptor genes, with a shared cytotoxicity and proliferation gradient. (A) Tc cells that score high on PC1 upregulate the Xcl1 chemotaxis factor, which facilitates the recruitment of DCs. Color bar: Xcl1 expression (LTQs). (B) Tc cells that score high on PC2 upregulate the co-stimulatory receptor Icos. (C) Each PC is driven by genes with shared function in TCR or NK lectin receptor signaling. Tc cells that score high on both PCs express genes involved in proliferation, cytotoxicity. (D) Tc cells can adopt two phenotypes, characterized by the expression of Cd8, Icos, Cd28, and genes from the TCR pathway or by the expression of Xcl1 and NK-style lectin receptor genes. Both phenotypes share a gradient of activation and transcriptional activity. Color bar: number of UMI per cell, a proxy for transcriptional activity. (E–G) DCs can adopt three phenotypes: cDC1, cDC2, and mregDC. (E) Genes with the strongest expression gradients associate specifically with one DC phenotype. (F) Expression patterns of DC markers identify the phenotypes of the three DC clusters. Expression of MHCI and MHCII represents the average expression of all genes that comprise the MHC class I and II, respectively. (G) MregDCs express significantly higher levels of MHCI compared to classical DCs. (H) The temporal dynamics of mregDC, cDC1, cDC2, and monocytes follow a stable colonization pattern of the TME. All data points are represented (dots), and the data are also summarized as mean (solid line) ± one standard deviation (error bars).

Figure 4

The temporal dynamics of Th cells, Tregs, and ILCs are characterized by a wave-like pattern (A–C) The transcriptional heterogeneity of Th, Tregs, and group 1 ILCs is described by an activation gradient. (A) Group 1 ILCs are a mix of ILC1s (low PC1) and NK cells (high PC1), with a shared activation/cytotoxicity gradient along PC2. (B) Th transcriptional heterogeneity is described by a Th1 response gradient along PC1, and an axis of matrisome and antigen presentation versus proliferation along PC2. All data points are represented (dots), and the data are also summarized as mean (solid line) ± one standard deviation (error bars). (C) A linear phenotypic continuum of proliferation – activation – matrisome, ECM phenotypes characterizes the transcriptional heterogeneity of Tregs. (D) Time dynamics of Tregs follow PC1. (E) Treg phenotypic heterogeneity shifts from low PC1 (day 11) to medium PC1 (day 18) to high PC1 (day 24), indicating a temporal shift from a proliferative phenotype initially, to an activated phenotype at intermediate times, to a matrisome and extracellular phenotype at the latest time points. The phenotypic density (y axis) represents the prevalence of Tregs with a specific phenotype (defined as PC1 score), so that the area under each curve equals the fraction of Tregs at this time point. (F) During tumor progression, the transcriptional heterogeneity of Th cells shifts from a proliferative phenotype to a matrisome and antigen presentation phenotype. (G) The temporal dynamics of activated Th cells, Tregs, and ILCs follow a wave-like pattern. (H) The multi-cellular phenotypic dynamics of the TME during the progression of a breast tumor model is summarized by three temporal patterns: progressive increase, stable colonization, and wave-like.