Single-Cell Analysis Unveils the Role of the Tumor Immune Microenvironment and Notch Signaling in Dormant Minimal Residual Disease Free

Tumor relapse and metastasis are a major cause of death in women with breast cancer. By the time of primary tumor diagnosis, cancer may have already metastasized, and patients have either developed metastatic tumors or dormant disseminated tumor cells (DTC). DTCs have been detected even at the early preinvasive stages of tumor evolution (1). Tumor dormancy is a specific stage in which residual cells remain inactive and therapy resistant, possibly due to quiescence, niche protection, or both (2). Even after five years of adjuvant endocrine therapy, breast cancer recurrence continues to occur steadily between 5 and 20 years and is strongly correlated with DTCs in lymph nodes and tumor grades (1).

The role of the microenvironment in regulating the survival and proliferation of residual cells following therapy is understudied and a very limited number of microenvironmental factors have been shown to promote dormancy. Microenvironment niches have been implicated in either a conversion to cells with characteristics of cancer stem cells (CSC) or better survival of CSCs as compared with non-CSCs at the metastatic sites (3–7). The sprouting microvascular endothelial niche has been implicated in dormancy and survival of breast cancer cells at the metastatic sites and that inhibition of Notch signaling was shown to reduce neovascular tips, and therefore to suppress tumor cell outgrowth (8, 9). The majority of these studies used established cell lines. Various studies suggest that neutrophils awaken tumor cells from dormancy and facilitate therapeutic resistance. Neutrophils have been found in conjunction with circulating tumor cells (CTC) in patient samples and exerted higher survival capabilities and induced proliferation of CTCs while in circulation (10). Neutrophils are a major component of the lung premetastatic niche (11) and awaken dormant cells in part by releasing neutrophil extracellular traps (NET) and mediating microenvironment extracellular matrix (ECM) remodeling (12).

A conditional mouse model for Her2/neu-driven mammary tumors has been used to demonstrate that following Her2 removal, Notch signaling becomes activated in a subset of residual cells, resulting in accelerated tumor recurrence (13). Furthermore, following Her2 removal, TNFα/NFκB signaling led to elevated CCL5 expression. CCL5 promoted tumor recurrence by recruiting CCR5-expressing macrophages, which contributed to collagen deposition in residual tumors (14). Inhibition of CCL5, however, was not sufficient to prolong the time to relapse (14). One limitation of these studies is that they were performed in immunocompromised hosts.

Here, we report generation of a Wnt/iFGFR^GFP (iFGFR^GFP) tumor model to study cell signaling and the tumor microenvironment during minimal residual disease (MRD) and recurrence. One advantage of this model is that these mice have an intact immune system and that the GFP reporter facilitates the isolation and characterization of the small number of residual cancer cells in a tolerized host. We generated four stages of tumorigenesis: primary, MRD (dormant), long-term dormant (long-term DT), and recurrent and performed single-cell RNA sequencing (scRNA-seq) and immunophenotyping to unveil the changes in microenvironment.

