Luminal epithelial cells integrate variable responses to aging into stereotypical changes that underlie breast cancer susceptibility

Abstract

Effects from aging in single cells are heterogenous, whereas at the organ- and tissue-levels aging phenotypes tend to appear as stereotypical changes. The mammary epithelium is a bilayer of two major phenotypically and functionally distinct cell lineages: luminal epithelial and myoepithelial cells. Mammary luminal epithelia exhibit substantial stereotypical changes with age that merit attention because these cells are the putative cells-of-origin for breast cancers. We hypothesize that effects from aging that impinge upon maintenance of lineage fidelity increase susceptibility to cancer initiation. We generated and analyzed transcriptomes from primary luminal epithelial and myoepithelial cells from younger <30 (y)ears old and older >55 y women. In addition to age-dependent directional changes in gene expression, we observed increased transcriptional variance with age that contributed to genome-wide loss of lineage fidelity. Age-dependent variant responses were common to both lineages, whereas directional changes were almost exclusively detected in luminal epithelia and involved altered regulation of chromatin and genome organizers such as SATB1. Epithelial expression variance of gap junction protein GJB6 increased with age, and modulation of GJB6 expression in heterochronous co-cultures revealed that it provided a communication conduit from myoepithelial cells that drove directional change in luminal cells. Age-dependent luminal transcriptomes comprised a prominent signal that could be detected in bulk tissue during aging and transition into cancers. A machine learning classifier based on luminal-specific aging distinguished normal from cancer tissue and was highly predictive of breast cancer subtype. We speculate that luminal epithelia are the ultimate site of integration of the variant responses to aging in their surrounding tissue, and that their emergent phenotype both endows cells with the ability to become cancer-cells-of-origin and represents a biosensor that presages cancer susceptibility.

Keywords: aging; breast cancer; cancer biology; cancer susceptibility; computational biology; gene expression variance; human; lineage fidelity; lumininal epithelia; systems biology.

Figure 1.. Genome-wide loss of lineage-specific expression with age.

(A) Immunofluorescence staining of normal breast tissue showing the mammary epithelium with an apical LEPs (KRT19) surrounded by basal MEPs (KRT14). (B) Representative FACS enrichment plot of HMECs stained with LEP-specific CD133 (PROM1) and MEP-specific CD271 (NGFR). (C) Boxplots of subject-level gene expression rlog values of canonical LEP-specific markers EPCAM and ERBB2 (top), and MEP-specific markers TP63 and EGFR (bottom) in FACS isolated LEPs and MEPs in all age groups. Lineage-specific DE adj. p-values annotated (*<0.05, **<0.01, ***<0.001, ****<0.0001). (D) DE LEP-specific and MEP-specific genes (adj. p<0.001, abs(lfc) ≥1) in younger <30 y (left) and older >55 y (right) women. Strata plot shows changes in lineage-specific DE with age, showing the number of LEP- and MEP-specific genes gained (cyan and magenta), lost (light and dark gray), and maintained (green and red) in older women. Number of subjects (n) and sample replicates (m) in each DE analysis annotated; number of DE genes (k) in each age group and DE status indicated. (E) Distribution of lfc in expression between LEPs and MEPs in younger and older subjects for either DE LEP-specific (top panel) or MEP-specific (bottom panel) genes. KS p-values for equality of distributions of lfc between younger and older women annotated. (F–G) Histogram of pairwise differences in lfc in expression between LEPs and MEPs in younger vs. older women for (F) all genes with lineage-specific expression in younger women or (G) only genes that maintain lineage-specific expression in older women. Genes with LEP-specific and MEP-specific expression are shown in the top and bottom panels respectively. The percent of genes with higher lfc in younger women (light blue) or higher lfc in older women (blue gray) are indicated. One-sided t-test p-values annotated.

Figure 1—figure supplement 1.. Isolation of luminal and myoepithelial lineages.

