Age-related gene expression in luminal epithelial cells is driven by a microenvironment made from myoepithelial cells

The relationship between promoter DNA methylation and transcription of key lineage-specific genes in vivo is preserved in primary culture

In order to tractably examine the impacts of age and microenvironment on lineage specificity, we developed LEP- and MEP-marker probe sets that enabled exploration of the relationships between promoter methylation and gene expression. These were used in functional cell-based experiments that require small numbers of cells and allow many replicates. To maintain consistency, all experiments used fourth passage (4p) pre-stasis finite lifespan HMEC from discarded reduction mammoplasty tissue for cell function studies [6]. The tissues were obtained from women, who do not have breast cancer, differing in chronological age at the time of surgery. The so called in vitro “replicative ages” are the same at 4p (Table S1).

Transcriptome-wide differential expression analysis (Illumina HumanHT-12 v4 BeadChips Set 1, n=24,965 gene probes, m=19,499 mapped genes) of FACS-enriched CD10+/CD227- MEP and CD10-/CD227+ LEP from 4 different HMEC strains were used to identify lineage specific genes. Genes selected for use as probes showed >3-fold differential expression (DE) between LEP and MEP from women <30y (Benjamini-Hochberg, BH, adj. p-val < 0.05) (Fig. 1A) while also having CpG islands in their 5’ region within 5kb of the ATG start codon. These genes did not exhibit culture adaptive expression changes. Probe sets were designed to facilitate qPCR analyses of gene expression and promoter methylation of DNA by McrBC methylation sensitive enzyme digestion. DKK3, COL7A1, IGFBP6 and TMP2 were selected as MEP marker genes, and KRT19, ELF5, RBM47 and COBL were selected as LEP marker genes. These genes showed lineage-specific expression (Fig. 1B) that was inversely correlated with promoter DNA methylation (Fig. 1C) in LEP and MEP from eight different women <30y, suggesting that transcription of these genes was at least partly regulated by DNA methylation. The LEP and MEP markers PROM1 and TP63, respectively, also were used as well-known mammary epithelial lineage markers that do not appear to be controlled by methylation in CpG islands within 5kb of the start codons. The probe sets showed excellent correlation between 4p HMEC and FACS-enriched LEP and MEP from uncultured human mammary epithelial organoids, both in terms of lineage-specific gene expression (Fig. 1D, r²=0.96) and DNA methylation (Fig. 1E, r²=0.93). Results from the probe-based assays were comparable to the methylation and expression data from uncultured mammary epithelia in the NIH Roadmap Epigenomics Project (Fig. S1) [19]. Thus these lineage-specific probe sets were validated both with in vivo samples and in publicly available data.

Myoepithelial cells maintain phenotypic fidelity of luminal cells

In vivo, LEP organize into an inner layer that borders a lumen, and is surrounded by MEP, which produce basement membrane components (Fig. 2A). To study LEP-MEP interactions, we established a biologically relevant LEP/MEP co-culture system that mimicked a LEP-microenvironment in mammary glands. LEP from <30y HMEC were cultured on top of completely confluent MEP feeder layers. When mammary epithelial cells are cultured they produce ECM that adsorb to the plastic and establish an apical (the side facing up) and a basal surface (side attached to ECM) [20]. We used this property of MEP to advantage and established a co-culture method that recapitulates a bilayered epithelia as shown in Figure 2A. Maintenance of CD227⁺/CD10^- LEP cellular phenotype was determined by FACS measurement of these surface protein markers over 10 days in culture (Fig. 2B). In LEP/MEP co-cultures, the proportion of LEP, defined as the fraction of cells maintaining CD227⁺ expression levels, increased steadily by day 10, showing a 15% increase compared to day 0, whereas the proportion of LEP declined steadily when cultured directly on tissue culture plastic (TCP) (Fig. 2B and 2C). The proportion of CD227⁺ LEP decreased at a rate similar to culture on TCP when co-cultured on isogenic breast stromal fibroblasts (Fig. 2C). Furthermore, LEP genes were down-regulated while IGFBP6, a MEP gene, increased in LEP on TCP and fibroblast co-cultures, whereas LEP genes remained high and MEP genes were suppressed in LEP grown on MEP feeders (Fig. 2E). Changes in expression of number of LEP- and MEP-specific genes in LEP that were cultured on TCP were anti-correlated with methylation changes (Fig. 2D). The other MEP genes also increased expression in LEP cultured in TCP and fibroblast microenvironments, but IGFBP6 proved to be particularly responsive, and thus most useful in subsequent experiments. To determine whether the ability to maintain LEP was dependent on isogenicity, co-cultures were established with LEP and MEP from three individuals aged 16-28y in the nine possible combinations. In all cases, maintenance of LEP on MEP co-cultures was superior to LEP on plastic, and was independent of the individual (Fig. 2F).

