Parity-induced changes to mammary epithelial cells control NKT cell expansion and mammary oncogenesis.

Pregnancy reprograms mammary epithelial cells (MECs) to control their responses to pregnancy hormone re-exposure and carcinoma progression. However, the influence of pregnancy on the mammary microenvironment is less clear. Here, we used single-cell RNA sequencing to profile the composition of epithelial and non-epithelial cells in mammary tissue from nulliparous and parous female mice. Our analysis indicates an expansion of γδ natural killer T-like immune cells (NKTs) following pregnancy and upregulation of immune signaling molecules in post-pregnancy MECs. We show that expansion of NKTs following pregnancy is due to elevated expression of the antigen-presenting molecule CD1d on MECs. Loss of CD1d expression on post-pregnancy MECs, or overall lack of activated NKTs, results in mammary oncogenesis. Collectively, our findings illustrate how pregnancy-induced changes modulate the communication between MECs and the immune microenvironment and establish a causal link between pregnancy, the immune microenvironment, and mammary oncogenesis.

INTRODUCTION

Changes to the functions of immune cells modulate both the mammary immune microenvironment and mammary epithelial cell (MEC) lineages during all stages of mammary development, with CD4+ T cells guiding lineage commitment and differentiation of MECs, while macrophages provide growth factors and assist in removal of cellular debris from apoptotic events (Dawson et al., 2020; Hitchcock et al., 2020; Plaks et al., 2015; Rahat et al., 2016; Stewart et al., 2019; Wang et al., 2020). Accordingly, changes that impact immune cell function and abundance can also influence the development and progression of mammary oncogenesis, particularly in tissue reconstruction during postpartum involution (Bach et al., 2021; Ibrahim et al., 2020; Lyons et al., 2011; Martinson et al., 2015; Freire-de-Lima et al., 2006; Guo et al., 2017; Fornetti et al., 2012; O’Brien et al., 2010).

Conversely, cell-autonomous processes in post-pregnancy MECs contribute to a lasting effect that decreases the risk of breast cancer by ~30% in rodents and humans (Medina et al., 2004; Britt et al., 2007; Terry et al., 2018). Epigenetic-mediated alterations of post-pregnant MECs have been shown to suppress mammary oncogenesis via oncogene-induced senescence (Feigman et al., 2020). Given that oncogene-induced senescence signals influence the immune system, a link between normal pregnancy-induced mammary development, the immune microenvironment, and oncogenesis needs to be addressed to fully understand the effects of pregnancy on breast cancer development.

In this study, we characterize the interactions between cell-autonomous (MECs) and non-cell-autonomous (immune cells) factors that are part of normal pregnancy-induced mammary development and are involved in inhibiting mammary oncogenesis. Our analysis identified that pregnancy induces the expansion of natural killer T-like cells (NKT) during the late stages of involution, which preferentially populates the fully involuted, post-pregnancy mammary tissue. Unlike typical NKTs that bear αβ T cell receptors (TCRs), pregnancy-induced NKTs express γδTCRs on their surface, indicating a role in specialized antigen recognition. NKT cell expansion was linked with increased expression of the antigen-presenting molecule, CD1d, on the surface of post-pregnancy MECs, which was associated with the stable gain of active transcription marks at the Cd1d locus and increased mRNA levels. Further analysis demonstrated that gain of CD1d expression on post-pregnancy MECs, and expansion of γδNKTs was observed in tissues that failed to undergo mammary oncogenesis in response to oncogenic signals, such as cMyc overexpression or Brca1 loss of function. Altogether, our findings elucidate how signals brought to MECs during pregnancy-induced development regulate epigenomic changes, gene expression, and immune surveillance, which together control mammary oncogenesis.

RESULTS

Identification of transcriptional programs and immune cellular heterogeneity in mammary tissue from parous female mice

The use of single cell RNA sequencing (scRNA-seq) has elucidated the dynamics of epithelial cell-lineage specification and differentiation across major mammary developmental stages (Bach et al., 2017; Chung et al., 2019; Li et al., 2020a; Pal et al., 2017, 2021). Previous studies have indicated that post-pregnancy epithelial cells bear an altered transcriptome and epigenome, thus suggesting that pregnancy stably alters the molecular state of MECs (Blakely et al., 2006; Feigman et al., 2020; Huh et al., 2015; dos Santos et al., 2015). However, it is unclear whether pregnancy leads to disproportionate changes in the transcriptome of specific mammary cell populations.

In order to characterize the effects of parity on the cellular composition and heterogeneity of mammary glands, we used scRNA-seq to compare the abundance, identity, and gene expression of mammary epithelial and non-epithelial cells from nulliparous (virgin, never pregnant) and parous (20 days gestation, 21 days lactation, 40 days post-weaning) female mice. scRNA-seq clustering defined 20 clusters (TCs), which were further classified into five main cell types: epithelial cells (Krt8+ and Krt5+), B lymphocytes (CD20+), and T lymphocytes (CD3e+) and two smaller clusters, encompassing fibroblast-like cells (Rsg5+) and myeloid-like cells (Itgax+), with similar cell-cycle states (Figures S1A-S1C).

To characterize the cellular heterogeneity across pre- and post-pregnancy MECs, we used a re-clustering approach that resolved 11 clusters of mammary epithelial cells (ECs) (Henry et al., 2021) (Figure 1A). Analysis of cellular abundance and lineage identity revealed that clusters EC7 (mature myoepithelial MEC), EC9 (luminal common progenitor-like MEC), EC10, and EC11 (bi-potential-like MECs) were evenly represented in pre- and post-pregnancy mammary tissue, thus demonstrating populations of cells that are mostly unchanged by a pregnancy cycle. We also identified clusters predominantly represented in pre-pregnancy mammary tissue (EC2, EC4, and EC8), and those biased toward a post-pregnancy state (EC1, EC3, EC5, and EC6), classified as luminal alveolar-like clusters (EC1, EC2, and EC6), myoepithelial progenitor-like clusters (EC3 and EC4), and luminal ductal-like clusters (EC5 and EC8) (Figures S1D-S1F). Comparative gene-expression analysis indicated that processes associated with immune cell communication were markedly enriched in luminal and myoepithelial cell clusters biased toward the post-pregnancy state (Figure 1B; Figures S1G and S1H; Table S1). This observation was supported by analysis of previously published pre- and post-pregnancy bulk RNA-seq data, which suggested an overall enrichment for immune communication signatures in epithelial cells after a full pregnancy cycle (Feigman et al., 2020) (Figure S1I; Table S2).

Figure 1.

Identification of transcriptional programs and immune cellular heterogeneity in mammary tissue from parous female mice

(A) UMAP of mammary epithelial cells from pre- and post-pregnancy mammary glands.

(B) mRNA levels of senescence-associated, immune communication genes Cxcl1, Ccl2, Il6, Cxcl5, Mhc-ii, and Cd1d in pre- and post-pregnancy MECs.

(C) UMAP of T cells (CD3e+ cells) from pre- and post-pregnancy mammary glands.

(D) Feature plots showing the expression of T cell markers Cd4, Cd8, Klrk1, and Gzma.

(E) Dendrogram clustering and dot plot showing the molecular signature and lineage identity of pre- and post-pregnancy mammary resident CD3+ immune cells.

See also Figures S1-S4.

Changes in the immune microenvironment are known to contribute to pregnancy-induced mammary development and cancer (Coussens and Pollard, 2011; Bach et al., 2021; Dawson et al., 2020; Saeki et al., 2021). Therefore, and in light of the potentially altered epithelial-immune cell communication identitied in post-pregnancy MECs suggested above, we set out to understand the effects of pregnancy on the mammary resident immune compartment using scRNA-seq. Transcriptional analysis of clusters representing B lymphocytes (CD20+) did not identify major differences between cells from pre- and post-pregnancy mammary glands, suggesting that B cells may not be significantly altered in fully involuted mammary tissue (Figure S2A). Re-clustering of CD3e+ T lymphocytes identified nine distinct immune cell clusters (IC) marked by the expression of immune lineage genes such as Cd4, Cd8, Klrk1, and Gzma (Figures 1C and 1D). Classification according to cell abundance and lineage identity of mammary resident lymphocytes revealed two cell clusters, IC1 (CD4+ memory-like T cells) and IC2 (CD8+ T cells), which were evenly represented across pre- and post-pregnancy mammary tissue (Figures S2B and S2C). Differential gene-expression analysis of clusters IC1 and IC2 identified minimal expression changes, suggesting that the transcriptional output of CD8+ T cells (IC2), and certain populations of CD4+ T cells (IC1) were not substantially altered by parity (Figures S2D and S2E).

Analysis of clusters biased toward a pre-pregnancy state identified several populations of CD4+ T lymphocytes, with gene identifiers supporting their identity as CD4+ Tregs (IC3), CD4+ naive T cells (IC7 and IC8), and CD4+ helper T cells (IC4), suggesting that pre-pregnancy mammary tissues are enriched for populations of CD4+ T cells (Figure 1E). Conversely, clusters enriched with post-pregnancy mammary immune cells (IC5, IC6, and IC9) were classified as NKT cells, a specialized population of T cells involved in immune recruitment and cytotoxic activity (Godfrey et al., 2004) (Figure 1E). These clusters expressed master regulators of NKT cell fate, including transcription factors (TFs) Tbx21 (Tbet) and Zbtb16 (Plzf) (Townsend et al., 2004; Savage et al., 2008).

While natural killer (NK) cells are known to play a role in mammary gland involution and parity-associated mammary tumorigenesis (Fornetti et al., 2012; Martinson et al., 2015), the role of NKT cells in this process has yet to be determined. Therefore, we analyzed clusters of immune cells expressing the common NK/NKT marker Nkg7 to further define the influence of pregnancy on the abundance and identity of NK and NKT cells. Deep-clustering analysis of Nkg7+ immune cells revealed six distinct cell clusters (NC1–6). Cells classified under cluster NC5, which includes cells from both the pre- and post-pregnancy mammary tissue, lacked expression of Cd3e and therefore represents the only cluster with an NK cell identity in our dataset (Figures S2F-S2H). Further gene-expression analysis confirmed that post-pregnancy mammary glands are enriched with a variety of NKTs, including those expressing markers of cell activation (Gzmb and Ccr5) and of a resting state (Bcl11b) (Figure S2H). In agreement, each of the post-pregnancy-biased NKT cell clusters was enriched with an array of immune-activation signatures, suggesting an altered state for these cell populations after pregnancy (Figure S2I).