(A) Representative FACS enrichment plot of HMECs stained with LEP-specific CD227 (MUC1) and MEP-specific CD10 (MME). (B–C) Individual FACS enrichment plots of HMEC strains from (B) younger and (C) older women stained with LEP-specific CD133 (PROM1) and MEP-specific CD271 (NGFR). (D) Unsupervised hierarchical clustering of all LEP and MEP samples based on genome-wide expression of 17,328 genes. Individual LEP and MEP strains are annotated, and replicate bridge samples enriched across both combinations of FACS markers and sequenced across RNA-seq batches are numerically indicated. Approximately unbiased (AU) p-values and bootstrap probability (BP), which assesses uncertainty in hierarchical clustering analysis, are annotated. Clusters with AU p≥0.95 are highlighted in red, with largest clusters in solid red. (E) Boxplots of subject-level gene expression rlog values of canonical LEP-specific markers KRT19, MUC1, CDH1, KIT, ELF5, KRT8, and KRT18 (top), and MEP-specific markers KRT14, THY1, ITGA6, KRT5, SERPINB5, CNN2, CALD1 (bottom) in FACS isolated LEPs and MEPs in all age groups. Lineage-specific DE adj. p-values annotated (*<0.05, **<0.01, ***<0.001, ****<0.0001). (F) FACS enrichment plots showing isolated CD133+/CD271- LEP cells staining positive for KRT19, and CD133-/CD271 +MEP cells are staining positive for KRT14. (G) Representative FACS enrichment plot of dissociated uncultured organoids stained with LEP-specific CD133 (PROM1) and MEP-specific CD271 (NGFR). (H) Boxplots of subject-level gene expression regularized log (rlog) values of canonical LEP-specific marker PROM1 (top), and MEP-specific marker NGFR (bottom) in FACS isolated LEPs and MEPs from uncultured organoids (left) and HMECs (right) in all age groups. Lineage-specific DE adj. p-values annotated (*<0.05, **<0.01, ***<0.001, ****<0.0001).

Figure 1—figure supplement 2.. Luminal and myoepithelial lineages in organoids and human mammary epithelial cell cultures.

(A–D) Pairwise comparison of subject-level rlog gene expression means between primary organoids and 4th passage in FACS-enriched (A, C) LEP or (B, D) MEP cells isolated from finite-lifespan HMECs derived from reduction mammoplasties of younger (A, B) or older (C, D) women. Linear regression line (dotted red) and standard error (red) shown with regression R², p-value, slope and y-intercept annotated. Distribution of residuals shown in the inset. Number of subjects (n) in each group annotated. Established lineage-specific markers from (Figure 1—figure supplement 1E) are shown as reference.

Figure 2.. Loss of lineage fidelity with age leads to disrupted lineage-specific signaling.

(A) Schematic illustrating cell-cell communication and dysregulated lineage-specific signaling with age either through loss of lineage-specific expression of the ligand (L) or its cognate receptor (R). (i) MEP_L→ LEP_R and LEP_L→ LEP_R lineage-specific signaling and active pathways (>>) associated with both ligands and receptors in <30 y young cells; (ii) LEP-specific expression of the receptor is lost with age leading to disrupted MEP_L→ LEP_R and LEP_L→ LEP_R lineage-specific signaling. LEP receptor pathway, as well as MEP ligand and LEP ligand pathways are detected as dysregulated by functional enrichment analysis methods; and (iii) Lineage-specific expression of the MEP ligand is lost with age and MEP_L→ LEP_R lineage-specific signaling is disrupted. MEP ligand pathway is dysregulated. MEP-directed LEP receptor pathway is also dysregulated as cell-cell signaling homeostasis is shifted and only LEP → LEP signaling is driving the LEP receptor pathway. (B) Interactome map of DE lineage-specific ligand-receptor pairs (LRPs) (adj. p<0.001, abs(lfc) ≥1) that show loss of lineage-specific expression of either ligands and/or their cognate receptors in older LEPs (light gray) or MEPs (dark gray). LRPs are connected by chord diagrams from the cell type expressing the ligand (L) to the cell type expressing the cognate receptor (R). Number of LRPs, and genes (k) in each category annotated. Summary of functionally enriched KEGG pathways (FDR p<0.001) associated with loss of lineage-specific DE in ligands and receptors are shown.

Figure 2—figure supplement 1.. Loss of lineage fidelity with age leads to disrupted lineage-specific signaling.

(A) Interactome map of the ligand-receptor pairs (LRPs) in younger women based on lineage-specific DE of ligands and their cognate receptors in LEPs (green) or MEPs (red) (adj. p<0.001, abs(lfc) ≥1). LRPs are connected by chord diagrams from the cell type expressing the ligand (L) to the cell type expressing the cognate receptor (R). Number of LRPs, and genes (k) in each category annotated. (B–C) Network functional enrichment of top KEGG pathways (FDR p<0.001) associated with lineage-specific DE of (B) ligands and/or (C) cognate receptors in LEPs and MEPs in younger women. (D–E) Network functional enrichment of top KEGG pathways (FDR p<0.001) associated with loss of lineage-specific DE of (D) ligands and/or (E) cognate receptors in LEPs and MEPs in older women.

