Comprehensive single-cell aging atlas of healthy mammary tissues reveals shared epigenomic and transcriptomic signatures of aging and cancer

Age is the greatest risk factor for breast cancer, with two-thirds of cancers occurring in women over 50 years of age¹. Understanding the cellular and molecular changes in mammary cells during aging can reveal novel insights into age-related cancer initiation. The mammary gland undergoes extensive dynamic remodeling at the cellular and molecular level throughout development, menstruation, pregnancy, lactation and involution^{2,3,4,5,6,7,8,9,10,11}. Mammary tissues are composed of duct-forming epithelial cells embedded in a stromal compartment containing fibroblasts (Fib), adipocytes and vascular, endothelial and immune cells. Two epithelial cell types are critical for mammary functions: luminal cells, which form the inner layer of ducts and lobules and from which most breast cancers are thought to originate, and myoepithelial, cells which surround the luminal layer and limit epithelial cell dissemination¹². Although studies have reported age-related changes in epithelial cell proportions and lineage fidelity in mouse and human breast tissues^{13,14,15,16,17,18,19,20}, knowledge of the other cell types and their roles during mammary aging remains limited^18,19,20,21. Furthermore, there is limited understanding of the underlying molecular drivers of these age-related changes in mammary cells and how they might relate to age-related breast cancer.

Changes in epigenetics, a hallmark of aging, could serve as a molecular driver of age-related changes in transcription and cell identity²². Indeed, widespread changes in chromatin accessibility have been observed across multiple aged human and mouse tissues^23,24. Profiling of younger mice revealed a distinct chromatin landscape of mammary epithelial subtypes^25,26; however, how aging impacts these profiles remains to be characterized. In human tissues, aged luminal cells exhibited distinct methylation patterns compared to younger individuals, impacting genes involved in lineage fidelity and cancer susceptibility^27,28. Finally, aging of the mammary microenvironment triggered DNA methylation and gene expression changes in human luminal cells in vitro¹⁴, highlighting the need to further investigate the complex role of epigenetics.

In this study, we leveraged single-cell RNA sequencing (scRNA-seq) and single-nucleus assay for transposase-accessible chromatin with sequencing (snATAC-seq) to comprehensively study gene expression and chromatin accessibility profiles during mammary gland aging in healthy mice. Our mouse aging mammary gland atlas captured cell compositional changes in epithelial, immune and stromal populations, along with transcriptomic and epigenomic changes with age. By integrating expression and chromatin accessibility data, we provide mechanistic insights into the transcriptional programs regulating mammary aging. With age, epithelial, immune and stromal clusters exhibited altered expression of cell identity marker genes, suggesting decreased lineage fidelity and increased cell plasticity. We also identified gene expression and chromatin accessibility changes in cancer-associated, senescence-associated and inflammation-associated genes. Furthermore, using spatial transcriptomics, we localized a subset of the age-related cell types and investigated predicted cell co-localization patterns. Finally, we identified age-related signatures of mammary cells that are found in human breast tumors, suggesting that these could be linked to pre-neoplasia.

Mammary glands undergo cell compositional changes with aging

To characterize the regulatory landscapes of aging mammary tissues, we isolated mammary glands from co-housed young adult (3 month) and older (18 month) virgin female C57BL/6J mice, which correspond to 20–30-year-old and more than 55-year-old humans, respectively. Mammary tissues were dissociated using a two-step lysis protocol to preserve sensitive immune cells while simultaneously capturing epithelial and stromal cells (Methods). Gene expression and chromatin accessibility in dissociated single cells were profiled using 10x Chromium scRNA-seq (n = 6 per age) and matched snATAC-seq from half of the samples (n = 3 per age) (Fig. 1a–c, Extended Data Fig. 1a,b and Supplementary Table 1a,b). Cell clustering of scRNA-seq and snATAC-seq data revealed three epithelial (luminal alveolar (AV), luminal hormone sensing (HS), myoepithelial), five immune (naive T cells, memory T and natural killer (NK) cells, plasma cells, dendritic cells (DCs) and macrophages, B cells) and three stromal (fibroblast, pericyte, vascular) clusters annotated using marker genes^7,18 (Fig. 1b,c, Extended Data Fig. 1c–f and Supplementary Table 2a). Cell compositional analysis revealed overall consistent changes between replicates in scRNA-seq and snATAC-seq for most cell types (Fig. 1d,e, Extended Data Fig. 1g and Supplementary Table 2b,c). In epithelial cells, two luminal subtypes—AV (Csn3⁺, Trf⁺, Mfge8⁺) and HS (Esr1⁺, Cited1⁺, Prlr⁺) cells—both decreased in proportion with age in snATAC-seq and scRNA-seq, respectively, whereas cells expressing myoepithelial markers (Krt17⁺, Acta2⁺, Myl9⁺) increased with age in both datasets (Fig. 1d,e and Extended Data Fig. 1g). For immune cells, aged animals had higher proportions of myeloid (Itgax⁺ or C1qa⁺), plasma (Jchain⁺, detected only in scRNA-seq) and memory T cells (Cd3d⁺, S100a4⁺) and decreased naive T cells (Cd3d⁺, Sell⁺) (Fig. 1d,e and Extended Data Fig. 1g). Finally, fibroblasts (Fn1⁺, Col1a1⁺) increased with age in the mammary stroma (Fig. 1d,e and Extended Data Fig. 1g).

a, Experimental approach using scRNA-seq and snATAC-seq on cells isolated from freshly dissociated mammary glands from 3-month-old (3M) and 18-month-old (18M) virgin female C57BL/6J mice. b,c, UMAP visualization of epithelial (luminal AV n = 2,843; luminal HS n = 1,494; and myoepithelial n = 2,971), immune (memory T and NK cells n = 6,684; naive T cells n = 13,771; B cells n = 12,044; plasma cells n = 116; and DCs and macrophages n = 3,485) and stromal (pericytes n = 417; vascular cells n = 686; and fibroblasts n = 3,130) clusters captured by scRNA-seq identified based on characteristic marker genes (b) and by snATAC-seq upon annotation transfer from scRNA-seq (c). d,e, Average proportions of epithelial, immune and stromal cells in 3M and 18M mice captured by scRNA-seq (n = 6, mean ± s.e.m.) (d) and by snATAC-seq (n = 3, mean ± s.e.m.) (e) (paired two-sided t-test; *P ≤ 0.05, **P ≤ 0.01, ***P ≤ 0.001,****P ≤ 0.0001; P values are listed in Supplementary Table 2b,c; replicates are shown in Extended Data Fig. 1g). With age, luminal HS cells decreased significantly, and plasma cells were detected and increased significantly in scRNA-seq (d). Luminal AV cells and pericytes decreased significantly in snATAC-seq (e). Myoepithelial cells, DCs and macrophages, memory T and NK cells and fibroblasts increased significantly in both analyses, and naive T cells decreased significantly in both analyses (d,e). See also Extended Data Fig. 1 and Supplementary Tables 1 and 2.

Full size image

Gene expression and chromatin changes in aged epithelial cells

In addition to age-related cell compositional changes, we detected changes in cell-intrinsic gene expression patterns across epithelial subtypes, with luminal cells displaying a bias toward upregulated genes, whereas myoepithelial cells exhibited more downregulated genes (Fig. 2a and Supplementary Table 3a). Overall, 80%, 71% and 82% of age-related differentially expressed (DE) genes were cell type specific in luminal AV, luminal HS and myoepithelial cells, respectively (Extended Data Fig. 2a). However, 29 genes were shared across all epithelial cell types, and several were involved in age-related pathways, such as the downregulation of ribosomal proteins and upregulation of interferon-gamma-related genes (Extended Data Fig. 2a and Supplementary Table 3b). Interestingly, aged luminal AV and HS cells exhibited increased expression of genes with tumor suppressor activity (Fig. 2b and Supplementary Table 3a). In addition, multiple cancer hallmark gene sets were differentially expressed with age across cell clusters, including epithelial-to-mesenchymal transition (EMT), mTORC1 signaling, estrogen responses, p53 signaling, MYC signaling, inflammatory response, hypoxia and TNF signaling (Extended Data Fig. 2b,c and Supplementary Table 3c). Although the same pathways were altered with age in multiple epithelial clusters, the specific genes were often distinct between clusters (Fig. 2b, Extended Data Fig. 2b,c and Supplementary Table 3c).

a,c, UMAP visualization of epithelial cell clusters captured by scRNA-seq (a) and snATAC-seq (c). The number of significant DE genes detected by single-cell and pseudobulk analysis with age (a) and DA peaks with age (c) is shown per cell cluster. b, Number of significant DE genes with tumor suppressor activity in luminal AV, luminal HS and myoepithelial cell clusters from 18M versus 3M mice (detected by single-cell and pseudobulk analysis). d, Differential TF activity score with age. Asterisks indicate significant differential motifs (Wilcoxon test; Bonferroni P_adj < 0.05). e, Top DE genes from 18M versus 3M mice across replicates from pseudobulk scRNA-seq data. f–k, Examples of DE genes (red in e) with DA peaks in luminal AV (f,g), luminal HS (h,i) or myoepithelial (j,k) clusters in 3M versus 18 Mmice. Left: log-normalized values are shown for individual cells (Welch two-sample t-test; **P_Fndc4 = 0.0033, ****P_Hp = 2.16 × 10⁻⁷, **P_Pygl = 0.0011, ****P_Epb41l3 = 9.46 × 10⁻⁵, **P_Fli1 = 0.0013, ****P_Prrx1 = 4.73 × 10⁻¹¹; cells with zero counts were excluded), along with a pie chart depicting the percentage of expressing cells versus non-expressing cells. Right: pseudobulk snATAC-seq tracks in 3M versus 18M mice are shown, along with gene structures and significant DA peaks with corresponding log₂fold change (FC) values. Predicted TF binding motifs from JASPAR are indicated within the DA peaks. See also Extended Data Figs. 2 and 3 and Supplementary Table 3.

Full size image