Collectively, our scRNA-seq analysis of fully involuted mammary tissue confirmed that pregnancy leads to a stable alteration of the transcriptional output of post-pregnancy MECs, including gene-expression signatures that suggest enhanced communication with the immune microenvironment. In addition, our study also indicates that mammary resident NKTs are present at higher levels in post-pregnancy glands, suggesting that pregnancy plays a role in inducing changes to the mammary immune microenvironment.

Pregnancy induces the expansion of a specific population of NKTs

The post-partum mammary gland involution is marked by an influx of infiltrating mast cells, macrophages, neutrophils, dendritic cells, and natural killer cells, which remove apoptotic epithelial cells and support the remodeling of the gland (Guo et al., 2017; Kordon and Coso, 2017; O’Brien et al., 2010; Schwertfeger et al., 2001). Since our scRNA-seq analyses suggested that post-pregnancy mammary glands are enriched for populations of NKT cells, we next used a series of flow cytometry analyses to validate this observation.

Analysis using the markers NK1.1 and CD3, which defines NKTs (NK1.1+CD3+), identified a 12-fold increase in the abundance of NKTs in post-pregnancy mammary tissue, consistent with the results of our scRNA-seq data (Figure 2A). Further analysis indicated a 2.3-fold higher abundance of NKT cells in recently involuted mammary tissue (15 days post-offspring weaning), compared to mammary glands from nulliparous mice, or those exposed to pregnancy hormones for 12 days (mid-pregnancy), suggesting that the expansion of NKTs takes place at the final stages of post-pregnancy mammary involution (Figure S3A). The selective expansion of NKTs was further supported by the analysis of markers that define mammary resident neutrophils (Ly6G+) and macrophages (CD206+), which were largely unchanged between pre- and post-pregnancy mammary tissue (Figures S3B and S3C). Immunofluorescence analysis of Cxcr6-GFP-KI mammary tissue, previously described to label NKTs (Germanov et al., 2008), demonstrated several GFP+ cells surrounding ductal structures, an observation that supports the presence of NKTs within the mammary tissue (Figure S3D). Moreover, analysis of bone marrow and spleen from nulliparous and parous mice showed no difference in the abundance of NKTs, suggesting that pregnancy-induced expansion of these cells is mammary specific (Figures S3E and S3F).

Figure 2.

Pregnancy induces the expansion of a specific population of NKTs

(A) Flow cytometry analysis of resident CD45+ NK1.1+CD3+ NKTs from pre- and post-pregnancy mammary tissue. n = 5 nulliparous and 5 parous female mice. *p = 0.0004.

(B) Flow cytometry analysis of the classic NKT cell markers T-bet, CD335, and IFN-γ in NKTs (CD45+NK1.1+CD3+) from pre- and post-pregnancy mammary tissue. For Tbet analysis, n = 4 nulliparous and 4 parous female mice. *p = 0.016. For CD335 analysis, n = 7 nulliparous and 7 parous female mice. *p = 0.03.

(C) Flow cytometry analysis of β and γδ TCRs of pre- and post-pregnancy mammary NKTs. n = 5 nulliparous and 5 parous female mice. *p = 0.005.

(D) Gene set enrichment analysis of differentially expressed genes in FACS-isolated NKTs from pre- and post-pregnancy mammary tissue.

(E) Venn diagram demonstrating the number of shared and exclusive ATAC-seq peaks of FACS-isolated NKTs from pre- (blue circle) and post-pregnancy (orange circle) mammary tissue.

(F) Genome browser tracks showing distribution of MACS-called, ATAC-seq peaks at the Pdk4, Maged1, and Lypla1 genomic loci from pre- and post-pregnancy NKTs. For all analyses, error bars indicate standard error of mean across samples of the same experimental group. Statistically significant differences were considered with Student’s t test p < 0.05.

See also Figures S5 and S6 and Table S3.

To further characterize the identity of the post-pregnancy, mammary resident NKTs, we combined cellsurface and intracellular staining to detect canonical NKT lineage markers, including the NKT master regulator Tbet, the NKT/T cell secreted factor interferon-γ (Ifn-γ), and the NKT lineage marker Nkp46 (CD335) (Yu et al., 2011). Pre- and post-pregnancy, mammary resident NK1.1+CD3+ cells expressed all three markers, supporting their NKT identity. However, we detected a 2-fold increase in the percentage of post-pregnancy cells expressing Tbet, Ifn-γ, and CD335, suggesting that specific populations of NKTs are expanded in post-involuted mammary tissue (Figure 2B).

We also investigated whether pregnancy-induced NKTs represented a specialized population of CD8+ T cells, a cytotoxic cell type reported to reside in mammary tissue (Wu et al., 2019). We found that a fraction of NKTs present in both pre- and post-pregnancy mammary tissue expressed CD8 on their surface, accounting for 41% and 35% of the total NKTs, respectively (Figure S3G). To determine whether the triple-positive (CD3+NK1.1+CD8+) cells contributed significantly to the expanded population of post-pregnancy NKTs, we analyzed the mammary tissue of nulliparous and parous RAG1 knockout (KO) mice, which lack mature CD8+ T cells (Mombaerts et al., 1992). We observed a 10-fold expansion of NKTs in RAG1KO post-pregnancy mammary tissue, suggesting that CD8-expressing cells do not comprise a significant fraction of pregnancy-induced NKTs (Figure S3H). These results are consistent with our scRNA-seq data and further validate the existence of specific NKT subtypes in mammary glands after a full pregnancy cycle.

NKTs have multiple roles, including tissue homeostasis, host protection, microbial pathogen clearance, and anti-cancer activity, mediated through their ability to recognize both foreign- and self-antigens via TCRs (Balato et al., 2009). Therefore, we next investigated changes to the TCR repertoire of mammary resident, post-pregnancy NKTs. We found that 17% of pre-pregnancy NKTs expressed γδTCRs, in marked contrast to post-pregnancy NKTs, which mostly expressed γδTCRs (44%) (Figure 2C, top panels). A pregnancy cycle did not alter TCR composition across all immune cells, given that mammary resident, pre- and post-pregnancy CD8+ T cells mostly express αβTCRs, suggesting that parity promotes expansion of subtypes of NKTs that bear a specific TCR repertoire (Figure 2C, bottom panels).

We next investigated the molecular signatures of fluorescence-activated cell sorting (FACS)-isolated, mammary resident, NKTs. Unbiased pathway analysis of bulk RNA-seq datasets revealed the enrichment of post-pregnancy NKTs for processes controlling overall NKT development and activation, such as Notch signaling, tumor necrosis factor-α (TnfNF-α) signaling, transforming growth factor-β (Tgf-β) signaling, response to estrogen, and cMYC targets (Oh et al., 2015; Almishri et al., 2016; Doisne et al., 2009; Huber, 2015; Mycko et al., 2009). Conversely, pre-pregnancy NKTs were mainly enriched for processes previously associated with reduced immune activation, such as Ifn-α response (Bochtler et al., 2008) (Figure 2D; Table S3).

The activation of specific processes in post-pregnancy NKTs was also evident from analysis of their accessible chromatin landscape. ATAC sequencing (ATAC-seq) profiles showed similar genomic distributions of accessible regions across pre- and post-pregnancy NKTs, with a 93% overlap of their total accessible chromatin regions, suggesting that parity-induced changes did not substantially alter the chromatin accessibility associated with NKT lineage (Figure 2E; Figure S4A). General TF motif analysis identified chromatin accessible regions bearing classic NKT regulator DNA binding motifs such as Tbet, Plzf, and Egr2, further supporting their NKT lineage identity (Seiler et al., 2012) (Figure S4B). Analysis of accessible chromatin exclusive to post-pregnancy NKTs showed an enrichment for terms/genes associated with regulation of the adaptive immune response, killer cell activation, and antigen presentation, such as Pdk4, Maged1, and Lypla1, all involved in enhanced immune activation (Na et al., 2020; Connaughton et al., 2010; Lee et al., 2016; Jehmlich et al., 2013) (Figure 2F; Figure S4C). DNA motif analysis at accessible regions exclusive to post-pregnancy NKTs identified enrichment of specific TF motifs, including those recognized by Maf, a factor associated with an activated NKT state, and previously predicted by our scRNA-seq data to be expressed in cell clusters with an NKT identity (Figure S4D).

Overall, our analyses confirmed that post-pregnancy mammary tissue has an altered γδNKT cell composition, which bears molecular and cellular signatures of activated and mature adaptive immune cells.

NKT expansion requires CD1d expression on post-pregnancy MECs

Classically, NKTs are subdivided based on their activating antigens, including the main antigen-presenting molecules MHC class I, MHC class II, and the non-classical class I molecule, CD1d, which can be expressed on the surface of macrophages and dendritic cells, as well on the surface of epithelial cells (Gapin et al., 2013; Rizvi et al., 2015; Thibeault et al., 2009). Therefore, we next analyzed whether the expression of antigen-presenting factors on the surface of mammary epithelial and non-epithelial cells could underlie NKT cell expansion after pregnancy.

Flow cytometry analysis detected a 5-fold increase in CD1d levels on the surface of post-pregnancy luminal and myoepithelial MECs, in marked contrast to the levels of MHC-I and MHC-II proteins, which were largely unchanged across pre- and post-pregnancy MECs (Figures 3A and 3B; Figures S5A and S5B). No difference in surface expression of CD1d on mammary CD45+ immune cells was detected, suggesting that signals provided by CD1d+ MECs could promote the post-pregnancy expansion of mammary NKT cells (Figure S5C).

Figure 3.

NKT expansion depends on CD1d expression on post-pregnancy MECs

(A and B) Flow cytometry analysis and quantification of CD1d+ MECs harvested from pre-pregnancy (black bars, n = 8) and post-pregnancy (pink bars, n = 10) mammary tissue. *p = 0.0036 for luminal MECs and **p = 0.0006 for myoepithelial MECs.

(C) Genome browser tracks showing MACS-called, H3K27ac ChIP-seq peaks at the Cd1d genomic locus in FACS-isolated, pre- and post-pregnancy luminal MECs.

(D) H&E-stained histological images and duct quantification from mammary glands harvested from nulliparous (top left, n = 6) and parous (bottom left, n = 7) CD1dWT female mice and nulliparous (top right, n = 6) and parous (bottom right, n = 7) CD1dKO female mice. p = 0.86 for pre-pregnancy glands and p = 0.78 for post-pregnancy glands. Scale: 7 mm. Zoom-in panels, scale 500 μm.

