Divergent aging of nulliparous and parous mammary glands reveals IL33+ hybrid epithelial cells

Domain: biology
Published: 2026-01-21
Editorial outcome: major_revision

Abstract
Aging increases breast cancer risk while an early first pregnancy reduces a woman’s life-long risk. Several studies have explored the effect of either aging or pregnancy on mammary stem/progenitor cells, however, the combined effect of both remains unclear. Here, we interrogate the functional and transcriptomic changes at single-cell resolution in the mammary gland of aged nulliparous and parous mice to discover that pregnancy normalizes age-related imbalances in lineage composition, while also inducing a differentiated cell state. Importantly, we uncover a minority population of Il33-expressing epithelial cells that express both luminal and basal markers (i.e. hybrid), which accumulate in aged nulliparous mice but are significantly reduced in aged parous mice. Functionally, IL33 treatment of mammary epithelial cells from young mice phenocopies aged nulliparous epithelial cells, induces proliferation and promotes formation of organoids with Trp53 knockdown. Collectively, our study demonstrates that pregnancy blocks the age-associated imbalances in lineage integrity in the basal layer, including a decrease in Il33+ hybrid cells, that could potentially contribute to pregnancy-induced breast cancer protection.

Introduction
A woman’s risk of breast cancer increases with age (median age at diagnosis = 62 years1), but an early first pregnancy (below the age of 30) significantly reduces lifetime risk2–8. Previous studies that have examined changes shortly after pregnancy suggest that post-pregnancy mammary epithelial cells (MECs) have reduced stem cell capacity9,10 and increased expression of differentiation markers11,12. On the other hand, the studies that examined changes with aging in nulliparous (no pregnancy) MECs found an accumulation of dysfunctional luminal progenitors13, lineage infidelity14, and altered differentiation programs15,16. Importantly, there is no clear understanding of how these processes combine and evolve with aging, i.e., how does pregnancy alter the process of aging in MECs?

Most of the previous studies focused on time points immediately post-pregnancy or 40 days post-involution (~3.3 human years), corresponding to a transient period when the mammary gland is at a higher risk of breast cancer17. However, pregnancy-induced protection does not take effect until almost 10 years post-birth4,17,18. To address this gap in knowledge, we simulated conditions that mimic an early first pregnancy (20–30 years in humans, 3–8 months in mice) and analyzed mice at the post-menopausal stage (>50 years in humans, ~18 months in mice). The long-term effect of pregnancy on the aging of MECs is important to determine because 75% of breast cancer diagnoses occur over the age of 5019, while most women in the United States have their first pregnancy between 20 and 33 years of age20,21. Using this 18-month time point in conjunction with 18-month nulliparous and 3-month (nulliparous and parous) mice, we interrogate the aging process of MECs with and without pregnancy while removing confounding factors affecting human samples22. Our study not only delineates the combined effects of aging and pregnancy but also uncovers a previously unknown Il33+ hybrid MEC population that accumulates with age and is reduced in mice that have undergone pregnancy.

Results
Parity induces lasting changes in the mammary gland with age
To understand whether pregnancy alters the aging of the mammary gland, we compared 18-month-old parous mice (18MP) that have undergone multiple pregnancies between 3 and 8 months of age with 3–4 month-old nulliparous (3MNP), 3–4 month-old parous that had undergone at least 2 pregnancies, are >3 weeks post-involution (3MP) and 18-month-old nulliparous (18MNP) (Fig. 1a). Hematoxylin & Eosin (H&E) staining and whole mount imaging showed no major differences in ductal density, morphology, or complexity of branching between 18MNP and 18MP (Supplementary Fig. 1a–d), indicating that parity does not confer long-lasting morphological changes to the mammary gland.

The mammary gland is a bilayered tree-like structure composed of two main epithelial cell lineages, luminal and basal cells23. Luminal cells can be further subdivided into two functionally distinct cell types: hormone receptor (HR)-high and -low cells24. To understand how these populations evolve with aging in parous mammary glands, we performed flow cytometry analysis using previously established markers of MEC lineages, basal (CD49fhigh/EPCAMmed/low) and luminal (CD49fmed/EPCAMhigh) (Fig. 1a, Supplementary Fig. 2). We found no differences in the overall epithelial or stromal compartments among the groups (Fig. 1b, c). In 18MNP mice, there was a significant increase in the basal population (% of lineage-/live cells) as compared to 3M NP mice (Fig. 1b, c), consistent with previous studies25. However, aged parous mice showed a decrease in the basal population as compared to the aged nulliparous, suggesting that pregnancy partially normalizes the aged-induced expansion of basal cells (Fig. 1b, c). Nevertheless, there is still an age-related increase in the basal population in mice that have undergone pregnancy (3MP versus 18MP), although it does not reach statistical significance (Fig. 1b, c).

To further interrogate stem/progenitor subpopulations within the mammary epithelium, we assessed expression of TSPAN8, which has been reported to mark a quiescent, hormone-responsive mammary stem cell population26. Our analysis revealed that TSPAN8+ MECs were more abundant in young parous glands compared to age-matched nulliparous glands (Supplementary Fig. 3), suggesting an early expansion of this population following pregnancy. Interestingly, aged nulliparous glands harbored higher levels of TSPAN8+ cells relative to aged parous glands (Supplementary Fig. 3), indicating that parity may restrict the age-associated accumulation of TSPAN8+ hormone-responsive stem-like population.

Similarly, within the luminal compartment, we observed a decrease in the luminal cells with aging (3MNP versus 18MNP), which is normalized by pregnancy (3MNP versus 18MP) (Fig. 1b, c). Moreover, the percentage of luminal cells does not change with age in mice that have undergone pregnancy (3MP versus 18MP) (Fig. 1b, c). Our analysis of the luminal sub-lineages showed that nulliparous mammary glands have a decreased population of HR-high (CD14-/CD133+) luminal cells as a consequence of age (3M NP vs. 18M NP, Fig. 1d, e), consistent with previous studies in 13–14 month mice24. In contrast, the HR-high population was not affected by age in parous mammary glands (3MP vs. 18MP, Fig. 1d, e). However, regardless of age, parity resulted in a robust reduction in the HR-high population, suggesting that the pregnancy-driven reduction in HR-high cells is maintained with age. Alternatively, this reduction of the HR-high luminal population in parous mice may represent a consequence of the expanded HR-low (CD14+/CD133-) luminal population, which was exclusive to parous mammary glands and unaffected by age (Fig. 1d, e). Given the observed upregulation of CD14 in luminal cells during involution27, our data suggest that while the age-induced expansion of basal cells is normalized by pregnancy, luminal cells retain a residual involution program, with an increased proportion of CD14+ HR-low luminal cells.

Terminal differentiation of stem and progenitor cells during pregnancy has been hypothesized to contribute to pregnancy-induced breast cancer protection6,8,9. To understand whether there are functional differences in the clonogenicity of basal or luminal cells between the 3MNP, 3MP, 18MNP, and 18MP mice, we performed in vitro organoid assays. In line with previous studies that have shown an age-induced increase in organoid formation25, we found that aged basal cells had the highest capacity to form organoids as compared to young basal cells, regardless of parity status (Fig. 1f, g). While basal cells from aged parous mice had a higher organoid formation capacity than young mice (3MP and NP vs. 18MP), the number of organoids was significantly lower compared to aged nulliparous mice (18MNP versus 18MP) (Fig. 1f, g). These data suggest that pregnancy normalizes some of the age-induced increase in the clonogenicity of basal cells. Luminal cells from 3MP mice had significantly reduced organoid formation capacity as compared to 3MNP (Fig. 1f, g). Interestingly, aging reduced organoid-forming capacity of luminal cells, and this was further reduced in aged parous mice (lowest organoid formation capacity) (Fig. 1f, g). These findings suggest that MECs from the aged parous glands are diminished in clonogenicity, supporting the hypothesis that pregnancy reduces regenerative potential.

Single-cell RNA-seq reveals rare hybrid MECs in aged nulliparous mice
Recent studies have demonstrated that minority populations of stem/progenitor cells evolve with age and are likely precursors of tumor initiation25,28. To investigate the cellular heterogeneity and aging-associated transcriptional changes that are altered with pregnancy at single-cell resolution, we performed single-cell RNA-sequencing (scRNA-seq) on mammary glands from 18MNP and 18MP mice. To increase the power of our analysis, we integrated single-cell transcriptomes from the Tabula Muris29 18M mammary gland dataset (see Methods). After filtering cells with low and high-expressing genes, high mitochondrial counts, and doublets (Supplementary Fig. 4, see Methods), we identified epithelial cells (basal, HR-high luminal, HR-low luminal), immune cells (B cells, T cells, dendritic, macrophages), endothelial cells, pericytes, and fibroblasts using the expression of well-established marker genes (Fig. 2a, b, Supplementary Fig. 5a). Intriguingly, we identified a minority population of cells co-expressing basal and luminal markers (hybrid) (Fig. 2c), that are present in 18MNP mice (~97% of cell cluster) but are significantly reduced in 18MP mice (~3% of cell cluster) (Fig. 2b). Moreover, we found no major differences in the proportion of B cells, T cells, macrophages, endothelial cells, dendritic cells, pericytes, fibroblasts, and HR-low luminal cells between 18MP and 18MNP (Fig. 2b). However, 18MP had a lower proportion of basal cells and HR-high luminal cells (Fig. 2b), consistent with our flow cytometry analyses (Fig. 1b–e). Moreover, we did not find differences in cell cycle between 18MNP and 18MP mice (Supplementary Fig. 5b).

To investigate the unknown hybrid population further, we calculated basal cell and luminal cell gene signature enrichment scores for all cells based on previously published transcriptomic analyses of adult MECs30. We found that the hybrid MEC population lies between the luminal and basal populations, implying that hybrid MECs display a mixed gene signature of both basal and luminal lineages (Fig. 2d). Pseudotime inference using PAGA31 and Slingshot32 also predicted the hybrid MECs lineage trajectory between basal and luminal populations (Fig. 2e, Supplementary Fig. 6a). To further understand whether these hybrid MECs represent a plastic/intermediate cell state, we applied CytoTRACE233 to agnostically predict cellular potency. The average cellular potency of hybrid MECs scored between basal and luminal cells (Fig. 2f), supporting our hypothesis that these cells represent a hybrid intermediate.

Interestingly, we found that the hybrid MECs are enriched in expression of Krt6a and Il33 (Fig. 2c, Supplementary Fig. 6b). KRT6A has been previously identified as a marker of bipotent luminal progenitor cells in the mammary gland34. KRT6A+ cells are normally found only in the luminal layer during development35 and do not express the basal marker Krt14. However, a previous report showed that a population of KRT14+/KRT6A+ cells expands during early stages of pregnancy35. IL33 has been previously implicated in activating pro-tumorigenic signaling pathways36–38, cancer stem cell maintenance39, and establishing an immunosuppressive microenvironment40. Although other genes such as Aldh3a1, Dmkn, Alcam, Ly6d, Dsc3, Ppp1r14c, S100a14, Lgals7, Dapl1, Atp1b1, Cldn10, and Anxa8 are enriched in the hybrid population (Fig. 2c), their expression is either not exclusive to this population or limited to a smaller fraction of cells in the cluster. We therefore used Krt6a and Il33 as markers to further understand the hybrid population.

Il33+/Krt6a+ cells are present in public datasets
To validate our findings in other published aging datasets, we used the recently published data from 3M and 18MNP datasets from Angarola et al.41 and the Smart-Seq2 data from Tabula Muris Senis (note: independent from the 10x dataset used for integration). Consistent with our dataset, we identified a minority Krt6a+/Il33+ population that appears as a distinct cluster (shown in pink) and was labeled as myoepithelial in the Angarola et al. dataset (Fig. 3a–c, 39 cells in the cluster) and basal in the Smart-Seq data (Fig. 3d–f, 103 cells in the cluster), although a small percentage of luminal cells are also Krt6a+/Il33+ (Fig. 3f). The proportion of Krt6a+/Il33+ cells is highest in the aged mammary glands as compared to young mammary glands, further confirming that an increase in the Krt6a+/Il33+ population occurs with aging (Fig. 3c, f).

To understand the dynamics of the hybrid population with pregnancy and aging, we integrated single-cell RNA-sequencing data from the mammary glands of 3-month-old nulliparous and 3-weeks post-involution mice from the Pal et al.42 dataset with our data. The highest proportion of Krt6a+/Il33+ cells were found in the 18M NP samples (~7%), with a much smaller proportion found in both 3MNP and 18MP samples (<1% in both, Fig. 3g). Intriguingly, there is a decrease in Krt6a+/Il33+ cells post-pregnancy (Fig. 3g, <0.3%), suggesting that the change in the proportion of cells occurs immediately after pregnancy but the accumulation in nulliparous mice is an age-dependent effect.

To understand whether the Krt6a+/Il33+ population is present during the different stages of mammary gland development, we analyzed the Pal et al.42 scRNA-seq dataset containing all mammary gland developmental stages. We found that Krt6a and Il33 expression peaks during early pubertal development, particularly within terminal end buds (TEBs)42 (Fig. 3h, Supplementary Fig. 6c). Taken together, these data suggest that aging coincides with an aberrant activation of Krt6a and Il33 expression, likely conferring a cell state similar to TEBs, but this is diminished by pregnancy.

An early pregnancy primarily protects women from ER+ breast cancer43; hence, we wanted to understand the relevance of our hybrid population found in aged nulliparous mice to this disease in the pre-malignant phase. Therefore, we used the recently published single-cell dataset generated from mouse mammary glands carrying mutations in Esr1 that arise after developing resistance to anti-estrogen therapies44 (Fig. 3i). Remarkably, we find that the cells that emerge after estrogen deprivation or inhibition in mutant Esr1 mammary glands, labeled as endocrine-suppressed mutant-enriched (ESME), are highly enriched in the proportion of Krt6a+/Il33+ cells (Fig. 3j) and expression of Krt6a and Il33 (Fig. 3k). These findings suggest that aging results in the emergence of a hybrid cell population that shares features of Esr1 mutant cells with high cellular plasticity.