Figure 3.. The luminal lineage is a hotspot for age-dependent directional changes.

(A) Models of loss of lineage fidelity illustrate hypothesized mechanisms leading to loss of lineage fidelity: (i) Age-dependent DE shifts in gene expression in LEPs and/or MEPs of older relative to younger cells; or (ii) An increase in gene expression variance in older LEPs and/or MEPs that lead to loss of detection of lineage-specific DE between LEPs and MEPs with age. (B) Number of DE genes (adj. p<0.05) between younger and older LEPs or MEPs. Number of subjects (n) and sample replicates (m) in DE analysis annotated. (C–D) Hierarchical clustering of all LEP and MEP samples based on sample-level expression of age-dependent (C) DE genes in LEPs (adj. p<0.05) and (D) DE transcription factors in LEPs (adj. p<0.05). Number of DE genes (k) indicated. Gene expression scaled regularized log (rlog) values are represented in the heatmap; clustering performed using Euclidean distances and Ward agglomerative method. (E) Distribution of lfc in expression between LEPs and MEPs in younger and older women for LEP-specific (top panel) or MEP-specific (bottom panel) genes that are lost with age (DE adj. p<0.001, abs(lfc) ≥1) and that have at least a two-fold age-dependent directional change in the older cohort. KS p-values annotated. (F) Number of DV genes (adj. p<0.05) between younger and older LEPs or MEPs. Number of subjects (n) in DV analysis annotated. (G) Distribution of lfc in expression between LEPs and MEPs in younger and older women for LEP-specific (top panel) or MEP-specific (bottom panel) genes that are lost with age (DE adj. p<0.001, abs(lfc) ≥1) and that have at least a two-fold age-dependent increase in variance in the older cohort. KS p-values annotated. (H–I) MSigDB Hallmark gene sets identified by GSEA to be enriched (adj. p<0.05) in younger and older (H) LEPs and (I) MEPs based on age-dependent DE (top) or DV (bottom).

Figure 3—figure supplement 1.. The luminal lineage is a hotspot for age-dependent directional changes.

(A) Boxplots of subject-level gene expression rlog values of five age-dependent DE genes identified in both epithelial lineages: (i) LRRC4, (ii) PSORS1C1, (iii) SCNN1B, (iv) ZNF518B, and (v) ZNF521 in younger and older LEPs and MEPs with age-dependent DE adj. p-values. Number of subjects (n) and sample replicates (m) in DE analysis annotated. (B) Hierarchical clustering of all LEP and MEP samples based on sample-level expression of age-dependent DE genes in MEPs (adj. p<0.05). Number of DE genes (k) indicated. Gene expression scaled regularized log (rlog) values are represented in the heatmap; clustering performed using Euclidean distances and Ward agglomerative method. (C) Hierarchical clustering of isolated LEP and MEP samples from uncultured organoids based on expression of 495 age-dependent DE genes identified in 4th passage LEPs and MEPs. Individual LEP and MEP organoid strains are annotated. Approximately unbiased (AU) p-values and bootstrap probability (BP), which assesses uncertainty in hierarchical clustering analysis, are annotated. Clusters with AU p≥0.95 were highlighted in red, with largest clusters in solid red. (D) Boxplots of subject-level gene expression rlog values of TFs: (i) GRHL2, (ii) SGSM2, (iii) HES4, (iv) ZNF827, (v) SATB1, (vi) SP3, and (vii) ZNF503 in younger and older LEPs with age-dependent DE adj. p-values. Number of subjects (n) and sample replicates (m) in DE analysis annotated.

Figure 3—figure supplement 2.. Comparisons of age-dependent directional and variant responses in the luminal and myoepithelial lineages with obesity and parity signatures.