(E) Flow cytometry analysis of mammary resident NKTs from pre- and post-pregnancy CD1dKO mammary tissue. n = 4 nulliparous and n = 4 parous female mice. *p = 0.3.

(F) Flow cytometry analysis of α and γδ TCRs of mammary resident NKTs from pre- (n = 3) and post-pregnancy (n = 3) CD1dKO mammary tissue. *p = 0.5.

For all analyses, error bars indicate standard error of mean across samples of the same experimental group. Statistically significant differences were considered with Student’s t test p < 0.05. See also Figure S7.

Gene-expression analysis of scRNA-seq datasets and qPCR quantification of FACS-isolated epithelial cells confirmed that post-pregnancy MECs express higher levels of Cd1d mRNA, supporting that pregnancy-induced molecular alterations may represent the basis for the observed increase in percentage of CD1d+ MECs (Figure 1D; Figure S5D). In agreement, we observed increased levels of the active transcription marker histone H3 lysine 27 acetylation (H3K27ac) at the Cd1d genomic locus in FACS-isolated post-pregnancy MECs, suggesting that increased mRNA levels could be associated with parity-induced epigenetic changes at the Cd1d locus (Figure 3C). These observations were confirmed in organoid systems that mimic the transcription and epigenetic alterations brought to MECs by pregnancy signals (Ciccone et al., 2020), where pregnancy hormones induced upregulation of Cd1d mRNA levels and increased H3K27ac levels at the Cd1d locus (Figures S5E and S5F). Thus, pregnancy-associated signals may induce epigenetic alterations that subsequently increase Cd1d mRNA and CD1d protein levels in post-pregnancy MECs.

To investigate whether CD1d expression is required for the expansion of NKTs after parity, we analyzed mammary glands from CD1dKO mice, which bear reduced levels of activated NKTs (Faunce et al., 2005; Macho-Fernandez and Brigl, 2015; Mantell et al., 2011). Mammary glands from nulliparous and parous CD1dKO mice displayed similar numbers of ductal structures as CD1d wild-type (WT) female mice, suggesting that loss of CD1d does not majorly alter mammary tissue homeostasis (Figure 3D). Further flow cytometry analysis indicated no statistically significant changes in the percentage of NKTs in mammary glands of nulliparous CD1dKO mice (2.2% ± 0.8), compared to nulliparous CD1dWT mice (3% ± 1.6) (Figures 2A, left panel, and 3E, left panel). Conversely, we found a 7-fold decrease in the percentage of NKTs in mammary tissue from fully involuted, parous CD1dKO female mice (3% ± 1.5) compared to parous CD1dWT mammary tissue (26% ± 4), supporting the role of CD1d in regulating NKT activation after pregnancy (Figures 2A, right panel, and 3E, right panel). Moreover, we found no difference in the abundance of NKTs in glands from pre- and post-pregnancy CD1dKO female mice, consistent with lack of CD1d expression reducing the activation of NKTs (Figure 3E). The analysis of an additional mouse strain that is deficient in mature/activated NKTs, due to the deletion of the histone-demethylase Kdm6 (UtxKO mouse model), failed to detect an expansion of NKTs post-pregnancy, thus supporting that pregnancy induces the expansion of mature/active subtypes of NKTs (Beyaz et al., 2017) (Figure S5G). Moreover, NKTs observed in post-pregnancy CD1dKO mammary tissue mainly expressed αβTCR on their surface, in contrast to the γδNKTs observed in CD1dWT post-pregnancy glands, further confirming that loss of CD1d expression affects the expansion and activation of specific populations of NKTs in post-pregnancy mammary tissue (Figure 3F).

Collectively, our studies identified pregnancy-induced epigenetic changes that may control the increased expression of Cd1d mRNA in post-pregnancy MECs and elucidated a role for CD1d in mediating communication between MECs and the γδNKTs, unique to post-pregnancy mammary glands.

Lack of mammary oncogenesis is marked by NKT expansion and CD1d+ MECs

Parity resulted in the expansion of γδNKTs in the mammary gland in response to the upregulation of CD1d on MECs, thus pointing to a mechanistic connection between pregnancy-induced changes to MECs and immune cell biology. Pregnancy-induced molecular modifications to MECs have also been associated with an oncogene-induced senescence response to cMyc overexpression, and suppression of mammary oncogenesis (Feigman et al., 2020). Therefore, we next investigated whether pregnancy-induced mammary cancer protection was associated with the expansion of NKTs.

Flow cytometry analysis of pre- and post-pregnancy mammary tissue from cMyc-overexpressing female mice (DOX-treated, CAGMYC model) demonstrated a 1.5-fold increase in the abundance of total CD3+ T cells (Figure S6A). CD3+ T cell expansion was also observed in mammary tissue transplanted with CAGMYC post-pregnancy MECs and in organoid cultures derived from post-pregnancy CAGMYC MECs; both conditions previously shown to lack mammary oncogenesis, thus further suggesting a link between pregnancy-induced tumorigenic inhibition and specific changes to the adaptive immune system (Figures S6B and S6C). This selective expansion of CD3+ T cells was further supported by the analysis of markers that define mammary resident neutrophils (Ly6G+) and macrophages (CD206+), which were largely unchanged in mammary tissue transplanted with either pre- or post-pregnancy CAGMYC MECs (Figure S6B).

Further flow cytometry analysis identified a 6-fold increase in the percentage of NKTs in mammary tissue from parous CAGMYC female mice, which predominantly expressed γδTCRs (Figure 4A; Figure S6D). No changes in the abundance of CD8+ T cells or CD4+ T cells was observed between mammary tissue from nulliparous and parous CAGMYC female mice, supporting the parity-induced expansion of γδNKTs (Figures S6E and S6F) and suggesting that specific constituents of the mammary immune microenvironment may control tumorigenesis. In agreement, we also found a 5-fold higher percentage of CD1d+ luminal MECs in post-pregnancy mammary tissue, thus linking gain of CD1d expression and the expansion of γδNKTs, which may collectively play a role in blocking tumorigenesis (Figure 4B).

Figure 4.

Lack of mammary oncogenesis is marked by NKT expansion and CD1d+ MECs

(A) Flow cytometry analysis of mammary resident NKTs (CD45+NK1.1+CD3+) from DOX-treated, nulliparous (left panel, n = 5) and parous (right panel, n = 5) CAGMYC female mice. *p = 0.002.

(B) Flow cytometry quantification of CD1d+ luminal and myoepithelial MECs from DOX-treated, nulliparous (left panel, n = 16) and parous (right panel, n = 11) CAGMYC female mice. *p = 0.02.

(C) Mammary tumor-free survival plot of nulliparous (black line, n = 5) and parous (pink line, n = 5) Brca1KO female mice.

(D) H&E-stained histological images from mammary tissue and tumors from nulliparous (top panels) and parous (bottom panels) Brca1KO female mice. Scale: 5 mm. Zoom-in panels, scale: 500 μm.

(E) Flow cytometry quantification of CD1d+CD24high luminal MECs from Brca1KO pre-pregnancy mammary tumors (black bar, n = 3), Brca1KO post-pregnancy healthy mammary tissue (pink bar, n = 4), and Brca1KO post-pregnancy mammary tumor (blue bar, n = 1). *p = 0.02.

(F) Flow cytometry analysis of NKTs in normal mammary tissue from nulliparous, tumor-bearing, Brca1KO female mice (left panel, n = 4) and normal mammary tissue from healthy parous Brca1KO female mice (right panel, n = 4). *p = 0.003.

(G) Quantification of γδNKTs in normal mammary tissue from nulliparous, tumor-bearing, Brca1KO female mice (black bar panel, n = 4), in mammary tumor tissue from nulliparous Brca1KO female mice (blue bar, n = 3), and in normal mammary tissue from healthy parous Brca1KO female mice (black bar panel, n = 2). *p = 0.023 and **p = 0.008.

For all analyses, error bars indicate standard error of mean across samples of the same experimental group. Statistically significant differences were considered with Student’s t test p < 0.05. See also Figures S8-S12.

cMYC overexpression is present in approximately 60% of basal-like breast cancers, with cMYC gain of function commonly found in BRCA1-mutated breast cancers (Chen and Olopade, 2008; Grushko et al., 2004). Interestingly, women harboring BRCA1 mutations with a full-term pregnancy before the age of 25 benefit from pregnancy-induced breast cancer protection (Medina et al., 2004; Terry et al., 2018). Therefore, we developed an inducible mouse model of Brca1 loss of function, for the purpose of investigating how pregnancy-induced changes influence Brca1-null mammary tumor development. In this model, tamoxifen (TAM) induces homozygous loss of Brca1 function in cells that express the cytokeratin 5 gene (Krt5+ cells), which include MECs(dos Santos et al., 2013), cells from gastrointestinal tract (Sulahian et al., 2015), reproductive organs (Ricciardelli et al., 2017), and additional epithelial tissue (Castillo-Martin et al., 2010; Majumdar et al., 2012), in a p53 heterozygous background (Krt5CRE-ERT2Brca1fl/flp53−/+, hereafter referred as Brca1KO mouse).

Nulliparous Brca1KO mice exhibited signs of mammary hyperplasia approximately 12 weeks post-TAM treatment, which gradually progressed into mammary tumors at around 20 weeks after Brca1 deletion (Figures S6G and S6H). Brca1KO mammary tumors display cellular and molecular features similar to those previously described in human breast tissue from BRCA1 mutation carriers and animal models of Brca1 loss of function, including high EGFR and KRT17 protein levels and altered copy-number variation marked by gains and losses of genomic regions (Annunziato et al., 2019) (Figures S6I and S6J).

To investigate the effects of pregnancy on the mammary immune microenvironment and oncogenesis, age-matched, TAM-treated, Brca1KO nulliparous and parous female mice were monitored for tumor development (Figure S7A). Our study demonstrated that only 20% of the parous Brca1KO female mice developed mammary tumors (one out of five), compared to 100% of nulliparous Brca1KO female mice with mammary tumors (five out of five mice), thus indicating that a full pregnancy cycle decreases the frequency of Brca1KO mammary tumors by 80% (Figures 4C and 4D).