We also detected age-related changes in chromatin accessibility (Fig. 2c and Supplementary Table 3d). Similar to DE profiles, most changes were cell type specific, with only 257 genes with differentially accessible (DA) peaks shared across clusters (Extended Data Fig. 2d and Supplementary Table 3e). Roughly 25% of DA peaks were detected in promoter regions, and remaining peaks mapped within exons, introns, untranslated regions (UTRs) and intergenic regions (Extended Data Fig. 2e). DA peaks were associated with cancer hallmark gene sets, including estrogen responses, hypoxia, inflammatory response and TNF signaling (Extended Data Fig. 2f and Supplementary Table 3f). To define putative transcriptional regulators, we conducted chromVar analyses to infer transcription factor (TF) activity based on chromatin accessibility associated with TF binding sites in young and older samples (Fig. 2d). Luminal AV cells displayed increased activity of JUN and NFκB family members (RELA/B) with age, which are involved in regulating pro-inflammatory responses. Activity of TP53 and family members increased in luminal HS cells with age, whereas AP-1 factor (FOS and JUN) activity decreased (Fig. 2d and Supplementary Table 3g). In addition to changes in TF activity, changes in TF expression included downregulation of JUN family members in luminal HS and AV cells (Extended Data Fig. 2g and Supplementary Table 3a).

We identified the top DE genes with age (Fig. 2e and Supplementary Table 3a), noting several mammary-gland-associated and breast-cancer-associated genes with concordant peak accessibility changes (Fig. 2f–k and Extended Data Fig. 3a–c). For example, Fndc4, which encodes an extracellular matrix (ECM) protein with anti-inflammatory activity, showed decreased expression in aged luminal AV cells and closing of a DA peak in its promoter region (Fig. 2f). Conversely, both the expression and promoter accessibility of Hp (haptoglobin), a gene implicated in metabolic reprogramming and breast cancer, increased in luminal AV cells with age (Fig. 2g). In luminal HS cells, Pygl and Epb41l3 were upregulated with age, accompanied by the opening of DA peaks in their respective promoters (Fig. 2h,i). Pygl is linked to metabolism in normal and breast cancer cells²⁹, whereas Epb41l3 is a tumor suppressor gene (TSG) often demethylated in breast tumors³⁰. Aged myoepithelial cells upregulated genes encoding TFs Fli1 and Prxx1 and displayed opening of DA peaks in their promoter regions (Fig. 2j,k). The proto-oncogene Fli1 plays a role in breast cancer cell proliferation³¹, whereas Prxx1 activates EMT and promotes drug resistance³². Furthermore, DA peaks in Fndc4, Hp, Pygl, Epb41l3, Fli1 and Prxx1 contained putative binding motifs for TFs with differential activity in old versus young (Fig. 2d), including BACH2 and FOX in luminal AV cells, FOS and TFAP2 in luminal HS cells and IRF and LBX1 in myoepithelial cells (Fig. 2f–k). Additional genes that exhibited both age-related expression and chromatin accessibility changes included Pdk4, encoding a kinase implicated in glucose metabolism; Rspo1, encoding a Wnt signaling agonist important in stem cell regulation; Alox12e, encoding a protein involved in lipid metabolism; and others (Extended Data Fig. 3a–c). Finally, several DE genes in mouse luminal cells were also detected in aged human luminal cells in vitro³³ (Extended Data Fig. 3d and Supplementary Table 3h), including the downregulation of Elf5, encoding a TF and marker of aged human luminal cells³⁴.

Changes in cancer-associated genes and cell identity markers

To further deconvolute age-related changes in epithelial cells, we performed unsupervised subclustering of the three epithelial cell populations, uncovering three luminal HS (Epi-C1–C3), four luminal AV (Epi-C4–C7) and four myoepithelial (Epi-C8–C11) subclusters (Fig. 3a–c, Extended Data Fig. 4a–c and Supplementary Table 3i). Although these subclusters shared similar expression patterns, we detected subcluster-specific expression patterns, potentially reflecting changes in cell states (Fig. 3d–h and Supplementary Table 3j).

a, UMAP visualization of epithelial subclusters Epi-C1–C3 (luminal HS), Epi-C4–C7 (luminal AV) and Epi-C8–C11 (myoepithelial) captured by scRNA-seq. b,c, Differences in cell numbers with age are shown by UMAP visualization of subclusters Epi-C1–C11 captured by scRNA-seq shown per age (b), along with cell number ratio between 18M and 3M mice (c) (n = 6; two-sided t-test; ***P_Epi-C8 = 0.00026, **P_Epi-C6 = 0.0084, *P_Epi-C2 = 0.048, *P_Epi-C2 = 0.023; NS, not significant; boxplots: center line, median; box limits, upper and lower quartiles; whiskers, 1.5× interquartile range; points, replicates). d–f, Expression of canonical marker genes in scRNA-seq luminal subclusters (d). For each luminal subcluster, the proportions of cells from 3M and 18M mice are shown below. Feature plots for selected marker genes are shown for luminal HS Epi-C1–C3 (e) and luminal AV Epi-C4–C7 (f) subclusters. g,h, Expression of canonical marker genes in scRNA-seq myoepithelial subclusters (g). For each myoepithelial subcluster, the proportions of cells from 3M and 18M mice are shown on the right. Feature plots for selected marker genes are shown for myoepithelial Epi-C8–C11 subclusters (h). Epi-C11 cells originated mostly from one replicate (Supplementary Table 3i), and, therefore, their pro-inflammatory population should be further investigated. See also Extended Data Fig. 4 and Supplementary Table 3. Lum., luminal; Sub., subcluster.

Full size image

In luminal HS subclusters, Epi-C1 expressed the TSG Rcan1 and decreased with age (Fig. 3d,e). Epi-C2 expressed ductal cell marker genes Fxyd2 and Tph1 (encoding proteins implicated in mammary expansion and milk production) and increased with age (Fig. 3d,e). Epi-C2 also expressed Fam3c, encoding a cytokine expressed in proliferative tissues that activates MAPK signaling and is associated with metastasis and poor survival³⁵ (Fig. 3d). Cells from Epi-C3 expressed marker genes of luminal HS, luminal AV and progenitor cells (reminiscent of a luminal HS-AV cluster in mice^7,18) and showed a decreasing trend with age (Fig. 3c–e) but did not reach significance, possibly due to low cell numbers (n = 97).

In luminal AV subclusters, both Epi-C4 and Epi-C5 highly expressed Hey1, encoding a downstream effector of Notch signaling, and the TSG Txnip (Fig. 3d). Epi-C4 expressed high levels of luminal AV marker genes (Csn3, Mfge8, Cst3, Igfbp5) and trended toward depletion in older animals (Fig. 3d,f), suggesting a partial loss of cell identity with age. Epi-C6, enriched with age, expressed Alox15 and Alox12e, encoding regulators of metabolism and inflammation, and Palmd, encoding a p53 target (Fig. 3d,f). Epi-C5, Epi-C6 and Epi-C7 all expressed Rspo1 (Fig. 3d,f). Finally, Epi-C7 was a small subcluster (n = 108) and expressed cycling markers (Mki67, Cdk1, Stmn1), similar to cells detected in adult mouse mammary gland⁷ and human breast²⁰ (Fig. 3d,f and Supplementary Table 3i).

In myoepithelial subclusters, Epi-C8 cells were more abundant with age and downregulated Krt17 and Krt5 (Fig. 3g,h), suggesting a loss of cell identity markers with age. Furthermore, cells from Epi-C9, which did not significantly change in number but decreased in proportion with age, expressed subcluster-specific (Tagln, Postn, Actg2) and myoepithelial marker (Krt17, Krt5) genes (Fig. 3c,g,h). Epi-C11 expressed an inflammatory and interferon-gamma signature (Ccl2, Cxcl10, Irgm1, Stat1) (Fig. 3g,h). Finally, differential expression analyses revealed aging-associated or cancer-associated genes changing with age across epithelial subclusters (Extended Data Fig. 4c and Supplementary Table 3k).

Fibroblasts express ECM and senescence markers with age

In the stroma, fibroblasts, pericytes and vascular cells showed both age-related gene expression and chromatin accessibility changes (Fig. 4a,b and Supplementary Table 4a,c). We further focused on fibroblasts due to their robust cell compositional, expression and chromatin accessibility changes with age and their critical role in the mammary gland, including ECM deposition and remodeling, paracrine signaling and interactions with epithelial cells³⁶. DE changes revealed increases in TNF signaling (Fos and Jun) and senescence-associated proteins (Cdkn1a and Cdkn2a) and decreases in EMT-related genes (including collagens) and translation-related genes (Extended Data Fig. 5a and Supplementary Table 4b). At the chromatin level, approximately 20% of DA peaks were found in promoter regions (Extended Data Fig. 5b). chromVar analyses suggested decreased activity of fibrosis and EMT-related TFs Twist1, Tcf12 and Nfya, along with increased activity of Hand2 with age (Extended Data Fig. 5c and Supplementary Table 4e). Genes associated with opening DA peaks were enriched in TNF signaling, pro-inflammatory and breast-cancer-related gene sets (Extended Data Fig. 5d and Supplementary Table 4d), similar to pathways enriched in age-related fibroblast DE genes. We detected expression and chromatin accessibility changes in genes related to fibroblast function and adiposity, such as upregulation of Fapb4, which encodes a protein that is critical for fatty acid transport and has been shown to promote breast tumorigenesis and metastasis³⁷ (Fig. 4c) and other genes (Extended Data Fig. 5e).