Upset plots showing intersection of (A) age-dependent DE genes in LEPs (adj. p<0.05), (B) age-dependent DE genes in MEPs (adj. p<0.05), (C) age-dependent DV genes in LEPs (adj. p<0.05), and (D) age-dependent DV genes in MEPs (adj. p<0.05) identified in this study and genes previously shown to be associated with obesity (Burkholder et al., 2020), and parity and time since full time pregnancy (FTP) Santucci-Pereira et al., 2019 in normal breast tissues. (A–D) 41 of 70 genes identified to be associated with obesity found in the LEP/MEP dataset (bottom panel 2nd row); 1 gene found to be DE with age in LEP, and 1 gene found to be DV in LEP annotated (top panel, 1st column); (A–D) 23 of 43 genes identified to be associated with parity found in the LEP/MEP dataset (bottom panel 3rd row); 1 gene found to be DE with age in LEP annotated (top panel, 2nd column); and (A–D) 185 of 287 genes identified to be associated with time since FTP found in the LEP/MEP dataset (bottom panel 4th row); 13 genes found to be DE with age in LEP are listed, 1 gene found to be DV in LEP, and 2 genes found to be DV in MEP annotated (top panel, 3rd column).

Figure 3—figure supplement 3.. Aging-associated increase in variance contributes to loss of lineage fidelity.

(A–D) Gene expression (A–B) means and (C–D) variances of LEPs and MEPs from younger women are categorized into quantile levels: very low, low, intermediate, high, very high. Corresponding categories of gene expression means and variances in older cells, as defined by threshold values in younger cells, are fractionally represented. (E–F) TF with significant increase in variances in older cells (DV adj. p<0.05) are shown with subject-level rlog gene expression (E) means in rank order from lowest to highest expressed in the younger cohort and (F) corresponding variances. Number of subjects (n) in analysis annotated. (G) Boxplots of subject-level gene expression rlog values of TFs: (i) EHF, (ii) KDM2B, (iii) HES4, (iv) MYCL, (v) GLI1, and (vi) DMRTA1 in younger and older LEPs, and (vii) HES6 in younger and older MEPs with age-dependent DV adj. p-values. Number of subjects (n) in DV analysis annotated.

Figure 4.. Age-dependent gene expression changes detected at the single cell level.

(A–C) UMAP projections of scRNA-seq data from (A) 13 non-tumorigenic breast tissue samples (19-69y) from reduction mammoplasties via (Pal et al., 2021), (B) 11 non-tumorigenic non-carrier (24-50y) breast tissue samples from reduction mammoplasties, prophylactic mastectomies, and contralateral to DCIS/tumor via (Nee et al., 2023), and (C) 28 healthy reduction mammoplasty tissue samples (19-39y) via (Murrow et al., 2022). Reanalysis of scRNA-seq data focused on three identified epithelial cells types, with UMAP projections showing single cells identified as Luminal1 (luminal adaptive secretory precursors), Luminal2 (luminal hormone-sensing), or Basal (basal-myoepithelial) (top left panel), and normalized expression levels of LEP-specific markers PROM1 and MUC1 (top middle and right panels), ER + LEP-specific marker AR (bottom left panel), and MEP-specific markers NGFR and MME (bottom middle and right panels). (D–F) Distribution of scGSEA signature scores capturing enrichment of age-dependent DE gene sets across cell types (boxplots) and clinical conditions (split violin plots). Distribution of signature scores showing enrichment of age-dependent DE genes (D) upregulated (left) or downregulated in old LEPs across cell types and between premenopausal and menopausal women via (Pal et al., 2021); (E) upregulated in old LEPs between women with cancer who received chemotherapy and those that did not (left), or downregulated in old LEPs between non-tumorigenic BRCA1 mutation carriers and non-carriers (right) via (Nee et al., 2023); and (F) upregulated in old LEPs between parous and nulliparous (left), or obese and nonobese (right) individuals via (Murrow et al., 2022). Standardized mean difference was used to compare mean signatures between cell types and biological groups, and to indicate significance levels.

Figure 4—figure supplement 1.. Lineage-specific gene expression changes detected at the single cell level.