Histopathological analysis suggested that pre-pregnancy mammary tumors were quite diverse, as previously reported for tumors from Brca1KO mice (Brodie et al., 2001). These included poorly differentiated tumors, such as micro-lobular carcinomas with squamous trans-differentiation (Figure 4D, top rows, far-left panel), medullary-like carcinomas (Figure 4D, top rows, right panel), and solid carcinomas resembling high-grade invasive ductal carcinoma (IDC) (Figure 4D, top rows, left and far-right panels). Accordingly, the only tumor-bearing parous Brca1KO female mouse developed a poorly differentiated carcinoma with extensive squamous trans-differentiation and extensive necrosis, also previously reported for tumors from Brca1KO mice (Figure 4D, bottom rows, far-right panels). Additional histopathological analysis confirmed that mammary tissues from the remaining parous Brca1KO female mice (four out of five) were largely normal (Figure 4D, bottom rows, far-left, left and right panels; Figure S7B). Immunofluorescence analysis confirmed that both pre-pregnancy mammary tumors and post-pregnancy normal mammary tissue were indeed deficient for Brca1+ epithelial cells, indicating that the lack of mammary tumors in parous female mice was not due to inefficient Brca1 deletion (Figure S7C).

Flow cytometry analysis of Brca1KO MECs demonstrated a progressive loss of CD24midCD29high myoepithelial cells in tumor tissue from nulliparous (2.5-fold) and parous (2-fold) Brca1KO female mice, and a marked increase in the percentage of CD24highCD29low luminal-like MECs (Figure S7D). These results suggest that tumor progression in this model is accompanied by changes to the population of CD24high MECs, which has been associated with poor clinical outcomes in patients with triple-negative breast cancer (Chan et al., 2019). Further cellular analysis indicated a 2.7-fold increase in the percentage of CD24high/luminal CD1d+ cells in healthy, post-pregnancy Brca1KO mammary tissue compared to tissue from tumor-bearing nulliparous and parous Brca1KO mice, supporting that parity-induced expression of CD1d at the surface of MECs associates with inhibition of mammary oncogenesis (Figure 4E).

Given the increased levels of CD1d expression, we next investigated the presence of NKTs in mammary tissue from nulliparous and parous Brca1KO female mice. Flow cytometry analysis demonstrated a 3.8-fold increase in the percentage of NKTs in healthy, post-pregnancy Brca1KO mammary tissue compared to non-affected normal mammary tissue from tumor-bearing nulliparous Brca1KO mice, and mammary tumors from parous Brca1KO mice (Figure 4F; Figure S7E). Additional flow cytometry analysis demonstrated that approximately 70% of total NKTs from healthy, post-pregnancy Brca1KO mammary tissue expressed γδTCR, in marked contrast to NKTs from healthy (2.7%) and tumor-bearing (8.6%) mammary tissue from nulliparous Brca1KO mice (Figure 4G).

Collectively, our findings show that pregnancy-induced gain of CD1d expression at the surface of MECs and expansion of γδNKTs associates with lack of mammary oncogenesis in response to cMyc overexpression or Brca1 loss of function. These results support the link between pregnancy-induced molecular changes, mammary tissue immune alteration, and inhibition of mammary tumorigenesis in clinically relevant mouse models of breast cancer.

Functionally active NKTs are required to block malignant progression of post-pregnancy MECs

Given that we demonstrated that pregnancy-induced changes block mammary oncogenesis in two distinct models (Figure 4), and that cMyc gain of function is commonly found in Brca1-mutated breast cancers, we utilized the cMyc overexpression mouse model to further characterize the effects of the immune microenvironment on the malignant development of post-pregnancy MECs. Analysis of fat-pad transplantations into severely immune-deficient NOD/SCID female mice, which lack T cells, B cells, NK, and NKTs, indicated that 100% of mammary tissue injected with pre-pregnancy (n = 5) or post-pregnancy (n = 5) CAGMYC MECs developed adeno-squamous-like carcinomas with acellular lamellar keratin, high levels of cell proliferation (Ki67 staining), and increased collagen deposition (Trichrome blue staining) (Figures S8A-S8C). Therefore, NKTs, or associated adaptive immune cells, are required for the parity-associated protection from oncogenesis in the CAGMYC model.

Bulk RNA-seq analysis demonstrated that post-pregnancy CAGMYC MECs transplanted into the fat pad of NOD/SCID female mice were less effective at activating the expression of canonical cMyc targets and estrogen response genes, compared to transplanted pre-pregnancy CAGMYC MECs, in agreement with the previously reported transcriptional state of post-pregnancy CAGMYC MECs (Feigman et al., 2020) (Figure S8D). We also found that organoid cultures derived from post-pregnancy CAGMYC MECs transplanted into NOD/SCID female mice retained a senescent-like state, characterized by reduced p300 protein levels and moderately increased p53 protein levels, in agreement with the previously reported senescent state of post-pregnancy CAGMYC MECs (Feigman et al., 2020) (Figure S8E). Together, these findings indicate that oncogenic progression of post-pregnancy CAGMYC MECs is associated with the immune-deficient mammary microenvironment of NOD/SCID mice.

While our investigation of post-pregnancy CAGMYC MECs that were transplanted into the mammary tissue of immunosuppressed animals alluded to the importance of a robust immune system in blocking mammary tumorigenesis, it did not uncouple whether functionally active NKTs, or CD1d expression at the surface of MECs, act to block oncogenesis in post-pregnancy mammary tissue. Therefore, to determine whether signaling between CD1d+ MECs and NKTs is critical for the development of mammary oncogenesis after pregnancy, we developed a double-transgenic mouse model by crossing the DOX-inducible CAGMYC mice into a CD1dKO background (CAGMYC-CD1dKO).

Histology analysis indicated that mammary tissue from DOX-treated, nulliparous, and parous CAGMYC-CD1dKO female mice showed signs of hyperplasia with atypia and abnormal ductal structures (Figure 5A, left and far-right panels; Figure S9A). Conversely, mammary tissue from DOX-treated, CAGMYC-CD1dWT parous female mice lacked malignant lesions in response to cMyc overexpression, thus suggesting that CD1d expression is required to inhibit the development of malignant lesions in post-pregnancy mammary gland (Figure 5A, right panels; Figure S9A). Flow cytometry analysis showed a lack of NKTs in mammary tissue from both nulliparous and parous CAGMYC-CD1dKO female mice, in marked contrast to the observed expansion of γδNKTs in healthy post-pregnancy CAGMYC-CD1dWT mammary glands that lacked tissue hyperplasia, supporting that CD1d expression may control pregnancy-induced expansion/activation of NKTs, and thus block mammary tumorigenesis (Figure S9B; Figure 4A). To further determine whether loss of CD1d expression underlies the malignant transformation of post-pregnancy MECs, we performed mammary transplantation assays of CAGMYC-CD1dKO MECs into the fat pad of syngeneic animals (CD1dWT female mice). We found that 100% of mammary tissue injected with pre-pregnancy CAGMYC-CD1dKO MECs and 70% of mammary glands injected with post-pregnancy CAGMYC-CD1dKO MECs developed signs of malignant lesions, supporting that loss of CD1d expression impacts with pregnancy-induced breast cancer protection (Figure 5B, black font; Figures S9C and S9D). This last observation was in marked contrast to the finding in glands transplanted with post-pregnancy CAGMYC-CD1dWT MECs, which, as previously reported, did not present signs of malignant transformation (Feigman et al., 2020) (Figure 5B, blue font; Figures S9E and S9F).

Figure 5.

Functionally active NKTs are required to block malignant progression of post-pregnancy MECs

(A) H&E-stained images of mammary tissue harvested from DOX-treated (DD5), nulliparous CAGMYC-CD1dWT (far-left panels), nulliparous CAGMYC-CD1dKO (left panels), parous CAGMYC-CD1dWT (right panels), and parous CAGMYC-CD1dKO (far right panels) female mice. Green arrows indicate signs of malignant lesions/mammary hyperplasia. Green asterisks indicate normal-like ductal structures. Scale: 1 mm.

(B) H&E-stained images of DOX-treated, CD1dWT mammary tissue transplanted with pre-pregnancy CAGMYC-CD1dWT MECs (blue font, top far left panel), pre-pregnancy CAGMYC-CD1dKO MECs (black font, top panel), post-pregnancy CAGMYC-CD1dWT (blue font, bottom far left panel), or post-pregnancy CAGMYC-CD1dKO MECs (black font, bottom panel). Green arrows indicate signs of malignant lesions/mammary hyperplasia. Green asterisks indicate normal-like ductal structures. Scale: 500 μm. See also Figure S13.

Altogether, these results suggest that loss of CD1d, with concomitant loss of pregnancy-induced expansion of NKTs, supports the development of mammary malignant lesions, independently of parity. Moreover, our study elucidates that parity blocks the malignant transformation of MECs, both by inducing cell-autonomous, epigenetic alterations within the MECs, and non-autonomous communication between CD1d+ MECs and NKTs in the mammary gland.

DISCUSSION

In mammals, reprogramming of the immune system is initiated after birth and continues throughout the lifespan of an individual due to exposure to pathogens, hormonal fluctuations, and aging. This dynamic reprogramming is part of an immune surveillance system that detects abnormal cells across many tissues, helping to prevent cancer. Here, we characterized a population of NKT-like immune cells (NKTs) in post-pregnancy mammary tissue, and their role in inhibiting mammary oncogenesis.

Our findings suggest that post-pregnancy mammary homeostasis does not rely on the presence of γδNKTs, given the normal histology of mammary tissue in mice deficient for this cell type. It is possible that NKTs expand in response to the re-setting of whole-body immunity post-partum, with the child-bearing event providing signals that alter antigens across all maternal tissues as well as expanding specific immune cell populations. γδNKTs have been found in the pregnant uterus across many mammalian species, linking NKT specialization and the pregnancy cycle (Mincheva-Nilsson, 2003). Our results support that the expansion of NKTs was predominantly observed in post-involution tissue, thus suggesting that the immune reprogramming of mammary tissue takes place after lactation.

Several other immune subtypes have been described to be enriched in mammary tissue during gestation, lactation, and involution stages of mammary gland development. These studies identified alterations in leukocyte interaction with mammary ductal structures, as well as specific transcriptional changes, suggesting that cell interaction and cellular identity of mammary resident cells are affected by pregnancy-induced development (Dawson et al., 2020; Hitchcock et al., 2020). Our analysis of leukocytes, specifically macrophages and neutrophils, did not show alterations to such cell populations in healthy parous murine mammary tissue or in post-pregnant CAGMYC mammary tissue lacking malignant lesions. However, given that leukocytes have been implicated in the activation of NKTs (Macho-Fernandez and Brigl, 2015; Rizvi et al., 2015), it is possible that molecular alterations, rather than changes to cellular abundance or antigen presentation of leukocytes, could play a role in inducing or sustaining the population of NKTs in post-pregnancy mammary tissue.