a,b, UMAP visualization of stromal cell clusters captured by scRNA-seq (a) and snATAC-seq (b). The number of significant DE genes detected by single-cell and pseudobulk analysis with age (a) and DA peaks with age (b) is shown per cell cluster. c, Example of DE gene with DA peak in fibroblast clusters in 3M versus 18M mice. Left: log-normalized values are shown for individual cells (Welch two-sample t-test; ***P = 1.32 × 10⁻⁴; cells with zero counts were excluded), along with a pie chart depicting the percentage of expressing cells versus non-expressing cells. Right: pseudobulk snATAC-seq tracks in 3M versus 8 M mice are shown, along with gene structures and significant DA peaks with corresponding log₂FC values. d,e, UMAP visualization of fibroblast subclusters Fib-C0–C5 captured by scRNA-seq (d), along with expression of canonical marker genes (e). The proportions of cells from 3M and 18M mice are shown on the right. f, Differences in cell number ratios with age (18M/3M) per fibroblast subclusters captured by scRNA-seq (n = 6; two-sided t-test; **P_Fib-C2 = 0.0028, **P_Fib-C0 = 0.0081, *P_Fib-C4 = 0.043; NS, not significant; boxplots: center line, median; box limits, upper and lower quartiles; whiskers, 1.5× interquartile range; points, replicates). g,h, Expression of marker genes in fibroblast subclusters captured by scRNA-seq. i, log-normalized counts of Cdkn1a in fibroblasts with age (two-sided t-test; ****P = 1.2 × 10⁻¹³; cells with zero counts were excluded). j, log-normalized counts of Cdkn1a in fibroblast subclusters with age with feature axis on log scale (two-sided t-test; ****P_Fib-C0 = 9.1 × 10⁻⁸, ****P_Fib-C2 = 3.8 × 10⁻⁹, *P_Fib-C3 = 0.012; NS, not significant; cells with zero counts were excluded). k, Senescence gene signatures in 3M versus 18M mice across fibroblast subclusters (one-sided t-test; ****P_Fib-C0 = 2.7 × 10⁻¹², ****P_Fib-C2 = 5.1 × 10⁻²⁰, ****P_Fib-C3 = 3.5 × 10⁻¹⁴, ****P_Fib-C4 = 1.4 × 10⁻²⁰, ****P_Fib-C5 = 3.5 × 10⁻¹⁴). l, Top DE genes with age across replicates from pseudobulk scRNA-seq data for indicated fibroblast subclusters (with >100 cells per age). See also Extended Data Fig. 5 and Supplementary Table 4.

Full size image

To dissect the expression patterns of senescence-related and cancer-related genes in the stroma, we resolved 11 subclusters (Fib-C0–C5 and Str-C6–C10) (Extended Data Fig. 5f–h and Supplementary Table 4f,g). Among these, Fib-C0–C5 expressed fibroblast marker genes (Col1a1⁺, Pdgfra⁺) and could be further subclassified into universal Col15a1⁺ fibroblasts that secrete basement membrane proteins (Fib-C0–C3) and Pi16⁺ fibroblasts that may develop into specialized fibroblasts³⁸ (Fib-C4) (Fig. 4d–g and Extended Data Fig. 5i). The expression patterns of these subclusters aligned with fibroblasts previously identified in younger animals³⁹: ECM-remodeling (Fib-C0 and Fib-C1), high-adipogenic-capacity (Fib-C2), adipo-regulatory (Fib-C3) or Dpp4⁺ (Fib-C4) fibroblasts. Fib-C2, Fib-C0 and Fib-C4 showed the greatest increases in number with age, with a striking approximately six-fold increase in Fib-C2 (Fig. 4f). Interestingly, Fib-C2 expressed a distinct repertoire of ECM proteins, including reduced expression of classical marker gene fibronectin (Fn1) and increased expression of collagens Col6a3, Col5a3, Col4a1 and Col4a2 (Fig. 4e,g). This suggested a change in cell identity in Fib-C2 with age compared to other subclusters, especially Fib-C4, which highly expressed Fn1 (Fig. 4e,g). Fib-C2 also expressed Fap, a marker of cancer-associated fibroblasts (CAFs)⁴⁰, and genes related to lipid metabolism (Lpl, Fabp4, Pparg), suggesting an adipogenesis commitment, but lacked adipocyte markers (Adipoq, Plin1) (Fig. 4e,h). Furthermore, aged mammary glands have more stromal cells expressing senescence markers Cdkn2a (p21) and Cdkn1a (p16) and a higher senescence gene signature compared to younger tissues (Fig. 4i–k, Extended Data Fig. 5j and Supplementary Table 4a,i). Several subclusters showed differential scoring using the myofibroblastic (myCAF), inflammatory (iCAF) or antigen-presenting (apCAF) CAF signatures, with Fib-C0, Fib-C2, Fib-C3 and Fib-C4 scoring higher for apCAF with age (Extended Data Fig. 5k).

Finally, differential expression analysis of fibroblast subclusters identified age-related differences in cancer-related or age-related genes in Fib-C0, Fib-C2, Fib-C3 and Fib-C4 (Fig. 4l). For example, Cdkn1a expression increased in older cells in Fib-C0 and Fib-C2, whereas expression of Ptn, encoding a secretory growth factor associated with angiogenesis and breast cancer progression⁴¹, decreased (Fig. 4l and Supplementary Table 4e). Fib-C3 decreased expression of Meg3, encoding a long non-coding RNA and tumor suppressor in breast cancer⁴² (Fig. 4l and Supplementary Table 4e). Finally, aged Fib-C4 cells downregulated five out of 10 genes (Gpx3, Spon1, Plac8, Ctsh, Col3a1) implicated in estrogen response of Dpp4⁺ fibroblasts³⁹, suggesting a decline in hormone responsiveness with age (Supplementary Table 4e).

Cytotoxic T cells expand and exhibit exhaustion markers with age

To further dissect changes in T and NK clusters with age, we resolved 10 scRNA-seq subclusters corresponding to distinct populations of naive (Ccr7⁺) and memory (S100a4⁺) CD4⁺ and CD8⁺ T, γδ and MAIT (Trdc⁺, Zbtb16⁺) as well NK (Ncr1⁺) cells (Fig. 5a,b, Extended Data Fig. 6a,b and Supplementary Table 5a–e). With age, naive CD4⁺ and CD8⁺ T cells declined in numbers in mammary tissues, whereas Gzmk⁺ T, γδ T and MAIT, memory CD4 and NK cells expanded with age (Fig. 5a). Gzmk⁺ T cells (encompassing CD8⁺ and CD4⁺) expanded approximately 30-fold with age and expressed Tox, Eomes, Pdcd1 and Lag3, in line with an exhaustion phenotype (Fig. 5a, Extended Data Fig. 6a,c,d and Supplementary Table 5b). Memory CD4⁺ T cells, which included regulatory T cells (Tregs) (Foxp3⁺, Il2ra⁺, Ctla4⁺), expanded approximately 1.8-fold with age (Fig. 5a, Extended Data Fig. 6e and Supplementary Table 5f). Interestingly, a fraction of Tregs expressed Itgae (Cd103), which mediates cell migration and lymphocyte homing through interaction with epithelial cells, suggesting a tissue-resident phenotype (Extended Data Fig. 6a,f). Finally, γδ T and MAIT cells, which also expressed Itgae, expanded by approximately three-fold with age (Fig. 5a,b and Extended Data Fig. 6b), similar to studies in aged lung and liver⁴³. γδ T cells play an important role in cancer cell surveillance and have been linked to favorable prognosis in solid tumors⁴⁴.

a, UMAP visualization of based T and NK cell subclusters (left) along with differences in cell number ratios with age (18M/3M, right) captured by scRNA-seq (n = 6; two-sided t-test; ***P_Gzmk+ = 1.29 × 10⁻⁴, ***P_{γδ T & MAIT} = 9.00 × 10⁻⁴, **P_{Memory CD4} = 2.89 × 10⁻³, *P_NK = 2.90 × 10⁻², **P_{Naive CD8-2} = 2.69 × 10⁻³, ***P_{Naive CD4-1} = 2.37 × 10⁻⁴, *P_{Naive CD8-1} = 1.03 × 10⁻², **P_{Naive CD4-2} = 4.10 × 10⁻³; NS, not significant; boxplots: center line, median; box limits, upper and lower quartiles; whiskers, 1.5× interquartile range; points, replicates). b, Expression of marker genes in scRNA-seq subclusters. c, UMAP visualization of based T and NK cell subclusters (left) along with differences in cell number ratios with age (18M/3M, right) captured by snATAC-seq (n = 3; two-sided t-test; **P_Gzmk+ = 9.99 × 10⁻³, **P_{Naive CD8-2} = 6.69 × 10⁻³, *P_{Naive CD4-1} = 1.66 × 10⁻², **P_{Naive CD8-1} = 1.92 × 10⁻³, **P_{Naive CD4-2} = 3.95 × 10⁻³; NS, not significant; boxplots: center line, median; box limits, upper and lower quartiles; whiskers, 1.5× interquartile range; points, replicates). d, Examples of marker genes that display chromatin accessibility signatures shown as pseudobulk snATAC-seq tracks per cell subcluster. The Gzmk promoter was the most accessible in Gzmk⁺ T cells, whereas the Gzmm promoter was the most accessible in Gzmm⁺ CD8⁺ T cells. Pdcd1 and Ctla4 promoters were highly accessible in Gzmk⁺ T and memory CD4⁺ subsets. The loci around Ccl5 and T-Box TF Eomes were highly accessible in Gzmk⁺ T cells and Gzmm⁺ CD8⁺ T cells. e, DE genes with age across replicates from pseudobulk scRNA-seq data related to cell function of memory CD4, CD8 Gzmk⁺ and CD8 Gzmm⁺ immune subclusters. f–h, Examples of DE genes with DA peaks in memory CD4, CD8 Gzmk⁺ and CD8 Gzmm⁺ immune clusters (red in e) in 3M versus 18M mice. Left: log-normalized values are shown for individual cells (Welch two-sample t-test; memory CD4: **P_S100a4 = 0.001, ****P_S100a6 = 2.52 × 10⁻²¹, ****P_Pdcd1 = 4.42 × 10⁻⁵; Gzmk⁺: **P_S100a6 = 0.0042, **P_Pdcd1 = 0.0063; CD8 Gzmm⁺: ****P_S100a6 = 1.24 × 10⁻¹²⁰; NS, not significant; cells with zero counts were excluded), along with a pie chart depicting the percentage of expressing cells versus non-expressing cells. Right: pseudobulk snATAC-seq tracks in 3M versus 18M mice are shown, along with gene structures and DA peaks. i, Examples of DE genes in γδ T and MAIT cell subclusters in 3M versus 18M mice. Left: log-normalized values are shown for individual cells (Welch two-sample t-test; **P_Ccl5 = 0.0013, ****P_Tnf = 5.26 × 10⁻⁷, ***P_Itgae = 0.0002, *P_Jag1 = 0.0435; cells with zero counts were excluded), along with a pie chart depicting the percentage of expressing cells versus non-expressing cells. See also Extended Data Fig. 7 and Supplementary Table 5.