scRNA-seq analysis of (A, D) 13 non-tumorigenic breast tissue samples (19-69y) from reduction mammoplasties via (Pal et al., 2021), (B, E) 11 non-tumorigenic non-carrier (24-50y) breast tissue samples from reduction mammoplasties, prophylactic mastectomies, and contralateral to DCIS/tumor via (Nee et al., 2023), and (C, F) 28 healthy reduction mammoplasty tissue samples (19-39y) via (Murrow et al., 2022). (A–C) Violin plots show distribution of quality control variables: ‘n-Feature’ which shows the number of detected genes (left); ‘n-Count’ which shows the number of detected unique molecular identifier (UMI) (middle); and ‘percent.mito’ which shows the percentage of mitochondrial genes (right) for every cell. (D–F) Violin and boxplots of scGSEA signature scores capturing enrichment of lineage-specific DE gene sets, with LEP-specific (left) and MEP-specific (right) enrichment scores shown across the three epithelial cell types Luminal1 (inclusive of luminal progenitors), Luminal2 (luminal hormone-sensing), or Basal (basal-myoepithelial). (G) UMAP projection of Nee et al., 2023 dataset by subject ID (only noncarriers shown for simplicity). (H) Fractional distribution of annotated cell states from Nee et al., 2023 by cell type and subject ID.

Figure 5.. Age-dependent directional changes in the luminal lineage are indicators of aging breast tissue.

(A–B) GSEA enrichments plots from age-dependent DE analysis of bulk tissue showing gene ranks based on DE test statistics and gene set enrichment scores. Age-dependent enrichment of two gene sets composed of (A) differentially upregulated genes in younger <30 y LEPs and (B) differentially upregulated genes in older >55 y LEPs in bulk tissue are shown. Negative enrichment scores indicate upregulation of specified gene set in tissue from younger <30 y women, while positive enrichment scores indicated upregulation in tissue from older >55 y women. GSEA enrichment BH adj. p-values are annotated. (C) Boxplots of SATB1 gene expression: (i) subject-level rlog values in LEPs and MEPs of younger and older women; (ii) log₂ values in normal breast tissue (GSE102088); and (iii) log₂ FPKM values in the TCGA breast cancer cohort by PAM50 subtype in cancers and in matched normal tissue. Age-dependent DE adj. p-values in LEPs and normal breast tissue, and lineage-specific DE adj. p-values in LEPs and MEPs are indicated (i–ii). Kruskal-Wallis (KW) test performed across PAM50 breast cancer subtypes as well as matched normal tissue; p-values adjusted across all age-dependent DE and DV genes identified in LEPs found in the TCGA dataset (iii). Post hoc pairwise Wilcoxon test adj. p-value significance levels annotated (*<0.05, **<0.01, ***<0.001, ****<0.0001); p-values adjusted across all pairwise comparisons and across all DE and DV genes identified in LEPs found in the TCGA dataset (iii). Number of subjects (n) and sample replicates (m) in each analysis annotated. (D–E) Venn diagram of genes with age-dependent DE in LEPs (adj. p<0.05) and at least nominal DE (unadj. p<0.05) in normal primary breast tissue. Genes commonly (D) upregulated and (E) downregulated in LEPs and bulk tissue with age are listed. (F) PPI network of common age-dependent DE genes in LEPs (adj. p<0.05) and bulk tissue (unadj. p<0.05) with TFs annotated in bold. Seven gene communities identified; corresponding network functional enrichment (FDR p<0.05) of selected processes annotated.

Figure 5—figure supplement 1.. Age-dependent directional changes in the luminal lineage are indicators of aging breast tissue.

(A) Boxplot of SATB1 gene expression: log₂ FPKM values in the TCGA breast cancer cohort by age at diagnosis in matched normal samples. Wilcoxon p-values between groups annotated. (B) Boxplots of log₂ FPKM values of SATB1 in paired matched normal and primary tumor across each PAM50 subtype in TCGA. Paired Wilcoxon test adj. p-value significance annotated; p-values adjusted across subtypes. Number of subjects (n) in each analysis annotated. (C) Heatmap of gene expression correlation matrix of SATB1 and published SATB1 gene targets identified in the MDA-MB-231 breast cancer cell line (Han et al., 2008). Clustering performed using complete agglomerative method with 1 – Pearson correlation as distances. (D) Venn diagram showing intersection of age-dependent DE genes in LEPs (adj. p<0.05) identified in this study and 515 of 778 published SATB1 gene targets found in the LEP/MEP dataset. Twenty-six genes were found to be DE with age in LEP with 15 genes annotated as SATB1-activated, and 16 genes as SATB1-repressed (5 genes identified as both SATB1-activated and repressed). (E) Boxplots of ssGSEA signatures scores capturing enrichment of age-dependent DE SATB1-activated (left) and SATB1-repressed (right) gene sets in LEPs of younger and older women. Wilcoxon adjusted p-values between groups annotated.