Our studies also provide evidence linking pregnancy-induced immune changes with the inhibition of mammary oncogenesis. Our previous research focused on how post-pregnancy MECs assume a senescence-like state in response to cMyc overexpression, an oncogene-induced response that activates the immune system via the expression of senescence-associated genes (Braig and Schmitt, 2006). Here, we found that CD1d expression at the surface of post-pregnancy MECs and the presence of NKTs were linked with the inhibition of mammary oncogenesis in two independent models of breast cancer, illustrating how epithelial and immune cells communicate to support pregnancy-induced mammary cancer prevention. Given that NKTs were previously shown to interact with senescent cells, it is possible that pregnancy-induced activation of CD1d expression and NKTs expansion represent additional responses to oncogene-induced cellular senescence (Kale et al., 2020).

Women completing a full-term pregnancy before the age of 25 have an approximate one-third reduction of the risk of breast cancer (Medina et al., 2004). This benefit applies to the risk of all breast cancer subtypes, including those from women harboring BRCA1 mutations (Terry et al., 2018). Thus, our findings supporting a role for pregnancy in inhibiting the development of Brca1KO mammary tumors lends a clinical relevance to our studies. Interestingly, the mammary tumor from parous Brca1KO female mouse was associated with low abundance of γδNKTs and CD1d+ MECs, suggesting that loss of the pregnancy-induced epithelial to immune microenvironment communication may be part of cellular changes that support mammary tumorigenesis. In fact, the genetically engineered loss of CD1d expression, with a consequent deficiency in activated NKTs, supported the malignant progression of cMYC-overexpressing MECs, further illustrating a link between epithelial and immune cells in supporting pregnancy-induced mammary cancer prevention.

Our findings are based on studies performed in mice that became pregnant at a young age (~8 weeks old), which reinforced pregnancy-induced changes to epithelial cells, and their effect on immune recruitment and oncogenesis inhibition. However, it remains unclear why such strong, pregnancy-induced changes do not fully prevent the development of breast cancer (Nichols et al., 2019). It has been suggested that specific mammary epithelial clones with oncogenic properties reside within the mammary tissue after pregnancy and may give rise to late-onset mammary oncogenesis in aged mice (Li et al., 2020b). It is possible that such populations of rare MECs lose some of their pregnancy-induced molecular signatures over time, thereby bypassing oncogene-induced senescence and immune recognition, and ultimately developing into mammary tumors. Moreover, given that pregnancy-induced breast cancer protection becomes apparent ~5–8 years after pregnancy, it is possible that additional immune reprogramming influenced by genetic makeup, age at pregnancy, and/or overall post-partum health may further modify breast tissue and erase pregnancy-induced changes that inhibit breast cancer development.

Nonetheless, the connection between pregnancy, immunity, and oncogenesis could be used to develop therapies to block cancer development. Indeed, a series of preclinical models have been developed to optimize the delivery of CD1d stimulatory factors, such as αGalcer and KRN7000, and induce expansion of NKTs (Zhang et al., 2019). Such strategies are mostly side-effect free and, if proven to support the expansion of pregnancy-induced NKT cells, could be used in cases of high breast cancer risk, including those with genetic alterations and/or family histories of breast cancer. Additionally, the characterization of specific, pregnancy-induced TCR rearrangements may be leveraged in CAR-NKT immunotherapy, for example, which could also efficiently target disease that has already developed. Collectively, such strategies could improve breast health and decrease cancer risk in women who experience their first pregnancy after 35 years of age, when they are at a greater risk to develop breast cancer.

Limitations of the study

The majority of existing transgenic and knockout models of breast cancer utilize mammary gland-specific promoters to control oncogene activation, such as MMTV, BLG, and WAP, which are enhanced/activated by signals present during pregnancy and lactation, thus potentially confounding the analysis of the molecular basis of pregnancy and mammary cancer risk. Therefore, the development of new model systems of mammary tumorigenesis, that do not rely on pregnancy-induced promoters, will allow us to further understand the effect of pregnancy on oncogenesis across all breast cancer subtypes.

STAR★METHODS

RESOURCE AVAILABILITY

Lead contact

Further information for resources and reagents should be directed to and will be fulfilled by the lead contact, Camila dos Santos (dossanto@cshl.edu).

Materials availability

All unique/stable reagents generated in this study are available from the lead contact upon request.

Data and code availability