Full size image

Epigenomic analyses of T and NK cells recapitulated compositional changes revealed by scRNA-seq (Fig. 5c). Chromatin around the promoters of cell-type-specific marker genes was more accessible in the relevant cell types (Fig. 5d). chromVar analyses revealed that cytotoxic cells, including age-related Gzmk⁺ T cells, had increased activity for TFs EOMES and TBX1-6 (Extended Data Fig. 6g and Supplementary Table 5h). Eomes was also expressed in aged Gzmk⁺ T cells, in line with its role as a transcriptional regulator of these cells in other tissues⁴³ (Extended Data Fig. 6d). In memory CD4⁺ cells, JUN and FOS family members had increased activity (Extended Data Fig. 6g), consistent with reports of AP-1 binding motifs in open chromatin regions of exhausted T cells⁴⁵.

CD8⁺ Gzmm⁺, Gzmk⁺ and memory CD4⁺ T cells exhibited the most age-related cell-intrinsic transcriptional changes (Extended Data Fig. 6h and Supplementary Table 5f). Gene signatures associated with cytotoxicity were upregulated in Gzmk⁺ and Gzmm⁺ T cells with age (Extended Data Fig. 6i and Supplementary Table 5j), suggesting that, in addition to the increase in abundance of cytotoxic cells, their intrinsic cytotoxicity also increased with age. In addition, genes encoding pro-inflammatory alarmins S100a4 and S100a6 were upregulated with age in all three subclusters and displayed open chromatin regions around these genes (Fig. 5e–g), suggesting an age-related shift from a naive to a more effector cell state. In all three subclusters, Ccl5, encoding a chemokine associated with breast cancer metastasis⁴⁶, was upregulated with age, suggesting that these cells switch to a more inflammatory and migratory state with age (Extended Data Fig. 6j). We also observed upregulation of checkpoint inhibitors Pdcd1, Ctla4 and Lag3 with age in Gzmk⁺ and memory CD4⁺ populations (Fig. 5h and Extended Data Fig. 6j). Downregulation of the protein synthesis machinery, including ribosomal genes, was observed across immune populations with age (Supplementary Table 5k). Finally, a senescence-associated gene signature increased in CD8⁺ Isg15⁺, Gzmk⁺ T and memory CD4⁺ T cells with age (Extended Data Fig. 6k).

Finally, in addition to cell number changes, we also detected 148 age-related DE genes in γδ T and MAIT cells, including the upregulation genes encoding (1) chemokine Ccl5 and chemokine receptor Cxcr3; (2) Itgae; (3) Tnf, a cytokine key for immune surveillance; and (4) Jag1, a Notch ligand (Fig. 5e,i and Supplementary Table 5f). Furthermore, chromVar analyses revealed that γδT and MAIT cells had significant TF activity for tumor suppressor TP53 and RAR-related orphan receptors (Extended Data Fig. 6g).

Myeloid cell subsets expand with age

Further clustering of myeloid cells revealed seven major subclusters (Mye-C1–C7) captured by scRNA-seq and snATAC-seq (Fig. 6a,c). These subclusters expressed marker genes for monocytes (S100a8, S100a9), M1-like macrophages (Cd86, Cd38), M2-like macrophages (Cd163, Mrc1) and conventional DCs, including markers for conventional type 1 DCs (cDC1) (Clec9a, Xcr1) and mature DCs enriched in immunoregulatory molecules (mregDCs) (Ccr7, Fscn1)⁴⁷ (Fig. 6b,d, Extended Data Fig. 7a and Supplementary Table 6a–e).

a, UMAP visualization of myeloid subclusters (left) along with differences in cell number ratios in 18M versus 3M mice (right) captured by scRNA-seq (n = 6; two-sided t-test; *P_Mye-C3 = 0.0244, **P_Mye-C5 = 0.0018, ****P_Mye-C1 = 3.33 × 10⁻⁵, *P_Mye-C6 = 0.0315; NS not significant; boxplots: center line, median; box limits, upper and lower quartiles; whiskers, 1.5× interquartile range; points, replicates). b, Expression of canonical marker genes in scRNA-seq DC and macrophage subclusters. c, UMAP visualization of myeloid subclusters (left), along with differences in cell number ratios in 18M versus 3M mice (right) captured by snATAC-seq (n = 3; two-sided t-test; **P_Mye-C3 = 0.0021, **P_Mye-C5 = 0.0010; NS, not significant; boxplots: center line, median; box limits, upper and lower quartiles; whiskers, 1.5× interquartile range; points, replicates). d, Examples of cell cluster marker genes that display chromatin accessibility signatures shown as pseudobulk snATAC-seq tracks per cell cluster. e, Expression of selected markers genes in scRNA-seq DC and macrophage subclusters. f, Averaged gene expression of the top 10 DE genes in every subcluster versus every other scRNA-seq subcluster. g, Differential TF activity score per subcluster. Significant differential motifs compared to every other cluster are colored; non-significant are in gray (Wilcoxon test; Bonferroni-corrected P_adj < 0.05; P_adj values are indicated in Supplementary Table 6e). See also Extended Data Fig. 7 and Supplementary Table 6.

Full size image

Mye-C3, Mye-C5, Mye-C1 and Mye-C6 increased in proportion with age (Fig. 6a,c). Mye-C3 expressed M2-like macrophage markers and expanded by approximately four-fold with age, similar to human breast aging^19,48, and expressed chemokines Ccl7 and Ccl8, suggesting an ability to recruit other immune cells (Fig. 6a–f and Extended Data Fig. 7a). M2-like macrophages mediate tissue repair, resolve inflammation and share similarities with tumor-associated macrophages⁴⁹. On the other hand, cells expressing M1-like macrophage markers (Mye-C2) did not significantly change with age (Fig. 6a,c,e,f). Cluster Mye-C5, which expanded by approximately 2.5-fold with age, expressed interferon-stimulated genes (Ifitm3, Ifi27l2a) and Cebpb, encoding a TF that regulates immune response genes, and showed increased variation in CEBP binding sites compared to other clusters (Fig. 6a,c,f,g, Extended Data Fig. 7b and Supplementary Table 6c,e). Although not changing with age, Mye-C4 co-expressed immuno-regulatory genes (Cd274), maturation genes (Ccr7, Cd40) and Fscn1, involved in cell migration and cellular interactions (Fig. 6b,e and Extended Data Fig. 7a). These marker genes are consistent with a recently described mregDC population in lung in mice and human⁴⁷. Overall, age-related transcriptional changes in myeloid cells were restricted to a few molecules (Supplementary Table 6f), suggesting that there are more cell compositional changes with age than cell-intrinsic changes.

Spatial profiling of cell interactions in aged mammary gland

To investigate the co-occurrence of epithelial, immune and fibroblast cells in situ in aged tissue, we generated spatial transcriptomic (ST) data from two mammary glands from 18M (ST1 and ST2). Using marker gene expression, we identified six ST clusters: one epithelial-enriched, one immune-enriched and four stromal-enriched ST clusters (Fig. 7a,b and Extended Data Fig. 8a,b). ST captured similar clusters to our single-cell data, with the addition of adipocytes, which were largely excluded in our 10x approach due to the lysis protocol, adipocyte buoyancy and size constraints of cell sorting. The ST annotations largely matched histopathology annotations from hematoxylin and eosin (H&E) staining: immune-enriched spots corresponded to the lymph node, whereas epithelial-enriched spots matched with regions containing mammary ducts and lobules, and stromal clusters matched adipose and connective tissue-dense regions.