Figure 5—figure supplement 2.. Age-dependent directional changes in the luminal lineage across breast cancer subtypes and matched normal tissue in TCGA.

(A) Boxplots of log₂ FPKM values of (i) BCL11A, (ii) MLLT3, (iii) ZNF521, (iv) PHC1 and (v) PCGF3 in the TCGA breast cancer cohort. Kruskal-Wallis (KW) test performed across PAM50 breast cancer subtypes as well as matched normal tissue; p-values adjusted across all age-dependent DE and DV genes identified in LEPs found in the TCGA dataset. Post hoc pairwise Wilcoxon test adj. p-value significance levels annotated (*<0.05, **<0.01, ***<0.001, ****<0.0001); p-values adjusted across all pairwise comparisons and across all DE and DV genes identified in LEPs found in the TCGA dataset. (B) Boxplots of log₂ FPKM values of (i) BCL11A, (ii) MLLT3, (iii) ZNF521, (iv) PHC1 and (v) PCGF3 in paired matched normal and primary tumor across each PAM50 subtype in TCGA. Paired Wilcoxon test adj. p-value significance annotated; p-values adjusted across subtypes. Number of subjects (n) in each analysis annotated.

Figure 6.. Gap Junction protein GJB6 is a mediator of the non-cell autonomous mechanism of aging in breast epithelia co-cultures.

(A) Age-dependent modulation of apical junction-associated genes in LEPs and MEPs (Lepage test p<0.05). Number of age-modulated genes (k) indicated. Genes annotated with their respective HUGO Gene Nomenclature Committee (HGNC) gene family and functional role in adherens junctions, cell adhesion molecules, desmosomes, tight junctions, or gap junctions. (B) Boxplot of GJB6 subject-level rlog expression values in LEPs and MEPs in younger and older women. Lepage test p-values are indicated. Number of subjects (n) in analysis annotated. (C) Schematic illustrating integration of directional and variant responses in older epithelial cells. Different genes are dysregulated in LEPs and MEPs of older individuals leading to an increase in variance in expression across aged cells. Through cell-cell signaling, variant responses in MEPs (gene A or gene B) lead to variant responses in LEPs (gene A’ or gene B’). Where these variant responses integrate and affect common downstream genes in LEPs (gene C’) lead to detectable age-dependent directional changes (***) that are seen as stereotyped responses in the lineage. (D) Relative expression of ELF5 in younger LEPs co-cultured with either younger (Y/Y) or older (Y/O) MEPs. Two-tailed t-test p-value indicated. (E) Schematic of co-culture methodology with HMEC cells from younger and older women enriched by FACS for LEPs and MEPs; GJB6 knock-down in older MEP feeder layer by shRNA; younger LEPs are co-cultured on top of the older MEP feeder layer for 10 days; LEPs separated from MEPs by FACS; and LEP expression levels measured by qPCR. (F) Relative expression of ELF5 in younger LEPs co-cultured with either shControl or shGJB6 older MEPs. Two-tailed t-test p-value indicated. Number of subjects (n) and technical replicates (l) annotated.

Figure 6—figure supplement 1.. Genes encoding for cell-cell junction proteins are dysregulated in aging epithelia.

(A) Boxplots of subject-level rlog expression values of (i) DSC3, (ii) DSG3 (iii) CLND10, and (iv) CLDN11, in LEPs and MEPs from younger and older women. Age-dependent DE adj. p-values are indicated. Number of subjects (n) and sample replicates (m) in DE analysis annotated. (B) Linear regression of GJB6 rlog expression in MEPs vs. LEPs across both age groups. Pearson correlation coefficient (R) and regression p-value are indicated. (C) GJB6 ChIP-seq (Cistromics) binding signal within +/-10 from the TSS in mammary gland (SPP Ominer database). (D) Boxplots of subject-level gene expression rlog values of canonical LEP-specific marker ELF5 in FACS isolated LEPs and MEPs from uncultured organoids (left) and HMECs (right) in both age groups. Lineage-specific DE adj. p-values annotated (****<0.0001). (E) Relative expression of GJB6 in either shControl or shGJB6 older MEPs in two individuals. Two-tailed t-test p-value indicated. Number of subjects (n) and technical replicates (l) annotated. (F) Relative expression of ELF5 in younger LEPs co-cultured with older MEPs treated with increasing concentrations of 18aGA compared to Y/Y and Y/O treated with DMSO. Number of technical replicates (l) annotated.