To understand what type of physical interactions between LEP and MEP were required, LEP were cultured in a number of defined conditions. LEP from a 19y woman were cultured on (I) plastic, (II) co-cultured with isogenic MEP, (III) plastic with MEP-conditioned medium, (IV) mono-cultured with LEP attached to Boyden chamber inserts or (V) co-cultured with isogenic MEP in which LEP were physically separated by a Boyden chamber insert (Fig. 2G). CD227 and CD10 expression was monitored every two days by FACS (Fig. 2H). Co-cultures II and V, regardless of direct or no-direct contact of LEP with MEP, showed increased proportions of CD227+ cells over time, whereas LEP declined in mono-culture conditions (I, III and IV) (Fig. 2H). Although the CD227+ cells population increased in both co-culture conditions II and V, the no-direct contact co-culture (V) showed reduced expression of PROM1 and ELF5 and increased IGFBP6 expression compared to direct cell-cell contact (II)(Fig. 2I), suggesting that only a partial LEP phenotype was maintained in co-cultures where MEP and LEP were physically separated (in this case by 900 μm). These experiments revealed that TCP does not provide an optimal microenvironment for maintaining the LEP phenotype even in presence of MEP-conditioned media. Indeed, LEP on TCP are eventually reprogrammed to more MEP-like states, and senesce more quickly than MEP (not shown), which serves as an extreme example of how microenvironment can negatively impact the LEP phenotype. These data show that direct LEP-MEP contact provided crucial microenvironment components for maintenance of the LEP phenotype. That stromal fibroblasts – a cell type that is ostensibly at least one cell-diameter removed from LEP in vivo, did not provide a satisfactory microenvironment for maintenance of LEP underscores the specificity of the LEP-MEP interactions for maintaining the LEP phenotype.

Loss of luminal-fidelity with age coincides with disruption of the relationship between promoter methylation and gene expression

To determine the consequences of age on fidelity of lineage-specific gene expression, we interrogated the corresponding expression (Illumina HumanHT-12 v4 BeadChip Set 1) and DNA methylation (Infinium Human Methylation 450K BeadChip) patterns of the selected lineage-specific probeset genes in LEP and MEP cells from 9 different HMEC strains representing <30y and >55y age groups. Differential expression (DA)analysis focused on our probeset genes (n=19 Illumina probes) shows that there is a shift towards smaller fold-differences between LEP and MEP genes in the >55y group compared to the <30y group (Fig. 3A). That is, fidelity with which lineage-specific genes are expressed decreased with age, LEP and MEP looked more alike.

Differential methylation (DM) analysis of M-values (the log2 ratio of methylated to unmethylated probes) at CpG sites associated with the probeset genes (n=247 Illumina CpG sites) revealed a gross inverse relationship between expression and DNA methylation; LEP-specific genes had higher methylation levels at these loci in MEP, and vice versa (Fig. 3B). The distribution of fold changes in methylation decreases in magnitude of lineage-specific differences in the >55y group compared to the <30y group (Fig. 3B). The effect of age is most evident in the LEP-specific probeset genes. Viewing the methylation architecture of these probeset genes at the resolution of annotated gene regions, we find age-dependent changes in methylation mainly occur in the genes’ regulatory regions, e.g. transcription start sites (TSS), 5’-untranslated region (UTR), and some in the 1^st exon (Fig. S2A).