a,b, Identification and annotation of ST spots within aged mammary gland tissues (n = 2, 18M ST1 and 18M ST2). Each ST spot (shown as 1.7× size) (a) is annotated as epithelial enriched, immune enriched or stromal enriched, using the expression of selected marker genes (b) and colored accordingly. c, Expression of indicated immune marker genes in ST spots in each tissue after excluding the lymph node. Dot plot shows the fraction of ST spots expressing the marker gene per ST cluster (epithelial enriched, immune enriched or stromal enriched) colored by mean expression. Zoomed-in images show examples of ST spots (1× size, colored by ST cluster annotation described in (a) that expressed indicated immune marker genes and are located near epithelial ducts as identified based on H&E staining. d,g,j, Ligand–receptor interactions inferred from scRNA-seq using CellPhoneDB between epithelial or fibroblast versus Gzmk⁺ clusters (d), Gzmk⁺ versus epithelial or fibroblast clusters (g) and γδ T cells versus epithelial or fibroblast clusters (j). Dot size represents P value scaled to −log₁₀ value (calculated using CellPhoneDB permutation test; P values available in Supplementary Table 7), and color represents the mean of the average expression of the first interacting molecule in the first cluster and second interacting molecule in the second cluster. e,h,k, Expression of indicated ligand–receptor pairs in ST spots (shown as 3× size) in each tissue after excluding the lymph node, colored by ST cluster annotations described in (a). Dot plot shows the fraction of ST spots expressing the gene pair per ST cluster (epithelial enriched, immune enriched or stromal enriched), colored by mean expression normalized in each ligand–receptor pair. f,i,l, Co-localization of ST spots (shown as 3× size) expressing indicated immune marker gene (yellow), ligand–receptor pair (blue) or both (green) in each tissue (n = 2, 18M ST1 and 18M ST2) after excluding the lymph node. Zoomed-in images show examples of co-occurring (green) or directly adjacent ST spots (1× size) located near epithelial ducts as identified based on H&E staining. See also Extended Data Figs. 8 and 9 and Supplementary Table 7.

Full size image

We leveraged the ST data to determine whether the age-related T cells identified by scRNA-seq were found in situ near mammary ducts and lobules. As predicted, most spots expressing Gzmk⁺ T markers (Cd3⁺ or Cd8⁺ and Gzmk⁺) were classified as immune spots and co-localized with the lymph node (Extended Data Fig. 8c); however, Gzmk⁺ T spots located outside of the lymph node were epithelial enriched (Fig. 7c). Similarly, spots expressing Pdcd1⁺ T markers (Cd3⁺ or Cd8⁺ and Gzmk⁺) and γδ T markers (Cd3⁺ and Il17a⁺ or Trdc⁺) were epithelial enriched after excluding the lymph node (Fig. 7c and Extended Data Fig. 8c). We also generated ST data on two mammary glands from 3-month-old mice and found that they had fewer total ST spots expressing Gzmk⁺, Pdcd1⁺ and γδ T markers, and a smaller proportion of these spots were epithelial enriched compared to older animals (Extended Data Fig. 9). Thus, the ST data supported the presence of immune cells within epithelial-enriched ST spots in situ in aged mice, suggesting that these epithelial–immune interactions might occur more frequently in aging tissues.

To further investigate cellular interactions, we inferred putative ligand–receptor interactions from our scRNA-seq data using CellPhoneDB (Supplementary Table 7). Predicted interactions between a4b1 integrins on Gzmk⁺ cells and molecules (including Plaur and Vcam-1) on epithelial cells and fibroblasts increased with age (Fig. 7d and Supplementary Table 7). Age-related increases in integrin-mediated interactions of Gzmk⁺ cells were previously predicted in other aging mouse tissues⁴³. ST data identified multiple epithelial ST spots that co-expressed both receptor–ligand pairs (that is, integrins and Plaur or Vcam-1) (Fig. 7e and Extended Data Fig. 8d). A subset of these ST spots co-localized with Gzmk⁺ T spots in aged mammary glands (Fig. 7f), further supporting that these predicted interactions occur in our cells of interest in situ. Furthermore, we detected interactions between the immune-inhibitory receptor Pdcd1 and Fam3c (Fig. 7g and Extended Data Fig. 8e). Fam3c was upregulated in luminal HS cells and served as a marker gene for Epi-C2 that expanded with age (Fig. 3d). Multiple ST spots co-expressed Fam3c and Pdcd1, with many being epithelial enriched, including a subset that co-localized with Gzmk⁺ T spots (Fig. 7h,i). Additionally, CellPhoneDB revealed increased interactions with age between chemokine ligands expressed by immune cells and receptors expressed by epithelial cells, including Ccl5-Ackr4, Ccl3-Ide and Ccl4-Slc7a1 (Fig. 7g), and multiple ST spots co-express these chemokine ligand–receptor pairs (Fig. 7h,i and Extended Data Fig. 8f).

Finally, CellPhoneDB suggested increased cellular interactions between γδ T cells and luminal HS, myoepithelial and fibroblasts cells through Jag1 and Notch proteins (Fig. 7j). Multiple Notch family members were upregulated with age, including Notch1 in myoepithelial cells and Notch3 in fibroblasts (Supplementary Tables 3 and 4). The Jag–Notch interaction was supported by the co-occurrence in situ of ST spots co-expressing Notch3 and Jag1, with most spots being epithelial enriched (Fig. 7k,l).

Conserved signatures of murine aging and human breast cancer

Among the luminal age-related DE genes from healthy nulliparous cancer-free mice, several genes were previously associated with cancer, including (1) upregulation of epidermal growth factor receptor (Egfr) and TSG Trp53 and downregulation the genes encoding apoptotic factor Bcl-2 and epithelial cellular adhesion molecule Epcam in luminal AV cells; (2) upregulation of TSG Socs2 and stemness gene Mex3a and downregulation of TSG Rbms2 and oncogene Myc in luminal HS cells; and (3) upregulation of the genes encoding Notch1 and Fgfr3 receptors and downregulation of TSG Slit2 and cyclin Ccnd2 in myoepithelial cells (Supplementary Table 3a).

To further reveal which age-related changes might be linked to breast cancer, we compared DE genes found in aged epithelial cells from healthy nulliparous cancer-free mice with those found in human breast tumors versus normal breast tissues from The Cancer Genome Atlas (TCGA). We focused on luminal A and B tumors, which are thought to originate from luminal cells and affect older patients⁵⁰ (Extended Data Fig. 10a,b). Comparing aged luminal AV cells from healthy nulliparous cancer-free mice and human luminal tumors, we identified 242 genes upregulated and 137 genes downregulated in both (Fig. 8a, Extended Data Fig. 10c and Supplementary Table 8b). For example, both aged cells and tumors upregulated NKD2, encoding a component of Wnt signaling; CXCL17, encoding a chemokine associated with poor survival in patients with breast cancer; and CRIP1, with both tumor suppressor and oncogene properties; whereas they decreased expression of genes encoding secreted anti-inflammatory factor FNDC4 and SEMA6A, a semaphorin with tumor suppressor activity (Fig. 8a,b). Similarly, 143 genes were upregulated in both aged mouse luminal HS cells and human luminal tumors as well as 92 downregulated in both (Fig. 8c, Extended Data Fig. 10d and Supplementary Table 8a). Among the top 25 DE genes, both aged cells and tumors upregulated genes encoding insulin-like growth factor binding protein IGFALS and stemness regulator MEX3A and downregulated TSG RCAN1 and TF genes ARID5B and SOX9 (Fig. 8c,d). Although a subset of these shared expression changes was found exclusively in luminal tumors (137 in luminal AV and 122 in luminal HS), many were also found in other breast tumor subtypes, including in basal tumors (Extended Data Fig. 10c–e and Supplementary Table 8a,b). These age-related changes identified in healthy nulliparous cancer-free aging mice and TCGA tumors were statistically significantly compared to random shared gene sets obtained via random permutations (Extended Data Fig. 10f).

a,c, Overlap between DE genes in luminal AV (a) or luminal HS (c) cells from 18M versus 3M mouse (n = 6 per age) and DE genes in TCGA human luminal A (LumA, n = 547) or luminal B (LumB, n = 207) tumors versus normal breast tissues (n = 112) (Wald test and Benjamini–Hochberg correction; log₂FC > |0.5|, P_adj < 0.05). Only the top 25 significant upregulated or downregulated genes changing in the same direction are shown. P_adj values are indicated. b,d, Examples of overlapping DE genes from a and c (red). Left: normalized values are shown from single-cell analysis of luminal AV (b) or luminal HS (d) cells from 3M versus 18M mice (two-sided t-test), along with a pie chart depicting the percentage of expressing cells versus non-expressing cells. Right: normalized values are shown for normal tissue (n = 112) and TCGA luminal A (n = 547) and luminal B (n = 207) tumors (two-sided t-test). P values are indicated. e,f, Examples of DE fibroblast (e) and T cell (f) marker genes in 3M versus 18M mice and DE genes from TCGA human luminal A (n = 547) or luminal B (n = 207) tumors versus normal breast tissues (n = 112). Left: normalized values are shown from single-cell analysis of fibroblasts or T cells from 3M versus 18M mice (two-sided t-test), along with a pie chart depicting the percentage of expressing cells versus non-expressing cells. Right: normalized values are shown for normal tissue (n = 112) and TCGA luminal A (n = 547) and luminal B tumors (n = 207) (two-sided t-test). P values are indicated. See also Extended Data Fig. 10 and Supplementary Table 8.

Full size image

Although tumors are enriched with epithelial cells, stromal cells (including fibroblasts and immune cells) infiltrate the cancerous tissue and support tumor growth. Therefore, we looked for age-related fibroblast changes in TCGA human breast tumors and identified 330 genes upregulated and 108 genes downregulated in both aged fibroblasts from healthy nulliparous cancer-free mice and human luminal tumors (Fig. 8e and Supplementary Table 8d). For example, fibroblast markers Fap and Pi16 increased and decreased, respectively, in expression in tumors and aged fibroblasts (Fig. 8e). Similarly, several pro-inflammatory and checkpoint inhibitor-related genes were upregulated with age and in tumors, including Pdcd1, Gzmk and Lag3 (Figs. 5 and 8f and Extended Data Fig. 8).

Aging is associated with widespread alterations in epigenetic, transcriptional and post-transcriptional programs across multiple cell and tissue types²². To define how mammary tissues are affected by age and how these age-related changes relate to breast cancer, we generated an aging atlas capturing both expression and chromatin accessibility changes at the single-cell level in mammary glands from younger and older healthy nulliparous cancer-free mice (https://mga.jax.org/). Our single-cell analysis revealed shifts in cell states across multiple cell types, highlighting common age-related changes, including dysregulation of cell function and identity, increased cell plasticity, decreased ribosomal expression and increased markers of inflammation and senescence. Our findings revealed novel insights into the underlying epigenetic mechanisms that might be driving some of these cellular and molecular alterations. Furthermore, using spatial transcriptomics, we inferred cell co-localization patterns predicted by single-cell analyses.

In epithelial populations, we found that cells expressing myoepithelial markers became more abundant with age in C57BL/6J mice, similar to studies using other mouse strains and ages^18,21, and reminiscent of age-dependent increases of luminal cells expressing myoepithelial markers in human breast^14,15,17. Similarly, decreases in the ratio of luminal HS versus luminal AV cells with age were consistent with recent mouse studies^18,21. Although luminal cells were shown to modestly increase with age in humans²⁰, single-cell studies uncovered that these changes might be driven by specific subsets, potentially highlighting a higher complexity of human tissue^13,17,19,20. Critically, data from aged human and mouse^{13,17,18,19,20}, including our study, converge upon misregulated expression of classical epithelial marker genes, which could drive a loss of normal cellular plasticity required for tissue repair and a gain of abnormal plasticity that could ultimately lead to cancer⁵¹.

In addition to epithelial cell compositional changes, we detected age-related transcriptional changes consistent with breast tumors. For example, aged luminal HS cells from healthy nulliparous cancer-free mice and human breast tumors increased expression of Tmprss6, Mex3a and Tph1 and decreased expression of Rcan1, Cav1 and Dmbt1. Similarly, increased Crip1, Ecm1 and Tk1 expression and decreased Fndc4, Masp1 and Sema6a expression were observed in aged luminal AV cells and tumors. Interestingly, several TSGs increased in expression and accessibility with age, suggesting that aged cells might activate tumor suppressor programs to prevent cancer⁵². This was particularly evident in aged luminal cells, which decreased in numbers with age, suggesting a possible decline in cell populations with an anti-tumor or protective role in the aging breast. Therefore, we speculate that age is associated with cellular dysregulation, leading to the expansion of cells that might contribute to cancer development. However, further studies using samples from younger patients are needed to refine the mechanistic links between breast aging and cancer, as the analysis of the TCGA dataset is limited by the number of normal samples and older median age of patients.

Our study revealed expression changes in fibroblast marker genes and upregulated senescence and CAF markers in aged fibroblasts. Although total fibroblast numbers increased in aged animals, Pi16⁺ Fib-C4 became less abundant in proportion to other fibroblast subclusters. Interestingly, Fap⁺ and Col15a1⁺ Fib-C2 increased with age, expressed genes related to adipogenicity and reduced expression of fibroblast marker Fn1. A recent study found FAP⁺ fibroblasts to be enriched in human lung tumors and PI16⁺ fibroblasts in surrounding tissue⁵³, suggesting that Fap⁺ fibroblast cells may contribute to a pro-tumorigenic environment. Our ST data further validated that Col15a1⁺ fibroblasts (Fib-C1–C3) are enriched in epithelial spots and Pi16⁺ fibroblasts (Fib-C4) in Stroma2 spots, suggesting that aging Col15a1⁺ fibroblasts might directly influence the aged microenvironment surrounding epithelial cells. As mice in our study were healthy and cancer free, it is unlikely that we detected true CAFs; however, given the links among adipose tissue, inflammation and altered ECM in cancer, these aged fibroblasts could contribute to a cancer-promoting environment. Future studies could assess whether these aged mammary fibroblasts could support tumor growth by secreting ECM as described for senescent CAFs⁵⁴ and/or directly impact cancer progression and metastasis as described in melanoma⁵⁵.

Our data point to an expansion with age of Gzmk⁺ cells, memory CD4⁺ T cells (including Tregs) and γδ T and MAIT cells in mammary tissues. This Gzmk⁺ T cell expansion is in line with age-related increases seen in other mouse tissues and in human blood⁴³. Mammary Gzmk⁺ T cells expressed markers of inflammation (Ccl5) and of exhaustion (Pdcd1). Furthermore, our ST analysis and a recent human breast atlas detected tissue-resident T cells near ductal and lobular regions¹⁹, suggesting a functional interaction with epithelial cells and changes to the aged microenvironment. The interaction between cell types may be functionally important as changes in the aged microenvironment can impact responses to immunotherapies. For example, in triple-negative breast cancer (TNBC), the aged tumor microenvironment is unable to generate a proper anti-tumor response to immune checkpoint blockade⁵⁶. In addition, loss of BRCA1 or BRCA2, which is frequent in TNBC and associated with accelerated aging^17,57, results in differential responses to checkpoint blockade, partially due to microenvironment changes⁵⁸. Moreover, a recent study showed increased lymphocytes, particularly pro-inflammatory T cells, in mammoplasties from BRCA1 and BRCA2 healthy carriers²⁰, furthering the link between immune microenvironment and age-related phenotypes. Therefore, further in-depth studies to define the therapeutic impact of age on clinical efficacy and toxicity of checkpoint blockade targeting are greatly needed, especially in luminal tumors.

Finally, we uncovered age-related epigenetic changes that can be linked to changes in cell identity and gene expression across several cell types. For example, in myeloid and T cells, open chromatin regions may contribute to the expression of critical marker genes and lineage commitment of age-related cell populations. Previous studies reported epigenetic regulation of specific lineage genes and age-related epigenetic dysregulation in blood and bone marrow⁵⁹, thus supporting our hypothesis that changes in chromatin accessibility may be, at least in part, driving the age-related population shifts and cellular dysfunction in mammary glands. Furthermore, we speculate that changes in chromatin accessibility may help to poise cells to respond to insults. For example, promoters of checkpoint inhibitors Pdcd1 and Ctla4 were highly accessible in naive T cells even though they are not frequently or highly expressed in these cells, which could be a poised epigenetic signature of checkpoint inhibitors to facilitate rapid immune responses. In addition, we identified gene expression changes that could be mechanistically explained by age-related changes in chromatin accessibility. For example, opening of the chromatin at the Pdk4 promoter with age provides a putative mechanism to explain its upregulation in older luminal AV cells detected in our study and in another mice model¹⁸. Lastly, snATAC-seq allowed us to calculate age-related TF activity scores, revealing that JUN and FOS family TFs, which regulate inflammatory responses and were previously associated with aging and cancer, had increased activity in mouse memory CD4⁺ in our study, similar to findings from Tabula Senis²¹. Future studies are needed to uncover how age-driven epigenetic changes contribute to tumor initiation and to untangle the contribution of other gene expression regulatory programs, including post-transcriptional, translational and post-translational, and their impact on the biology of aged tissues.

In sum, our mouse mammary atlas reveals similarities in cellular composition between mouse and human breast profiled in single-cell atlases^{17,19,20,60,61,62} as detailed above. Further work is needed to understand how hormonal and environmental changes impact aging of the mammary gland in model systems. In this study, we modeled breast aging in mice to enable longitudinal and well-controlled studies that are not possible in humans. Although there are structural differences between human and mouse mammary glands, both are composed of ducts comprising luminal epithelial and myoepithelial cells surrounding a lumen, suggesting structural and functional similarities between both species⁶³. Furthermore, mice exhibit age-dependent transcriptional and epigenetics changes in other tissues and are frequently used to model breast cancer. However, a limitation of the mouse model is that female mice do not experiences a notable decline in estrogen with age, similar to human menopause⁶⁴. Future studies are needed to investigate how aging mammary glands are impacted by hormonal decline using, for example, mouse or rat models of estrogen withdrawal^65,66.

Furthermore, our understanding of how parity, which affects post-partum and lifetime breast cancer risk in humans and rodents⁶⁷, might impact age-related changes in breast remains incomplete. Currently, we know that pregnancy changes the cellular composition of human breast and reprograms the chromatin landscape of mammary epithelial cells and blocks the development of premalignant lesions in mice after pregnancy^17,20,68. Moreover, post-pregnant mammary tissues exhibit an expansion of γδ NK T cells and an upregulation of immune signaling molecules (CD1d) by mammary epithelial cells, offering a protective mechanism against oncogenesis⁶⁹. However, it is unclear how parity-related molecular and cellular changes, including the changes listed above, are impacted by aging; in fact, changes in cellular abundance might go in opposite directions when comparing aging and parity²⁰. Furthermore, it remains to be determined if parity-related changes impact any of the aging-associated and cancer-associated events described in this study. Decoupling the individual and concomitant contributions of menopause, parity and aging is challenging, but future studies, in combination with the available atlases, will help elucidate the impact of these changes on the mammary gland. Toward this goal, our resource will help to guide which age-related molecular and cellular events detected in women can be effectively modeled in aging C57BL/6J mice.

Experimental model

All animal work was performed in accordance with protocols approved by the Institutional Animal Care and Use Committee at The Jackson Laboratory (no. 06005). Virgin female young (3–4 months old) and older (17–18 months old) C57BL/6J mice (JAX no. 000664), which correspond to 20–30-year-old and more than 55-year-old humans⁷⁰, respectively, were used in this study. All mice were bred and maintained in-house under regular conditions, with a 12-h light/dark cycle and food (LabDiet, 5K0G) and water ad libitum. Young and older animals were housed together for at least 1 month before tissue collection to eliminate environmental differences and to synchronize their estrous cycles.

Tissue dissociation for single-cell analysis

Each sample (n) was prepared from mammary glands pooled from three mice. To isolate mammary cells from young and aged mice, we adapted protocols from refs. ^71,72,73. Fresh mammary glands were surgically excised, finely minced and incubated in a filter-sterilized digestion solution containing DMEM/F-12 (Gibco, 11320-033), 4% heat-inactivated FBS (Seradigm, 1500-500), 1.5 mg ml⁻¹ collagenase IV (Sigma-Aldrich, C5138), 0.2% trypsin (Corning, 25-054), 5 µg ml⁻¹ insulin (Sigma-Aldrich, I-1882) and 2.5 µg ml⁻¹ gentamycin (Gibco, 15750) for 15 min at 37 °C with gentle manual agitation. The digestion was briefly centrifuged (2 min at 600g) to separate aqueous (easy-to-dissociate, sensitive cells such as immune cells) and undigested (stromal cells and epithelial cells) fractions. The aqueous fraction was pelleted (7 min at 600g), resuspended in DMEM/F-12 + 10% FBS and stored on ice. Undigested tissues (floating tissue and pelleted tissue) were collected and further digested at 37 °C for 20 min before a second centrifugation (7 min) at 600g. Pellets from both fractions were combined and further digested using a solution of DMEM/F-12 (Gibco, 11320-033) and 4 U ml⁻¹ DNAse I (Invitrogen, 18068) for 2–5 min at room temperature with gentle manual shaking. Cells were pulse centrifuged three times at 520g as previously described^71,72, and cells found in the pellet and supernatant fractions at each step were separated and collected. Epithelial cells (enriched in the pellet fractions) were further digested for 10 min at 37 °C in a solution containing TrypLE (Gibco, 12604) and 4 U ml⁻¹ DNAse I (Invitrogen, 18068) and monitored for viability and digestion using a microscope. Cells found in the supernatant fraction were subjected to red blood cell lysis in ACK buffer (Gibco, A10492) for 3–4 min. Dissociated cells from both fractions were filtered through a 70-μm cell strainer and stained with propidium iodide (PI) and calcein (Invitrogen, C1430, and BD, 556463). Live cells (calcein⁺ and PI⁻) from both fractions were isolated separately using flow cytometry (nozzle 130, flow rate 1) and collected in DMEM/F-12 with 10% FBS. Fractions were combined at a 1:1 ratio for downstream library preparations.

scRNA-seq and snATAC-seq library preparation and sequencing

Dissociated cells isolated from young and aged mice as described above were counted and assessed for viability on a Countess II automated cell counter (Thermo Fisher Scientific), and 12,000 viable cells were loaded into one lane of a 10x Chromium microfluidic chip for scRNA-seq for a targeted cell recovery of 6,000 cells per lane. From the remaining cell suspensions, nuclei were isolated according to 10x Genomics protocol (CG000212, protocol 1.2).

For scRNA-seq, single-cell capture, barcoding and library preparation were performed using 10x Chromium v.3 chemistry according to the manufacturer’s protocol (CG00183) by The Jackson Laboratory Single Cell Biology core service. cDNA and single-cell libraries for RNA-seq were checked for quality on an Agilent 4200 TapeStation, quantified by KAPA qPCR and sequenced on a NovaSeq 6000 instrument (Illumina) to an average depth of 50,000 reads per cell. For snATAC-seq, nuclei suspensions were incubated in a transposase-containing mix; nuclei were counted; and 9,250 nuclei were loaded into one lane of a 10x Chromium microfluidic chip. Single-nuclei capture, barcoding and library preparation were performed using 10x Chromium v.1 chemistry according to the manufacturer’s protocol (CG000168) by The Jackson Laboratory Single Cell Biology core service. Libraries for snATAC-seq were checked for quality on an Agilent 2100 Bioanalyzer, quantified by KAPA qPCR and sequenced on a NovaSeq 6000 instrument (Illumina) to an average depth of 25,000 read pairs per nucleus by The Jackson Laboratory Genome Technologies core service.

scRNA-seq data processing

The base call files from Illumina were demultiplexed and converted to FASTQ files using bcl2fastq (version 2.20.0.422, Illumina). The Cell Ranger (version 3.1.0, 10x Genomics) pipeline was used to align the sequence reads against the mm10 reference genome, deduplicate reads, call cells and generate cell × gene count matrices for each library. Paired-end reads were processed and mapped to the mm10 mouse genome using Cell Ranger pipeline version 4.0.0. We performed preliminary filtering of low-quality cells (gene counts <200). Cell doublets were estimated using Scrublet version 0.2.1 (ref. ⁷⁴) and DoubletDecon version 1.1 (ref. ⁷⁵). Additional filtering was applied using the Seurat package to eliminate cells with (1) gene counts less than 500 and (2) mitochondrial gene ratio higher than 10%, yielding 48,180 cells from 12 samples (n = 6 per age; Supplementary Table 1). Filtered data matrices were then analyzed using Seurat version 4 (ref. ⁷⁶). NormalizeData() was used for normalization (log normalization method). FindVariableFeatures() was used to identify highly variable features, including all genes as features. ScaleData() was used to scale data, using all genes as features. To reduce the dimensionality of the data, we ran principal component analysis (PCA) using the RunPCA() function. We used 10 principal components to define 11 clusters⁷⁷ that are annotated using annotated marker genes: B cell (Blnk, Cd79a, Cd79b), plasma cell (Jchain), T cell (Cd3d, Cd3e, Cd3g), macrophage (Itgax), DC (Fcgr2b, Cd209a, Itgam), luminal AV (Mfge8, Trf, Csn3, Wfdc18, Elf5, Ltf), luminal HS (Prlr, Cited1, Pgr, Prom1, Esr1), myoepithelial (Krt17, Krt14, Krt5, Acta2, Myl9, Mylk, Myh11), fibroblast (Col1a1, Col1a2, Col3a1, Fn1), vascular (Pecam1, Cdh5, Eng, Pdgfra, Pdgfrb, Fap) and pericyte (Rgs5, Des, Notch3)¹⁸. Neighborhood graph computing was performed with FindNeighbors(), and clusters were determined using FindClusters(). Harmony version 0.1 (ref. ⁷⁸) was used for batch correction between the two batches (batch 1 replicates 1–2–3; batch 2 replicates 4–5–6). The neighborhood graph was visualized as a uniform manifold approximation and projection (UMAP) using RunUMAP(). A doublet cluster expressing B and T cell marker genes was excluded from downstream analysis (n = 539), resulting in 47,641 cells.

snATAC-seq data processing

Paired-end reads were processed and mapped to the mm10 mouse genome using the Cell Ranger ATAC pipeline version 1.2.0. Doublets were removed using AMULET version 1.0 (ref. ⁷⁹). Simple repeats, segmental duplications and repeat masker and blacklisted regions obtained from the UCSC Genome Browser and the ENCODE project were filtered. We kept cells based on the following criteria: (1) total number of fragments in peaks >1,000 or <100,000; (2) fraction of fragments (percent reads) in peaks >40; (3) blacklist ratio <0.01; (4) nucleosome signal <4; and (5) transcription start site enrichment score >2. Filtered reads were analyzed using Signac version 1.1 (ref. ⁸⁰). The combined filtered data for downstream analysis yielded 173,699 chromatin accessibility sites (that is, peaks) associated with 22,018 genes detected in 22,842 cells (n = 3 per age; Supplementary Table 1b). Data normalization and dimensionality reduction were performed using Signac with latent semantic indexing (LSI), consisting of term frequency–inverse document frequency (TF-IDF) normalization and singular vector decomposition for dimensionality reduction. Clustering was performed using the smart local moving algorithm, and the anchors were transferred from scRNAseq to snATAC using the canonical correlation analysis (CCA) reduction method. Genome browsers were generated using igvtools version 2.5.3 (ref. ⁸¹).

Cell compositional changes with age

To quantify age-related cell proportion changes (Fig. 1 and Extended Data Fig. 1), we used a two-sided paired t-test. To calculate age-related cell ratio changes via fold change (FC) enrichment (subclustering analyses), we calculated the log-transformed ratio of cell number in old mice versus cell number in young mice and determined significance with a t-test, using 0 FC as the reference group. For the FC enrichment analysis of epithelial and fibroblast subclusters, we added 1 to every cluster to prevent division by 0.

Epithelial cell subset analyses

For scRNA-seq, epithelial cells (luminal HS, luminal AV and myoepithelial clusters; n = 7,308 cells) were subset and reprocessed using Seurat version 4. We used all genes as features in the PCA analysis. We used Harmony to correct for batch effects (batch 1 replicates 1–2–3; batch 2 replicates 4–5–6), 14 principal components and resolution 0.8 to resolve subclusters. Suspected doublets (expressing immune marker genes) and low-quality cells were filtered before final clustering, resulting in 12 subclusters (n = 6,953 cells) remaining for downstream analyses.

For snATAC-seq, epithelial cells were subset (luminal HS, luminal AV and myoepithelial clusters; n = 6,963 cells) and reprocessed using Signac. We performed LSI (n = 50) and created the UMAP using the first 20 LSI components. We excluded the first LSI component as it had a strong correlation with the total number of cell counts. These cells were assigned to cell types using the annotation from scRNA-seq using the FindTransferAnchors() function in Signac.

Stromal cell subset analysis

For scRNA-seq, stromal cells (fibroblast, pericyte and vascular clusters; n = 4,233 cells) were subset and reprocessed using Seurat version 4. We used all genes as features in the PCA analysis. We used Harmony to correct for batch effects (batch 1 replicates 1–2–3; batch 2 replicates 4–5–6), 14 principal components and resolution 0.4 to resolve subclusters. Suspected doublets (expressing immune marker genes), cells expressing putative muscle marker genes and low-quality cells were filtered before final clustering, resulting in 11 subclusters (n = 3,832 cells) remaining for downstream analyses.

For snATAC-seq, stromal cells (fibroblast, pericyte and vascular clusters; n = 2,029 cells) were subset and reprocessed using Signac. Subsequent analysis proceeded as describe above (epithelial cell subset analysis).

T cell subset analyses

For scRNA-seq, T cells (naive T cell and memory T and NK cell clusters; n = 20,455 cells) were subset and reprocessed using Seurat version 4. We calculated the top 2,000 features and ran PCA for these variable features. We used Harmony to correct for batch effects (batch 1 replicates 1–2–3; batch 2 replicates 4–5–6), 30 principal components and resolution 0.5 to resolve subclusters. The analyses yielded 11 clusters annotated using cell markers. One subcluster was filtered out from further downstream analysis because it included markers for different cell types (that is, potential doublets), resulting in 10 subclusters (n = 19,693 cells) remaining for downstream analyses.

For snATAC-seq, T cells (naive T cell and memory T and NK cell clusters; n = 6,443 cells) were subset and reprocessed using Signac. We performed LSI (n = 50) and clustered the data using the first 10 LSI components (excluding first LSI component) at resolution 0.5. This analysis generated 11 subclusters that were assigned to cell types using the annotation from scRNA-seq using the FindTransferAnchors() function in Signac. We similarly filtered out a doublet cluster, resulting in nine subclusters (n = 6,078 cells) for downstream analyses. For robust epigenomic comparisons, more than 100 cells are needed; therefore, cell-intrinsic epigenomic comparisons were not conducted if cell numbers in either age group were less than 100 (for example, Gzmk⁺ cells).

Myeloid cell subset analysis

For scRNA-seq, myeloid cells (DC and macrophage cluster; n = 3,485 cells) were subset and reprocessed using Seurat version 4. We calculated the top 2,000 features and ran PCA for these variable features. We used Harmony to correct for batch effects (batch 1 replicates 1–2–3; batch 2 replicates 4–5–6), 50 principal components and resolution 0.5 to resolve subclusters. This resulted in 13 subclusters annotated using cell markers. One of these subclusters was filtered out from further downstream analysis as potential doublets, and several of them were merged due to similar transcriptional profiles, resulting in seven subclusters (n = 3,407 cells) for downstream analyses.

For snATAC-seq, myeloid cells (DC and macrophage cluster; n = 1,509 cells) were subset and reprocessed using Signac. We performed LSI (n = 50) and clustered the data using the first 10 LSI components (excluding first LSI component) at resolution 0.2. The remaining seven subclusters were assigned to cell types using the annotations from scRNA-seq using the FindTransferAnchors() function in Signac.

Differential gene expression between old and young

For epithelial cells and fibroblasts (at the cluster and subcluster level), we used single-cell and pseudobulk differential analysis pipelines. DE analysis at the single-cell level was performed using logistic regression (LR) in the Seurat FindMarkers() function⁷⁶. Pseudobulk DE analysis was conducted as follows: https://hbctraining.github.io/scRNA-seq/lessons/pseudobulk_DESeq2_scrnaseq.html (DESeq2 version 1.34.0 (ref. ⁸²)). DESeq2 internal normalization was performed to correct for library size and RNA composition bias. We used a false discovery rate (FDR) < 0.05 cutoff to determine significance. The union of the DE genes from single-cell and pseudobulk analyses was used for downstream analysis. To functionally annotate the DE genes, hypergeometric gene set enrichment testing was conducted using cinaR version 0.2.4 (ref. ⁸³) or the online Enrichr interface⁸⁴ and databases including Kyoto Encyclopedia of Genes and Genomes (KEGG), WikiPathways, Gene Ontology (GO) and Molecular Signatures Database (MSigDB) Hallmark. To identify tumor suppressors, we used TSGene 2.0 (https://bioinfo.uth.edu/TSGene/). For violin plots comparing 18 M versus 3 M expression (for example, Fig. 2f–k), we plotted the log-normalized counts store in [[“RNA”]]@data in the Seurat object, and cells with zero counts were filtered. To visualize cells with zero counts, a pie chart was included depicting cells expressing and not expressing the gene.

For T cell and myeloid cluster-specific DE analysis, we conducted one comparison using a Wilcoxon rank-sum test (logfc.threshold = 0.25 and min.pct = 0.1). Hypergeometric gene set enrichment testing of age-specific and cell-specific DE genes was performed using phyper and WikiPathways and immune system–related modules^23,85 from CinaRgenesets⁸³. The enrichment results were adjusted using the Benjamini–Hochberg FDR adjustment method (FDR = 10%).

Differential peak analysis between old and young

For epithelial cells and fibroblasts, cell-type-specific alignment files were obtained per sample using the sinto package version 0.8.0 (https://github.com/timoast/sinto). MACS2 version 2.2.7.1 (ref. ⁸⁶) was used for peak calling for each sample per cell type using the BAMPE option. We used the Diffbind version 3.4.11 (ref. ⁸⁷) package to generate the consensus sequences per cell type, which were then used for differential analyses between young and old using DESeq2 with absolute log₂FC ≥ 1 in the cinaR R package⁸³. Differential peaks were annotated using ChIPseeker version 1.30.3 (ref. ⁸⁸). In cases where ChIPseeker was unable to annotate the peaks to the nearest gene, we used HOMER version 4.11 (ref. ⁸⁹).

For single-cell differential peak analysis in immune cell subclusters, we used the Signac FindMarkers() function to calculate age-related differential accessibility. We used a Wilcoxon rank-sum test with a minimum of 10% cells accessible to the peak as our cutoff. Peaks with P_adjusted > 0.05 were filtered out from downstream analysis.

ChromVar analyses

We added motif information to the snATAC-seq object using the AddMotifs() function in Signac for the mm10 genome using the JASPAR 2020 database⁹⁰. We calculated a per-cell motif activity score using chromVAR⁹¹ and added this information to the snATAC-seq object. We used these motif activity scores to conduct differential analysis using the FindMarkers() function in Signac between old and young mice using the Wilcoxon rank-sum test with no cutoffs being used for FC or minimum percentage of cells expressing the motif.

Senescence and CAF scoring

We calculated the sum of normalized counts (slot = ‘data’) for genes in module and divided by the total number of transcripts per cell. We then scaled this ratio between 0 and 1 and plotted the scaled scores using the VlnPlot() function. We scored fibroblasts for the CAF signatures⁹² and T cells for the SenMayo senescence signature⁹³. For fibroblasts, we scored the senescence signature using SenePy⁹⁴ and senescence marker genes identified in fibroblasts of several mouse tissues (n = 939 genes; Supplementary Table 4i).

Ligand–receptor interactions

CellPhoneDB⁹⁵ was used to calculate ligand–receptor interactions between different cell types in young and old mice. The normalized cell counts and metadata were extracted from the scRNA-seq object. The results were obtained by running the statistical_analysis() function available in CellPhoneDB. The analysis was run separately for old and young mouse cells.

TCGA analyses

To uncover DE genes from human TCGA data, we downloaded tumor and normal tissue samples from TCGAbiolinks^96,97,98 and tumors that passed quality testing as described⁹⁹ (normal tissue n = 112; luminal A tumors n = 547; luminal B tumors n = 207; normal-like tumors n = 39; basal tumors n = 183; Her2 tumors n = 77). DE analysis was performed using DESeq2 (ref. ⁸²). A Wald test was used to calculate P values, and the Benjamini–Hochberg procedure was used to calculate corrected P values. Differential genes were selected based on P_adj < 0.05 and log₂ FC > |0.5|.

To assess the specificity of the associations between age and cancer, we shuffled the tumor and normal labels of the TCGA-BRCA samples and newly defined their groups (10,000 times). From the permutated data, we identified significant genes (P_adj < 0.05 and log₂FC > |0.5|) repeatedly. We then compared these gene lists with the mouse aging-related DE genes to capture shared gene sets (also considering FC direction). From these 10,000 shared gene sets, we demonstrated that our real association set contains significantly more shared genes compared to the number of shared genes of 10,000 random trials. To verify robustness of the top 10 enriched pathways of our real associations, we selected 1,000 shared gene lists and conducted a GO test for each gene list. Using these 1,000 results of the GO test, we assumed empirical distribution of corrected P values for each enriched term and calculated empirical P values our observed event.

ST profiling of aged mammary tissues

ST experiments were performed using the Visium Platform (10x Genomics) according to the manufacturer’s protocols. Fresh mammary tissues from 3-month-old or 18-month-old mice were formalin-fixed and paraffin-embedded (FFPE). Two 15-µm sections from each tissue block were used for total RNA extraction with an Rneasy Micro Kit (Qiagen), and the percentage of RNA fragments larger than 200 base pairs (bp), as determined by an Agilent Bioanalyzer system (DV200 score), was used as a measure of RNA quality. Tissue blocks with DV200 scores were above 50%.

Sections from each tissue block were placed on microscopy slides (Thermo Fisher Scientific, ColorFrost Plus) and subsequently deparaffinized, H&E stained and imaged in bright-field using a NanoZoomer SQ (Hamamatsu) slide scanner. Each slide was incubated with mouse-specific probe sets provided by the manufacturer for subsequent mRNA labeling, probe transfer using CytAssist (10x Genomics) onto a Visium CytAssist slide and subsequent library generation per the manufacturer’s protocol (10x Genomics, CG000495). Library concentration was quantified using a TapeStation High Sensitivity DNA ScreenTape (Agilent) and fluorometry (Thermo Fisher Scientific, Qubit) and verified using KAPA qPCR. Libraries were pooled for sequencing on an Illumina NovaSeq 6000 200-cycle S4 flow cell using a 28-10-10-90 read configuration, targeting 100,000 read pairs per spot covered by tissue.

Illumina base call files for all libraries were converted to FASTQs using bcl2fastq version 2.20.0.422 (Illumina). For each tissue section and corresponding library, the whole-slide bright-field image and the CytAssist image were aligned manually using the Loupe Browser (version 6.4.1) with landmark registration. Each whole-slide image was uploaded to a local OMERO server where a rectangular region of interest (ROI) containing just the tissue was drawn using OMERO.web, and OMETIFF images of each ROI were generated using the OMERO Python API. FASTQ files, the image registration JSON file and associated OMETIFF corresponding to high-resolution bright-field images were used for further processing, including alignment to the GRChm38 mm10-specific filtered probe set (10x Genomics Mouse Probe Set version 1.0.0) using the version 2.1.0 Space Ranger count pipeline (10x Genomics).

ST data analysis

Sequencing reads from the Visium Spatial Gene Expression slide were pre-processed with Scanpy version 1.9.3 (ref. ¹⁰⁰). Spots with fewer than 245 genes and more than 10% of mitochondrial gene expression (defined using MitoCarta2 version 2.0 (ref. ¹⁰¹)) were filtered. Additionally, we filtered out rarely expressed genes (<3 reads per spot). As a result, filtered data from older samples contained 6,428 ST spots (18 M ST1) and 8,662 ST spots (18 M ST2), and younger samples contained 4,383 ST spots (3 M ST1) and 4,011 ST spots (3 M ST2), with 19,203 genes.

Raw read counts in each sample were normalized using the Pearson residual method¹⁰². We used BBKNN version1.5.1 (ref. ¹⁰³) to remove batch effects and identified clusters using the Leiden clustering method (resolution 0.15). Young samples had different neighbor distances due to differences in gene expression compared to the old samples; therefore, we applied a different Leiden clustering (resolution 0.1). We defined the dominant cell type in each cluster by identifying marker genes in each cluster compared to other clusters using log-normalized counts.

To test co-localization of specific cell types and ligand–receptor pairs in the ST data, we transformed the read counts of marker genes (for cell type and ligand–receptor pairs) into binary values (positive and negative). Consequently, each spot was assigned a value of 1 (positive signal) if it had at least one read for the marker genes. Using this binary count matrix, we defined Cd3d or Cd3g or Cd3e or Cd247 or Cd8a or Cd8b1 as Cd8⁺ cells and Cd3d or Cd3g or Cd3e or Cd247 as Cd3⁺ cells. Based on these definitions, we considered a spot to be Cd3⁺ or Cd8⁺ and Gzmk⁺ if it had a positive signal for both Cd3⁺ or Cd8⁺ markers and Gzmk. We applied the same process to other target cells, such as Cd3⁺ or Cd8⁺ and Pdcd1 and γδ T cells. After that, the number of target cell spots was normalized by the total count of spots in each cluster to compare their enrichment between different clusters. To test putative ligand–receptor interactions, we considered a spot to be a specific ligand–receptor ‘positive’ if it had a positive signal for both ligand and receptor genes. We applied the same normalization method as used for specific cell type spots to compare their enrichment between different clusters. Additionally, using this ligand–receptor spot information, we measured the co-localization of ligand–receptor pairs and specific cell type spot. In that case, we considered a spot to be co-localized if it had a positive signal for both the ligand–receptor pair and the specific cell type or was close to both signals.

Statistics and reproducibility

No statistical method was used to predetermine sample size, but our sample sizes are similar to those reported in previous publications^18,21. No data were excluded from the analyses. The experiments were not randomized. The investigators were not blinded to allocation during experiments and outcome assessment. Data distribution was assumed to be normal, but this was not formally tested.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.