Figure 7.. Age-dependent dysregulation in LEPs shape predictors of normal breast tissue and PAM50 subtypes.

(A) Unsupervised hierarchical clustering of TCGA samples from matched normal and primary tumor tissue based on expression of age-dependent DE and DV genes identified in LEPs (adj. p<0.05). PAM50 intrinsic subtypes and patient age at diagnosis are annotated. Gene expression scaled log₂ FPKM values are represented in the heatmap; clustering performed using Euclidean distances and Ward agglomerative method. (Note: extreme outlier values are set to either the minimum or maximum value of the scale bar). Number of age-dependent DE and DV genes (k) included in analysis annotated. (B) Number of individuals in each ML class in the training, test and independent test sets. (C–D) Multi-class classification model performance in predicting normal tissue and breast cancer breast cancer subtypes in the (C) TCGA test set and (D) GTEx/GSE81540 independent test set. Macro AUC, micro AUC, and AUC of each group vs. rest are shown.

Figure 7—figure supplement 1.. Breast cancer expression levels of aging biomarkers identified in normal luminal epithelia are associated with PAM50 subtypes in TCGA.

Median expression levels of age-dependent DE and DV genes from LEPs were assessed in the TCGA cohort across PAM50 LumA, LumB, Her2 and Basal intrinsic subtypes. (A) Number and fractional distribution of genes by PAM50 subtype and by highest and lowest median expression levels for: DE genes upregulated in young LEP (top left); DE genes upregulated in old LEP (top right); DV genes with higher variance in young LEP (bottom left); and DV genes with higher variance old LEP (bottom right). Fisher’s exact test performed on each of the contingency tables; p-values were adjusted using the Bonferroni method.

Figure 8.. Relative contribution of top aging-associated biomarkers in predictive models of normal breast tissue and PAM50 subtypes in TCGA.

(A) Gene predictors with scaled variable importance ≥25% in prediction of at least one class: normal breast tissue, PAM50 LumA, PAM50 LumB, PAM50 Her2 or PAM50 Basal, in TCGA. Rank ordered heatmap of DE and DV lfc in LEPs with (+) lfc higher in older and (-) lfc higher in younger LEPs (left); scaled variable importance of each gene in each TCGA class (right). Number of gene predictors (k) annotated.

Figure 8—figure supplement 1.. Breast cancer expression levels of top aging-associated predictors of normal breast tissue and PAM50 subtypes in TCGA.

(A–E) Boxplots of log₂ FPKM values in TCGA matched normal and PAM50 subtypes for the top five genes with the highest scaled variable importance contribution in each class in the predictive ML model: (A) normal tissue, (B) PAM50 LumA, (C) PAM50 LumB, (D) PAM50 Her2, and (E) PAM50 Basal. Kruskal-Wallis (KW) test performed across PAM50 breast cancer subtypes as well as matched normal tissue; p-values adjusted across all age-dependent DE and DV genes identified in LEPs found in the TCGA dataset. Post hoc pairwise Wilcoxon test adj. p-value significance levels annotated (*<0.05, **<0.01, ***<0.001, ****<0.0001); p-values adjusted across all pairwise comparisons and across all DE and DV genes identified in LEPs found in the TCGA dataset. Significance level shown only for pairwise comparisons to the reference class. (Note: 3 extreme outlier points in ZNRF3 outside the plotted range in (C) was excluded from visualization only). Number of subjects (n) in each analysis annotated.

Figure 8—figure supplement 2.. Association of top aging-associated predictors of PAM50 subtypes with survival outcomes in TCGA.

(A–C) Multivariate Cox proportional hazard models were run to assess the simultaneous effect of top ML predictors of each subtype along with age at diagnosis and cancer stage (early: Stage I and II vs. late: Stage III and IV) on overall survival (OS, top) and progression-free interval (PFI, bottom) in (A) PAM50 LumA, (B) PAM50 LumB, and (C) PAM50 Basal subtypes. Top ML predictors for each subtype (k_{LumA_top _pred}=47, k_{LumB_top _pred}=37, and k_{LumB_top _pred}=14 genes) were defined as the subset of DE and DV genes with scaled variable importance contribution ≥25% to predictive models. Scaled patient-level gene expression values from primary tumors were used to build the models. Table and forest plots show the number of subjects, the Cox proportional hazard ratios with 95% confidence intervals, and the Wald p-values for each covariate.