Consistent with a loss of lineage fidelity, but not specificity, qPCR results showed KRT19 and ELF5 significantly decreased expression (Fig. S2C), while MEP-specific genes modestly increased in LEP from the >55y women (Fig. S2B). Age-dependent gene expression patterns were similarly associated with inverse changes in promoter DNA methylation when evaluating McrBC digested gDNA with qPCR (Fig. S2B and S2C). The gene expression and DNA methylation patterns in MEP were stable during aging, with the exception of RBM47 and DKK3, which show modest increased and decreased in expression with age, respectively. When the expression/methylation relationships of all the lineage-specific genes were viewed together using the qPCR probeset, the relationship between age, lineage, and lineage-specific gene regulation became most apparent. Both younger MEP and LEP, as well as older MEP, show tightly regulated lineage-specific gene expression that were inversely correlated with DNA methylation (Fig. 3C and 3D). However, in older LEP, the LEP-specific pattern of gene expression and methylation were significantly disrupted (Fig. 3D).

We next asked whether artificial disruption of the lineage-specific gene expression/methylation patterns in normal HMEC could recapitulate the biochemical loss of lineage fidelity in phenotypes that arise with age. After HMEC from a 19y woman were treated for 48h with 5-aza-2’dC (5aza), the inhibitor of DNA methyltranferases, the expression of otherwise suppressed lineage-specific genes were observed in both MEP (Fig. 3E) and LEP (Fig. 3F). HMEC from women <30y typically exhibit distinctive CD10+/CD227- MEP and CD10-/CD227+ LEP populations by FACS analysis, and the divisions between the two populations are typically become less obvious in HMEC from older women (representative examples of <30y and >55y 4p HMEC in Fig. 3G). Treatment of the 19y HMEC with 5aza phenocopied a lineage distribution more characteristic of HMEC from a much older person (Fig. 3H). Although the effects of using 5aza are non-specific, the approach showed that disruption of the DNA methylation-gene expression relationship during aging may partly explain loss of lineage fidelity. Taken together we speculate that changes in methylation of LEP-specific genes is fundamental to loss of fidelity with age in this lineage.

Age-related lineage fidelity in luminal cells is a switchable state that is determined by the age of myoepithelial cells

Because contact with MEP is important for the maintenance of LEP lineage specificity, we hypothesized that MEP from women of different age groups could impose age-dependent gene expression patterns in LEP. LEP from women <30y were co-cultured on isogenic MEP ((Y)oung/Y) or MEP from women >55y (Y/(O)lder). LEP- and MEP-marker gene expression were measured in LEP after 10 days of co-culture. LEP cultured on MEP from younger women (Y/Y) maintained robust LEP gene expression and suppressed MEP genes. However, in LEP cultured on MEP from older women (Y/O), expression of PROM1 and ELF5 decreased and IGFBP6 increased (Fig. 4A). In non-isogenic co-culture experiments, LEP from 3 young HMEC strains, cultured on MEP from 3 younger or 3 older women consistently showed that the age of MEP determined ELF5 and IGFBP6 levels in LEP (Fig. 5A and 5B).

We next determined whether LEP from older women could be rejuvenated in a younger microenvironment. On isogenic MEP (O/O), older LEP exhibited high IGFBP6 expression and low levels of ELF5 consistent with older primary LEP (Fig. 4B). In contrast, when cultured on younger MEP (O/Y), ELF5 expression was increased 7-fold and IGFBP6 decreased to background levels compared to LEP in O/O co-culture (Fig. 4B). The age-dependent MEP-imposed gene expression in LEP was consistent with age-dependent expression in FACS-enriched LEP from both 4p HMEC strains (Fig. 5C) and directly from breast tissue (Fig. 5D).

We next examined transciptome-wide expression of <30y LEP cells in either Y/Y (n=3) or Y/O (n=3) co-cultures in parallel with 4p LEP and MEP from both <30y (n=5) and >55y (n=4) strains isogenic to the MEP strains used in co-culture (Illumina HumanHT-12 v4 BeadChips Set 2). Unsupervised hierarchical clustering based on transciptome-wide gene expression levels (n=26,599 gene probes, m=20,577 mapped genes) show that the co-cultured LEP expression profile is determined by the age-group of the MEP used for co-culture. Clustering of <30y LEP in Y/Y co-culture with <30y 4p LEP is highly supported (Approximately Unbiased, AU, p-val > 0.95), and separate from >55y 4p LEP and the corresponding isogenic <30y LEP in Y/O co-culture (Figure 4C).