KRT6A+ cells are present in aged nulliparous but not parous mammary glands
To expand our analysis from mRNA to protein, we performed immunofluorescence staining on 3MNP, 3MP, 18MNP, and 18MP mammary glands for KRT6A, together with KRT8 (luminal marker) and KRT5 (basal marker) (Fig. 4a). We found a 4.5-fold increase in the percentage of KRT6A+ cells in the 18MNP (6.1%) mice relative to 3M NP (1.34%) (Fig. 4b), suggesting that these KRT6A+ MECs accumulate over the course of aging. However, 18MP mice had a small percentage of KRT6A+ cells (<1%), less than the 3MNP (Fig. 4b), consistent with our single-cell analysis (Fig. 3g). We could not detect more than a couple of KRT6A+ cells in 3MP mice, suggesting that there is a decrease in the population at 3 weeks post-involution (Fig. 4b), again consistent with our single-cell analysis (Fig. 3g).

To further understand the localization of the KRT6A+ cells, we quantified the percentage of KRT6A+ cells in the luminal and basal layer (i.e., KRT6A+/KRT8+ or KRT6A+/KRT5+). We found that 95% of KRT6A+ cells in 3MNP glands were positive for KRT8 and resided in the luminal layer (Fig. 4b), consistent with previous studies34. In contrast, over 50% of KRT6A+ cells in the 18MNP gland were positive for KRT5 and localized to the basal layer. Interestingly, in the 18MP mice, the localization of the KRT6A+ cells was similar to 3MNP mice (~10%, Fig. 4b). This suggests that the KRT6A+ MECs that accumulate with age are likely due to the loss of lineage integrity of basal cells, but this integrity is maintained in mice that have undergone pregnancy. Moreover, we found an increased percentage of KRT5+ basal cells and a decreased percentage of KRT8+ luminal cells in 18M NP mice (Fig. 4c), consistent with our flow cytometry data (Fig. 1b, d), but 18MP cells were more similar to 3MNP cells.

To determine whether the phenotype is maintained ex vivo, we performed immunofluorescence staining on 3MNP, 3MP, 18MNP, and 18MP organoids derived from basal or luminal cells (hereafter referred to as basal organoids or luminal organoids) for KRT6A, KRT8, and KRT5. Consistent with the staining in the mammary gland, we found that basal organoids from 18MNP mice had the highest percentage of KRT6A+ cells (Fig. 4d, e). Moreover, we observed very few KRT6A+ cells in organoids derived from luminal cells (Fig. 4f, g). Flow cytometry analysis of the organoids using markers in Fig. 1 showed trends in basal and luminal cell proportions comparable to immunofluorescence staining experiments (Supplementary Fig. 7). Of note, the proportions of epithelial populations appear to be mostly dependent on age in basal organoids, but parity status in luminal organoids. These data suggest that basal cells from 18MNP mice retain the aberrant differentiation program in short-term in vitro culture.

In vitro IL33 treatment mimics aged nulliparous profile
Cellular plasticity and hybrid cell states have been proposed as permissive cell states that allow for malignant transformation. Although IL33 has been studied in the context of breast cancer36,39, its role in the normal mammary gland, cellular plasticity, and aging is not defined. To determine whether IL33 treatment affects primary MECs in vitro, we cultured basal or luminal cells from 3-month old nulliparous (3MNP) mice with recombinant IL33 (5 or 10 ng/mL). The treatment of basal and luminal cells with IL33 resulted in an increased number of basal-derived organoids (Fig. 5a) but had no effect on organoids derived from luminal cells (Fig. 5b). Interestingly, IL33 treatment of basal organoids for 2 weeks inhibited spontaneous in vitro luminal cell differentiation, resulting in an increase in the proportion of basal cells (Fig. 5c, d) but no differences were observed in luminal cells (Fig. 5e). Moreover, immunofluorescence staining and 3D imaging of organoids demonstrated an increase in the percentage of KRT6A+ and KRT5+ cells upon treatment with IL33 in basal cells, and a slight but non-significant decrease in KRT8+ cells (Fig. 5f, g). These changes were not observed in organoids derived from luminal cells (Fig. 5g, bottom panel), suggesting that IL33 treatment of young nulliparous basal cells, but not luminal cells, phenocopies the increased proportion of basal cells and accumulation of KRT6A+ hybrid MECs observed in aged nulliparous mice (Fig. 4a, b).

To test whether IL33 treatment increases clonogenicity in cells with tumor-promoting mutations, we transduced 3M NP basal and luminal cells with lentiviral shRNA targeting Trp5345 (Supplementary Fig. 8) and exposed them to IL33. We found that Trp53 knockdown resulted in increased clonogenicity of basal cells when treated with IL33, but not luminal cells (Fig. 5h). Thus, our data suggest that the age-induced increase in Il33+ hybrid MECs could confer a survival advantage as cells acquire oncogenic mutations.

In vivo IL33 treatment drives proliferation and alters lineage identity
To determine whether IL33 treatment phenocopies an age-related change in vivo, we treated young nulliparous mice (3MNP) with IL33, or PBS for 5 days as previously described46 and analyzed mammary glands at 1 and 8 days post-treatment (dpt) (Fig. 6a). Whole mount analysis of the mammary glands from IL33-treated mice showed a significantly increased area of the lymph node and dramatic dilation of the ducts at 1dpt (Fig. 6b, c). This phenotype closely resembles mammary glands transplanted with MECs expressing constitutively active PIK3CA47 and in mammary glands where FGF3 expression is ectopically induced48. There were no differences in the length of the primary duct (Fig. 6d). H&E staining confirmed a significant increase in size of the lumen in the IL33-treated group at 1dpt (Fig. 6e, f) but did not reveal fluid accumulation within the ducts; however, we cannot entirely exclude the possibility of dilation secondary to edema. Aside from the area of the lymph node, analysis of mice at 8dpt showed no significant changes in any of the parameters above, suggesting the effects of short-term IL33 treatment were transient (Fig. 6b–f).

To further investigate whether IL33 treatment altered the cell-type composition in the mammary gland, we performed flow cytometry analysis at 1dpt and 8dpt for epithelial and immune cell types. We find that IL33 treatment reduces the total number of MECs, largely due to an increase in the percentage of stromal cells and a decrease in the percentage of luminal cells resulting in an overall increase in basal cells (Fig. 6g). To understand whether these changes are due to proliferation or apoptosis, we performed EdU labeling during and after treatment (Fig. 6a) as well as Annexin V staining. Analysis of mice at 1dpt showed that IL33 induces significant proliferation of both basal (~2-fold increase) and luminal cells (~4-fold increase) (Fig. 6h), while no differences were observed by 8dpt. At 1dpt and 8dpt, there are small but significant differences in apoptosis in the basal compartment (Supplementary Fig. 9). Specifically, we observed a decreased percentage of cells in early apoptosis (DAPI-/AnnexinV+) in basal cells at both 1 and 8dpt (Supplementary Fig. 9b, left). While insignificant at 1dpt (p=0.06), we also observed a slightly elevated live population (DAPI-/AnnexinV-) in IL33-treated mice (Supplementary Fig. 9b, left). Within the luminal population, we find a non-significant increase in apoptosis with IL33 treatment at 1dpt, which is reversed at the 8dpt (Supplementary Fig. 9b, right). To determine whether ductal dilation was due to an increase in cell size, we used forward scatter area (FSC-A) as a proxy for cell size. We observe a small decrease in median FSC-A in IL33-treated MECs at 1dpt (Supplementary Fig. 9c, d). However, this difference was lost when we excluded the EdU+ cells, suggesting that the ductal dilation is not primarily due to an increase in cell size (Supplementary Fig. 9e). Given the relatively minor changes in apoptosis and cell size, we speculate that the effects of IL33 are largely due to a dramatic increase in proliferation and not cell death. Moreover, we did not see significant changes in immune cell composition (Supplementary Fig. 10a, b) except an increase in ILC2 cells, the responder cell to IL3346,49,50, confirming the response of the treatment (Supplementary Fig. 10c, d).

We next tested whether IL33 treatment induces KRT6A expression in the mammary gland. Remarkably, IL33 increased the percentage of KRT6A+ hybrid cells in mice at 1dpt but not in mice at 8dpt (Fig. 6i, j). Moreover, while IL33 treatment in vitro induces KRT6A expression in basal cells, the in vivo treatment induces expression in luminal cells. This suggests that although IL33 induces the KRT6A+ cell state, the responding cell type to exogenous IL33 is distinct in vitro and in vivo in young mice. This discrepancy could be due to the short-term in vivo treatment, lower proliferation of basal cells in vivo, and/or limited differentiation into luminal cells in vivo as compared to organoid cultures. In addition, IL33 treatment in vivo increased the percentage of SMA+ basal cells and decreased the percentage of KRT8+ luminal cells (Fig. 6i, j) consistent with our flow cytometry analysis (Fig. 6g), suggesting the IL33 treatment can partially recapitulate some age-related phenotypes in vivo (Fig. 1).

IL33+ hybrid MECs are present in normal human breast tissue
We next sought to determine whether the Il33+ hybrid MEC population is present in human samples. To this end, we sub-clustered epithelial populations from the Human Breast Cell Atlas (HBCA)51 scRNA-seq dataset. We identified a population of cells (sub-cluster 5, pink) that were classified as HR low luminal cells, though they formed a distinct cluster separate from the major HR low luminal cluster (Fig. 7a). Cells in sub-cluster 5 co-expressed basal and luminal markers (Fig. 7b, Supplementary Fig. 11a), suggesting that this cluster represents an intermediate cell state between basal and luminal cells. Within sub-cluster 5, we discovered a minority population of cells enriched in the expression of genes present in our hybrid MEC population, as seen in Fig. 3c, such as IL33, LY6D, ALCAM, DMKN, PPP1R14C, and DAPL1, among others (Fig. 7c, Supplementary Fig. 11b). Interestingly, this hybrid population is enriched in expression of MKI67 and TOP2A (Supplementary Fig. 11b), potentially indicating a connection between IL33 and aberrant proliferation.

However, we were unable to identify robust expression of KRT6A (Supplementary Fig. 11c), likely reflecting differences in mouse and human keratin markers. Moreover, expression of these hybrid genes confined to a relatively small population of cells within the cluster; this may be due in part to technical limitations such as dropout events or low sequencing depth that specifically affects genes that are expressed at low levels52,53. Furthermore, we did not observe differences in age or parity within this subcluster (Supplementary Fig. 11d). We therefore explored whether the percentage of IL33+ cells change with age or parity. We found that there is a significant increase in IL33+ cells with age in nulliparous women but no difference between aged nulliparous and parous women (Supplementary Fig. 11e). However, since the HBCA51 dataset does not include information on age at first pregnancy and a late pregnancy may not have the same effect, we are unable to conclusively determine whether this population declines in older women who experienced early pregnancy. Similarly, we analyzed the iHBCA22 to test whether the percentage of IL33+ cells changes with age (Supplementary Fig. 12a). First, we were able to identify a population of hybrid cells that is enriched in expression of IL33 (Supplementary Fig. 12b). Specifically, this distinct hybrid cluster expresses both basal and luminal markers (Supplementary Fig. 12b, black dashed circle) and is enriched in the Kumar et al. dataset. Moreover, in the Twigger et al. dataset54, there is also a population of IL33+ cells that only expresses luminal markers (Supplementary Fig. 12b, pink dashed circle). As the Twigger et al. dataset was generated using cells isolated from human milk, we speculate that the expression of IL33 could reflect unique biology associated with lactation and immunity. Although a distinct IL33+ cluster could not be identified in the other datasets that are part of the iHBCA, we were able to identify IL33+ cells in all of them (Supplementary Fig. 12c). Therefore, we explored whether the percentage of IL33+ epithelial cells and/or expression of IL33 changes with age in the iHBCA. Overall, we find that the percentage of IL33+ cells increases with age in all datasets except the Twigger et al. dataset (Supplementary Fig. 12d). Like the HBCA51, we could not identify enough aged patients that had an early pregnancy and no mutations in DNA-repair pathways, precluding our analysis with early pregnancy. Despite these discrepancies between our analysis of mouse and human scRNA-seq datasets, it is apparent that the minority population of IL33-expressing hybrid MECs is present in the adult human breast, and the percentage of IL33+ cells increases with age.

IL33 induces basal differentiation and growth in HMECs
To functionally test whether IL33 impacts differentiation of primary human mammary epithelial cells (HMECs) as seen in mice, we treated cells with IL33 in 2D and 3D. IL33 treatment increased the organoid formation capacity of HMEC cells (Fig. 7d, e), suggesting that IL33 has similar functions in primary mouse and human mammary epithelial cells. Flow cytometry analysis of both 2D and 3D cultures treated with IL33 showed an increase in the percentage of CD49f+/EpCAM- (basal-like) cells and a decrease in the percentage of CD49f-/EpCAM+ (luminal-like) cells (Fig. 7f–h, Supplementary Fig. 13a). Furthermore, we find that HMECs treated with IL33 in both 2D and 3D culture had an increased percentage of CD44+ cells, a marker for immature stem-like cells in breast cancer (Fig. 7i–k, Supplementary Fig. 13b). Thus, our analysis reveals a minority population of IL33+ cells in the human mammary gland that expresses luminal and basal markers, and IL33 treatment induces expression of basal markers in HMECs.