Author response image 1.. CDH1 shows no age-dependent change in expression in luminal epithelial cells.

(A) Expression of CDH1 in young vs. old LEP with DE adj. p-values. (B-D) Linear regression of expression values in LEPs between (B) GRHL2 and SGSMS2, (C) GRHL1 and CDH1, (D) SGSM2 and CDH1.

Author response image 2.. Differential expression of SATB1 in breast cancer between PAM50 subtypes is consistent across age groups.

Gene expression levels of SATBI in the TCGA cohort stratified by age group: Young <30y (n=12), MiddleAged >30y <55y (n=481) , and Old >55y (n=712). Kruskal-Wallis (KW) test performed across PAM50 breast cancer subtypes; p-values adjusted across age groups. Post hoc pairwise Wilcoxon test adj. p-value significance levels annotated; p-values adjusted across age groups.

Author response image 3.. Expression of ELF5 and target genes show age specific changes.

Adapted from Miyano et al., 2017. (A) Schematic shows co-culture conditions with young LEP atop of young or old MEP, respectively. Bar graph shows differences in LEP-specific gene expression of LEP markers, and of IGFBP6, a MEP-specific gene, in co-cultured LEPs after 10 days culture on young or old MEP feeders. * p<O.05 and ** p<O.01 annotated. (B) (i) Schematic of experimental outline for extended co-cultures. LEP were separated by FACS after 10 days co-culture either with young MEP (YY) or old MEP (Y/O). LEP from Y/Y and Y/O were further co-cultured with young MEP (Y/Y/Y and Y/O/Y) or older MEP (Y/O/O) for 7 days. (ii) Bar graphs showing normalized expression levels of KRT19, ELF5 and IGFBP6 in cocultured LEP following extended co-culture experiments. (C) Unsupervised hierarchical clustering of <30y LEP in Y/Y (n=3) and Y/O (n=3) co-cultures in parallel with <30y (n=5) and >55y (n=4) 4p LEP and MEP isogenic to the MEP strains used in co-culture (k=20,577 mapped genes). Percent Approximately Unbiased (AU) p-values denoted in red, and percent Bootstrap Probability (BP) in green; AU p≥95% are highlighted in red. (D) Gene-gene correlation matrix of ELF5 and 92 ELF5-target genes found to have absolute correlation 2 0.5 with ELF5 in LEP from <30y and >55y age groups across 9 HMEC strains. (E-F) Hierarchical clustering based on expression levels of ELF5 and the anti-/correlated ELF5-target genes (E) in <30y and >55y 4p pre-stasis LEP (n=9) and (F) in Y/Y (n=3) and Y/O (n=3) co-cultures with <30y (n=5) and >55y (n=4) 4p LEP isogenic to the MEP strains used in co-culture. Percent AU p-values denoted in red, and percent BP in green.

Author response image 4.. HES4 and ZNF503 are differentially expressed in breast cancer between PAM50 subtypes.

Gene expression levels of (A) HES4 and (B) ZNF503 in the TCGA cohort. Kruskal-Wallis (KW) test performed across PAM50 breast cancer subtypes as well as matched normal tissue; pvalues adjusted across all age-dependent DE and DV genes identified in LEPs found in TCGA dataset. Post hoc pairwise Wilcoxon test adj. p-value significance levels annotated; p-values adjusted across all pairwise comparisons and acriss all DE and DV genes.

Author response image 5.. Patient BMI does not affect expression of early and late estrogen response pathway signatures in the studied reduction mammoplasty dataset.

(A) Overweight and obesity statistics from the 2017-2018 NHANES data via the NIDDKD. Percentage of US adults with overweight and obesity by sex (left); and prevalence of obesity among adults age 20 and over by age and sex (right). (B-C) BMI distribution in the GSEI 02088 reduction mammoplasty (RM) dataset across all samples (B) and by age group (C). (D) Early (left) and late (right) estrogen response GSEA signature scores in GSE102088 RM dataset by BMI group.

Author response image 6.. Knockdown of genes that show age-dependent upregulation or increased variability in old MEPs lead to increase expression of ELF5 in co-cultured young LEP.

Genes identified as DE or DV via RNA-seq and proteins identified as DE via mass spec are knocked-down in old MEP, followed by 10 day heterochronus culture with young LEP. ELF5 expression level in co-cultured young LEP are then assessed via qPCR.