Because down-regulation of KRT19 in LEP with age is observed by immunostaining in human breast tissue and in 4p cultured cells, and is one of the most prominent age-related phenotypes that we have identified so far [6], it was unexpected that reduction of KRT19 mRNA was not detected in mixed-age co-culture (Fig. 4A). We hypothesized that longer time scales might be necessary to observe broader phenotypic changes. Consequently, we examined the effect of protracted exposure to the aged microenvironments on the expression of KRT19, ELF5 and IGFBP6 in LEP (Fig. 4D). Young LEP were re-isolated after 10-days co-culture with MEP, then cultured on fresh (Y or O) MEP feeders for a further 7 days. Under prolonged Y/O/O exposures, ELF5 and eventually KRT19 were down-regulated, and IGFBP6 was up-regulated (Fig. 4E). The pattern of ELF5 and IGFBP6 expression in the Y/O/Y conditions followed the age of the MEP donor and KRT19 remained high (Fig. 4E). Thus, at least a subset of age-dependent gene expression patterns in LEP are malleable, able to be directed by the age of the MEP in the microenvironment. The extent of this age-dependent reprogramming of the LEP, i.e. the number of lineage-specific genes and their degree change in expression, was further dependent on the duration of exposure to microenvironments.

Microenvironments impose age-specific promoter DNA methylation states at ELF5 and IGFBP6 loci

We next determined if exposure of LEP to the different aged MEP was also associated with changes in DNA methylation in the promoter regions. ELF5 and IGFBP6 were the focus of these experiments because age-dependent changes in expression were reliably linked to changes in promoter methylation. DNA methylation was examined in young LEP co-cultured with young (Y/Y) or old (Y/O) MEP for 10 days. CDX1 and BCLAF1 were used as hyper- and hypo-methylation control genes, respectively, which did not change in any of our co-culture conditions (Fig. 6A). DNA methylation in LEP in the regulatory regions of ELF5 increased 30±9% in Y/O cultures and methylation at the IGFBP6 promoter decreased 13±0.7% compared to Y/Y (Fig. 6A). These changes in methylation measured with McrBC digests in co-culture were consistent with age-dependent changes in ELF5- and IGFBP6- DNA methylation at the gene regulatory regions (Infinium Human Methylation 450K BeadChip) of <30y vs. >55y FACS enriched 4p LEP (Figure 6B-6C). In 4p LEP, differences in corresponding methylation beta-values, which measure percentage of methylation, show large age-specific increase in methylation levels in ELF5 across numerous CpG sites in the gene regulatory region from hypo-methylated (β-val < 0.25) in <30y LEP to hemi- (0.25 < β-val < 0.75) and hyper-methylated (β-val > 0.75) in >55y LEP. Furthermore, these changes in ELF5 methylation levels in >55y LEP shift towards more MEP-like methylation levels (compare Fig. 6B, Fig. S3A), showing loss of lineage-specific methylation with age. A smaller, though significant, decrease in methylation levels in IGFPB6 regulatory region is observed with hyper-methylated sites in <30y LEP becoming hemi-methylated in >55y LEP (compare Fig. 6C, Fig. S3B). These differentially methylated regulatory regions in ELF5 and IGFBP6 correspond to either annotated enhancer element regions or DNaseI hypersensitivity sites, suggesting that transition between hypo-, hemi- and hyper-metylation states in these regions could have a significant effect on transcription levels (Fig. SA and SB). Differential methylation between <30y and >55y 4p LEP for the remaining probeset genes are shown in Fig. S4A-S4F.