scRNA-seq, RNA-seq, ATAC-seq datasets were deposited into BioProject database under number PRJNA708263 [https://www.ncbi.nlm.nih.gov/bioproject/PRJNA708263], and are publicly available as of the date of publication. All accession numbers are listed in the key resources table. Results shown in Figure 1 (pre-pregnancy scRNA-seq) were previously deposited into BioProject database number PRJNA677888 [https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA677888]. Results shown in Figure S2C (pre- and post-pregnancy RNA-seq), Figure 3C (pre- and post-pregnancy H3K27ac ChIP-seq) were previously deposited in the BioProject database under numbers PRJNA192515 [https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA192515] and PRJNA544746 [https://www.ncbi.nlm.nih.gov/bioproject/PRJNA544746]. Results shown on Figure S7F (H3K27ac Cut&Run of organoid cultures) were previously deposited in the BioProject database under number PRJNA656955 [https://www.ncbi.nlm.nih.gov/sra/?term=PRJNA656955]. This manuscript does not report original code. Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

KEY RESOURCES TABLE

Reagent or resource	Source	Identifier
Antibodies
Biotinylated anti-CD45	Thermo Fisher Scientific	Cat# 13-0451-85; RRID:AB_466447
Biotinylated anti-CD31	Thermo Fisher Scientific	Cat# 13-0311-85; RRID:AB_466421
Biotinylated anti-Ter119	Thermo Fisher Scientific	Cat# 13-5921-85; RRID:AB_466798
Biotinylated anti-CD34	Thermo Fisher Scientific	Cat# 13-0341-82; RRID:AB_466425
eFluor 450 conjugated anti-CD24	Thermo Fisher Scientific	Cat# 48-0242-82; RRID:AB_1311169
PE-Cy7 conjugated anti-CD29	BioLegend	Cat# 102222; RRID:AB_528790
7-AAD viability staining solution	BioLegend	Cat# 420404; RRID:SCR_020993
PerCP-Cy5.5 conjugated anti-CD1d	BioLegend	Cat# 123514; RRID:AB_2073523
PE conjugated anti-CD1d	BioLegend	Cat# 140805; RRID:AB_10643277
APC conjugated anti-CD45	BioLegend	Cat# 103112; RRID:AB_312977
FITC conjugated anti-CD3	BioLegend	Cat# 100204; RRID:AB_312661
Alexa Fluor 700 conjugated. anti-NK1.1	BioLegend	Cat# 108730; RRID:AB_2291262
APC/Cy7 conjugated anti-CD8	BioLegend	Cat# 100714; RRID:AB_312753
PE conjugated anti-TCR g/d	BioLegend	Cat# 118108; RRID:AB_313832
APC conjugated anti-TCR b	BioLegend	Cat# 109212; RRID:AB_313435
APC conjugated anti-H-2Kb	BioLegend	Cat# 116517; RRID:AB_10568693
Pacific Blue conjugated anti-I-Ab	BioLegend	Cat# 116421; RRID:AB_10613291
Brilliant Violet 421 conjugated anti-CD206	BioLegend	Cat# 141717; RRID:AB_2562232
Alexa Fluor 700 conjugated anti-Ly6G	BioLegend	Cat# 127621; RRID:AB_10640452
PE conjugated anti-IFNg	BioLegend	Cat# 505808; RRID:AB_315402
Pacific Blue conjugated anti-T-bet	BioLegend	Cat# 644807; RRID:AB_1595586
eFluor 450 conjugated mouse IgG	Thermo Fisher Scientific	Cat# 48-4015-82; RRID:AB_2574060
FITC conjugated rat IgG	Thermo Fisher Scientific	Cat# 11-4811-85; RRID:AB_465229
PE-Cy7 conjugated mouse IgG	BioLegend	Cat# 405315; RRID:AB_10662421
biotinylated anti-CD1d	BioLegend	Cat# 123505; RRID:AB_1236543
anti-p300 antibody	Santa Cruz Biotechnology	Cat# SC-585; RRID:AB_2231120
anti-Vinculin antibody	Abcam	Cat# ab129002; RRID:AB_11144129
anti-p53 antibody	Leica Biosystems	Cat# P53-CM5P; RRID:AB_2744683
goat anti-rabbit IgG HRP	Abcam	Cat# ab6721; RRID:AB_955447
goat anti-mouse IgG HRP	Abcam	Cat# ab97051; RRID:AB_10679369
anti-Cytokeratin 5 (KRT5)	BioLegend	Cat# 905501; RRID:AB_2565050
anti-Cytokeratin 7/17 (KRT7/17)	Santa Cruz Biotechnology	Cat# sc-8421; RRID:AB_627856
anti-EGFR	Santa Cruz Biotechnology	Cat# sc-373746; RRID:AB_10920395
anti-AR	Santa Cruz Biotechnology	Cat# sc-7305; RRID:AB_626671
anti-Ki67	Spring Bioscience	Cat# M3062; RRID:AB_11219741
Alexa Fluor 647 conjugated anti-Cytokeratin 5 (KRT5)	Abcam	Cat# AB193895; RRID:AB_2728796
unconjugated rabbit anti-BRCA1	Bioss	Cat# bs-0803R; RRID:AB_10858843
Alexa Fluor 568 conjugated goat anti-rabbit IgG	Thermo Fisher Scientific	Cat# A-11036; RRID:AB_10563566
Alexa Fluor 488 conjugated anti-GFP	BioLegend	Cat# 338007; RRID:AB_2563287
Alexa Fluor 405 conjugated anti-Cytokeratin 8 (KRT8)	Abcam	Cat# ab210139; RRID:AB_2890924
Chemicals, peptides, and recombinant proteins
DNase I	Sigma	Cat #D4263
Collagenase A, type IV solution	Sigma	Cat #C5138-1G
ITS (Insulin/Transferrin/Sodium selenite)	GIBCO	Cat #41400-045
FGF-2	PeproTech	Cat #450-33
Progesterone	Sigma	Cat #P8783
17-β-Estradiol	Sigma	Cat #E2758
Prolactin	Sigma	Cat #L4021
Doxycycline	Clontech	Cat# 631311
Collagenase/Hyaluronidase 10x solution	Stem Cell Technology	Cat #07912
Growth factor reduced matrigel solution	Corning	Cat #356230
Trilogy	Cell Marque	Cat# 920P-10
ProLong Glass Antifade Mountant	Invitrogen	Cat# P36980
17b-Estradiol (0.5 mg/pellet) + Progesterone (10 mg/pellet)	Innovative Research of America	Cat# HH-112
Luminata Crescendo Western HRP substrate	Millipore	Cat# WBLUR0100
TrypLE Express	Thermo Fisher Scientific	Cat #12604-013
Dispase	Stem Cell Technology	Cat #07913
Critical commercial assays
Ovation ultralow DR multiplex system	Nugen Technologies	Cat #0331-32
Nextera DNA sample Preparation kit	Illumina	Cat #FC-121-1031
Ovation RNA-seq system (V2)	Nugen Technologies	Cat #7102-32
DNeasy Blood & Tissue Kit	QIAgen	Cat# 69504
SuperScript III kit	Thermo Fisher Scientific	Cat #18080-051
Deposited data
ATAC-seq data	This paper	PRJNA708263
RNA-seq data	This paper	PRJNA708263
WGS data	This paper	PRJNA708263
scRNA-seq data, Figure 1 (pre-pregnancy)	Henry et al., 2021	PRJNA677888
RNA-seq (pre- and post-pregnancy)	dos Santos et al., 2013	PRJNA192515
H3K27ac ChIP-seq (pre- and post-pregnancy)	Feigman et al., 2020	PRJNA544746
H3K27ac Cut&Run, Figure S7F	Ciccone et al., 2020	PRJNA656955
Experimental models: Organisms/strains
Mouse: BALB/c	Charles River	https://www.criver.com/
Mouse: NOD/SCID	Jackson Laboratory	https://www.jax.org/
Mouse: CAGMYC	Feigman et al., 2020	N/A
Mouse: Cxcr6-GFP KI	Jackson Laboratory	https://www.jax.org/strain/005693
Mouse: RAG1 KO	Jackson Laboratory	https://www.jax.org/
Mouse: UTX KO	Beyaz et al., 2017	N/A
Mouse: CD1d KO	Jackson Laboratory	https://www.jax.org/
Mouse: Krt5CRE-ERT2Brca1fl/flp53het	This paper	N/A
Oligonucleotides
Cd1d qPCR FWD: 5′ TCC GGT GAC TCT TCC TTA CA 3′	This paper	N/A
Cd1d qPCR REV: 5′ CTG GCT GCT CTT CAC TTC TT 3′	This paper	N/A
b-actin qPCR FWD: 5′ TGT TAC CAA CTG GGA CGA CA 3′	This paper	N/A
b-actin qPCR REV: 5′ GGG GTG TTG AAG GTC TCA AA 3′	This paper	N/A
Software and algorithms
Fiji	ImageJ	Version 2.1.0
Zen lite software, Blue edition	ZEN Digital Imaging forLight Microscopy	Version 2.0.0.0
FlowJo	BD Biosciences	Version 10.0
Prism	Graphpad	Version 9.0
CellRanger	Zheng et al., 2017	Version 3.1.0
Seurat	Stuart et al., 2019	Version 3.1.1
GSEA	Broad Institute	Version 3.0
BD FACSDiva Software	BD Biosciences	Version 6.0
STAR	Dobin et al., 2013	Version 2.4.0
Bowtie2	Langmead et al., 2009	Version 2.4.2
MACS2	Zhang et al., 2008	Version 2.2.5
GREAT	McLean et al., 2010	Version 4.0.4
HOMER	Benner et al., 2017	Version 4.11
Bedtools	Quinlan and Hall, 2010	Version 2.28.0
UCSC Genome Browser	Dreszer et al., 2013	N/A
Hisat2	Kim et al., 2015	version 2.1.0
DNAcopy	Seshan and Olshen, 2014	version 1.50.1
DESeq	Anders and Huber, 2010	N/A

EXPERIMENTAL MODEL AND SUBJECT DETAILS

Animal studies

All experiments were performed in agreement with approved CSHL Institutional Animal Care and Use Committee (IACUC). All animals were housed at a 12 hour light/12 hour dark cycle, with a controlled temperature of 72°F and 40%–60% of humidity. Balb/C female mice were purchased from The Jackson Laboratory and Charles River. RAG1KO mice (B6.129S7-Rag1tm1Mom/J, IMSR Cat# JAX:002216, RRID:IMSR_JAX:002216) were purchased from The Jackson Laboratory. VavCre UTXKO were generated as previously described (Beyaz et al., 2017). CXCR6-KO-EGFP-KI mice (B6.129P2-Cxcr6tm1Litt/J, IMSR Cat# JAX:005693, RRID:IMSR_JAX:005693) were purchased from The Jackson Laboratory. CAGMYC transgenic mouse strain was generated as previously described (Feigman et al., 2020). CD1dKO CAGMYC transgenic mouse stain was generated by crossing CD1dKO (C.129S2-Cd1tm1Gru/J, IMSR Cat# JAX:003814, RRID:IMSR_JAX:003814) mice with CAGMYC mice. Krt5CRE-ERT2Brca1fl/flp53het (Brca1KO) transgenic mouse strain was generated by crossing BlgCREBrca1fl/flp53het transgenic mouse strain (Trp53tm1BrdBrca1tmAashTg(B-cre)74Acl/J, IMSR Cat# JAX:012620, RRID:IMSR_JAX:012620) with Krt5CRE-ERT2 transgenic mouse strain (B6N.129S6(Cg)-Krt5tm1.1(cre/ERT2)Blh/J, IMSR Cat# JAX:029155, RRID:IMSR_JAX:029155). Female mice ranging from 3 weeks old to 30 weeks old were utilized in the described research.

METHOD DETAILS

Antibodies

All antibodies were purchased from companies as indicated below and used without further purification. Antibodies for lineage depletion: biotinylated anti-CD45 (Thermo Fisher Scientific Cat# 13-0451-85, RRID:AB_466447), biotinylated anti-CD31 (Thermo Fisher Scientific Cat# 13-0311-85, RRID:AB_466421), biotinylated anti-Ter119 (Thermo Fisher Scientific Cat# 13-5921-85, RRID:AB_466798) and biotinylated anti-CD34 (Thermo Fisher Scientific Cat# 13-0341-82, RRID:AB_466425). Antibodies for cell surface flow cytometry: eFluor 450 conjugated anti-CD24 (Thermo Fisher Scientific Cat# 48-0242-82, RRID:AB_1311169), PE-Cy7 conjugated anti-CD29 (BioLegend Cat# 102222, RRID:AB_528790), 7-AAD viability staining solution (BioLegend Cat# 420404, RRID:SCR_020993), PerCP-Cy5.5 conjugated anti-CD1d (BioLegend Cat# 123514, RRID:AB_2073523), PE conjugated anti-CD1d (BioLegend Cat# 140805, RRID:AB_10643277), APC conjugated anti-CD45 (BioLegend Cat# 103112, RRID:AB_312977), FITC conjugated anti-CD3 (BioLegend Cat# 100204, RRID:AB_312661), Alexa Fluor 700 conjugated. anti-NK1.1 (BioLegend Cat# 108730, RRID:AB_2291262), APC/Cy7 conjugated anti-CD8 (BioLegend Cat# 100714, RRID:AB_312753), PE conjugated anti-TCR γ/δ (BioLegend Cat# 118108, RRID:AB_313832), APC conjugated anti-TCR β (BioLegend Cat# 109212, RRID:AB_313435), APC conjugated anti-H-2Kb (BioLegend Cat# 116517, RRID:AB_10568693), Pacific Blue conjugated anti-I-Ab (BioLegend Cat# 116421, RRID:AB_10613291), Brilliant Violet 421 conjugated anti-CD206 (BioLegend Cat# 141717, RRID:AB_2562232), Alexa Fluor 700 conjugated anti-Ly6G (BioLegend Cat# 127621, RRID:AB_10640452). Antibodies for intracellular flow cytometry: PE conjugated anti-IFNγ (BioLegend Cat# 505808, RRID:AB_315402), Pacific Blue conjugated anti-T-bet (BioLegend Cat# 644807, RRID:AB_1595586). Antibodies for negative controls: eFluor 450 conjugated mouse IgG (Thermo Fisher Scientific Cat# 48-4015-82, RRID:AB_2574060), FITC conjugated rat IgG (Thermo Fisher Scientific Cat# 11-4811-85, RRID:AB_465229), and PE-Cy7 conjugated mouse IgG (BioLegend Cat# 405315, RRID:AB_10662421). Antibody for MaSC enrichment: biotinylated anti-CD1d (BioLegend Cat# 123505, RRID:AB_1236543). Antibodies for Western Blot: anti-p300 antibody (Santa Cruz Biotechnology Cat# SC-585, RRID:AB_2231120), anti-Vinculin antibody (Abcam Cat# ab129002, RRID:AB_11144129), anti-p53 antibody (Leica Biosystems Cat# P53-CM5P, RRID:AB_2744683), goat anti-rabbit IgG HRP (Abcam Cat# ab6721, RRID:AB_955447) and goat anti-mouse IgG HRP (Abcam Cat# ab97051, RRID:AB_10679369). Antibodies for Immunohistochemistry (IHC) staining: anti-Cytokeratin 5 (KRT5) (BioLegend Cat# 905501, RRID:AB_2565050), anti-Cytokeratin 7/17 (KRT7/17) (Santa Cruz Biotechnology Cat# sc-8421, RRID:AB_627856), anti-EGFR (Santa Cruz Biotechnology Cat# sc-373746, RRID:AB_10920395), anti-AR (Santa Cruz Biotechnology Cat# sc-7305, RRID:AB_626671), and anti-Ki67 (Spring Bioscience Cat# M3062, RRID:AB_11219741). Antibodies for Immunofluorescence (IF) staining: Alexa Fluor 647 conjugated anti-Cytokeratin 5 (KRT5) (Abcam Cat# AB193895, RRID:AB_2728796), unconjugated rabbit anti-BRCA1 (Bioss Cat# bs-0803R, RRID:AB_10858843), Alexa Fluor 568 conjugated goat anti-rabbit IgG (Thermo Fisher Scientific Cat# A-11036, RRID:AB_10563566), Alexa Fluor 488 conjugated anti-GFP (BioLegend Cat# 338007, RRID:AB_2563287), Alexa Fluor 405 conjugated anti-Cytokeratin 8 (KRT8) (Abcam Cat# ab210139, RRID:AB_2890924).

Mammary gland isolation

Female mice classified as Pre-pregnancy (nulliparous, never pregnant), Post-pregnancy (parous, 1 gestation cycle, 21 days of lactation and 40 days of involution post offspring weaning), were housed together for 1-2 weeks to allow for estrous cycle synchronization prior to mammary gland isolation. For the experiments utilizing exposure to pregnancy hormones (EPH), never pregnant female mice (~8 weeks old) were implanted with 21 days-slow-release estrogen and progesterone pellets (17β-Estradiol (0.5 mg/pellet) + Progesterone (10 mg/pellet) – Innovative Research of America Cat# HH-112) prior to mammary gland isolation (at D12 post pellet implantation). Females classified as involution D15 had 1 gestation cycle, 21 days of lactation and 15 days of involution post offspring weaning. In all cases, mammary gland isolation was performed as previously described (dos Santos et al., 2013). In short, mammary glands (one to four pairs per mouse) were harvested, minced, and incubated for 2 hours with 1x Collagenase/Hyaluronidase (10x solution, Stem Cell Technology Cat# 07912) in RPMI 1640 GlutaMAX supplemented with 5% FBS. Digested mammary gland fragments were washed with cold HBSS (Thermo Fisher Scientific Cat# 14175103) supplemented with 5% FBS, followed by incubation with TrypLE Express (Thermo Fisher Scientific Cat# 12604-013) and an additional HBSS wash. Cells were incubated with 2 mL of Dispase (Stem Cell Technology Cat# 07913) supplemented with 40 μL DNase I (Sigma Cat# D4263) for 2 minutes and then filtered through a 100 μm Cell Strainer (BD Falcon Cat# c352360). The single cell suspension was incubated with lineage depletion antibodies and loaded onto a MACS magnetic column (Miltenyi Biotec Cat# 130-042-401). Lineage negative, flow-through cells (epithelial cells) were utilized for flow cytometry, and transcriptomic analysis. Lineage positive cells (immune cells) were eluted from column with 3ml of MACS buffer and utilized for flow cytometry, transcriptomic and epigenomic analysis. For cell analysis, Dual Fortessa II cell analyzer (BD Biosciences) was used. Data analysis was performed using BD FACSDiva Software (RRID:SCR_001456) or FlowJo (FlowJo, RRID:SCR_008520). Statistically significant differences were considered with Student’s t test p-value lower than 0.05 (p < 0.05).

Flow cytometry analysis

Mammary resident cells (epithelial and non-epithelial) were harvested from both top and bottom mammary glands, and analyzed according to the bellow indicated strategy. For all flow cytometry analysis an average of 300,000 cells live cells (7-AAD negative) were recorded. Gating strategy for all flow cytometry analysis is available in Methods S1.

Mammary organoid culture

Mammary tissue dissected was minced and digested for ~40 minutes in Collagenase A, type IV solution (Sigma, Cat# C5138-1G), following a series of centrifugations to enrich for mammary organoids. Freshly isolated mammary organoids were cultured with Essential medium (Advanced DMEM/F12, supplemented with ITS (Insulin/Transferrin/Sodium selenite, GIBCO Cat# 41400-045, and FGF-2 (PeproTech, Cat# 450-33)) prior to analysis. For experiments shown in Figures S5E and S5F, organoid cultures were derived from normal mammary tissue from pre- or post-pregnancy Balb/C female mice (RRID:IMSR_CRL:028), cultured in the presence of FGF-2 for 6 days, following FGF-2 withdrawal for 24 hr and then incubated with Complete medium (AdDF+++, supplemented with ITS (Final Concentration:1x, Insulin/Transferrin/Sodium Selenite, GIBCO Cat# 41400-045), 17-β-Estradiol (Final concentration: 40ng/mL, Sigma Cat# E2758), Progesterone (Final concentration: 120ng/mL, Sigma Cat# P8783), Prolactin (Final concentration: 120ng/mL, Sigma Cat# L4021), as previously described (Ciccone et al., 2020). For experiments shown in Figure S6C, organoids cultures were derived from pre- or post-pregnancy CAGMYC MECs, following treatment with doxycycline (DOX, 0.1mg/mL, Clontech Cat# 631311) for 2 days (DD2). For experiments shown in Figure S8E, organoid cultures were derived from NOD/SCID female mice, transplanted with either pre- or post-pregnancy CAGMYC MECs, following treatment with doxycycline (DOX, 0.1mg/mL) for 2 days (DD2).

RT-qPCR

Lineage depleted MECs or organoid cultures were washed with 0.5mL 1x PBS, following RNA extraction with Trizol (0.5mL, Thermo Fisher Scientific, Cat# 15596018). Reverse transcription was carried out using SuperScript III kit (Thermo Fisher Scientific, Cat# 18080-051). RT-qPCR was performed using a Quantstudio 6 with SYBR Green Master mix (Applied Biosystems, Cat# 4368577). Relative mRNA expression of target gene was calculated via the ΔΔCt method and normalized to β-actin mRNA levels.

Cd1d qPCR primers: FWD: 5′ TCC GGT GAC TCT TCC TTA CA 3′ and REV: 5′ CTG GCT GCT CTT CAC TTC TT 3′.

β-actin qPCR primers: FWD: 5′ TGT TAC CAA CTG GGA CGA CA 3′ and, REV: 5′ GGG GTG TTG AAG GTC TCA AA 3′.

Mammary fat pad transplantation

MaSCs-enrichment was performed as previously described (dos Santos et al., 2013). In short, lineage depleted MECs were incubated with biotinylated anti-CD1d antibody, to allow for MaSC enrichment. CD1d-enriched MEC fractions were resuspended with 50% growth factor reduced matrigel solution (Corning, Cat# 356230) and injected into the cleared fat-pad of the inguinal mammary gland (anterior part of the gland). For experiments presented on Figure S5B CD1d-enriched MECs fractions (~100K) were injected into the mammary fatpad of 12 weeks old CAG-only female mice, followed by DOX-treatment and histology analysis. For experiments presented on Figure S8 CD1d-enriched MECs fractions (~100K) were injected into the mammary fatpad of 12 weeks old NOD/SCID (RRID:IMSR_JAX:001303) female mice, followed by DOX-treatment and histology analysis. For experiments presented on Figure 5 and Figure S9, pre- or post-pregnancy CAGMYC-CD1dWT MECs (~10K) or CAGMYC-CD1dKO MECs (~10K) were injected into the mammary fatpad of 8-10 weeks old CD1d WT female mice, and allowed 3-days of tissue engraftment prior to DOX-treatment for 5 days.

Histological analysis

For histological analysis, the left inguinal mammary gland was harvested and fixed in 4% Paraformaldehyde overnight prior to paraffin embedding. For conventional histological analysis, mammary gland tissue slides were stained with Hematoxylin and Eosin (H&E). For ductal quantification, mammary gland H&E histological images were uploaded into Fiji (Fiji, RRID:SCR_002285), and ducts present in the posterior part of the gland were manually counted. Immunohistochemistry staining (IHC) was performed on a Roche Discovery Ultra Automated IHC/ISH stainer. For Masson’s trichrome staining, Leica Multistainer Stainer/Coverslipper Combo (ST5020-CV5030) was used to stain slides according to standard reagents and protocols. Images were acquired using Aperio ePathology (Leica Biosystems) slide scanner in 40X lenses.

Immunofluorescence analysis

Paraffin-embedded mammary gland sections were deparaffinized in Xylene (Sigma Cat# 534056) and rehydrated, followed by antigen retrieval in Trilogy (Cell Marque Cat#920P-10). Tissue was washed in 1x PBS (phosphate-buffered saline) for 1 min then blocked with blocking solution (10mM Tris-HCl pH 7.4, 100mM MgCl2, 0.5% Tween 20, 10% FBS, 5% goat serum) for 4 hours in a humidified chamber. Sections were stained with the appropriate conjugated primary antibodies in blocking solution for 16 hours at 4°C. After subsequent washings with 1x PBS and blocking solution, tissues were incubated with DAPI (Sigma Cat# 10236276001) for 10 minutes to stain nuclei, and slides were mounted in ProLong Glass Antifade Mountant (Invitrogen Cat# P36980). Cell visualization and image collection was performed on a Zeiss LSM780 confocal laser-scanning microscope utilizing Zen lite software, Blue edition (ZEN Digital Imaging for Light Microscopy, RRID:SCR_013672) version 2.0.0.0.

Doxycycline treatment

Doxycycline was purchased from Takara Bio USA, Inc. (Cat# 631311) and sucrose was purchased from Sigma (Cat# S7903). DOX drinking solution (1 mg/mL) was prepared using sterile 1% sucrose water.

Tamoxifen treatment

Tamoxifen USP grade was purchased from Sigma-Aldrich (Cat# 1643306) and sunflower seed oil (European Pharmacopoeia grade) was purchased from Sigma-Aldrich (Cat# 88921). To prepare the working solution, the Tamoxifen powder was weighed and dissolved in ethanol by vortexing. Heat sterilized sunflower oil was added at a ratio of 19:1 oil:ethanol mixture to a final concentration of 5mg/100ul (one dose), heated to 55°C and shaken vigorously to homogenize the mixture. Krt5CRE-ERT2Brca1fl/flp53het transgenic female mice received a total of three intraperitoneal doses of Tamoxifen warmed to 37°C on alternate days.

Monitoring tumor growth

3 week old Krt5CRE-ERT2Brca1fl/flp53−/+ female mice were treated with TAM. Half of TAM-treated female mice were housed together (pre-pregnancy/nulliparous group), and the other half were paired with a male (1 female and 1 male per breeding cage). Breeding TAM-treated females were allowed to give birth, nurse the offspring (21 days), and were considered post-pregnant (parous) after 40 days from offspring weaning. Both pre- and post-pregnancy mice were monitored for signs of tumor growth, and added to the Kaplan-Meier curve as soon as there was a palpable tumor. Mice with a tumor burden exceeding the limit of the animal’s well-being (> 2 cm), or mice showing signs of distress independently of tumor development were euthanized. At experimental end point, mammary tissue or mammary tumors were harvested for histological and flow cytometry analysis. Statistical analysis was performed with Logrank (Mantel-Cox) test.

Western blot

DOX-treated and control organoid cultures were homogenized in 1x Laemmli sample buffer (Bio-Rad, Cat# 1610747). Samples were loaded into homemade 10% SDS-PAGE gel and transferred overnight to PVDF membrane (Bio-Rad, Cat# 162-0177) using wettransfer apparatus. Membranes were blocked with 1% BSA solution and incubated overnight with a diluted solution of primary antibody, followed by incubation with HRP-conjugated antibody for 40 minutes. HRP signal was developed with Luminata Crescendo Western HRP substrate (Millipore, Cat# WBLUR0100) in autoradiography film (Lab Scientific, Cat# XARALF2025). Developed films were scanned on Epson Perfection 2450 photo scanner.

scRNA-seq data analysis

Single cell RNA-seq data (pre-pregnancy mammary glands = 3,439 cells from n = 2 biological replicates; post-pregnancy mammary glands = 4,412 cells from n = 2 biological replicates) were aligned to mm10 using CellRanger v.3.1.0 (10x Genomics) (Cell Ranger, RRID:SCR_017344) (Zheng et al., 2017), and downstream processing was performed using Seurat v3.1.1 (SEURAT, RRID:SCR_007322) (Stuart et al., 2019). Cells with fewer than 250 features or higher than 10% mitochondrial gene content were removed prior to further analysis. Genes with fewer than 3 cells expressing them were removed, and the data were then log-normalized. Post-filtering analysis was performed on 3,075 cells (pre-pregnancy) and 4,029 cells (post-pregnancy). Principal component analysis was performed using the top 2,000 variable genes. This analysis was used to identify the number of significant components before clustering. Clustering was performed by calculating a shared nearest neighbor graph, using a resolution of 0.6. Subsetting into different cell types was performed using known markers for MECs, T cells, Myeloid cells, B cells and NK cells. Epithelial cells for both datasets were defined by the expression of Epcam, Krt8, Krt18, Krt5 and Krt14(cluster average expression > 2). Non-epithelial were cells considered having low expression of Epcam, Krt8, Krt18, Krt5 and Krt14. Epithelial lineage identification and T cell lineage identification was performed utilizing a previously validated gene signature (Henry et al., 2021). Genes used to define each immune cluster (differentially expressed genes, DEGs) were determined using known cell type markers and using the FindAllMarkers function, which uses a Wilcoxon Rank Sum test to identify differentially expressed genes between all clusters in the dataset. Cell cycle scoring was performed with the CellCycleScoring function, using the default gene lists provided by Seurat. Cell dendrograms were generated using the BuildClusterTree function in Seurat, using default arguments. Diffusion mapping was performed using the DiffusionMap function from the “destiny” R package (Angerer et al., 2016). Gene Set Enrichment Analysis (GSEA, RRID:SCR_003199) (Subramanian et al., 2005) was used for global analyses of differentially expressed genes.

RNA-seq library preparation and analysis

FACS-isolated pre- and post-pregnancy NKTs were collected and homogenized in TRIzol LS (Thermo Fisher Scientific, Cat# 10296010) for RNA extraction. Double stranded cDNA synthesis and Illumina libraries were prepared utilizing the Ovation RNA-seq system (V2) (Nugen Technologies, Cat# 7102-32). RNA-seq libraries were prepared utilizing the Ovation ultralow DR multiplex system (Nugen Technologies, Cat# 0331-32). Each library (n = 2 biological replicates per experimental condition) was barcoded with Illumina TrueSeq adaptors to allow sample multiplexing, followed by sequencing on an Illumina NextSeq500, 76bp single-end run. Analyses were performed with command-line interfaced tools such as FastQC (FastQC, RRID:SCR_014583) (Andrews, 2015) for quality control and Trimmomatic (Trimmomatic, RRID:SCR_011848) (Bolger et al., 2014) for sequence trimming. We used STAR (STAR, RRID:SCR_004463) for mapping reads (Dobin et al., 2013), FeatureCounts (featureCounts, RRID:SCR_012919) for assigning reads to genomic features (Liao et al., 2014) and DESeq (DESeq, RRID:SCR_000154) to assess changes in expression levels simultaneously across multiple conditions and in multi-factor experimental designs, incorporating information from multiple replicates (Anders and Huber, 2010). Genes with a statistically significant pvalue of p < 0.05 were considered differentially expressed. Gene Set Enrichment Analysis (GSEA) (Gene Set Enrichment Analysis, RRID:SCR_003199) was used for global analyses of differentially expressed genes (Subramanian et al., 2005). GSEA terms with statistically significant pvalue of p < 0.05 were selected for data plotting and data interpretation. For experiments presented on Figure 2D, FACS-isolated, pre- and post-pregnancy CD45+NK1.1+CD3+ NKT cells (n = 2 females per experimental group, n = 4 pairs of mammary glands per female, n = 2 biological replicates per experimental group) were utilized. For experiments presented on Figure S8D, total mammary tissue isolated from DOX-treated, NOD/SCID female mice transplanted with either pre- or post-pregnancy CAGMYC MECs (n = 2 biological replicates per group) were utilized.

ChIP-seq library analysis

Previously published H3K27ac ChIP-seq datasets (Feigman et al., 2020) were mapped to the indexed mm9 genome using bowtie2 short-read aligner tool (Langmead et al., 2009), using default settings. MACS2 peak-calling program (MACS, RRID:SCR_013291) (Zhang et al., 2008) was used to identify enriched genomic regions in this data by comparing the pulldown ChIP data to the control (Input) data using a q-value cutoff of 1.00−3 . Identification of genes closest to these differentially called peaks was performed using Genomic Regions Enrichment of Annotations Tool (UCSC Genome Browser, RRID:SCR_005780) (McLean et al., 2010). Peak visualizations were generated using the UCSC Genome Browser (UCSC Genome Browser, RRID:SCR_005780) (Dreszer et al., 2013).

Cut&Run library analysis

Previously published H3K27ac Cut&Run datasets (Ciccone et al., 2020), were mapped to the indexed mm9 genome using bowtie2 short-read aligner tool (Langmead et al., 2009) using default settings. Sparse Enrichment Analysis for Cut&Run (SEACR) peak-calling program (Meers et al., 2019) was used to identify enriched genomic regions with an empirical threshold of n = 0.01, returning the top n fraction of peaks based on total signal within peaks. The stringent argument was implemented, which used the summit of each curve. Identification of genes closest to these differentially called peaks was performed using Genomic Regions Enrichment of Annotations Tool (UCSC Genome Browser, RRID:SCR_005780) (McLean et al., 2010). Peak visualizations were generated using the UCSC Genome Browser (UCSC Genome Browser, RRID:SCR_005780) (Dreszer et al., 2013).

ATAC-seq library preparation and analysis

Nuclei of FACS-isolated, pre- and post-pregnancy NKTs were isolated utilizing hypotonic lysis buffer and incubated with Tn5 enzyme from Nextera DNA sample Preparation kit (Illumina, Cat# FC-121-1031) for the preparation of ATAC libraries. Each library (n = 2 per experimental condition) was amplified and barcoded as previously described (Buenrostro et al., 2013), then pooled for sequencing on an Illumina Nextseq500, 76bp single-end run. ATACseq library reads (n = 2 per cell condition) were mapped to the indexed mm9 genome using Bowtie2 short read-aligner (Bowtie 2, RRID:SCR_016368) (Langmead et al., 2009) and replicate alignment files were merged. MACS2 (MACS, RRID:SCR_013291) (Zhang et al., 2008) was used to identify enriched genomic regions in both conditions using a tag size of 25bp and a q-value cutoff of 1.00−2. Peaks were annotated using Homer (HOMER, RRID:SCR_010881) (Benner et al., 2017) with standard mm9 genome reference. Location of peaks was then grouped into intergenic, promoter and genic (containing 5′UTR, Exons, Introns, Transcription Termination Sites, 3′UTR, ncRNA, miRNA, snoRNA, and rRNA) regions. The UCSC genome browser (UCSC Genome Browser, RRID:SCR_005780) (Dreszer et al., 2013) was used to analyze genomic regions for overlap, using the Bedtools intersect function (BEDTools, RRID:SCR_006646) (Quinlan and Hall, 2010) . Any base pair overlap was enough to consider two regions “shared” and regions where no overlap existed defined the regions as exclusively being in one condition. The comparison was made into a Venn-diagram using tool available at https://www.meta-chart.com/venn.

DNA motif analysis

Peaks from pre- and post-pregnancy NKTs ATAC-seq libraries were utilized as input for an unbiased transcription factor analyses using Analysis of Motif Enrichment (AME) (McLeay and Bailey, 2010) and Find Individual Motif Occurrences (FIMO) (MEME Suite - Motif-based sequence analysis tools, RRID:SCR_001783) (Grant et al., 2011) was used to computationally define DNA binding motif regions to identify sequences of known motifs, with a statistical threshold of 0.0001.

Genomic library preparation and copy-number variation analysis

Mammary normal tissue and tumor from nulliparous Brca1KO female mice were dissociated as above described. Lineage depleted tumor cells were utilized for DNA extraction using DNeasy Blood & Tissue Kit (QIAGEN Cat# 69504). Genomic DNA was sonicated to an average of 300bp using Covaris E220 Focused-ultrasonicator. For library preparation, fragmented DNA went through standard end-repair (NEB Cat# E6050), dA-tailing (NEB Cat# E6053), and sequencing adaptor ligation (NEB Cat# M2200) steps. Following universal adaptor ligation, eight cycles of PCR was performed for each sample. During the PCR step, a unique pair of Illumina TrueSeq i7 index and i5 index was added to each sample. The PCR library was purified with AMPure XP beads (Beckman Coulter Cat# A63881), and quantified using NanoDrop spectrophotometer and Agilent Technologies 2100 Bioanalyzer. Whole-genome-sequencing libraries with different combination of Illumina indexes were pooled together for one lane of Illumina MiSeq. 150 base pairs from both ends were sequenced along with two 8-bp indexes. For CNV analysis, Read 1 of the sequence data was mapped to the mm9 reference genome using Hisat2 version 2.1.0 in single read alignment mode (Kim et al., 2015). The reference genome was divided into 5,000 variable-length bins with equal mappability as previously described (Baslan et al., 2012). The ratio of mapped reads in the tumor sample to mapped reads in the diploid sample (normal tissue) was used to compute a fitted piecewise constant function (segmentation). This segmentation used DNAcopy version 1.50.1 implementation of the circular binary segmentation algorithm (Seshan and Olshen, 2014) and the copy number profiles were plotted using R version 3.4.4 (R Core Team, 2019).

QUANTIFICATION AND STATISTICAL ANALYSIS

Data represent results from three or more independent biological replicates, unless otherwise specified. Sequencing data are from two biological replicates from each condition. All statistical analyses were performed using GraphPad Prism V9 software. For all analyses, error bars indicate standard error of mean across samples of the same experimental group. Statistically significant differences were considered with p-values lower than 0.05 (p < 0.05) from unpaired Student’s t tests, or otherwise indicated, as described in the figure legends.