Methods
Mice
All mice used for this study were maintained at the UCSC Animal Facility/Vivarium in accordance with the guidelines of the Institutional Animal Care and Use Committee (Protocol #Sikas2311dn).
Wild-type C57BL/6 mice (young), retired breeders, and age-matched mice when appropriate were purchased from Charles River Laboratories and The Jackson Laboratory. Simultaneously, an aging colony was maintained in-house. Only female mice were used in the study, as male mice do not undergo pregnancy. The age of mice is indicated in the figures and figure legends. All mice were housed in climate-controlled rooms with a minimum air exchange rate of eight times per hour. The ambient temperature was maintained between 20 and 24°C, and the relative humidity ranged from 45 to 65%. An automated 12:12 h light-dark cycle was implemented in all animal holding rooms. The mice were recorded for estrous, and 18M mice were found to be non-cycling. Mice are euthanized by CO2 followed by cervical dislocation.

H&E staining and quantification
H&E staining was performed as previously described75. Briefly, mammary glands were dissected from parous and nulliparous mice, fixed in 4% PFA, and paraffin-embedded for histology. Sections underwent a standard staining protocol, with 5-min incubation in 100% xylene and ethanol (100%, 70%) followed by DI water for deparaffinization. Sections were stained with Hematoxylin for 30s, then bluing solution (NH4OH + MilliQ H2O, 2 min) and Eosin for 1 min, followed by dehydration and mounting (Permount). Visualization of mammary gland ducts was performed using the Leica Widefield microscope, and ductal quantification was conducted using Fiji.

Mammary gland whole mounts and branching analysis
Inguinal mammary glands were removed, adhered to a Superfrost Plus Microscope slide (Fisherbrand Cat. No. 1255015), and incubated in Carnoy’s Fixative overnight at room temperature. Mammary glands were then incubated in 70% ethanol, 50% ethanol, and then DI water for 15 min each. Fixed mammary glands were stained overnight with Carmine Alum at room temperature, then dehydrated in 15-min incubations in increasing ethanol grades, followed by clearing in xylene and mounted in Permount.
Images of the entire mammary gland whole mount for both parous and nulliparous samples were imported as a TIFF file onto the Fiji/ImageJ software, where a three-by-three tiled ROI was chosen from a consistent distance from the lymph node. The ROI was processed to remove background noise. The image was skeletonized, and the branches dilated to improve clarity. A Sholl analysis (neuroanatomy) plugin was used on the skeletonized image, and a vertical radius from the side of the image was set as a consistent measurement throughout all samples to account for decreased branching normally seen at the edges of an image. After setting the starting radius as zero microns and selecting the best-fitting Sholl methods, the Sholl analysis program was run. This generated a “Sholl Log Plot,” which indicated branching intersections over an area, as well as a “Sholl Regression Coefficient (k)”, which described branching complexity. A k value closer to zero indicated a higher branching complexity. A “Sholl Mask (heatmap)” was also generated to visually depict regions of higher branching complexity, where the red areas indicated maximal branching.

Immunofluorescence staining of mammary gland paraffin sections
Mammary glands were fixed in 4% paraformaldehyde (PFA) in PBS and embedded in paraffin for immunostaining. 5 μm sections were deparaffinized, dehydrated, and autoclaved for 15 min in Tris-EDTA buffer [10 mM Tris 1 mM EDTA (pH 9.0)] for antigen retrieval. Tissue sections were incubated overnight at 4°C with primary antibodies diluted in tris-buffered saline (TBS) + 5% bovine serum albumin (BSA) (antibodies are listed in Supplementary Table 2). Samples were subsequently washed (2X) with TBS + 0.05% Tween for 10 min and were incubated with donkey anti-rat Alexa Flour 647 (1:400), donkey anti-chicken Alexa Flour 488 (1:400), donkey anti-rabbit Alexa Flour 594 (1:400) conjugated secondary antibodies (Jackson ImmunoResearch Laboratories) in TBS + 5% BSA 0.1% Tween for 1 h at room temperature (RT). Samples were subsequently washed (3×) with TBS + 0.1% Tween and were incubated with DAPI in TBS + 5% BSA + 0.1% Tween (1:10,000) for 10 min. All the immunofluorescence sections and cells were mounted in Fluoromount-G (Genesee). Images were acquired by a Solamere Spinning Disk Confocal microscope. Images were processed using Fiji.

Immunofluorescence staining of primary MEC organoids
Organoids were collected at day 12–14 in culture and recovered from the growth factor-reduced Matrigel by incubating with Cell Recovery Solution (Corning) for 1 h on ice with gentle shaking. 5 mL collection tubes and pipette tips were precoated in PBS + 5% BSA. Organoids were washed with cold PBS (2×) before being fixed with 4% PFA in PBS for 45 min on ice. Fixed organoids were then washed in PBS (1×) and incubated with 0.2% glycine in PBS for 20 min at room temperature. After a final wash, fixed organoids were resuspended in 0.05% sodium azide in PBS and stored at 4°C.
Organoids were permeabilized in cold 100% methanol for 10 min on ice, washed in PBS (1×), and incubated in blocking buffer (0.1% BSA, 0.3% Triton-X, 5% normal donkey serum) for 3 h at room temperature with light shaking, followed by overnight incubation with primary antibodies in blocking buffer at 4°C with light shaking. The next day, they were washed (3×) in PBS + 0.3% Triton-X and incubated overnight with secondary antibodies in a blocking buffer at 4°C with light shaking in the dark. Organoids were again washed (3×) in PBS + 0.3% Triton-X and incubated with DAPI for 10 min at room temperature. After a final round of washes (3X), organoids were transferred to µ-Slide 8 Well Glass Bottom slide (Ibidi) pre-coated in poly-L-lysine (Sigma-Aldrich) and imaged on a Spinning Disk Confocal microscope. Images were processed using Fiji.

Organoid assays
Primary MECs were sorted into complete organoid media as previously described75,76 (Advanced DMEM F/12, 10% FBS, 1% PSA, 50 ng/mL EGF, 100 ng/mL Noggin, 250 ng/mL R-Spondin-1, 1X N2, 1X B27, 1X GlutaMAX, 10 mM HEPES) supplemented with Y-27632 (10 uM). 1000 cells were plated in 96-well ultra-low attachment plates seeded with a 50 µL mixture of growth factor-reduced Matrigel and irradiated L-Wnt3a-secreting-3T3 feeder cells (11,000 feeder cells per 50 uL growth factor-reduced Matrigel). Organoids were cultured at 37°C, 5% CO2 with added humidity. On day 4, Y-27632 supplemented media was removed and replaced with complete organoid media. Fresh media was added every other day thereafter. Basal and luminal organoids were imaged at day 6 and 10 in culture, respectively, on a Zeiss Live Cell microscope and analyzed on Biodoc.AI to collect data on organoid count and average size. For IL33 treatment experiments, organoids were cultured as described above with the addition of IL33 (0 ng/mL, 5 ng/mL, 10 ng/mL).

Mammary gland digestion and processing
L2-5 and R2-5 mammary glands were harvested, minced, and chemically digested overnight in Advanced DMEM F/12 with 1% PSA, collagenase/hyaluronidase, and DNAse I at 37°C, 5% CO2, and added humidity as previously described77. Briefly, partially digested glands were then mechanically digested by pipetting with a serological pipette until no tissue pieces were visible. Digested glands were washed with staining buffer (Hank’s Balanced Salt Solution, 2% Bovine Calf Serum, 1% PSA) and centrifuged (300×g) at 4°C for 5 min. Red blood cells were lysed with 5 mL of ACK Lysis buffer for 5 min, and cells were washed with 15 mL of staining buffer. Cells were treated with 0.25% Trypsin with EDTA and gently pipetted continuously for 2–3 min to digest the basement membrane. Cells were then treated with DNAse I and Dispase and pipetted continuously for 2–3 min to prevent clumping. The single-cell suspension was then filtered through a 40 µm mesh strainer and pelleted via centrifugation (300×g, 4°C, 5 min). Cells were then resuspended in staining buffer and transferred to FACS tubes for staining.

Fluorescence-activated cell sorting (FACS)
All antibodies used in the study have been previously validated by the manufacturer or by previous studies. Cells were stained with antibodies listed in Supplementary Table 2 for 15 min at room temperature, as previously described75,76. Stained cells were then washed with a staining buffer, resuspended with DAPI (1:10,000), filtered, and analyzed on a BD Biosciences FACSAria cell sorter. See Supplementary Fig. 2 for FACS gating strategies. Data were analyzed using FlowJo software (10.10.0).

Lentiviral Transduction
For Trp53 knockdown experiments, pSicoR-GFP-shTrp53 was a gift from Tyler Jacks (Addgene plasmid #12090; http://n2t.net/addgene:12090; RRID:Addgene_1209045). Viruses were produced in HEK293T cells using the second-generation lentiviral system and transfection using Lipofectamine 2000 (Life Technologies) as previously described77. Supernatants were collected at 48 h, filtered with a 0.45 µm filter, and precipitated with lentivirus precipitation solution (Alstem LLC) per the manufacturer’s instructions. Viral titers were determined by flow cytometry analyses of HEK293T cells infected with serial dilutions of concentrated virus.

IL-33 treatment of 3-month nulliparous mammary glands
Three-month-old nulliparous C57BL/6 mice were purchased from The Jackson Laboratory. Mice were injected peritoneally daily with 1 µg of mouse recombinant IL-33 (Catalog number 580504, R&D Systems) or PBS for five consecutive days as previously described46. Mice were humanely euthanized in accordance with the guidelines of the Institutional Animal Care and Use Committee (Protocol #Sikas2311dn), and mammary glands were collected 1 and 8 days after treatment for analysis.

EdU labeling and detection by flow cytometry
Wild-type 12-week-old C57BL/6 female mice were purchased from The Jackson Laboratory (#000664). An 8 mg/mL stock EdU solution was prepared by dissolving EdU (Life Technologies, #A10044) in sterile PBS (Life Technologies, #14190250). 10 mg/kg EdU was delivered via intraperitoneal (I.P.) injection into each mouse daily, as illustrated in Fig. 6a. After processing the mammary gland tissue to a single cell suspension, EdU was detected according to the manufacturer’s instructions (Life Technologies, Click-iT Plus EdU Alexa Fluor 488 Flow Cytometry Assay Kit #C10632). Since EdU detection requires fixed cells, samples were stained with Zombie UV (Biolegend #423107) according to the manufacturer’s instructions prior to fixation in order to assess viability. Samples were filtered and analyzed on a BD Biosciences FACSAria cell sorter.

Annexin-V cell death flow cytometry analysis
Following antibody labeling and washes, cell pellets were resuspended in Annexin-V Binding Buffer (Biolegend #422201) to reach a concentration of 1×10^6 cells/mL. 5 µL of Annexin V-PE (Biolegend, #640908) was added for every 100 uL of single cell suspension. Samples were incubated at room temperature away from light for 15 min before washing with Annexin-V Binding Buffer. Phases of cell death were determined by DAPI and Annexin V signal (see Supplementary Fig. 10).

Primary human mammary epithelial cell (HMEC) culture
HMECs were purchased from ATCC (PCS-600-010), who authenticates and validates the cells, and cultured according to manufacturer instructions. All experiments involving HMECs were performed using low-passage cells (<P6). For 3D HMEC culture, 1000 HMECs were seeded in a 96-well ultra-low attachment plate on top of a layer of Corning™ Matrigel™ GFR Basement Membrane Matrix (Catalog number CB40230C). Colonies were grown in Mammary Epithelial Cell Basal Medium (ATCC PCS-600-030) supplemented with Mammary Epithelial Cell Growth Kit (ATCC PCS-600-040) and 2% GFR Matrigel. 2D and 3D cell cultures were maintained for seven days before analysis on a BD Biosciences FACSAria cell sorter.

Single-cell RNA-seq cell preparation and sequencing
Libraries were generated using the 10× Genomics Chromium Next GEM Single Cell 5’ Reagent Kits v2 (Dual Index) (cat#1000264) following the manufacturer’s protocol. Single-cell suspensions were prepared from mammary glands as described above and encapsulated using the Chromium Controller to generate Gel Beads-in-Emulsion (GEMs), allowing for barcoding of individual cells. Post-GEM generation, reverse transcription (GEM-RT) was performed within each GEM, followed by cleanup and cDNA amplification. The cDNA underwent fragmentation, end repair, and A-tailing. Dual index adapters were ligated to the cDNA fragments, which were then amplified using PCR. The libraries were quantified using a Qubit Flex Fluorometer (Invitrogen Q33327), and their size distribution was assessed using an Agilent TapeStation 4200. Sequencing libraries were constructed to target 200 million paired-end reads (400 million total) to achieve a coverage of 20,000 reads per cell from a targeted capture of 10,000 cells. Sequencing was performed on an Illumina NovaSeq 6000 following the manufacturer’s recommendations.

Single-cell data preprocessing and quality control
Sequencing data were processed using the 10× Genomics Cell Ranger pipeline, with reads aligned to the mouse reference genome (mm10). Quality control steps included removing: (1) genes detected in fewer than 3 cells, (2) cells with <600 or >8000 genes, (3) cells with total UMI counts <2000 or >12,000, and (4) cells with >1.5% mitochondrial gene content. The mitochondrial gene ratio was calculated as the percentage of UMIs mapped to mitochondrial genes relative to total cellular UMIs. Doublets were identified and removed using Scrublet78, excluding 723 cells. All preprocessing and downstream analyses were performed using Scanpy (v1.11.0).

Single-cell data integration, dimensionality reduction, and annotation
We integrated our single-cell dataset with a pre-processed 18-month nulliparous dataset from Tabula Muris Senis29. UMI counts were normalized to a target sum of 1e4 per cell, log-transformed, with a pseudocount of 1. The combined dataset included 18,341 cells and 26,056 genes. Principal component analysis (PCA) was performed for dimensionality reduction, and batch effects were corrected using Harmony (v0.0.4)79. Neighborhood graphs were computed on the corrected data and used for Leiden80 clustering (resolution=0.80), identifying sixteen distinct clusters. Clusters were manually annotated based on established marker genes. Data visualization was performed using Uniform Manifold Approximation and Projection (sc.tl.umap)81.

Cell cycle analysis and epithelial gene scoring
Cell cycle phase classification was performed using the Regev lab’s curated cell cycle gene set, assigning each cell to its respective phase82. To classify epithelial cells, we used established marker genes for basal, HR-high luminal, and HR-low luminal cells based on transcriptomic profiles of adult MECs30. Enrichment scores for these markers were computed across four epithelial populations: basal, HR-high luminal, HR-low luminal, and hybrid MECs. Fibroblasts and T cells (gray) are included as controls. Both cell cycle and epithelial subtype scoring were carried out using the sc.tl.score_genes function from the Scanpy library83.

Pseudotime inference and cellular potency prediction
Pseudotime analysis was conducted on epithelial subtypes (basal, HR-low, HR-high, and hybrid). Partition-based graph abstraction (PAGA, v1.2) was performed using Scanpy (v1.9.6), with diffusion maps computed using basal cells as the root. To validate trajectory inference, pseudotime was also estimated using pyslingshot (v0.1.3), the Python implementation of Slingshot, again using basal cells as the starting node. Slingshot pseudotime values were scaled between 0 and 1 for comparison with PAGA results. Cellular potency was assessed using CytoTRACE2 (v1.0.0) in Python. Visualizations were generated with Scanpy, Seaborn, and Matplotlib. To evaluate differences in pseudotime between epithelial subtypes, we performed pairwise Wilcoxon Rank Sum tests, applying Benjamini-Hochberg correction to control the false discovery rate. Statistical tests were implemented using scipy.stats.mannwhitneyu (SciPy v1.11.1) and multipletests from statsmodels.stats.multitest.

Analysis of published single-cell RNA sequencing data
We investigated the presence of hybrid cells in previously published and publicly available single-cell datasets as outlined in Supplementary Table 1. We analyzed three mouse and two human datasets. We used previously published processed single-cell data where available. Else, single-cell RNA sequencing data were processed and visualized using Scanpy (version 1.11.0) in Python. Following preprocessing and normalization, principal component analysis was performed for dimensionality reduction, followed by the construction of a nearest neighbor graph. Uniform Manifold Approximation and Projection (UMAP) embeddings were generated using the sc.tl.umap function to visualize cellular populations in two dimensions. We performed Leiden clustering at resolutions [0.02, 0.05, 0.1, 0.5, 2.0] and annotated and filtered for epithelial cells (luminal and basal) for downstream analysis.
Hybrid proportions were computed for each Leiden cluster as the (number of cells expressing Il33 and Krt6a/total cells in Leiden cluster) and plotted using seaborn’s boxplot function. In the human dataset, cells with no expression of the gene are plotted on the UMAP in gray.

Statistics and reproducibility statement
All graphs display the average as central values, and error bars indicate ±SD unless otherwise indicated. P values are calculated using paired or unpaired t-test, ANOVA, Wilcoxon rank-sum test, and Mann–Whitney U test, as indicated in the figure legends. All P and Q values were calculated using Prism (10.2.2) or Python (3.11.5), unless otherwise stated. For animal studies, sample size was not predetermined to ensure adequate power to detect a prespecified effect size; no animals were excluded from analyses, experiments were not randomized, and investigators were not blinded to group allocation during experiments.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this study Olander and colleagues have investigated cellular and molecular changes induced by pregnancy in mice. They analyzed mammary glands of 3-, 8-, and 18-month-old virgin or multiparous mice by flow cytometry to assess the relative proportion of luminal/basal cells and also performed RNA-seq on FACS sorted basal and luminal cells. Subsequently, they also performed single cell RNA-seq on mammary glands from 18-month-old virgin and multiparous mice. Lastly, they tested the effect of IL33, a gene they identified as differentially expressed in mammary epithelial cells 18-month-old virgin and multiparous mice. Overall, the manuscript has limited data and novelty.\nSpecific points:\n1. Full term pregnancy in early adulthood only reduces the risk of postmenopausal ER+ breast cancer, while the risk of triple-negative breast cancer (TNBC) is actually increased. Mice do not get ER+ mammary tumors, thus, they may not be the best model to use to assess the effect of pregnancy on ER+ breast cancer risk.\n2. There have been prior related studies in mice that compared retired breeders (old parous) and old virgin mice and also virgin mice and different stages of pregnancy/lactation. For example: https://pubmed.ncbi.nlm.nih.gov/25619437/ and scRNA-seq profiling (cited by the authors): https://pubmed.ncbi.nlm.nih.gov/34187545/\n3. Similarly, the increase in hybrid basal-luminal cells with age has been descried by multiple groups, both in human and mice, and IL33 has been shown to promote cancer stem cells in breast cancer.\n4. The authors have to include 3-month-old parous mice in their flow cytometry and RNAseq analyses (Figure 1-3) and also 3-month-old virgin and parous to see how quickly the differences due to pregnancy are observed. This is important, because the authors argue that their study is unique, since they looked at old parous mice compared to prior studies that analyzed mice shortly after pregnancy. But, w/o showing data for 3-month-old virgin and parous this is not possible to conclude.\n5. Similarly, in Figure 4 the authors argue that IL33 treatment of organoids from young mice phenocopies aging. But they do not analyze organoids from old mice and also do not compare 3-month-old virgin vs. parous. These additional samples would be critical to include.", "reviewer_2": "This study by Olander et al. investigates the combined effects of aging and pregnancy on mouse mammary epithelial cells (MECs). The authors report that pregnancy reduces the age-induced expansion of basal cells and partially normalizes the basal:luminal cell ratio in aged mice. They also report the presence of a small population of Il33-expressing hybrid MECs, which accumulate in aged nulliparous mice, but reduced in aged parous mice. The authors show that pregnancy induces transcriptomic changes in both basal and luminal cells, promoting a differentiated state and reversing some age-associated transcriptional programs. Finally, the authors show that IL33 treatment of young basal cells mimics the aged nulliparous MECs, increasing the proportion of basal cells and inducing a KRT6a+ hybrid state, which may confer a survival advantage for cells with oncogenic mutations. Overall, the study suggests that pregnancy-induced changes in MECs contribute to long-term breast cancer protection by maintaining lineage integrity and reducing the accumulation of potentially tumorigenic hybrid cells. The study is well presented and written. Generally, there is a slight over interpretation of some of the data which should be substantiated with further evidence. Below are the main points this reviewer believes should be addressed.\n1) Title: I personally feel that the title might be slightly overreaching, as the paper does not explicitly demonstrate that pregnancy reduces the formation of the hybrid progenitor. The findings just show that the hybrid progenitor is identified in the old nulliparous samples but is less prevalent in the old parous samples.\n2) Figure 2 claims that pregnancy “also reverses aged-associated transcriptional programs\". The authors have done DGE analysis to directly search for genes specific to that group (as far as I understand - they should specify the exact comparisons they do to generate these DGEs. Ie. are the genes for 18m NP a result of the test: 18m NP vs (3m NP + 18m P)? I suspect this is the case and thus forces the result that they claim. For the claim that Parity is undoing the effects of ageing, I would instead want to see how the genes for a 3m NP vs 18m NP test compares to those of the 18m NP vs 18m P. And potentially how they might oppose each other.\n3) Fig 2f, Again I am not sure this GO analysis has really shown this? To me it is just descriptive of the DEGs with each state but does not show evidence of parity undoing ages effects?\n4) Fig 3, It would be interesting to see the distribution of hybrid cells across samples. Are they primarily identified in specific samples, or are they generally upregulated in all aged NP samples?\n5) Fig 3. On the tabula muris data, it would be good to include graphs showing the integration, along with more detailed information on the distribution of cells across the samples. With this if they are claiming to have found a completely new cell type I would want them to find it in other datasets or provide some explanation why no one else has seen these cells?\n6) Fig3, I really want to see doublet scores for each of the clusters/cells calculated (scDblFinder, scrublet, or any other comparable tool). Additionally, why are there no pericytes/vascular mural cells in their dataset?\n7) Fig3d, the authors should consider adding a control plot here for fibroblasts (for example) to see what this analysis would look like for an unrelated population. Additionally I would want to know if these hybrid cells have any unique gene markers. This would help convince me they are a distinct cell type not some doublet/conglomeration of the two.\n8) Fig 4, the results are interesting BUT I'd really like to see what this treatment would do in the mice - does it do the same in vivo with a fully functioning immune system present?\n9) Finally, the authors need to perform a cross species analysis utilising the large human breast cell atlas datasets now available to see if these findings in the mouse do translate to the human.", "reviewer_3": "In this manuscript, Dr. Shaheen Sikandar and colleagues use RNAseq and scRNA techniques as well as organoid and immunostaining to investigate the cell and molecular changes in aged parous mammary glands compared to aged nulliparous glands as well as young nulliparous glands. This is an interesting developmental study considering both pregnancy and ageing that has not been well studied previously. The significant findings include the identification of a hybrid cell population that expands with ageing in nulliparous mice but this expansion is prevented in the parous aged group. In addition, the finding of elevated IL33 expression in this hybrid population is also interesting since IL33 has been implicated in cancer risk and the authors showed convincing data that IL33 treatment can promote organoid growth and hybrid cell expansion, mimicking ageing effects. However, there are also a number of significant weaknesses including over-interpretation of IL33 experiments on basal vs. luminal cells, premature experiments in testing IL33 and cancer relevance, and lack of adequate controls for some of the experiments. In addition, while data from the IL33 treatment of organoids provide an exciting and mechanistic link between ageing and hybrid cells, an in vivo validation, which should not take too long, would provide a much stronger argument. Detailed comments are listed below:\n1. Fig 1: An increase in the ratio of basal cells vs. luminal cells does not prove an increase of the basal population. It could be due to the drop of the luminal population, which seems to be the case here (Fig 1a). This ratio also neglects the change in the overall epithelial content in the comparing groups. Therefore, besides presenting the plots showing the luminal basal ratio and ER+ vs. ER- luminal ratio, please plot the basal and luminal population sizes for all three groups of mice. Please also make sure that the FACS images are representative of the plotted data.\n2. Fig 1: Without these absolute population size comparisons, the authors cannot conclude “the age-induced expansion of basal cells is normalized by pregnancy, the luminal cells retain a residual involution program, with an increased proportion of CD14+ HR low luminal cells.”\n3. Fig 2: it is shocking and highly interesting that their gene expression data on FACS sorted basal cells suggest parity leads to production of milk associated genes. However, since FACS sorting cannot guarantee no luminal cell contamination, it is critical that the authors confirm this suggestion by performing spatial assays such as IF or IHC to RNAscope to demonstrate the milk gene expression in the basal layer in aged glands.\n4. Paragraph starting line 123: please reference 18m B in the statement “We found significant enrichment of genes involved in cytoskeletal remodeling, extracellular matrix, and contractile functions (Fig. 2e).” The authors conclude this paragraph with the statement that “These data suggest that pregnancy normalizes age-induced transcriptional changes in luminal and basal epithelial cells, while also inducing a differentiated state.” But the data presented in this paragraph does not logically lead to this conclusion.\n5. IL33 treatment of 3m NP basal cells in organoids led to blocked luminal differentiation based on flow cytometry (Fig 4b), but no change of the proportion of K8+ luminal cells (Fig S6F). These two data seem to be contradictory.\n6. While IL33 treatment did not affect 3m NP luminal cells in forming organoids, this data does not mean that IL33 treatment would not affect other features such as cell differentiation and expression of K6 or K8 or K5. Please analyze these luminal derived organoids in the same way as the basal cell derived organoids were characterized. If the luminal organoids do not expand the K6+ population, the data would strength the claim of basal origin of the hybrid cell population.\n7. Studying IL33 impact on cells with Tp53 knockdown helps increase the cancer implication of this developmental work. However, this experiment is highly preliminary. Only one shRNA was used, and it lacks the confirmation of Tp53 knockdown. Also, non-targeting shRNA controls are missing but would provide the critical baseline comparison to understand the specific impact of p53 loss. In the absence of these data, it is premature and misleading to suggest that IL33 has a p53-dependent impact.\n8. There are no adequate and decisive data to conclude that the K6+ hybrid cells originate in the basal population as the authors claimed in Discussion.\n9. While the IL33 treatment of organoids experiments provide an exciting and mechanistic link between ageing and hybrid cells, an in vivo validation, which should not take too long, would provide a much stronger argument."}, "author_response": "Point by point response.\nOriginal comments are in blue, authors replies are in black.\nChanges are highlighted in green in the revised manuscript.\nPlease note: Figures included in the revised manuscript paper are referenced in the point-by-point\nresponse but not included in this document. Only figures specific to the point-by-point response are\nincluded here.\nReviewer #1 (Remarks to the Author):\nIn this study Olander and colleagues have investigated cellular and molecular changes induced by\npregnancy in mice. They analyzed mammary glands of 3-, 8-, and 18-month-old virgin or multiparous\nmice by flow cytometry to assess the relative proportion of luminal/basal cells and also performed RNA-\nseq on FACS sorted basal and luminal cells. Subsequently, they also performed single cell RNA-seq on\nmammary glands from 18-month-old virgin and multiparous mice. Lastly, they tested the effect of IL33, a\ngene they identified as differentially expressed in mammary epithelial cells 18-month-old virgin and\nmultiparous mice. Overall, the manuscript has limited data and novelty.\nWe thank the reviewer for their comments and assessment of our manuscript. We have\nsignificantly expanded the scope of our work, adding several new figures in the main and\nextended data as per the suggestions from all the reviewers. We believe that our revised work\nprovides an analysis of aging and pregnancy that has not been done before.\nSpecific points:\n1. Full term pregnancy in early adulthood only reduces the risk of postmenopausal ER+ breast cancer,\nwhile the risk of triple-negative breast cancer (TNBC) is actually increased. Mice do not get ER+\nmammary tumors, thus, they may not be the best model to use to assess the effect of pregnancy on ER+\nbreast cancer risk.\nWe appreciate the reviewer's insightful comment and welcome the opportunity to clarify\nour rationale. The goal of our study was to investigate how pregnancy influences the aging of\nstem/progenitor cells in the mammary gland, specifically examining whether certain subsets are\nenriched or lost with aging and pregnancy. Previous studies have long utilized mice as a model\nsystem to study parity-induced protection1,2 and understand underlying mechanisms3–5. More\nrecently studies have used mice to study aging in mammary epithelial cells and their relation to\nbreast cancer risk6–8. Thus, we followed past precedent in using mice as a model system to\nunderstand the long-term impacts of pregnancy.\nAs the reviewer correctly noted, the risk of breast cancer increases in the first 10 years\nfollowing pregnancy in humans. The key limitation of previous studies is that they looked at the\npost-pregnancy mammary gland when there is increased risk and it is not post-menopausal (peri-\nmenopause in mice)3–5. This makes our study unique, as it specifically examines a time point\nbeyond this 10-year window9–12, aiming to better mimic conditions relevant to human biology—\nwhere women who experience an early first pregnancy gain protective benefits much later in life,\ntypically after the age of 50. Moreover, none of the aging studies in mammary epithelial cells\naccount for pregnancy6–8.\nWhile it is true that mice usually get basal like breast cancer there are mouse models of\nER+ breast cancer that have been used to study human disease. Some of these include ESR1,\nCCDN1, prolactin, and WNT1 overexpression, PIK3CA gain of function and loss of STAT1 –\ncomprehensive reviews are here13,14. As one recent example published in Nature Cancer in 2025,\na group from Genentech generated and utilized a mouse model to study human ESR1 mutations\narising due to resistance to ERα-targeted therapies in breast cancer15,16. They showed that the\ncells present in the mammary gland of mice carry the ESR1 mutation that are resistant to\nendocrine suppression have a distinct gene expression pattern that is between luminal and basal\ncells. We have analyzed the single cell RNA-sequencing data from this manuscript15 (Fig. 3i). We\nfind that the cells arising from ESR1 mutant mice after endocrine suppression (labelled as\nendocrine suppressed mutant enriched or ESME cells in the published manuscript) are enriched\nin expression of Krt6a/Il33 (Fig. 3j,k).\nMultiple studies have also employed mouse ER+ breast cancer cell lines, which were\npreviously characterized17 and recently used to understand metastatic dissemination of ER+\nbreast cancer18. Lastly, a recent study showed that the same mutation that generates ER- (basal-\nlike) tumors in young mice gives rise to ER+ tumors in aged mice19. Since many of the earlier\nstudies introducing mutations to generate breast cancer have been done with young mice, it is\npossible that if these experiments were repeated with aged mice, they would give rise to ER+\nbreast cancer.\nFinally, as requested by Reviewer #2, we have extended our findings to human datasets\nand primary mammary epithelial cells. These results are now included as a new figure (Fig. 7) in\nthe main manuscript. Briefly, our findings identify a minority population of cells with a similar gene\nexpression pattern as seen in nulliparous aged mice (Fig. 7a-c, Extended Data Fig. 10). We also\ntreated primary human mammary epithelial cells (HMECs) with IL33 and observed results\nconsistent with basal epithelial cells from the mouse mammary gland. Specifically, IL33 enhances\norganoid formation capacity, restricts luminal differentiation, and promotes the expression of the\ncancer stem cell marker CD4420 (Fig. 7d-k). In summary, we believe that our model system is\nhighly relevant to studying the effects of aging and pregnancy on stem/progenitor cells in the\nmammary gland.\n2. There have been prior related studies in mice that compared retired breeders (old parous) and old\nvirgin mice and also virgin mice and different stages of pregnancy/lactation. For example:\nhttps://pubmed.ncbi.nlm.nih.gov/25619437/ and scRNA-seq profiling (cited by the authors):\nhttps://pubmed.ncbi.nlm.nih.gov/34187545/\nWe thank the reviewer for the comment. As mentioned above, these studies were done at\nan early age (28-36 weeks, 9 months in mice, ~30-38 human years), which is not post-\nmenopausal. Moreover, the first study cited by the reviewer was done using bulk-RNA sequencing\nobscuring cellular heterogeneity21 and second one was done immediately post-involution where\nthere is an increase in breast cancer risk22. Therefore, while these studies have informed our\nanalysis, they do not address our main objective: to characterize the molecular and cellular effects\nof pregnancy on the aged mammary gland. A more detailed description of both studies is given\nbelow.\nHuh et al.21 interrogated three experimental groups: virgin (3, 9, 24wk), pregnant (D10, 16,\n19), and retired breeder (28-36wk). Their analysis incorporated bulk DNA methylation and RNA\nsequencing, obscuring cellular heterogeneity and identification of minority cell populations.\nMoreover, they discovered “that pregnancy had the most significant effects on\nCD24(+)CD61(+)CD29(hi) and CD24(+)CD61(+)CD29(lo) cells, inducing distinct epigenetic\nstates that were maintained through life”. This conclusion is incomplete, as their oldest point\nsampled is 36wk/9m, which translates to ~30-38 human years. Given that most breast cancer\npatients receive a diagnosis over a decade later, we feel that our time point of 18m (~56 years in\nhumans) is a more comprehensive approach to address our objective.\nThe second study Pal B et al.22 explore the developing mammary gland using scRNA-seq\nof 8 key stages spanning E18.5 to 10wk virgin and 3wk post-involution. This narrow time frame\nprecludes any interrogations into the effects of aging on the mammary gland, which is a key\ncomponent and source of novelty for our study. In our study we have used this dataset to show\nthat Il33+ and Krt6a+ cells are present in the terminal end bud during development (Fig. 3h), but\nthe proportion of these cells decreases in the adult and post involution (Fig. 3h).\nAs requested by the reviewer in point #4, we also integrated the data from Pal B, et al.22\nwith our data to directly compare all the groups. We find that while there is a decrease in the\nproportion of Il33+/Krt6a+ cells right after pregnancy, the population is very small in young mice\nto begin with, it increases with age and is reduced in aged parous mice (Fig. 3g, Fig. 4a, b).\n3.Similarly, the increase in hybrid basal-luminal cells with age has been described by multiple groups,\nboth in human and mice, and IL33 has been shown to promote cancer stem cells in breast cancer.\nWe thank the reviewer for their comment. While the effects of aging on loss of lineage\nintegrity in luminal cells have been described23, no study has assessed how early pregnancy\naffects the aging process or the mechanisms leading to loss of lineage integrity with age.\nWe for the first time show that age induced accumulation of hybrid cells located in the\nbasal layer and demonstrate that a previous pregnancy reduces this accumulation. We also show\nthat IL33 mimics the age-induced loss of lineage integrity in young basal cells. Moreover, in a\nprevious study23 Vatter et al. propose that luminal cells acquire basal features whereas our study\nsuggests that basal cells could acquire luminal features. This is an important distinction because\nwhile it has been thought that luminal cells are the cell of origin for breast cancers recent studies\nhave showed that basal cells can be the cell of origin24,25.\nRegarding IL33, its role in normal mammary stem/progenitor cells, cellular plasticity and\ntumor initiation have never been studied. Previous studies have primarily examined the role of\nIL33 from the microenvironment on tamoxifen resistance and metastasis26,27. Therefore, our\nfinding that 1. in a minority population of cells, IL33 expression increases with age 2. is reduced\nin aged parous mice and 3. regulates growth and differentiation of normal mammary epithelial\ncells to our knowledge have not been described before.\n4. The authors have to include 3-month-old parous mice in their flow cytometry and RNAseq analyses\n(Figure 1-3) and also 3-month-old virgin and parous to see how quickly the differences due to pregnancy\nare observed. This is important, because the authors argue that their study is unique, since they looked\nat old parous mice compared to prior studies that analyzed mice shortly after pregnancy. But, w/o showing\ndata for 3-month-old virgin and parous this is not possible to conclude.\nWe thank the reviewer for this comment. As per the reviewer’s suggestion, we have now\nrevised Fig. 1 & 3 to include a 3-month-old parous mice (mice underwent a minimum of two\npregnancies and were analyzed at least 3 weeks post-involution).\nIn brief, we repeated all the experiments in Fig. 1 with a 3-month-old parous (3M P) mouse\nside by side with 3M NP, 18M NP and 18M NP. Our experiments characterize the mammary gland\nat these different stages through flowcytometry and organoid assays, which has never been done\nbefore. We find that parity has distinct effects on the aging of luminal and basal cells – 1.\nLuminal cells from parous mice do not undergo age related loss that is seen in nulliparous mice 2.\nParity restricts basal cell expansion that is seen with aging, but it is not a complete rescue. We also\nperformed organoid assays with the 4 groups to find that while aging increase organoid formation\ncapacity of basal cells, it reduces the organoid forming capacity of luminal cells.\nFor Fig 3., we were able to compute the fraction of Krt6a+/Il33+ cells across the 4 different\ngroups as described above and included in the revised manuscript (Fig. 3g). We find that while\nthere is a decrease in the Krt6a+/Il33+ cells immediately post-pregnancy, this percentage is\nsignificantly lower in 3M NP and 3M P mice (Fig 3g). Thus, our data suggests that accumulation\nof the Krt6a+/Il33+ cells occurs with aging in nulliparous mice. Moreover, we have confirmed this\nanalysis using immunofluorescence (Fig. 4a-c).\nFor the previous Fig. 2, we integrated the Pal et al. dataset22 with ours and performed\npseudo bulk analysis across all the groups (young nulliparous, young parous, aged nulliparous\nand aged parous). We have included this analysis below as part of the point-by-point response\n(Fig. R1).\nFirst, we confirmed successful isolation of basal and luminal cell populations from the\npseudobulk analysis by expression of Krt14 and Acta2 in basal samples and Krt8 and Krt18 in\nluminal samples (Fig. R1a). To identify transcriptomic differences between nulliparous and parous\nMECs as they age, we performed differential gene expression analyses on aged basal or luminal\ncells versus young basal or luminal cells from either parous or nulliparous mice (i.e., 18M NP vs.\n3M NP and 18M P vs. 3M P) and identified differentially expressed genes with opposing log2 fold\nchanges (Fig. R1b). Basal cells from parous mice displayed increased expression of the basal\ncell lineage marker Lipg28 and ECM-related genes (Ltbp2, Col4a1, Fig. R1c, d) with age, whereas\nbasal cells from nulliparous mice showed a decrease in expression with age. In agreement with\nprevious studies29, basal cells from nulliparous mice had an age-dependent increase in\nexpression of Arg1. However, we find that basal cells from parous mice showed an opposite trend,\nwith Arg1 being highly expressed in basal cells from young parous mice and returning to levels\nsimilar to young nulliparous basal cells in aged parous basal cells (Fig. R1c, d). Similarly, the\nexpression of Tspan8, which has been previously reported to mark a population of quiescent\nhormone responsive stem-cells30, increased with age in nulliparous basal cells but this was\nreversed in parous (Fig. R1c, d). We had independently analyzed and confirmed the age- and\nparity-dependent trend in Tspan8 expression at the protein level with flow cytometry analysis\n(Extended Data Fig. 3).\nThe number of age-specific differentially expressed genes (DEGs) with opposing log2 fold\nchanges between parous and nulliparous mice was far greater in luminal cells compared to basal\ncells (Fig. R1e), underscoring the dramatic impact of pregnancy on luminal cell transcriptomes.\nAs was observed in basal cells, Arg1 expression increased with age in nulliparous mice, but\ndecreased with age in parous mice (Fig. R1e, f). Other genes showing similar expression patterns\nwith age and parity include Kitl, Sparc, Serpina3n, and Igfbp3, which are often associated with\nthe basal cell lineage28, as well as the basal progenitor cell marker Cd3631 (Fig. R1e, f).\nOur findings support the hypothesis that MECs display a loss of lineage integrity over the\ncourse of aging and suggest that previous pregnancies attenuate this age-dependent phenotype\nwhile instilling a persistent transcriptomic memory of pregnancy. However, we feel this analysis\nis largely descriptive as pointed out by Reviewer 2 and tangential to the main point of the paper.\nWe have therefore removed it from the main manuscript.\n[FIGURE REDACTED]\nFigure R1: Parous and Nulliparous MECs display distinct transcriptomic patterns with\naging. (a) Expression of basal (Krt14, Acta2) and luminal (Krt18, Krt8) lineage markers in bulk\nRNA-seq samples collected from FACS-sorted basal and luminal cells of 18M NP (red) and 18M\nP (blue) mammary glands. 18M bulk RNA-seq samples were integrated with 3M NP (grey) and\n3M P (green) scRNA-seq datasets (GSE103275) for pseudo-bulk RNA-seq analyses. (b)\nSchematic diagram describing the strategy used to identify differentially expressed genes with\naging (18M NP vs. 3M NP and 18M P vs. 3M P) and exhibit opposing log2 fold changes (Log2FC).\n(c) Differentially expressed genes (DEGs) in 18M NP basal vs. 3M NP basal (pink) and 18M P\nbasal vs. 3M P basal (blue) with opposing Log2FC values. DEGs were filtered for p_adj. < 0.05\nand |Log2FC(NP) – Log2FC(P)| > 1.5. (d) Expression values of select genes identified in (c)\nacross comparison groups. (e) Differentially expressed genes (DEGs) in 18M NP luminal vs. 3M\nNP luminal (pink) and 18M P luminal vs. 3M P luminal (blue) with opposing Log2FC values. DEGs\nwere filtered for p_adj. < 0.05 and |Log2FC(NP) – Log2FC(P)| > 1.5. (f) Expression values of\nselect genes identified in (d) across comparison groups. Gene expression is represented as\ntranscript per million (TPM) × 103 (a) or normalized values from DESeq2 (d, f). n = 3 mice per\ngroup\n5. Similarly, in Figure 4 the authors argue that IL33 treatment of organoids from young mice phenocopies\naging. But they do not analyze organoids from old mice and also do not compare 3-month-old virgin vs.\nparous. These additional samples would be critical to include.\nWe thank the reviewer for their comment. As suggested by the reviewer we have included\nall the 4 groups using both the mammary glands from mice (Fig. 4a-c) and on organoids derived\nfrom them (Fig. 4d-g).\nConsistent with the single cell RNA-sequencing data, we find that the 18M NP mice have\nthe highest percentage of KRT6A+ cells and that this percentage is reduced in 18M P mice.\nMoreover, we find that there is a slight but significant reduction of KRT6A+ cells in the young\nparous mice (Fig. 4a, b). For the in vitro culture, we observe the increase in KRT6A+ cells in the\norganoids derived from basal cells (Fig. 4d, e) but do not see any differences in organoids derived\nfrom luminal cells (Fig. 4f, g) in 18M NP mice.\nReviewer #2 (Remarks to the Author):\nThis study by Olander et al. investigates the combined effects of aging and pregnancy on mouse\nmammary epithelial cells (MECs). The authors report that pregnancy reduces the age-induced expansion\nof basal cells and partially normalizes the basal:luminal cell ratio in aged mice. They also report the\npresence of a small population of Il33-expressing hybrid MECs, which accumulate in aged nulliparous\nmice, but reduced in aged parous mice. The authors show that pregnancy induces transcriptomic\nchanges in both basal and luminal cells, promoting a differentiated state and reversing some age-\nassociated transcriptional programs. Finally, the authors show that IL33 treatment of young basal cells\nmimics the aged nulliparous MECs, increasing the proportion of basal cells and inducing a KRT6a+ hybrid\nstate, which may confer a survival advantage for cells with oncogenic mutations. Overall, the study\nsuggests that pregnancy-induced changes in MECs contribute to long-term breast cancer protection by\nmaintaining lineage integrity and reducing the accumulation of potentially tumorigenic hybrid cells. The\nstudy is well presented and written. Generally, there is a slight over interpretation of some of the data\nwhich should be substantiated with further evidence. Below are the main points this reviewer believes\nshould be addressed.\nWe thank the reviewer for their assessment of our manuscript and constructive feedback.\nWe have addressed the concerns of the reviewer below by including new data.\n1) Title: I personally feel that the title might be slightly overreaching, as the paper does not explicitly\ndemonstrate that pregnancy reduces the formation of the hybrid progenitor. The findings just show that\nthe hybrid progenitor is identified in the old nulliparous samples but is less prevalent in the old parous\nsamples.\nWe appreciate the reviewer for raising this concern. In view of the reviewer’s concern, we\nhave revised our title to read, “Characterization of the mammary gland with pregnancy and aging\nidentifies a population of IL33+ hybrid progenitors."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have responded to each of the points raised by the reviewers and included additional experiments to strengthen the study. The revised manuscript is improved, but still there are some remaining concerns.\nThe single cell RNA-seq analyses are still suboptimal, there are much more sophisticated data analysis and presentation tools than UMAP and bar plots.\nMany bar plots lack statistical tests and p-values (like Fig 3 g,h), also it’s better to indicate actual p-values instead of asterisks for significance.\nThe in vivo experiment with IL33 (new Figure 6) is useful, but it raises several new issues:\n(1) IL33 acts on many other cell types, not just epithelial cells, so are the observed changes through IL33 directly acting on epithelial cells or changing the microenvironment?\n(2) The authors only treated 3 mice and only for 5 days – this is too few mice to ensure reproducibility and significant and they need an explanation for how only 5 days of treatment has such dramatic effects on mammary ductal morphology. Was this due to proliferation or apoptosis? The authors should test this.\n(3) similarly, in the scRNAseq data (Fig 7 1-b) it would be important to plot the Ki67+ cells, since the authors claim that IL33 induces aberrant proliferation.", "reviewer_2": "I thank the authors for addressing the questions raised. The manuscript is much improved.\nTwo minor comments:\n1) The authors should use the integrated HBCA (Reed et al 2024 Nat. Genetics) to assess the expression of IL33. Its a much larger dataset which integrated the 7 main human scRNAseq studies (2.1million cells). They should also show if these cells are represented more or less in any of the datasets.\n2) The discussion starts with \"Our study resolves and comprehensively maps the complex relationship between aging and pregnancy in the mammary gland, such as basal-cell bias and differentiation.\" I don't believe the current study fully resolves or comprehensively maps this relationship. I think this needs to be toned down to reflect the data presented.", "reviewer_3": "The authors have been responsive to the critiques and this revision has improved many aspects of the previous data presented. Figs 1-5 are excellent and present a great coherent story.\nIn response to reviewers asking for in vivo confirmation of IL33 effect, and for human validation, they added figures 6 and 7.\nHowever, these two new figures caused more concerns than solving problems.\nIL-33-induced changes in mammary glands in mice at day 1 post completion of 5 days of treatment are not consistent with the cell culture and organoid data presented before this new figure 6 - KRT6A expression was induced in luminal cells but not basal cells but the reserve was observed in the in vitro experiments. More convincing explanations should be provided for this inconsistency. Furthermore, the authors did not detect an impact of IL-33 treatment at day 8, calling into question about the critical role of IL-33 in mediating the ageing effect on expanding the hybrid cell population. The authors should explain why the lack of detecting a more long-term change did not weaken their central conclusions.\nHuman scRNA data did not detect a hybrid cell population that is consistent with their mouse studies, and they did not detect any impact of ageing or parity on this cell population, calling into the question of whether their analysis was appropriate or proper datasets were used. In the absence of more convincing data, this figure probably should be dropped."}, "author_response": "Point by point response to reviewer comments.\nOriginal comments are in blue, authors replies are in black.\nReviewer #1 (Remarks to the Author):\nThe authors have responded to each of the points raised by the reviewers and included additional\nexperiments to strengthen the study. The revised manuscript is improved, but still there are some\nremaining concerns.\nWe thank the reviewer for their thoughtful feedback, which has helped us further strengthen the\nmanuscript.\nThe single cell RNA-seq analyses are still suboptimal, there are much more sophisticated data analysis\nand presentation tools than UMAP and bar plots.\nWe appreciate this comment and note that we have applied several sophisticated and state-of-\nthe-art bioinformatics tools to analyze both newly generated and publicly available scRNA-seq\ndata in this work. These include pseudotime trajectory inference methods such as PAGA and\nSlingshot, as well as the recently developed potency prediction tool, CytoTRACE 2 (PMID:\n38562882), which leverages interpretable deep learning and is in press at Nature Methods (A.\nNewman, personal communication). Additionally, by employing UMAPs and bar plots (as well as\nviolin plots and dot plots) for summarizing and visualizing complex scRNA-seq data, we followed\nbest practice in the field, as exemplified by guidelines and display/analysis options in Scanpy and\nSeurat. However, we are happy to implement specific analyses if the reviewer has suggestions.\nMany bar plots lack statistical tests and p-values (like Fig 3 g,h), also it’s better to indicate actual p-values\ninstead of asterisks for significance.\nWe thank the reviewer for raising this point – we have included the statistics in the Fig. 3g, h, i\nand have indicated p-values on all figures instead of asterisks for significance. For Fig. 3g, h, i,\nwe used the Chi square test as the we have represented the fraction of IL33+/KRT6A+ cells over\nthe total cells in the dataset. We had to use this method as some of the published datasets\ndeposited in GEO do not have data on replicate samples. Moreover, we have performed\nimmunofluorescence staining in independent samples for KRT6A (see Fig. 4a, b). We show that\nthere is an accumulation of KRT6A+ cells in the aged nulliparous mice, confirming the single cell\nfindings in Fig. 3g at the protein level.\nThe in vivo experiment with IL33 (new Figure 6) is useful, but it raises several new issues:\n(1) IL33 acts on many other cell types, not just epithelial cells, so are the observed changes through IL33\ndirectly acting on epithelial cells or changing the microenvironment?\nWe thank the reviewer for raising this important point. Indeed, we do see an increase in the\nstromal population and increase in ST2+ immune cells with IL33 treatment (Extended Data\nFigure 9c, d). It is entirely plausible that these IL33-dependent changes within the\nmicroenvironment can in turn impact on the epithelial populations. However, our in vitro studies\n(Fig. 5) with mammary epithelial cells demonstrates that IL33 increases colony formation and\ninduces the proliferation and expression of KRT6A in basal epithelial cells. Moreover, we see an\nincrease in EdU+ cells in basal organoids treated with IL33, further suggesting that IL33 is directly\naffecting basal cell proliferation (Fig. R1). Therefore, our data support the notion that the effect of\nIL33 seen in vivo is due to its direct action on epithelial cells. Nevertheless, we cannot rule out\nadditional cross talk between the microenvironment and epithelial cells in vivo. We have now\nadded this part to our discussion (line 387-389).\nFigure R1: IL33 treatment directly enhances proliferation in basal organoids, but not luminal organoids. (a)\nRepresentative immunofluorescence images of basal organoids treated with (top panel) or without (bottom panel) IL33\n(10ng/mL). Organoids were generated from either luminal or basal cells sorted from 8-week-old nulliparous mice.\nOrganoids were cultured for 9 days before the addition of EdU (10uM). After a 24hr chase, organoids were harvested.\n(b). Quantification of EdU+ cells in basal and luminal organoids exposed to IL33. Data represent the total number of\nEdU+ nuclei per total nuclei in a field of view (FOV). Statistical significance was determined by performing multiple\nunpaired t-test with Holm-Šídák correction.\n(2) The authors only treated 3 mice and only for 5 days – this is too few mice to ensure reproducibility\nand significant and they need an explanation for how only 5 days of treatment has such dramatic effects\non mammary ductal morphology. Was this due to proliferation or apoptosis? The authors should test this.\nWe thank the reviewer for raising this point. We have now expanded our analysis to include\nadditional mice (n= 6, see revised Fig. 6). We see the same consistent changes in ductal\nmorphology at 1 day post treatment (1dpt). To determine whether IL33 treatment induces\nproliferation or apoptosis in vivo, we performed EdU and Annexin V analysis (Revised Fig. 6,\nExtended Data Fig. 10, and Fig. R2). We find that IL33 treatment induces proliferation of both\nluminal and basal cells at 1 dpt, but proliferation is normalized at 8dpt (Fig. 6h, Fig. R2f, g).\nAt 1dpt and 8dpt there are small but significant differences in apoptosis in the basal compartment\n(Extended Data Fig. 10). Specifically, we observed a decreased percentage of cells in early\napoptosis (DAPI-/AnnexinV+) in basal cells at both 1 and 8dpt. While insignificant at 1dpt (p =\n0.06), we also observed a slightly elevated live population (DAPI-/AnnexinV-) in IL33-treated\nmice. These results could suggest a pro-survival program enacted by IL33 in basal cells. Despite\nobserving significant differences in cell death within the luminal population between IL33- and\nPBS-treated mice, the differences were inconsistent and often conflicted between timepoints.\nGiven the relatively minor changes in apoptosis, we speculate that the effects of IL33 are largely\ndue to changes in proliferation and not cell death.\nFrankly, we were also surprised by the dramatic effects of IL33 on the mammary gland in such a\nshort time. We therefore performed an additional timepoint where we treated mice with IL33 or\nPBS for only 2 days, administered a single dose of EdU on the second day, and analyzed on the\nthird day (Fig. R2a, short exposure or “S.E.”). We find that even short treatment with IL33 induces\nsignificant changes in ductal morphology (Fig. R2b-e) and proliferation of mammary epithelial\ncells (Fig. R2f-g), consistent with the 1dpt cohort.\nFigure R2: Two day-exposure to IL33 is sufficient to induce altered ductal morphology and epithelial\nproliferation. (a) Schematic of experimental approach. Mice in the short exposure (S.E.) cohort were treated with IL33\nfor two consecutive days, pulsed with EdU on day 2, and sacrificed for analysis at day 3. (b) Representative images of\nCarmine-Alum stained mammary gland wholemounts across timepoints and treatment cohorts. (c) Quantification of\nwidth (μm) of the main epithelial duct in PBS- and rIL-33-injected mammary glands from S.E., 1dpt and 8dpt timepoints\n(n = 3 mice/condition, 5-10 measurements per mouse). (d) Representative images of PBS- and rIL-33-injected H&E\nstains from S.E., 1dpt and 8dpt timepoints. (e) Lumen area (μm2 x 103) in H&E stained sections of PBS- and rIL-33-\ninjected mammary glands from S.E., 1dpt and 8dpt timepoints (n = 3 mice/condition, 20 measurements per mouse). (f)\nRepresentative flow plots depicting EdU+ gates for basal (left two columns) and luminal (right two columns) cells across\ntreatment groups (PBS and IL33) and timepoints (S.E., 1dpt, and 8dpt). (g) Quantification of EdU+ basal (top) and\nluminal (bottom) cells across treatment groups (PBS and IL33) and timepoints (S.E., 1dpt, and 8dpt). Statistical\nsignificance was determined by performing 2-way ANOVA with Šídák’s multiple comparisons test.\nWhile these results are surprising and demonstrate a rapid response to IL33, we speculate that\nthese observations reflect IL33’s function as an “alarmin” to promote proliferation of epithelial cells\ninitiated by cellular stress or damage, as observed in Wu et al., 2021 (PMID: 34343135). Indeed,\nWu et al. found that, as basal epithelial stem cells activate expression of Il33 to proliferate in\nresponse to lung damage, Krt6a promoter activity is enhanced—an observation that was not\nrecapitulated in Il33-/- mice. These findings are in line with our work and illustrate an important,\nfunctional relationship that ties IL33, epithelial proliferation, and Krt6a expression together. Our\nfuture work will focus on understanding how IL33 induces these changes in the mammary gland\nepithelium.\n(3) similarly, in the scRNAseq data (Fig 7 1-b) it would be important to plot the Ki67+ cells, since the\nauthors claim that IL33 induces aberrant proliferation.\nWe have included UMAPs visualizing MKI67 and TOP2A expression (Fig. R3a and Extended\nData Fig. 11b). Interestingly, we do find that Cluster 5 has a high proportion of MKI67+ cells and\nTOP2A+ cells relative to all other clusters. However, it is important to note that in our aged\nnulliparous mice the hybrid cluster has no differences in cell cycle stage when compared with\nother epithelial clusters (Extended Data Fig. 5b).\nFigure R3: Expression of proliferation-related genes in IL33+ human breast epithelial cells. (a) UMAPS showing\nexpression of IL33, MKI67, and TOP2A in the HBCA (Kumar et. al) dataset.\nReviewer #2 (Remarks to the Author):\nI thank the authors for addressing the questions raised. The manuscript is much improved.\nWe are delighted to hear the reviewers kind support and welcome their continued feedback.\nTwo minor comments:\n1) The authors should use the integrated HBCA (Reed et al 2024 Nat. Genetics) to assess the expression\nof IL33. Its a much larger dataset which integrated the 7 main human scRNAseq studies (2.1million cells).\nThey should also show if these cells are represented more or less in any of the datasets.\nWe thank the reviewer for suggesting this dataset. We have now analyzed 889,281 epithelial cells\nwithin the iHBCA dataset (Fig. R4a) (Reed et al., E-MTAB-13664). Overall, we find that there is\na population of hybrid cells that is enriched in expression of IL33 (Fig. R4b). Specifically, this\ndistinct hybrid cluster is enriched in the Kumar et al. dataset that expresses both basal and luminal\nmarkers (Fig. R4b, black dashed circle).\nFigure R4: Analysis of iHBCA. (a) UMAP of iHBCA annotated by dataset (left) and cell type (right). (b) Expression of\nKRT18, KET5, ELF5 and IL33 in iHBCA. IL33+ cluster is highlighted with a dashed circle. (c) Percentage of IL33+ cells\nacross all datasets (see also Table 1). (d) Percentage of IL33+ cells in each dataset split by age (see also Table 1). (e)\nPercentage of IL33+ cells in the Kumar et al. dataset split by age and parity. Please note that age at first pregnancy,\nBRCA1/2 status are not available in the metadata. Aged >50 years, Mid <50 years but >35 years, Young <35 years.\nMoreover, in the Twigger dataset, there is also a population of IL33+ cells that only expresses\nluminal markers (Fig. R4b, pink dashed circle). As the Twigger dataset was generated using cells\nisolated from human milk, we speculate that the expression of IL33 could reflect unique biology\nassociated with lactation and immunity. Although a clear cluster wasn't visible in the other\ndatasets, we were able to identify IL33+ cells in all of them (Fig. R4c). Therefore, we explored\nwhether the percentage of IL33+ epithelial cells and/or expression of IL33 changes with age in\nthe iHBCA. Overall, we find that the percentage of IL33+ cells increases with age in all datasets\nexcept the Twigger dataset (Fig. R4d).\nThere are several important caveats to note regarding our analysis with human samples. Our\nmanuscript explores the impact of early first pregnancies on aging. The available human single\ncell RNA-sequencing datasets do not include the age at first pregnancy (Table 1), which can have\nan impact on the hybrid cell population. In the only dataset that includes age at first pregnancy\n(Reed et al.), there is only one patient with an early pregnancy that does not have mutations in\nBRCA1/2 (Table 1). Therefore, we cannot assess the combined effects of early pregnancy and\naging in any of the datasets in a rigorous manner. Nevertheless, we used the Kumar et al. dataset\nto understand whether the percentage of IL33+ cells changes with aging and pregnancy as this\ndataset has the greatest number of samples in the aged (>50 years of age) group (Table 1),\nalthough we do not have information of age at first pregnancy or BRCA1/2 status. We find that\nthere is a significant increase in the percentage of IL33+ cells with age (p value = 0.03), but there\nis no significant difference in the aged nulliparous and parous samples (Fig. R4e).\nTable 1. Metadata for the seven datasets within iHBCA.\n# of parous patients # of nulliparous patients\nFresh/ # of % of Age at first\nDataset >50 with no reported >50 with no reported\nFrozen cells IL33+ pregnancy\nmutations in BRCA1/2 mutations in BRCA1/2\nKumar Fresh 230298 1.49% unknown 12 9\nPal Fresh 75239 2.6% unknown 0 2\nGray Fresh 30677 1.5% unknown 5 1\nReed Frozen 343561 0.9% known 6 1\nNee Frozen 108184 0.4% unknown 2 0\nMurrow Frozen 53071 0.2% unknown 0 0\nTwigger Frozen 48251 9.7% unknown 3 0\nImportantly, we do not yet understand what other factors may affect the accumulation of this\nhybrid population. These could include race, ethnicity, BMI and mutation in BRCA1/2. Indeed, our\nanalysis presented in Extended Data Figure 11d demonstrates an over-representation of IL33+\nhybrid MECs in Asian and African American donors. Preliminary analysis of BRCA1 mutant breast\ntumors also shows an increase in the IL33+ cells (data not shown). These factors/variables also\ninfluence breast cancer risk (American Cancer Society. Breast Cancer Facts & Figures 2024-\n2025. Atlanta: American Cancer Society; 2024). It is therefore possible that our findings may not\ntranslate into humans within certain contexts and/or populations.\nFinally, although the iHBCA has generated an immensely useful resource for the community, four\nof the seven datasets in the iHBCA used frozen cells for sequencing (Table 1). Moreover, there\nwere diverse methods used for tissue digestion, cell prep and sequencing (e.g., Reed et.al dataset\nused frozen, partially digested tissues). This can significantly affect the representation of rare cell\ntypes in the population. We have now included a paragraph highlighting the limitations of our\nanalysis in human samples in the discussion section (lines 369-383).\n2) The discussion starts with \"Our study resolves and comprehensively maps the complex relationship\nbetween aging and pregnancy in the mammary gland, such as basal-cell bias and differentiation.\" I don't\nbelieve the current study fully resolves or comprehensively maps this relationship. I think this needs to\nbe toned down to reflect the data presented.\nWe sincerely appreciate the reviewer’s suggestions. We have now modified the sentence to read,\n“Our study lays the groundwork for understanding the complex relationship between aging and\npregnancy in the mammary gland, such as basal-cell bias and differentiation”. We have also\nrevised the discussion section to more accurately reflect the findings, scope, and limitations of\nour study.\nReviewer #3 (Remarks to the Author):\nThe authors have been responsive to the critiques and this revision has improved many aspects of the\nprevious data presented. Figs 1-5 are excellent and present a great coherent story.\nWe thank the reviewer for their enthusiasm towards our study and welcome their thoughtful\nfeedback.\nIn response to reviewers asking for in vivo confirmation of IL33 effect, and for human validation, they\nadded figures 6 and 7. However, these two new figures caused more concerns than solving problems.\nIL-33-induced changes in mammary glands in mice at day 1 post completion of 5 days of treatment are\nnot consistent with the cell culture and organoid data presented before this new figure 6 - KRT6A\nexpression was induced in luminal cells but not basal cells but the reserve was observed in the in vitro\nexperiments. More convincing explanations should be provided for this inconsistency.\nWe thank the reviewer for raising this important point. We acknowledge that the data from aged\nm\naged mice and in vitro data do not completely translate to our in vivo experiments in young mice.\nWe hypothesize that this difference is because mammary basal cells in vivo are largely non-\nproliferative and do not give rise to luminal cells under homeostatic conditions (PMID: 26511661,\nPMID: 30271991"}, {"round": 3, "reviewer_comments": {"reviewer_1": "The authors have responded to each of the remaining comments of the reviewers by revising the manuscript, including adding some additional experiments.\nHowever, the dramatic increase in duct size after just 1 day of IL33 administration is not resolved conclusively (Fig 1b, f). The authors checked proliferation and apoptosis and found that IL33 mainly increases proliferation. However, to see such a dramatic (at least 2-fold) increase in duct size in 24 hrs would mean that virtually all cells in the duct double in 24 hrs – which is not matching their EdU labeling data showing that 10-20% of the basal and luminal cells have incorporated EdU in ~6 days. There is clearly an increase in proliferation, but the magnitude and time course of this does not match the dramatic phenotypic change that occurs in a day. Have the authors considered alternative explanations like swelling of the ducts – checked cell size and potential fluid accumulation?", "reviewer_2": "Thank you for addressing my concerns. The analysis of the iHBCA presented in the rebuttal is informative and should be part of the main paper.\nI have no further comments.", "reviewer_3": "The authors have been responsible to the reviewer’s comments. The revised manuscript now include additional data and revised presentation to address the inconsistency between in vitro finding and in vivo validation and bioinformatic studies of published human scRNAseq data. Overall, the manuscript has strengthened, but the inconsistency moderates the significance of the study. Additional concerns are detailed below:\n1. Fig 2b: the proportion plot shows the 18M P group has a lower proportion of all cell types. If the same total cell numbers were analyzed in 18M P vs 18M NP, reduction of some cell groups should be accompanied by decrease of others. In this study, nearly twice as many cells were used in the P group compared to NP group (NP, red, 3,327 cells vs P, blue, 7,305 cells), to this reviewer, it is mathematically impossible that the 18M P has lower counts for all cell subsets. As a result, all the statements regarding this fig are a bit shaky. To help readers, the authors should ensure that the same number of cells were used for this comparison in Fig 2A and 2B. Also since 3 samples are used, intra-sample variations should be shown.\n2. Fig 2a/b: The integration with Tabula Muris (red, n = 2,993 cells) is confusing. If more cells were indeed needed, the authors should sequence more cells from their own experiments. Integrated analysis should be secondary not primary, and readers should be convinced that the public data by others was comparable.\n3. Fig 4a/b: this figure is used to show KRT6A+ cells are switched from the luminal layer to the basal cell layer with ageing. While the 18M NP panel shows convincing basal location of KRT6a+ cells, the 3M NP panel does not show convincing luminal cell staining for KRT6a.\n4. Fig 5H: the organoid culture of Tp53 WT control basal cells does not show response to IL33. This seems to be inconsistent with Fig 5A.\n5. Line 262-264: The authors state that “Wholemount analysis of the mammary glands from IL33-treated mice showed a significantly increased lymph node area and dramatic dilation of the ducts at 1dpt (Fig. 6b, c).” But the phrase “increase lymph node area” is confusing. Do the authors mean increased ducts around the lymph node area? If so, an explanation is needed for this difference around the LN area compared to other parts of the mammary glands.\n6. Line 266-268: “Aside from LN area, analysis of mice at 8dpt showed no significant changes in any of the parameters above, suggesting the effects of short term IL33 treatment were transient (Fig. 6b-f).” This sentence does not make sense. Do the author mean changes are seen around the LN area, but not other parts of glands? If so, explanations are again needed.\n7. Fig 6I and J: IL33 treatment for 1 day (but now 8 days) led to a change of the ratio of basal vs. luminal cells. But this potentially important observation is not stated in the results section nor discussed."}, "author_response": "Point by point response to reviewer comments.\nOriginal comments are in blue, authors replies are in black.\nReviewer #1 (Remarks to the Author):\nThe authors response to the remaining comment is fine. Could they also add the statement that expansion of the ducts due to swelling cannot be excluded to the manuscript also when they are describing the data? the enlargement of the ducts is a very dramatic phenotype (even if it was 48 hrs not 24 hrs treatment) so I'm sure others will also question this.\nWe thank the reviewer for the feedback. As requested, we have included the statement, “H&E staining confirmed that there was a significant increase in size of the lumen in IL33 group at 1dpt (Fig. 6e, f) but did not reveal fluid accumulation within the ducts; however, we cannot entirely exclude the possibility of dilation secondary to edema” to the manuscript (lines 267-269).\nWe have also included the cell size analysis as part of Extended Data Fig. 9c-e and in the manuscript (lines 287-290).\n***\nThe responses to reviewer 3 and revision seem OK to me, although I do not see the figures to make sure resolution is publishable quality- but this you can check during QC.\nWe thank the reviewer for their comments. We have now uploaded high resolution Illustrator files that are publishable quality."}, {"round": 4, "reviewer_comments": {"reviewer_1": "The authors response to the remaining comment is fine. Could they also add the statement that expansion of the ducts due to swelling cannot be excluded to the manuscript also when they are describing the data? the enlargement of the ducts is a very dramatic phenotype (even if it was 48 hrs not 24 hrs treatment) so I'm sure others will also question this.\n***\nThe responses to reviewer 3 and revision seem OK to me, although I do not see the figures to make sure resolution is publishable quality- but this you can check during QC."}, "author_response": "Point by point response.\nOriginal comments are in blue, authors replies are in black.\nChanges are highlighted in green in the revised manuscript.\nPlease note: Figures included in the revised manuscript paper are referenced in the point-by-point\nresponse but not included in this document. Only figures specific to the point-by-point response are\nincluded here.\nReviewer #1 (Remarks to the Author):\nIn this study Olander and colleagues have investigated cellular and molecular changes induced by\npregnancy in mice. They analyzed mammary glands of 3-, 8-, and 18-month-old virgin or multiparous\nmice by flow cytometry to assess the relative proportion of luminal/basal cells and also performed RNA-\nseq on FACS sorted basal and luminal cells. Subsequently, they also performed single cell RNA-seq on\nmammary glands from 18-month-old virgin and multiparous mice. Lastly, they tested the effect of IL33, a\ngene they identified as differentially expressed in mammary epithelial cells 18-month-old virgin and\nmultiparous mice. Overall, the manuscript has limited data and novelty.\nWe thank the reviewer for their comments and assessment of our manuscript. We have\nsignificantly expanded the scope of our work, adding several new figures in the main and\nextended data as per the suggestions from all the reviewers. We believe that our revised work\nprovides an analysis of aging and pregnancy that has not been done before.\nSpecific points:\n1. Full term pregnancy in early adulthood only reduces the risk of postmenopausal ER+ breast cancer,\nwhile the risk of triple-negative breast cancer (TNBC) is actually increased. Mice do not get ER+\nmammary tumors, thus, they may not be the best model to use to assess the effect of pregnancy on ER+\nbreast cancer risk.\nWe appreciate the reviewer's insightful comment and welcome the opportunity to clarify\nour rationale. The goal of our study was to investigate how pregnancy influences the aging of\nstem/progenitor cells in the mammary gland, specifically examining whether certain subsets are\nenriched or lost with aging and pregnancy. Previous studies have long utilized mice as a model\nsystem to study parity-induced protection1,2 and understand underlying mechanisms3–5. More\nrecently studies have used mice to study aging in mammary epithelial cells and their relation to\nbreast cancer risk6–8. Thus, we followed past precedent in using mice as a model system to\nunderstand the long-term impacts of pregnancy.\nAs the reviewer correctly noted, the risk of breast cancer increases in the first 10 years\nfollowing pregnancy in humans. The key limitation of previous studies is that they looked at the\npost-pregnancy mammary gland when there is increased risk and it is not post-menopausal (peri-\nmenopause in mice)3–5. This makes our study unique, as it specifically examines a time point\nbeyond this 10-year window9–12, aiming to better mimic conditions relevant to human biology—\nwhere women who experience an early first pregnancy gain protective benefits much later in life,\ntypically after the age of 50. Moreover, none of the aging studies in mammary epithelial cells\naccount for pregnancy6–8.\nWhile it is true that mice usually get basal like breast cancer there are mouse models of\nER+ breast cancer that have been used to study human disease. Some of these include ESR1,\nCCDN1, prolactin, and WNT1 overexpression, PIK3CA gain of function and loss of STAT1 –\ncomprehensive reviews are here13,14. As one recent example published in Nature Cancer in 2025,\na group from Genentech generated and utilized a mouse model to study human ESR1 mutations\narising due to resistance to ERα-targeted therapies in breast cancer15,16. They showed that the\ncells present in the mammary gland of mice carry the ESR1 mutation that are resistant to\nendocrine suppression have a distinct gene expression pattern that is between luminal and basal\ncells. We have analyzed the single cell RNA-sequencing data from this manuscript15 (Fig. 3i). We\nfind that the cells arising from ESR1 mutant mice after endocrine suppression (labelled as\nendocrine suppressed mutant enriched or ESME cells in the published manuscript) are enriched\nin expression of Krt6a/Il33 (Fig. 3j,k).\nMultiple studies have also employed mouse ER+ breast cancer cell lines, which were\npreviously characterized17 and recently used to understand metastatic dissemination of ER+\nbreast cancer18. Lastly, a recent study showed that the same mutation that generates ER- (basal-\nlike) tumors in young mice gives rise to ER+ tumors in aged mice19. Since many of the earlier\nstudies introducing mutations to generate breast cancer have been done with young mice, it is\npossible that if these experiments were repeated with aged mice, they would give rise to ER+\nbreast cancer.\nFinally, as requested by Reviewer #2, we have extended our findings to human datasets\nand primary mammary epithelial cells. These results are now included as a new figure (Fig. 7) in\nthe main manuscript. Briefly, our findings identify a minority population of cells with a similar gene\nexpression pattern as seen in nulliparous aged mice (Fig. 7a-c, Extended Data Fig. 10). We also\ntreated primary human mammary epithelial cells (HMECs) with IL33 and observed results\nconsistent with basal epithelial cells from the mouse mammary gland. Specifically, IL33 enhances\norganoid formation capacity, restricts luminal differentiation, and promotes the expression of the\ncancer stem cell marker CD4420 (Fig. 7d-k). In summary, we believe that our model system is\nhighly relevant to studying the effects of aging and pregnancy on stem/progenitor cells in the\nmammary gland.\n2. There have been prior related studies in mice that compared retired breeders (old parous) and old\nvirgin mice and also virgin mice and different stages of pregnancy/lactation. For example:\nhttps://pubmed.ncbi.nlm.nih.gov/25619437/ and scRNA-seq profiling (cited by the authors):\nhttps://pubmed.ncbi.nlm.nih.gov/34187545/\nWe thank the reviewer for the comment. As mentioned above, these studies were done at\nan early age (28-36 weeks, 9 months in mice, ~30-38 human years), which is not post-\nmenopausal. Moreover, the first study cited by the reviewer was done using bulk-RNA sequencing\nobscuring cellular heterogeneity21 and second one was done immediately post-involution where\nthere is an increase in breast cancer risk22. Therefore, while these studies have informed our\nanalysis, they do not address our main objective: to characterize the molecular and cellular effects\nof pregnancy on the aged mammary gland. A more detailed description of both studies is given\nbelow.\nHuh et al.21 interrogated three experimental groups: virgin (3, 9, 24wk), pregnant (D10, 16,\n19), and retired breeder (28-36wk). Their analysis incorporated bulk DNA methylation and RNA\nsequencing, obscuring cellular heterogeneity and identification of minority cell populations.\nMoreover, they discovered “that pregnancy had the most significant effects on\nCD24(+)CD61(+)CD29(hi) and CD24(+)CD61(+)CD29(lo) cells, inducing distinct epigenetic\nstates that were maintained through life”. This conclusion is incomplete, as their oldest point\nsampled is 36wk/9m, which translates to ~30-38 human years. Given that most breast cancer\npatients receive a diagnosis over a decade later, we feel that our time point of 18m (~56 years in\nhumans) is a more comprehensive approach to address our objective.\nThe second study Pal B et al.22 explore the developing mammary gland using scRNA-seq\nof 8 key stages spanning E18.5 to 10wk virgin and 3wk post-involution. This narrow time frame\nprecludes any interrogations into the effects of aging on the mammary gland, which is a key\ncomponent and source of novelty for our study. In our study we have used this dataset to show\nthat Il33+ and Krt6a+ cells are present in the terminal end bud during development (Fig. 3h), but\nthe proportion of these cells decreases in the adult and post involution (Fig. 3h).\nAs requested by the reviewer in point #4, we also integrated the data from Pal B, et al.22\nwith our data to directly compare all the groups. We find that while there is a decrease in the\nproportion of Il33+/Krt6a+ cells right after pregnancy, the population is very small in young mice\nto begin with, it increases with age and is reduced in aged parous mice (Fig. 3g, Fig. 4a, b).\n3.Similarly, the increase in hybrid basal-luminal cells with age has been described by multiple groups,\nboth in human and mice, and IL33 has been shown to promote cancer stem cells in breast cancer.\nWe thank the reviewer for their comment. While the effects of aging on loss of lineage\nintegrity in luminal cells have been described23, no study has assessed how early pregnancy\naffects the aging process or the mechanisms leading to loss of lineage integrity with age.\nWe for the first time show that age induced accumulation of hybrid cells located in the\nbasal layer and demonstrate that a previous pregnancy reduces this accumulation. We also show\nthat IL33 mimics the age-induced loss of lineage integrity in young basal cells. Moreover, in a\nprevious study23 Vatter et al. propose that luminal cells acquire basal features whereas our study\nsuggests that basal cells could acquire luminal features. This is an important distinction because\nwhile it has been thought that luminal cells are the cell of origin for breast cancers recent studies\nhave showed that basal cells can be the cell of origin24,25.\nRegarding IL33, its role in normal mammary stem/progenitor cells, cellular plasticity and\ntumor initiation have never been studied. Previous studies have primarily examined the role of\nIL33 from the microenvironment on tamoxifen resistance and metastasis26,27. Therefore, our\nfinding that 1. in a minority population of cells, IL33 expression increases with age 2. is reduced\nin aged parous mice and 3. regulates growth and differentiation of normal mammary epithelial\ncells to our knowledge have not been described before.\n4. The authors have to include 3-month-old parous mice in their flow cytometry and RNAseq analyses\n(Figure 1-3) and also 3-month-old virgin and parous to see how quickly the differences due to pregnancy\nare observed. This is important, because the authors argue that their study is unique, since they looked\nat old parous mice compared to prior studies that analyzed mice shortly after pregnancy. But, w/o showing\ndata for 3-month-old virgin and parous this is not possible to conclude.\nWe thank the reviewer for this comment. As per the reviewer’s suggestion, we have now\nrevised Fig. 1 & 3 to include a 3-month-old parous mice (mice underwent a minimum of two\npregnancies and were analyzed at least 3 weeks post-involution).\nIn brief, we repeated all the experiments in Fig. 1 with a 3-month-old parous (3M P) mouse\nside by side with 3M NP, 18M NP and 18M NP. Our experiments characterize the mammary gland\nat these different stages through flowcytometry and organoid assays, which has never been done\nbefore. We find that parity has distinct effects on the aging of luminal and basal cells – 1.\nLuminal cells from parous mice do not undergo age related loss that is seen in nulliparous mice 2.\nParity restricts basal cell expansion that is seen with aging, but it is not a complete rescue. We also\nperformed organoid assays with the 4 groups to find that while aging increase organoid formation\ncapacity of basal cells, it reduces the organoid forming capacity of luminal cells.\nFor Fig 3., we were able to compute the fraction of Krt6a+/Il33+ cells across the 4 different\ngroups as described above and included in the revised manuscript (Fig. 3g). We find that while\nthere is a decrease in the Krt6a+/Il33+ cells immediately post-pregnancy, this percentage is\nsignificantly lower in 3M NP and 3M P mice (Fig 3g). Thus, our data suggests that accumulation\nof the Krt6a+/Il33+ cells occurs with aging in nulliparous mice. Moreover, we have confirmed this\nanalysis using immunofluorescence (Fig. 4a-c).\nFor the previous Fig. 2, we integrated the Pal et al. dataset22 with ours and performed\npseudo bulk analysis across all the groups (young nulliparous, young parous, aged nulliparous\nand aged parous). We have included this analysis below as part of the point-by-point response\n(Fig. R1).\nFirst, we confirmed successful isolation of basal and luminal cell populations from the\npseudobulk analysis by expression of Krt14 and Acta2 in basal samples and Krt8 and Krt18 in\nluminal samples (Fig. R1a). To identify transcriptomic differences between nulliparous and parous\nMECs as they age, we performed differential gene expression analyses on aged basal or luminal\ncells versus young basal or luminal cells from either parous or nulliparous mice (i.e., 18M NP vs.\n3M NP and 18M P vs. 3M P) and identified differentially expressed genes with opposing log2 fold\nchanges (Fig. R1b). Basal cells from parous mice displayed increased expression of the basal\ncell lineage marker Lipg28 and ECM-related genes (Ltbp2, Col4a1, Fig. R1c, d) with age, whereas\nbasal cells from nulliparous mice showed a decrease in expression with age. In agreement with\nprevious studies29, basal cells from nulliparous mice had an age-dependent increase in\nexpression of Arg1. However, we find that basal cells from parous mice showed an opposite trend,\nwith Arg1 being highly expressed in basal cells from young parous mice and returning to levels\nsimilar to young nulliparous basal cells in aged parous basal cells (Fig. R1c, d). Similarly, the\nexpression of Tspan8, which has been previously reported to mark a population of quiescent\nhormone responsive stem-cells30, increased with age in nulliparous basal cells but this was\nreversed in parous (Fig. R1c, d). We had independently analyzed and confirmed the age- and\nparity-dependent trend in Tspan8 expression at the protein level with flow cytometry analysis\n(Extended Data Fig. 3).\nThe number of age-specific differentially expressed genes (DEGs) with opposing log2 fold\nchanges between parous and nulliparous mice was far greater in luminal cells compared to basal\ncells (Fig. R1e), underscoring the dramatic impact of pregnancy on luminal cell transcriptomes.\nAs was observed in basal cells, Arg1 expression increased with age in nulliparous mice, but\ndecreased with age in parous mice (Fig. R1e, f). Other genes showing similar expression patterns\nwith age and parity include Kitl, Sparc, Serpina3n, and Igfbp3, which are often associated with\nthe basal cell lineage28, as well as the basal progenitor cell marker Cd3631 (Fig. R1e, f).\nOur findings support the hypothesis that MECs display a loss of lineage integrity over the\ncourse of aging and suggest that previous pregnancies attenuate this age-dependent phenotype\nwhile instilling a persistent transcriptomic memory of pregnancy. However, we feel this analysis\nis largely descriptive as pointed out by Reviewer 2 and tangential to the main point of the paper.\nWe have therefore removed it from the main manuscript.\n[FIGURE REDACTED]\nFigure R1: Parous and Nulliparous MECs display distinct transcriptomic patterns with\naging. (a) Expression of basal (Krt14, Acta2) and luminal (Krt18, Krt8) lineage markers in bulk\nRNA-seq samples collected from FACS-sorted basal and luminal cells of 18M NP (red) and 18M\nP (blue) mammary glands. 18M bulk RNA-seq samples were integrated with 3M NP (grey) and\n3M P (green) scRNA-seq datasets (GSE103275) for pseudo-bulk RNA-seq analyses. (b)\nSchematic diagram describing the strategy used to identify differentially expressed genes with\naging (18M NP vs. 3M NP and 18M P vs. 3M P) and exhibit opposing log2 fold changes (Log2FC).\n(c) Differentially expressed genes (DEGs) in 18M NP basal vs. 3M NP basal (pink) and 18M P\nbasal vs. 3M P basal (blue) with opposing Log2FC values. DEGs were filtered for p_adj. < 0.05\nand |Log2FC(NP) – Log2FC(P)| > 1.5. (d) Expression values of select genes identified in (c)\nacross comparison groups. (e) Differentially expressed genes (DEGs) in 18M NP luminal vs. 3M\nNP luminal (pink) and 18M P luminal vs. 3M P luminal (blue) with opposing Log2FC values. DEGs\nwere filtered for p_adj. < 0.05 and |Log2FC(NP) – Log2FC(P)| > 1.5. (f) Expression values of\nselect genes identified in (d) across comparison groups. Gene expression is represented as\ntranscript per million (TPM) × 103 (a) or normalized values from DESeq2 (d, f). n = 3 mice per\ngroup\n5. Similarly, in Figure 4 the authors argue that IL33 treatment of organoids from young mice phenocopies\naging. But they do not analyze organoids from old mice and also do not compare 3-month-old virgin vs.\nparous. These additional samples would be critical to include.\nWe thank the reviewer for their comment. As suggested by the reviewer we have included\nall the 4 groups using both the mammary glands from mice (Fig. 4a-c) and on organoids derived\nfrom them (Fig. 4d-g).\nConsistent with the single cell RNA-sequencing data, we find that the 18M NP mice have\nthe highest percentage of KRT6A+ cells and that this percentage is reduced in 18M P mice.\nMoreover, we find that there is a slight but significant reduction of KRT6A+ cells in the young\nparous mice (Fig. 4a, b). For the in vitro culture, we observe the increase in KRT6A+ cells in the\norganoids derived from basal cells (Fig. 4d, e) but do not see any differences in organoids derived\nfrom luminal cells (Fig. 4f, g) in 18M NP mice.\nReviewer #2 (Remarks to the Author):\nThis study by Olander et al. investigates the combined effects of aging and pregnancy on mouse\nmammary epithelial cells (MECs). The authors report that pregnancy reduces the age-induced expansion\nof basal cells and partially normalizes the basal:luminal cell ratio in aged mice. They also report the\npresence of a small population of Il33-expressing hybrid MECs, which accumulate in aged nulliparous\nmice, but reduced in aged parous mice. The authors show that pregnancy induces transcriptomic\nchanges in both basal and luminal cells, promoting a differentiated state and reversing some age-\nassociated transcriptional programs. Finally, the authors show that IL33 treatment of young basal cells\nmimics the aged nulliparous MECs, increasing the proportion of basal cells and inducing a KRT6a+ hybrid\nstate, which may confer a survival advantage for cells with oncogenic mutations. Overall, the study\nsuggests that pregnancy-induced changes in MECs contribute to long-term breast cancer protection by\nmaintaining lineage integrity and reducing the accumulation of potentially tumorigenic hybrid cells. The\nstudy is well presented and written. Generally, there is a slight over interpretation of some of the data\nwhich should be substantiated with further evidence. Below are the main points this reviewer believes\nshould be addressed.\nWe thank the reviewer for their assessment of our manuscript and constructive feedback.\nWe have addressed the concerns of the reviewer below by including new data.\n1) Title: I personally feel that the title might be slightly overreaching, as the paper does not explicitly\ndemonstrate that pregnancy reduces the formation of the hybrid progenitor. The findings just show that\nthe hybrid progenitor is identified in the old nulliparous samples but is less prevalent in the old parous\nsamples.\nWe appreciate the reviewer for raising this concern. In view of the reviewer’s concern, we\nhave revised our title to read, “Characterization of the mammary gland with pregnancy and aging\nidentifies a population of IL33+ hybrid progenitors."}]

Decision letter