Methylation of the promoters should facilitate stabilization of microenvironment-imposed gene expression. Furthermore, when the regulated gene product is a lineage-specific transcription factor, we should be able to measure changes in downstream gene networks. As a transcription factor, ELF5 warranted further examination and we asked whether age-related decline of ELF5 alters expression of target genes. ELF5 target genes were previously identified by ELF5-ChIP in the breast cell line T47D (499 probes, 429 unique genes at FDR < 0.05, t = 48HR) [21]. Of this ELF5-ChIP gene set, 323 of the prospective ELF5-target genes were in our Illumina HumanHT-12 v4 BeadChips (Set 1) corresponding to 528 gene probes. Gene-gene correlation analysis showed 92 target genes corresponding to 105 gene probes to be correlated or anti-correlated with ELF5 expression in 4p LEP from 9 different HMEC strains from both <30y and >55y age groups (absolute correlation ≥ 0.5, Figure 7A). Pre-stasis 4p LEP clustered according to age based on ELF5 and ELF5-target gene expression (Figure 7B). Furthermore, young LEP co-cultured with young MEP (Y/Y) clustered with <30y while two out of three young LEP co-cultured with old MEP (Y/O) clustered with a majority of >55y LEP (Figure 7C). This suggests that age-related changes in ELF5 expression impact known target genes downstream of ELF5 in LEP in an age-dependent manner, and that these changes can be driven by cues from the MEP microenvironment. Taken together, these data converge on the conclusion that the age of MEP-generated microenvironments is a determinant of promoter DNA methylation states, which may drive, and stabilize some of the age-related gene expression phenotypes in LEP.

Cell culture

Primary HMEC strains were generated and maintained as described previously [45,46]. HMECs were grown in M87A medium containing cholera toxin and oxytocin at 0.5 ng/ml and 0.1nM, respectively. HMEC strains used in this study were listed in Table S1. In co-culture study, FACS-enriched MEP from 4^th passage HMEC strains or fibroblast cells were re-plated to 6-well plates and cultured until the cells were confluent. The cells were treated with mitomycin C (Santa Cruz Biotechnology) at 10µg/ml for 2.5h. Then, FACS-enriched LEP were seeded directly on the mitomycin C-treated cells or on a cell insert (BD) for separation from MEP cells. LEP in co-culutre were separated by FACS after 7-10 days for gene expression and DNA methylation analyses. 5-aza-2’-deoxycytidine (5aza) was added at 15 µM in culture medium. Cells were cultured for 48hr with culture medium containing 5aza and harvested for further analysis.

Flow cytometry

Cells dissociated from organoids or primary HMEC strains were stained with anti-human CD227-FITC (BD Biosciences, clone HMPV) and anti-human CD10-PE (BioLegnend, clone HI10a) by following standard flow cytometry protocol. Co-cultured LEP was stained with anti-human CD227-FITC and anti-human CD10-PE and separated by FACS. Cells were sorted by FACSVantage SE or analyzed by FACSCalibur (BD Biosciences). Data were analyzed with FlowJo software.

Quantitative gene expression analysis

Total RNA was isolated from FACS-enriched cells using Quick-RNA MicroPrep (ZYMO Research). cDNA synthesis was performed by iScript (BioRad) according to manufacturer’s manual. Gene expression was measured by LightCycler 480 (Roche) with iTaq universal SYBR Green supermix (BioRad). Data were normalized to RPS18 or GAPDH by relative standard curve method or normalized by Vandesompele method [47]. ANOVA was used to test for statistical significance. Primers are listed in Table S2.

DNA methylation analysis

Genomic DNA was purified using Quick-gDNA MicroPrep (ZYMO Research). Genomic DNA was digested with McrBC (New England BioLabs) and EcoRI (New England BioLabs), or EcoRI only as a control. DNA methylation was measured by real-time PCR using LightCycler 480. Amount of DNA methylation was normalized by internal primer control that targeted the DNA not containing CG dinucleotide. DNA methylation by McrBC method shows % of cells with methylated DNA.

Statistical analysis

We examined the gene expression correlation between 4^th passaged HMECs and uncultured breast tissue cells using Pearson correlation. Figures 1D and E show the correlation of each gene's expression by colors. Yellow indicates that the trend of one gene expression across other genes is correlated to both HMECs and uncultured cells. Blue color indicates there are no correlations.

Gene expression and DNA methylation analysis using Illumina beadchip

Total RNAs and genomic DNAs from LEP and MEP were isolated using Quick-RNA MicroPrep and Quick-DNA MicroPrep (ZYMO Rseach), respectively, after enrichment by FACS. Sample preparations for Illumina HumanHT-12 v4 Expression BeadChip array and HumanMethylation450 BeadChip array were performed in UCLA Neuroscience Genomics Core (UNGC).