Sterile surgical procedure and mouse handling were in accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the NIH (Bethesda, MD). The protocol was approved by the Association for Assessment and Accreditation of Laboratory Animal Care–accredited Baylor College of Medicine Institutional Animal Care and Use Committee (IACUC #AN-504).

MMTV-Wnt1/iFGRR1 (Wnt-iR1) mice were generated and characterized previously (15). The GFP mice [FVB.Cg-Tg(CAG-EGFP)B5Nagy/J- Stock No:003516], kindly provided by Dr. Michael Lewis (Baylor College of Medicine, Houston, TX), were bred with Wnt-iR1 bigenic mice to generate GFP⁺ tumors and tumor chunks were cryopreserved in 90% FBS + 10% DMSO for future experiments. The GFP⁺ tumors were then transplanted into tolerized Glow-Head (GH) recipients (FVB.rGH-Luc-GFP) kindly provided by Dr. Chi-Ping Day (NIH/NCI, Bethesda, MD; ref. 16).

Tumor pieces 1–2 mm³ in size were surgically transplanted into the cleared fat pad of 5-week-old GH mice. The fat-clearing and tumor transplantation technique were performed by modifying a procedure described previously (17). Briefly, instead of a traditional Y-shaped incision, a 4-mm medial semicircle incision was made around the fourth MG teat; for example, for the left (viewer's) gland cut the right side of the nipple. The skin was detached by a sterile cotton swab and the fat pad was pulled out by a narrow-head forceps (WPI-500451, World Precision Instruments) holding the tissue close to the lymph node. This part of tissue contained the epithelial cells that were to be removed. As we cauterized the major vessels and finally the tissue containing the lymph node, we retracted out the gland. While still holding on to the tissue, the forceps was placed on the body outside the cut. A very small tumor piece was placed into the fat pad. The tissue holding on to the forceps was finally cut and the wound was closed using 7-mm clips. Pre- and post-surgery procedures were performed in accordance with the Guide for the Care and Use of Laboratory Animals.

One day after tumor transplantation, mice were injected intraperitoneally (i.p) with the B/B homodimerizer (dimerizer, Clontech) every 3 days until tumors reached about 4–5 mm in diameter. At that point, tumors grew independently of the dimerizer. Once tumors reached about 7–9 mm in diameter, mice were treated with the FGFR inhibitor NVP-BGJ398 (BGJ; ref. 18) for 14 days by oral gavage (50–810–47, Thermo Fisher Scientific) and then drug treatment was stopped for another 14 days. At this timepoint, depending on the experiment, mice were: (i) sacrificed for residual tissue collection, or (ii) left untreated for recurrence, or (iii) treated with the α-Jagged-1 (Jag1) blocking antibodies kindly provided by Dr. Christian Siebel (Genentech), or α-IgG2b (BE0086, BioXcell).

Mice were treated with α-Jag1 or IgG control (15 mg/kg) once a week for 5 weeks and then treatment was stopped. Tumors were measured by a digital caliper and once they reached 5 mm in diameter, they were designated as “recurrent”.

Before tissue harvest, each mouse was perfused by 20 mL of heparin/PBS solution (10 U/mL heparin in PBS) using the infusion pump (14–831–200, Thermo Fisher Scientific) set to 4 mL/minute. Small pieces of tumor, lung, liver, and spleen were stored in 10% formalin for histopathology and the remainder of the tumors was immediately prepped for single-cell isolation using the gentleMACS Dissociators and mouse Tumor Dissociation Kit (130–096–730, Miltenyi Biotec) for 30 minutes. Red blood cell lysis was performed for 5 minutes (555899, BD Pharm Lyse) and after washing cells in PBS + 2% FBS, cells were frozen in 90% FBS + 10% DMSO for future experiments.

Live single cells that were frozen in 90% FBS + 10% DMSO described in the tissue harvest section were thawed and pooled to perform the Mouse Cytokine Antibody Array (ab133993, Abcam), according to the manufacturer's protocol. Single cells were lysed and 250 μg protein were used on each membrane. This assay only measured the intracellular cytokines as the excreted cytokines were washed off during the single cell prep. Membranes were developed by SuperSignal West Femto (34095, Thermo Fisher Scientific) and scanned on Amersham Imager 600 imagers (29083463, GE). Image was quantified by ImageJ.

Mouse tissues were collected and fixed in 10% formalin-neutral buffer. Paraffin embedding and hematoxylin and eosin (H&E) staining were performed at Baylor College of Medicine Breast Center Histopathology Core. The following antibodies were used: GFP (1:500, ab13970, Abcam), Ki67 (1:200, ab15580, Abcam), p-P27 (1:100, 71–7700, Zymed), Jagged-1 (1:100, ab7771), 2 different clones of S100A8 (1:5000, MAB3059, R&D Systems) used for mouse tissue and MRP8 (S100A8; 1:400, EPR3554, Abcam) used for human samples, CD4 (1:200, EPR19514, Abcam), CD8 (1:200, EPR21769, Abcam), F4–80 (1:200, CI:A3–1, Invitrogen), FoxP3 (1:200, FJK-16S, eBioscience), CD3 (1:400, A045229–2, Agilent), CD45 (1:200, 103102, BioLegend), and DAPI (R37606, Invitrogen). H&E and IHC images were scanned and quantified by Aperio ImageScope 11.2.0.780 (Aperio Technologies). IF images were scanned on an Axioscan and analyzed by using Zen Blue 3.1 (Zeiss). Human patient samples were purchased from Oregon Health and Science University (SR0494).

Live single cells from primary tumors that were frozen 90% FBS + 10% DMSO described in the tissue harvest section were thawed and washed with PBS + 2% FBS (FACS buffer) for flow cytometry analysis. Cells (2 × 10⁶) were blocked in Fc block (553142, BD Biosciences) and True-Stain Monocyte Blocker (426102, BioLegend) and stained for live (Live/Dead Aqua, L34957, Invitrogen) on ice for 30 minutes. After washing, cells were stained with antibodies on ice for 30 minutes in the dark: myeloid panel (Cd45, Cd11c, Cd11b, F4–80, Ly6c, Ly6g, PDCA-1, MHC2), lymphoid panel (CD45, CD3e, γσTCR, CD4, CD8α, CXCR3, granzyme B, PD-1, CD19, CD49b).

For the myeloid panel, cells were fixed in BD Cytofix (BDB554655) and saved in the FACS buffer. For the lymphoid panel, cells were stained for FoxP3 by using the Foxp3/Transcription Factor Staining Buffer Set (00–5523–00 and 00–8333, eBioscience). After staining, cells were kept in FACS buffer at 4°C to be analyzed on a BD LSRFortessa the next day. Data were analyzed by FlowJo (version 10.6.2).

Single-cell 3′ RNA sequencing was performed at the Single Cell Genomics Core at Baylor College of Medicine. After thawing frozen cells (recovery of about 40%–50% of cells), live single cells (DAPI- R37606, Thermo Fisher Scientific) were sorted on BD FACSAria I to obtain cells with the highest purity. A single-cell suspension was prepared at 1,000 cells per μL in DMEM + 10% FBS. The single-cell 3′ RNA-seq library was prepared according to the Chromium Single Cell 3′ Reagent Kit version 2 (120234 and 120235, 10x Genomics) for Fig. 2 and version 3 for Fig. 6. Briefly, single-cell suspension was loaded on a Chromium controller (10x Genomics) to generate single-cell GEMS (Gel Beads-In-Emulsions) where full-length cDNA was synthesized and barcoded. Subsequently the GEMs were broken and cDNAs from each single cell were pooled. Following cleanup using Dynabeads MyOne Silane Beads (370020, Thermo Fisher Scientific), cDNA was amplified by PCR. The amplified product was fragmented to optimal size before end-repair, A-tailing, and adaptor ligation. Finally, PCR was carried out to produce the library.

Sequencing was performed at Novogene using HiSeq (HiSeq 4000) and Novaseq (NovaSeq 6000, S4 lane). Sequencing strategy for the Novaseq S4 and HiSeq platform on these runs was PE150 (paired-end 150 bp) with an approximate data output of 800 Gb (roughly 2.6B reads per lane) and 110 Gb (roughly 360–400M reads per lane), respectively.

Raw sequencing reads were aligned to mouse mm10 reference genome following Cell Ranger v.3.1.0 pipeline. To ensure the quality, only the cells that express at least 200 genes and less than 10% mitochondrial genes and the genes expressed in less than three cells were retained for downstream analysis. In addition, suspected doublets with the cells expressing more than 4,000 genes (3,000 for α-IgG) were also filtered out. In each sample, cells were annotated by comparison to the reference datasets of pure mouse immune cells from the Immunologic Genome Project using SingleR package (19, 20). The GFP sequence was added to the reference gene (21). We noticed that some cancer cells that were annotated as “cancer cells”, lost their GFP expression. Therefore, for these analyses we used the population that was annotated as “cancer cells”, as opposed to GFP⁺ alone. The identity of the cells from closely related cell types, for example, T cells and natural killer (NK) T cells, and NK cells and ILC cells, was refined based on their expression profiles of canonical markers (22–24). Due to the low number, the cancer cells in the long-term DT were only included in a limited analysis.

Because of the lack of an apparent batch effect, the first four samples (primary, dormant, long-term DT, and recurrent) and the second two samples (α-IgG and α-Jag1–treated samples) were directly merged for downstream analysis, respectively, using Seurat v.3 (25). To unveil the underlying cell transition within both cancer cells and immune cells, we performed single-cell trajectory analysis for cancer cells, monocyte/macrophages, and T cells from the first four samples, respectively, using Monocle v. 2.8.0 (26). The high-dimensional data were reduced into a two-dimensional subspace using DDRTree algorithm, and the pseudotime and state were defined by orderCells function. Genes whose expression was significantly changed along the trajectory were identified using differentialGeneTest function with an q-value < 0.001. GO analysis was implemented on each cluster of differentially expressed genes (DEG) using clusterProfiler (27). GO terms were considered as significantly enriched if P_adj < 0.01.

Publicly available microarray data for patients (n = 1911) from eleven human primary breast cancer datasets were downloaded from NCBI GEO (GSE9893; GSE7390; GSE22226; GSE4922; GSE1378; GSE1379; GSE5327; GSE2740; GSE21653; GSE11121; GSE2990; GSE2034) along with the corresponding clinical annotations. Data were first log₂ transformed. Principal component analysis (PCA) was performed on the gene expression signature, and then k-means method was applied to cluster the samples into two categories. The association between sample category and relapse-free survival was estimated using the HR from Cox proportional hazards regression model.

Graphs were generated using GraphPad Prism 8 and R. To calculate recurrence, days were counted from the day when the treatment stopped to the day when recurrent tumors reached 5 mm in diameter and median survival was measured. The log-rank test was performed to measure the significance level for the difference in tumor recurrence between α-Jag1 and α-Jag1 treatment.

scRNA-seq gene expression data reported in this paper have been deposited in the Gene Expression Omnibus (GEO) database (www.ncbi.nlm.nih.gov/geo, accession no. GSE171464).

Previously, bigenic iFGFR transgenic mice developed in our laboratory and treated with a FGFR inhibitor provided a new model to study MRD and recurrence (18). In this syngeneic model, tumor regression induced tumor stroma remodeling and a decrease in myeloid-derived suppressor cells (MDSC). When tumors recurred EGFR signaling was upregulated and MDSCs increased (18). However, using this model it was not feasible to isolate and characterize the small number of residual cancer cells because they lacked a reporter gene. To facilitate their isolation and to begin to elucidate the mechanisms involved in MRD, we generated spontaneous GFP+ tumors from iFGFR^GFP mice. The GFP⁺ tumors were then transplanted into tolerized Glow-Head (GH) recipients (Fig. 1A; ref. 16). Briefly, as described previously (18), GFP⁺ tumor chunks were transplanted into the cleared fat pad of the GH mice and treated initially with a dimerizer to activate FGFR signaling until they were able to grow independently of the dimerizer. Mice bearing these large tumors then were treated with the FGFR inhibitor NVP-BGJ398 (BGJ; ref. 18) for 14 days, and subsequently drug treatment was stopped for another two weeks. We hypothesized that under these conditions the residual cells would not be affected by the immediate cytotoxic effect of BGJ treatment. At this timepoint, depending on the experiment, mice were either: (i) sacrificed for MRD studies (Fig. 1A, i and ii), or (ii) left untreated to monitor for recurrence. Tumor recurrence was measured once tumors reached 5 mm in diameter (Fig. 1B and C). Median survival was 90 days (1–4 months) consistent with previously published results using non-GFP-labeled cells (18). In three instances, there were no detectable recurrences and these samples were designated as long-term DT (Fig. 1B and C, in red). We characterized four stages of tumorigenesis: primary, MRD (dormant), long-term DT, and recurrent. Analysis of primary and recurrent tumors indicated that they had a similar histology, while the regressed lesions had increased stroma and immune cell infiltration (Fig. 1D). Interestingly, in the majority of the cases when tumors regressed, the epithelial cells formed duct-shape structures, which either lacked detectable or expressed low levels of Ki67 (Fig. 1E), but expressed a high level of phosphorylated-p27 (p-p27; Fig. 1F). GFP⁺ cells were undetectable by IHC in long-term DT due to very low number of remaining cells (Fig. 1E). Thus, we utilized this improved iFGFR^GFP tumor model to elucidate the properties of MRD that exhibited decreased proliferation and increased immune cell infiltration and stroma.

We next performed scRNA-seq analysis to investigate the transcriptional profiles of MRD (Supplementary Fig. S1A and S1B) and assess if the cells become dormant or are proliferating, as well as to identify any changes in their microenvironment. To ensure that we captured and characterized the dormant stage, we only selected lesions for further analysis that were negative for Ki67 in GFP⁺ tumor cells. Similar to IHC analysis of GFP (Fig. 1E), GFP⁺ cancer cells decreased when regressed and more dramatically in long-term DT (Fig. 2A). As indicated by the analysis of gene expression related to growth inhibition [MODULE_488, Gene Set Enrichment Analysis (GSEA)] and growth suppression (MODULE_407, GSEA; ref. 28), we observed that both MRD and long-term DT cells represented a quiescent state (Fig. 2B, i and ii). Furthermore, we also detected the elevated RNA expression of dormancy markers including Bhlhe41 (Dec2), Cdkn1b (p27), Zfp281, Nr2f1, Gas1, and Gas6 (29–33) in both MRD and/or long-term DT cells (Supplementary Fig. S1C). However, we did not detect any RNA expression of β-galactosidase, a senescence marker, in either the MRD or long-term DT cells. Moreover, senescence is by definition irreversible, leading to cell death and clearance by innate immune cells, whereas dormancy is reversible and can lead to recurrence (34). In addition, dormant cells activate variable degrees of pluripotency programs (34). GFP⁺ cells are present in all regressed lesions and eventually recur within 1–4 months, which argues against senescence, that is, a reversible process. In addition, comparing dormant to long-term DT, we found that these two states express different markers of quiescence and dormancy, suggesting that these two stages were a result of dormancy, but not senescence, hence these were designated as dormant.

We next applied trajectory analysis to determine the different status between dormant cells and the primary and recurrent cells (Fig. 2C). Pseudotemporal gene expression analysis identified four clusters of genes that were differentially expressed (Fig. 2D). Cluster1, which consists of the genes more highly expressed at the early stage of trajectory, was enriched in the dormant cells (Fig. 2E). Cluster3 and 4, which were highly expressed at the far end of the trajectory, were enriched in the primary and recurrent cells (Fig. 2E). Gene Ontology (GO) enrichment analysis showed that Cluster1 genes are involved in remodeling of stroma and immune cells and recruitment of both myeloid and lymphoid cells (Fig. 2F). Genes from Cluster2 predominantly function in DNA replication and the cell cycle, while Clusters 3 and 4 are enriched for innate immunity and response to virus, and metabolic pathways, respectively (Supplementary Table S1).

We also noticed that regulation of the Notch, Hippo, and Wnt signaling pathways are significantly enriched in Cluster1 (Fig. 2D and F) but we focused on Notch signaling pathway. Previously, we showed that MDSCs increase when tumors recur (18) and that MDSCs directly regulate tumor-initiating cells through regulating the Notch pathway (35). In addition, tumors can induce Jagged ligands in MDSCs through NFκB-p65 (36). Inhibiting MDSCs by using an anti-Jagged immunotherapy overcomes T-cell tolerance and enhances T-cell immunotherapy (36). Interestingly, we found that during dormancy, Jag1 not only participates in the Notch signaling pathway, but also regulates a variety of other biological processes including hemopoiesis, myeloid cell differentiation, stem cell (SC) development, mesenchymal cell differentiation, and epithelial-to-mesenchymal transition (EMT; Fig. 2F, annotated by #). Both Jagged-1 and -2 were expressed in dormant and long-term dormant (Fig. 2G). Moreover, different Notch signaling genes were found to be involved in the regulation of the four stages of tumorigenesis (Fig. 2G). For example, Notch2 is most likely involved in the regulation of primary and recurrent cells, while Notch3 and 4 may regulate dormant cells, and Notch1 may regulate long-term DT cells. The different subunits of the γ-secretase complex, presenilin-1 (Psen1), nicastrin (Ncstn), anterior pharynx-defective 1 (APH-1), and presenilin enhancer 2 (Psenen), were also shown to be differentially expressed at these different stages (Fig. 2G). Altogether, these results suggest that when tumor cells become quiescent, therefore dormant, they have enhanced Notch signaling.

Murine mammary gland (MG) stem/progenitor cells have been identified previously at different stages of mammary gland development. These include fetal mammary SCs that peak in late embryogenesis and their subsequent differentiation to adult basal and luminal SCs (37). Because the dormant cells were quiescent, we assessed any similarities to these different classes of MG stem/progenitor cells. Interestingly, primary and recurrent- compared to dormant and long-term dormant- cells were found to be enriched for the features of luminal progenitors (Supplementary Fig. S1D, i), whereas dormant cells were enriched for the signatures observed in fetal and basal progenitors as well as mature luminal cells (Supplementary Fig. S1D, ii, iii, and iv). Moreover, to further investigate response differences among the subgroups of cancer cells to the BGJ treatment, we performed unsupervised clustering on cancer cells and identified eight clusters in the combined samples (Supplementary Fig. S1E, CC1–8). There was a substantial change in the composition of the different clusters among the primary, dormant and recurrent stages (Supplementary Fig. S1F). Primary tumors mainly consist of CC2, 3, 4, and 7, which were related to metabolic pathways (2), remodeling of stroma and T-cell activation (3), DNA replication and cell cycle (4), and regulation of both myeloid and lymphoid cells (Supplementary Table S2; ref. 7). The major clusters observed in dormant lesions were CC3, 4, 5, 6, and 8, where CC6 and 8 seemed to appear de novo and did not exist in the primary tumor or were presumably present at an extremely low abundance. The GO analysis showed that CC5 and 8 participate in recruitment of immune cells both myeloid and lymphoid, and gland regeneration, respectively (Supplementary Table S2). CC6 was of a great interest that it exhibited the most diverse functions including remodeling of stroma and immune cells and recruitment of both myeloid and lymphoid cells, metabolic pathways, DNA and RNA biosynthetic processes, regulation of SC differentiation, and gland regeneration (Supplementary Table S2). In addition, Jag1 regulated many pathways in this cluster. Recurrent tumors primarily consisted of CC1, 2, 4, and 7, where the feature genes in CC1 are enriched for regulation of myeloid cells and gland development (Supplementary Table S2).

These results suggest that there either may be rare cell subpopulations or perhaps a transition into a de novo population that becomes the major repopulating cluster after cells become dormant. This is contrary to the concept that all the clusters regress identically with treatment, resulting in a similar probability of recurrence.

To determine whether the gene signatures of the dormant stage were associated with human breast residual disease, we compared the expression profiles of dormant-overexpressed genes in residual breast cancer following neoadjuvant targeted therapy from two public microarray datasets [GSE10281 (38) and GSE21974 (39)] to their corresponding pretreated samples. We found that signature genes of residual cells were more highly expressed in patient tumors with residual disease (Supplementary Fig. S2A and S2B). For example, CCL5, which can promote breast cancer recurrence through macrophage recruitment in residual tumors (14) was shown to be upregulated in the dormant stage as compared with primary and recurrent cells, and was also more highly expressed in the residual tumor samples. However, we did not find any association of the dormant gene signature with recurrence-free survival in the primary tumors of patients with breast cancer (Supplementary Fig. S2C).

Because MDSCs increase during recurrence and may regulate Notch signaling pathway through Jagged-1, we sought to identify Jagged-1 expressing MDSCs in patient samples. We performed Immunofluorescent (IF) staining of Jag-1 and S100A8 (a MDSC marker) in infiltrating ductal carcinoma (IDC) and adjacent areas with fibrocystic changes (FCC) and found that in both sites, MDSCs express Jag-1(Supplementary Fig. S2D). In addition, some cancer cells in IDC also express Jag-1 (Supplementary Fig. S2D). Interestingly, we found S100A8 expression in the cluster4, which consists of primarily recurrent cells, as well as some primary and dormant cells (Fig. 2D).

Because this mouse model has an intact immune system, we also were able to investigate the transcriptional changes in the stroma and immune cells during dormancy (Fig. 3A and B). Both primary and recurrent tumors have very similar microenvironments with a small decrease in T cells and a minor increase in fibroblasts. In contrast, the dormant cells showed substantial differences in the abundance of stroma and immune cell populations, including an increase in T cells, NK cells, B cells, monocytes, dendritic cells (DC), fibroblasts, and endothelial cells and a decrease in macrophages and neutrophils (Fig. 3A). In long-term DT cells, however, we found a decrease in T cells, NK cells, and B cells, as compared with dormant cells, and an increase in DCs, monocytes, macrophages, stromal cells, and fibroblasts (Fig. 3A), suggesting that the microenvironment during progression from dormant to long-term DT might have become more immunosuppressive. These results obtained from scRNA-seq were confirmed by immunofluorescence (IF) detection of some of the known markers for these immune cells (Fig. 3C and D). Interestingly, in the dormant stage, we found small aggregates of T cells (CD4⁺ and CD8⁺; Fig. 3C and D, dotted line circles). Consistent with the scRNA-seq data, we found macrophages (F4/80) decreased in the dormant stage as compared with other stages (Fig. 3D). To validate these results, we performed flow cytometry analysis of six frozen primary tumors using canonical markers for different immune cells grouped into myeloid or lymphoid panels (Supplementary Fig. S3A and S3B). Using functional markers, we identified a population of T regulatory cells (Tregs), CD4 T cells that are not Th1 activated CXCR3⁺, and CD8 T cells that are not cytotoxic (granzyme B+; Supplementary Fig. S3C; refs. 40, 41). Moreover, we measured the cytokines expressed in cells (intracellular) from different samples and found that CCL5 and IL2 were expressed in the dormant stage (Supplementary Fig. S3D), again suggesting an increased recruitment and proliferation of immune cells (42, 43).

In addition, we identified a population of CD34⁺ cells (Fig. 3B), which exhibits transcriptional signatures primarily of fibroblasts and endothelial cells. CD34⁺ fibroblast-like stromal cells can be recruited to the tumor-reactive stroma and are involved in wound healing, antigen presenting, and the secretion of various types of collagen and cytokines (44–46). These cells also had increased expression of Jag1 as compared with other populations (Fig. 3B). Together, these results show that numerous immune and stroma cells are being remodeled during dormancy and various cytokines are expressed at this stage that cause recruitment and proliferation of immune cells.

To obtain a better understanding of the functions of myeloid cells during dormancy, we first performed a trajectory analysis on monocytes and macrophages. This identified five distinct states of macrophages (MC1–5), where MC1 and MC5 are more prominent in the dormant and long-term DT samples (Fig. 4A). Interestingly, we observed that monocytes are more prominent in the dormant and long-term DT stages (Fig. 4B), suggesting an increase of monocyte recruitment. Macrophages display different polarization phenotypes often referring to M1 and M2, where M2 is immune-suppressive (47). To characterize the functions of the five distinct states of macrophages, we profiled the expression of the feature genes of M1, M2, breast cancer–associated macrophages (Br-TAM), and tumor-educated monocytes (TEMo; Fig. 4C–F; ref. 48). In general, a mixed phenotype was observed in all macrophage states. Dormant and long-term DT-prominent MC1 grouped together with monocytes at the far left of the trajectory (the early stage of pseudotime), suggesting that MC1 may be less differentiated. In addition, MC1 showed increased expression of genes characteristic of M2 and TEMo (Fig. 4E and F) but not Br-TAM (Fig. 4C). This suggests that dormant and long-term DT might maintain a higher level of immunosuppressive innate immunity. Comparatively, the primary- and recurrent-prominent macrophages, MC3 and MC4, showed a Br-TAM phenotype, and were more characteristic of M1 than M2 macrophages (Fig. 4C–E). However, although they express the Br-TAM gene signature, they did not express the two suggested TAM markers SIGLEC1 and CCL8 (48), but instead expressed CCL3 and CCL4- known to regulate TAMs (Fig. 4C; refs. 49, 50). The MC5 state was the only dormant- and long-term DT-prominent population that showed Br-TAM features and expressed SIGLEC1 and CCL8 (Fig. 4C). MC2 showed some expression of all phenotypes but more of M2, TEMo, and Br-TAM, in that order (Fig. 4C–F), suggesting a transient state. Furthermore, monocytes did not display TEMo features (Fig. 4F). Thus, these results demonstrate that there is a spectrum of macrophage and monocyte phenotypes rather than canonical M1-M2 polarization during these four stages of tumorigenesis.

MDSCs are a heterogenous population of neutrophil- and monocyte-like myeloid cells that play immunosuppressive roles (denoted as G- and M-MDSCs, respectively) (51). MDSCs promote tumor progression by various mechanisms including suppressing immune cells, supporting tumor-initiating cells (TIC), and enhancing angiogenesis and metastasis (52). G-MDSCs characterized in breast cancer models have a distinct gene signature as compared to normal neutrophils (51). In our studies, we found that neutrophils display high expression of the signature genes of breast cancer–associated MDSCs (Br-MDSCs; combining the G- and M-MDSC signature genes; Supplementary Fig. S4A). In contrast, monocytes did not express many Br-MDSC genes (Supplementary Fig. S4A). Furthermore, CD84, a suggested G-MDSCs marker (51), was not expressed in these neutrophils, but was detected on some macrophage states (Supplementary Fig. S4A).

Our previous studies showed that MDSCs decreased during dormancy and increased during recurrence (18) and that MDSCs regulated TICs in part through activating Notch (35). Here, we found a similar expression of tumor resident neutrophils as compared with Br-MDSCs (Supplementary Fig. S4A). In addition, we demonstrated that MDSCs expressed Jag1 that increased during long-term DT and recurrence (Supplementary Fig. S4B). Comparing the expression of Jag1 in circulating monocytes and neutrophils from primary and recurrent to non–tumor bearing mice showed that both monocytes and neutrophils in recurrent tumors expressed Jag1 (Supplementary Fig. S4C).

Utilizing IF staining and scRNA-seq, we also observed an increase of T cells during dormancy (Fig. 3A, C, and D). Trajectory analysis for T cells revealed that T cells are elevated in the dormant site as compared to primary, long-term DT, and recurrent stages (Fig. 5A). Moreover, the composition of the subclusters of CD4 and CD8 T cells changed during each stage of tumorigenesis (Fig. 5B and C). Specifically, Tregs decreased during dormancy and long-term DT as compared with primary and recurrent tumors (Fig. 5C), while total CD8 T cells increased (Fig. 5C). According to the expression profiles of canonical markers (53), the majority of these CD8 T cells (CD8+_1) were not cytotoxic (granzyme A⁻/B⁻ and perforin⁻; Fig. 5D). The CD8+_2 population, which was decreased during dormancy, showed the most activated and nonexhausted phenotype (Fig. 5C and D). The CD8+_3 population that was increased during dormancy resembled exhausted T cells because of the expression of Lag3, Cd244a, Cd274 (PD-L1), Pdcd1 (PD-1), and Cd160 (Fig. 5C and D; refs. 53, 54). Similarly, we observed that CD4 T cells also exhibit potentially different functions among samples (Fig. 5C and D). CD4 T cells during dormancy (CD4+_1) have Th1 activity, whereas CD4 T cells in primary tumors (CD4+_2) exhibit Th2 activity with increased expression of Lag3, Havcr2 (Tim3), Cd244a, Cd274, and Pdcd1 (Fig. 5D). The majority of T cells in long-term DT are CD4⁺/C8⁺ T cells (Fig. 5D). These results suggest that overall T-cell functions may favor a more antitumorigenic state in dormant samples because of a decrease in Tregs and Th2 cells and an increase in CD8 T cells and Th1 cells. However, as most of CD8 T cells appear not to be not fully activated/cytotoxic (granzyme A⁻/B⁻ and perforin⁻) or are already exhausted, cancer cells may have escaped from T cell immunity.

Innate lymphoid cells (ILC) are innate counterparts of T cells and are classified into ILC1, ILC2, ILC3, and NK cells. However, they have not been well characterized in cancers. Using a previously published gene lists for ILCs (22), we found various clusters of ILCs that resembled ILC1, 2, 3, and NK cells (Fig. 5E). Although the total abundance of ILCs were similar, the proportion of these different types of ILCs changed among the four samples (Fig. 5F). NK cells increased during dormancy compared with the primary, long-term DT, and recurrent stages, while ILC1 decreased during both dormancy and long-term DT (Fig. 5F). ILC2 decreased during dormancy but drastically increased during long-term DT and ILC3 increased during both dormancy and long-term DT as compared with primary and recurrent (Fig. 5F). Increased NK cells and decreased ILC2 during dormancy suggest an antitumor function of ILCs. The function of ILC1 and ILC3 remain unclear because of their dual roles in tumors.

The scRNA-seq data suggest a role of Jagged-Notch signaling pathway in dormancy. Notch signaling has been shown to accelerate tumor recurrence (13) and regulate endothelial–dormant cell interaction (8). Notch signaling is also implicated in MDSC-mediated regulation of TICs (35). Given that MDSCs increase during long-term DT and recurrence, we next investigated the effect of an anti-Jag1 (α-Jag1) blockade on tumor dormancy. We induced tumor dormancy using a protocol identical to the cohort described previously (Fig. 6A). Two weeks after BGJ treatment was stopped, mice were divided into α-Jag1 or control (α-IgG) groups. Mice were treated with α-Jag1 once a week for a total of five weeks, and then treatment was stopped. Recurrence was reported once tumors reached 5 mm in diameter (Fig. 6A). Compared with the control group, the α-Jag1 blockade markedly delayed tumor recurrence (median survival = 222 days in α-Jag1 vs. 52.5 days in α-IgG; Fig. 6B).

We repeated this experiment in a third cohort, collected samples one week after the cessation of the treatment (Fig. 6A), and performed scRNA-seq on α-IgG and α-Jag1–treated samples (Supplementary Fig. S5A–S5C). Similar to the first cohort, we also observed substantial changes in immune and stroma cells in α-Jag1 and α-IgG–treated samples (Fig. 6C). Specifically, macrophages, NK cells, T cells, and γσT cells increased in the α-Jag1–treated group compared with α-IgG, while monocytes, neutrophils, and B cells decreased (Fig. 6C). The expression of some of these immune cells was validated by IF staining of tissues for canonical markers of macrophages, MDSCs, CD4 or CD8 T cells, and Tregs in α-Jag1- and α-IgG–treated samples (Fig. 6D).

Additionally, we identified three subclusters of CD4 and CD8 T cells and a cluster of Treg cells in both sample types (Supplementary Fig. S5D). Treg cells were found to be slightly decreased while CD4 and CD8 T cells were increased overall in the α-Jag1 as compared with the α-IgG treatment group (Supplementary Fig. S5D). Two of the CD8 subclusters (CD8+_2 and CD8+_3) showed a cytotoxic phenotype (granzyme B⁺), with CD8+_3 showing some exhaustion features (Supplementary Fig. S5E). The proportion of CD4+_2, which has a Th1 phenotype, slightly increased following α-Jag1 treatment (Supplementary Fig. S5D and S5E).

Altogether, these results show that the effect of α-Jag1 treatment on delaying recurrence is most likely due to a pleotropic effect on immune cells and possibly cancer and endothelial cells. α-Jag1 treatment is associated with the increase of NK cells, macrophages, and CD4 and CD8 T cells and the decrease of monocytes, neutrophils, B cells, and Tregs.

This is the first time, to our knowledge, that a comprehensive characterization of stroma and immune cells has been analyzed during dormant MRD as summarized in Fig. 7. To understand the complex ecosystem that promotes dormant cell survival and remodeling of microenvironment, we utilized a novel transgenic model, iFGFR1^GFP with an intact immune system, to study MRD and relapse following treatment with an FGFR inhibitor. We demonstrated that Ki67⁻ MRD have gene expression profiles similar to growth-arrested cells with dormancy markers including Cdkn1b, Zfp281, Nr2f1, and Gas6. Dormancy is accompanied by remodeling of stroma including endothelial cells and fibroblasts and immune cells including recruitment of both myeloid and lymphoid cells, the majority of which show immunosuppressive activity or are not fully functional and cytotoxic (Fig. 7).

Interestingly, different groups of MG progenitor cells were detected during dormancy and recurrence perhaps as a function of either de novo reprogramming or by transition from other cell types. Thus, it does not appear that during dormancy the FGFR inhibitor treatment results only in tumor de-bulking and that similar cell populations exist throughout these different stages of tumorigenesis. Given that different MG progenitors appear during dormancy, we hypothesize that the CC6 cells identified by scRNA analyses may be the major group of progenitor cells that give rise to the recurrent tumors, but this will have to be established by lineage tracing experiments. Interestingly, a recent study has also demonstrated that dormant residual tumor cells express a mammary SC signature, which is only capable of forming recurrent tumors but not primary tumors, suggesting that distinct TICs may give rise to primary or recurrent tumors (55).

Jag1 has been shown to regulate many aspects of tumorigenesis including SC development, EMT, and immune cells homeostasis during dormancy and may be involved in the CC6 cells. The pleiotropic effects of Jag1 expression in cancer cells, fibroblasts, endothelial cells, and myeloid cells, may all account for the marked effect of α-Jag1 treatment to increase dormancy and delay recurrence.

This study focused on some myeloid cell types including monocytes, macrophages, neutrophils and NK cells. The increase of immune-suppressive myeloid cells may explain the lack of T-cell activity. Although we saw an increase of neutrophils (MDSCs) during long-term DT by IHC analysis, there was a discrepancy when we analyzed scRNA-seq. This can be explained by the fact that MDSCs decrease viability during cryopreservation (56, 57).

Although we saw an increase of T cells and the appearance of T-cell aggregates during dormancy, which appear to resemble a TLS (58), most of CD8 T cells seem not fully activated/cytotoxic (granzyme A⁻/B⁻ and perforin⁻) or are already exhausted. Cytotoxic NK cells also increased but dormant cancer cells may eventually find an opportunity to evade immune surveillance including NK cells (7) and T-cell immunity as suggested by the decrease of NK and CD8 T cells in long-term DT. Interestingly, a distinct population of CD8 T cells (CD39⁺PD-1⁺), were discovered to inhibit the outgrowth of lung mets and mediate the dormancy of DTCs in a mouse model but also was unable to fully eradicate all of the cancer cells due to expression of checkpoint molecules (59).

We treated mice once weekly for only five weeks with an α-Jag1 blockade antibody during dormancy, but notwithstanding we observed a significant sustained delay of tumor relapse. These results suggest that this brief window of treatment during dormancy was sufficient both for the survival of dormant cancer cells, and also to inhibit the changes observed in both the stroma and immune cells. Treatment with α-Jag1 had a pleotropic effect on immune cells. It induced an antitumor immunity by increasing macrophages, NK cells, and CD8 T cells (granzyme B⁺ and perforin⁺) and decreased myeloid cells including monocytes and neutrophils, as well as B cells and Tregs. The GO analysis of residual tumor cells suggests that there are other potential targets including those of Notch pathway signaling (Jagged-2, DLL1, DLL4, and etc.) and Hippo, and Wnt signaling pathways that can be evaluated to target dormant cells in the future.

Overall, our study demonstrates that there is significant remodeling of both the stroma and immune cells during dormant MRD. The ongoing immune surveillance process recruits a mixture of immune-active (T cells and NK cells) and suppressive (Tregs and myeloid cells) cells. In long-term dormancy however, cancer cells can evade the immune surveillance and the immune-suppressive cells become the prominent population. Cancer (stem) cells can also transition into a different state. The Notch signaling pathway plays an important role in regulating both the cancer and immune cells, and targeting this pathway provides a promising therapeutic strategy to inhibit tumor relapse. In future studies it will be important to determine if similar mechanisms are employed in the metastatic microenvironment. In addition, we did not find any association between our dormant gene signature and survival of patients with breast cancer. One caveat is that gene expression data obtained from all of these studies are from tumor chunks as opposed to single-cell analysis so there is no resolution to specifically assess the dormant signature in cancer cells. The other caveat is that gene expression data are only from primary tumors with their survival metadata, and not regressed tumors. As our study suggests, the CC6 is a rare or de novo population that appears during dormancy with the potential to contribute to tumor relapse. Without the availability of single-cell and spatial transcriptomic data from patient samples to specifically identify these populations during MRD, we are not able to interrogate patient survival by only investigating primary samples. Future studies using single-cell analysis of both human primary and MRD will be required to address the relevance of these mouse tumor models.

No disclosures were reported.

M. Janghorban: Conceptualization, resources, data curation, software, formal analysis, supervision, validation, investigation, visualization, methodology, writing–original draft, project administration, writing–review and editing. Y. Yang: Data curation, software, formal analysis, supervision, validation, visualization, writing–original draft, writing–review and editing. N. Zhao: Methodology. C. Hamor: Methodology. T.M. Nguyen: Methodology. X.H.-F. Zhang: Conceptualization, resources, data curation, software, formal analysis, supervision, visualization, writing–original draft, writing–review and editing. J.M. Rosen: Conceptualization, resources, supervision, funding acquisition, writing–original draft, project administration, writing–review and editing.

The authors thank Dr. Chris Siebel for kindly providing anti-Jagged-1 blockade antibody and helpful suggestions, Dr. Chi-Ping Day (NIH/NCI) for providing the Glow-Head mice and Dr. Michael Lewis for providing the GFP mouse, Dr. Yumei Li for prepping single cell libraries, Dr. Brian Gibson for training in mouse surgery techniques and protocol preparation, Dr. Alejandro Contreras for providing access to Aperio slide scanner at the Department of Translational Molecular Pathology-MD Anderson, Dr. Susan Hilsenbeck for helping with statistical power analysis, Dr. Rainer Lanz for helping with some data analysis, Dr. Brian Ruffell for helping with the flow cytometry, and Dr. Lisa Coussens for advice about some data analysis. The authors acknowledge the joint participation of the Diana and Adrienne Helis Malvin Medical Research Foundation and Baylor College of Medicine. This project was supported by grant CA016303 (to J.M. Rosen), the Quantitative Science Shared Resource in the Duncan Comprehensive Cancer Center (CA125123), the Cytometry and Cell Sorting Core at Baylor College of Medicine (BCM) with funding from the CPRIT Core Facility Support Award (CPRIT-RP180672), the NIH (CA125123, RR024574, S10 OD025251) and the assistance of Joel M. Sederstrom, the RNA In Situ Hybridization Core at BCM, which is, in part, supported by a Shared Instrumentation grant from the NIH (1S10OD016167), and the assistance of Dr. Cecilia Ljungberg. Single Cell 3′ RNAseq was performed at the Single Cell Genomics Core at BCM, which is partially supported by the NIH shared instrument grants (S10OD023469, S10OD025240) and P30EY002520 (to R. Chen). This work benefitted from data assembled by the ImmGen consortium (20).

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.