Gene expression and DNA methylation of ten LEP and ten MEP samples from nine individuals (four young subjects <30 years old, five older subjects > 55 years old) were analyzed across two Illumina HumanHT-12 v4 Expression (Set 1) and two Infinium 450K Methylation BeadChip arrays. BeadChips were designed to have balanced, well-stratified distribution of samples for factors of interest (cell type and subject age) with one reference individual (240L MEP and LEP samples, <30 years old) shared between the two chips. Additionally, gene expression from six LEP co-culture samples (three biological replicates of 240L <30y LEP co-cultured on 240L <30y MEP, and three biological replicates 240L <30y LEP co-cultured on 122L >55y MEP) along with 18 samples of 4p LEP and MEP isogenic to the co-culture MEPs (five biological replicates of 240L <30y LEP and MEP, and four biological replicates 122L >55y LEP and MEP) distributed evenly across five Illumina HumanHT-12 v4 Expression (Set 2) were also analyzed.

Raw gene expression data from Illumina HumanHT-12 v4 BeadChips were pre-processed with Bioconductor limma package neqc function which performs normexpr background correction using negative control probes, log2 transformation and quantile normalization between arrays [48,49]. Normalized data set were pre-filtered to remove gene probes with values less than negative control probes. This was done by calculating detection p-values using limma package detectionPValues function [48,49], and only gene probes with detection p-values < 0.05 for at least 3 samples were retained. Potential batch effects between chips were checked using Principal Component Analysis (PCA).

DNA methylation measured using Infinium 450K Methylation BeadChips were analyzed using a custom R script. The arrays were read and normalized using the lumi package which return methylation M-values, the calculated log2 ratio of the intensities of methylated probe versus unmethylated probe [50]. Color bias adjustment was performed using lumiMethylC and simple scaling normalization (SSN) was performed using lumiMethylN functions [51]. Methylation beta-values, which is an approximation of the percentage of methylation of a given CpG site, and is calculated as the ratio of the methylated probe intensity and the overall intensity (sum of methylated and unmethylated probe intensities), were computed from M-values using m2beta function [50].

Differential expression (DE) analysis was performed on the pre-filtered normalized data set as a one-step method for handling batch effects is preferred for DE analysis. DE analysis was performed using Bioconductor limma package which applies an empirical Bayes approach proposed to provide more stable inference when number of arrays is small [49,52]. Differential methylation analysis was performed on methylation M-values as M-values are found to be more statistically valid for differential analysis [50]. For both DE and DM analysis, a model was fitted with a coefficient for each of the four factor combinations (MEP/LEP cell type vs. young/old subject age group). Batch was added to the design matrix along with the factor combinations. Sample replicates, as well as paired nature of MEP/LEP samples were accounted for by blocking for individual in the duplicate Correlation function [53]. Array weights were also computed and included in the model. The four contrast terms and an interaction term were considered for analysis. Empirical Bayes moderation of computed statistics was applied. Benjamini-Hochberg (BH)-adjusted p-values (Benjamini and Hochberg’s method to control for false discovery rate) and log2 fold change statistics were calculated for each probe for both the lineage-specific (MEP vs. LEP) or age-specific (<30y vs. >55y) contrast terms.

For Exploratory Data Analysis (EDA) and visualization of gene expression data, normalized data were corrected for batch effects using ComBat [54] – an empirical Bayes approach for adjusting data for batch effects that is robust to outliers in small sample sizes as applied in the Bioconductor sva package [55], with no covariates included in the model for batch adjustment [56]. After batch adjustment, expression data were pre-filtered as described above using the detection p-values. PCA and hierarchical clustering were used pre- and post-ComBat treatment for visualization to illustrate removal of batch effects and the clustering of reference sample. All analyses using hierarchical clustering of samples were done using complete linkage of Euclidean distances in the pvclust package which calculates Approximately Unbiased (AU) p-values and Boostrap Probability (BP) values for each cluster [57].

For Exploratory Data Analysis (EDA) and visualization of DNA methylation data, methylation beta-values – which have a more direct biological interpretation, were plotted across CpG sites ordered by their chromosomal mapping. Beta values range from 0 to 1 and denote hypo- (β-val < 0.25), hemi- (0.25 < β-val < 0.75) and hyper-methylated (β-val > 0.75) methylation levels. Differential methylation was annotated based on DM analysis of methylation M-values.