The genomic landscape of estrogen receptor α binding sites in mouse mammary gland

Estrogen receptor α (ERα) is the major driving transcription factor in the mammary gland development as well as breast cancer initiation and progression. However, the genomic landscape of ERα binding sites in the normal mouse mammary gland has not been completely elucidated. Here, we mapped genome-wide ERα binding events by chromatin immunoprecipitation followed by high-throughput sequencing (ChIP-seq) in the mouse mammary gland in response to estradiol. We identified 6237 high confidence ERα binding sites in two biological replicates and showed that many of these were located at distal enhancer regions. Furthermore, we discovered 3686 unique genes in the mouse genome that recruit ER in response to estradiol. Interrogation of ER-DNA binding sites in ER-positive luminal epithelial cells showed that the ERE, PAX2, SF1, and AP1 motifs were highly enriched at distal enhancer regions. In addition, comprehensive transcriptome analysis by RNA-seq revealed that 493 genes are differentially regulated by acute treatment with estradiol in the mouse mammary gland in vivo. Through integration of RNA-seq and ERα ChIP-seq data, we uncovered a novel ERα targetome in mouse mammary epithelial cells. Taken together, our study has identified the genomic landscape of ERα binding events in mouse mammary epithelial cells. Furthermore, our study also highlights the cis-regulatory elements and cofactors that are involved in estrogen signaling and may contribute to ductal elongation in the normal mouse mammary gland.


Introduction
Estrogen is a key hormone for mammary epithelial proliferation during the mammary gland development, as well as during breast cancer initiation and progression [1][2][3][4] action of estrogen is mediated through two members of the nuclear receptor family, estrogen receptor α (ERα) and ERβ [5]. Of these two receptors, ERα plays a dominant role in mammary gland development, function, and tumorigenesis [6][7][8][9]. The classical pathway of estrogen action involves estrogen binding to the estrogen receptor, resulting in receptor dimerization and ligand-receptor complex entering into the nucleus where it then binds directly to the genomic DNA containing estrogen response elements (EREs) of the target genes. Various coregulators are then recruited to the gene promoter and/or distal enhancer region to modulate ER-mediated gene transcription [10,11]. In the non-classical pathway, ligand-receptor complexes bind indirectly to genomic regions by tethering to other transcription factors, including the activator protein 1 (AP-1) families of transcription factors (TFs), members of the specificity protein 1 (SP1), and nuclear factor kappa-light-chain-enhancer of activated B cells (NF-κB), which lead to recruitment of chromatin-modifying coregulator proteins and allow activation or repression of ER target genes to drive cell proliferation [12][13][14].
In the mammary gland, ERα is the major receptor operating during ductal morphogenesis, as determined by the analysis of knockout phenotypes in mice [15]. Studies have shown that deletion of ERα inhibits development of the rudimentary ductal structure, and signaling from estradiol through ERα during puberty is required for mammary epithelial proliferation, ductal elongation, bifurcation, and invasion throughout the mammary fat pad [8,16,17]. Further, it has been shown that only a subset of luminal epithelial cells express ERα and estrogen promotes mammary epithelial cell proliferation by inducing amphiregulin (Areg), an epidermal growth factor (EGF) receptor ligand produced in ERα positive cells. Areg acts through a paracrine mechanism on neighboring ER-negative mammary epithelial and stromal cells to mediate estrogen-induced cell proliferation that drives ductal elongation [18,19]. The processes of normal postnatal mammary gland development display many of the properties associated with breast cancer progression such as proliferation, invasion, angiogenesis, and resistance to apoptosis [20].
Several studies have focused on studying genome-wide ER binding sites in ER-positive breast cancer cell lines [21][22][23][24][25][26]. This approach has helped to identify the cis-regulatory elements and cofactors that are involved in mediating ER binding and ER target gene transcription in breast cancer cells. These studies have shown that the vast majority of ER binding sites are located at distal enhancer regions [21][22][23]. Furthermore, Brown and colleagues discovered that ER binding regions were enriched for the FOXA1 binding motif. Subsequently, FOXA1 was then identified as a pioneer factor for ER-chromatin interactions in ER-positive breast cancer cells [27,28]. Genome-wide mapping of FOXA1 binding demonstrated that more than 50% of FOXA1 binding sites overlapped with ER binding events in ER-positive breast cancer cells [29]. In addition to FOXA1, studies have implicated GATA3 and PBX1 as pioneer factors for ER-chromatin interactions, since depletion of these transcription factors leads to reduction of ER binding events and transcriptional activity in response to estrogen [30,31]. Furthermore, it has been shown that GATA3 acts upstream of FOXA1 in the ER binding [31].
Recently, Carroll and colleagues also mapped ER binding events in primary breast cancer samples from patients with different clinical outcomes [32]. The authors found that a subset of ER binding sites is maintained in good outcome, poor outcome, and metastatic breast tumors. Differential binding analysis revealed that ER binding events can discriminate between good and poor outcome tumors. Furthermore, motif analysis revealed that ERE motif is only enriched in good outcome ER-bound genomic regions, indicating that FOXA1 is not involved in these ER binding regions. On the other hand, in poor outcome tumors, ER binding events were associated with ERE and FOXA1 motifs. Nonetheless, it should be noted that genomewide ER binding studies have been limited to ER-positive breast cancer cell lines and primary breast cancer samples, not the normal mammary gland. Although one recent study has shown ER binding sites in the developing mouse mammary gland, this study focused only on control mammary glands which have not received any hormonal treatment [33].
Given the potential involvement of ERα in mammary gland development, as well as in breast cancer initiation and progression, we investigated for the first time, genome-wide ERα binding events by chromatin immunoprecipitation followed by high-throughput sequencing in mouse mammary gland under in vivo conditions of acute treatment with estradiol. Furthermore, we also used genome-wide transcriptome profiling (RNA-seq) to uncover the global transcriptome response to estradiol. This combined system approaches allowed us to identify the ER targetome in the mammary gland and this ER targetome consists of a subset of ER-regulated target genes whose acute regulation by estradiol is associated with direct binding to ER in luminal epithelial cells. Overall, our study provides a unique resource for the mechanisms underlying estrogen regulated gene expression in the mouse mammary gland.

Mice
Animal experiments were approved by the Institutional Animal Care and Use Committee (IACUC) at Baylor College of Medicine (Houston, TX). BALB/cJ mice were purchased from The Jackson Laboratory (Bar Harbor, ME). All mice were housed in the animal facility at Baylor College of Medicine and maintained in a conventional mouse facility with room temperature set at 22˚C with food and water provided ad libitum. The animal facility is accredited by the American Association of Laboratory Animal Care.
Hormone treatments. At 6 weeks of age, mice were ovariectomized and rested for 10 days and then mice were injected subcutaneously with sesame oil (50μl) or 17β-estradiol (100ng) for 2 hours. Sesame oil and 17β-estradiol were purchased from Sigma. After 2 hours, mice were sacrificed and both inguinal mammary glands (#4 mammary glands) were harvested from each mouse. The lymph node was removed from the # 4 mammary glands and used for RNA-seq and ChIP-seq analysis.
Total RNA isolation. For RNA-seq and quantitative Real-Time PCR (qPCR), total mammary gland RNA was isolated from #4 glands using the RNeasy Lipid Tissue Midi Kit according to the manufacturer's instructions (QIAGEN, Inc., Valencia, CA). To avoid potential contamination from muscle, only inguinal (rather than thoracic) glands were used in this study. Physical integrity of the RNA was assessed using the Agilent 2100 Bioanalyzer (Agilent Technologies, Inc., Santa Clara, CA). RNA quantitation was performed using the Nanodrop ND1000 spectrophotometer (Nanodrop Technologies, Wilmington, DE). For qPCR and RNA-seq analysis, total RNA pooled from a set of four to five mice per treatment group was analyzed. To ensure statistical significance, three separate sets of mice (12-15 mice in total) per treatment group were used in qPCR experiment.
RNA-seq library preparation and sequencing. Total RNA samples with RNA integrity number (RIN) �8 were used for transcriptome sequencing. Total RNA (10ng) was used for amplified double-stranded cDNA using the Ovation RNA-Seq System (NuGEN, San Carlos, CA). Double-stranded DNA was sheared to 200-300bp using the Covaris S2 sonicator (Covaris, Woburn, MA) and ligated to Illumina paired-end adaptors using the Illumina TruSeq DNA library preparation kit according to the manufacturer's instructions (Illumina, San Diego, CA). PCR amplification was performed to obtain the final cDNA library. Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA) analysis was used to verify fragment size after amplification, library size and concentration before clustering. A total of 10pM of the library was then used for paired-end sequencing on the HiSeq 2500 at the Genomic and RNA Profiling Core in Baylor College of Medicine.
Binary Analysis Results (BAR) files. Peak locations were determined using the Model-based Analysis of ChIP-seq Algorithm (MACS) with a p-value cutoff of 1E-8.
Validation of ERα -binding sites by ChIP-quantitative PCR. ERα ChIP was performed as described above in the ChIP-seq analysis. ChIP DNA and input DNA were analyzed by qPCR using SYBR Green Master Mix. Primers were generated corresponding to the regions identified by our ChIP-seq. The enrichment of ERα-binding found in each sample was normalized to input values. For primer sequences, see S1 Table.
Accession numbers. RNA-seq and ERα ChIP-seq data were deposited in the Gene Expression Omnibus (GEO, Accession number GSE130032).
Statistical analysis. Data are presented as mean±SEM. The significance of the differences between groups was determined by Student's t-test. Values were considered statistically significant at p < 0.05.

Genome-wide analysis of ERα binding sites
We mapped ERα binding events using chromatin immunoprecipitation followed by deep sequencing on mammary glands collected from ovariectomized mice that had been treated with estradiol for 2 hours. To ensure a robust and precise representation of the ERα binding sites in the mammary gland, we employed two independent biological replicates for ERα ChIP-seq analysis. Binding events were called using the MACS model based peak finding algorithm. A p-value cutoff 1E-8 identified 9607 binding regions in replicate 1 and 6819 binding regions in replicate 2. The total number of reads and uniquely mapped reads for each sample is shown in S2 Table. By using the web-based application Galaxy/Cistrome, intersection of the two biological replicates showed 78% of binding site concordance. This analysis yielded 6237 common ER binding events in mammary epithelial cells mapped near 3686 unique genes after acute treatment with estradiol ( Fig 1A). A heatmap was generated for visualization of the 6237 common ERα peaks found in both replicates ( Fig 1B). We used Pearson's correlation method to compare two biological replicates. As shown in Fig 1C, Pearson's correlation was very high (r = 0.95) between the two biological replicates suggesting that good reproducibility between the two replicate samples. Comparison of the sequences from ER intervals within various placental mammalian genomes showed a strong level of conservation in the regions of ERα binding only and not in surrounding regions ( Fig 1D). Furthermore, we also performed crosscorrelation analysis in ERα ChIP replicate 1 and 2. We observed two cross-correlation peaks one corresponding to the read length (Phantom peak) and other one to the average fragment length of library suggesting that our ChIP-seq is robust (S1 Fig). We next analyzed ERα binding site distribution in this data set. ERα binding sites were divided into three categories; distal upstream (-100kb to -1kb), proximal (-1kb to +400 bp), and distal downstream (+400 bp to +100kb), based on the peak locations relative to RefSeq genes. In this analysis, binding sites were assigned to the nearest transcription start site (TSS). As shown in Fig 1E, only 5% of ERα binding sites were distributed in the proximal region. However, majority of ERα binding sites were located in the distal upstream (40%) and downstream (55%) regions. These results indicate that the vast majority of ERα binding sites are located at the distal regions of ER target genes in the mouse mammary gland which is consistent with ER binding sites in ER-positive breast cancer cells [21][22][23]. Further, it has been shown that distal ERα binding sites are anchored at gene promoters through long-range chromatin interactions, suggesting that ERα functions by bringing genes together for coordinated transcriptional regulation by chromatin looping [41]

Motif analysis
To discover the network of transcription factors linked with genomic binding of ERα, we employed the Cistrome tool SeqPos and identified enriched binding motifs within the proximal and distal regulatory regions of ERα binding sites. In the distal regions (up and down) of ERα binding sites contained a canonical ERE and other motifs including PAX2, ESRRB, SF1, and AP1 motifs as the most highly enriched (p-value <10 −30 ) cis-elements. Additionally, there was substantial enrichment of FOXA1 motif in the distal ERα binding regions which is consistent with previously reported data that FOXA1 acts as an important pioneer factor for ER chromatin interaction in luminal breast cancer cells [27][28][29]. The proximal region (-1kb to +-400bp) of ERα binding sites contained canonical ERE as the most significantly enriched motif, along with other significantly enriched motifs such as ESRRB, SF1 and TAL (Fig 1F). The complete motifs enrichment is included in S3 Table. Furthermore, we also performed motif enrichment analysis by using HOMER in the overlapped ER binding sites (regardless of binding site distribution). As expected, most of the motifs identified by Cistrome were confirmed by HOMER. The complete HOMER motifs enrichment is also included in S4 Table.

Validation of ERα binding by ChIP qPCR
We next interrogated several known estrogen-regulated genes for their ability to recruit ERα in an estrogen-dependent manner. To test this, mammary gland chromatin was isolated from ovariectomized mice that had been treated with vehicle or estradiol for 2 hours and then subjected to ERα ChIP followed by qPCR as described in the Materials and Methods. Our ChIPseq data have identified several ERα binding sites in the distal upstream region of the Greb1 gene (Fig 2A). Of these, we verified two of the ERα binding sites (-3281 and -33380) by ChIP qPCR. As shown in Fig 2A, ERα was highly enriched in these regions in response to estradiol treatment, confirming direct estradiol-dependent ERα recruitment to Greb1. Our ERα ChIPseq analysis revealed peak of ERα binding of progesterone receptor (Pgr) in the distal downstream region (+70663) which was independently verified by ChIP qPCR (Fig 2B). Furthermore, our ChIP-seq data also discovered ERα binding in the distal downstream regions of Fos, Gata3, and Areg genes (Fig 2C-2E) and these binding sites were validated by ChIP qPCR in an estradiol-dependent manner (Fig 2C-2E). Collectively, our ERα ChIP-seq was independently confirmed by ChIP qPCR.

Genome-wide transcriptome analysis of ER-target genes
To disclose the global gene signatures that are acutely regulated by estradiol, we next performed genome-wide transcriptome analysis in the mammary gland of ovariectomized mice treated with vehicle or estradiol for 2 hours. We used biological replicate pools of mammary glands from 5 mice per pool for total RNA isolation and then subjected to global gene expression analysis by RNA-seq as described in the Material and Methods. The reads were mapped to the mouse Ensembl reference set of 82,508 transcripts (annotated to 31034 genes) using TopHat. Relative transcript abundance was identified using the RSEM algorithm [42]. The edgeR package was used to identify differentially regulated gene expression [43]. The total number of reads and uniquely mapped reads for each sample is shown in S2 Table. Data quality was assessed using principal component analysis (PCA) and clustering of RNA-seq samples using Euclidean distance. Our PCA plot of the four RNA-seq samples showed that the control and estradiol samples are well separated along the principal component 1 (PC1) which explained the 75% of the total variance (S2A Fig). Further, we also used the Euclidean sample distances to show that the samples from the same conditions are well clustered (S2B Fig). We identified 493 genes differentially regulated due to an acute treatment with estradiol. We observed that 220 genes (45%) were upregulated and 273 (55%) genes were downregulated at FDR <0.05 (S5 Table). A heatmap was also generated to show the differentially regulated genes (Fig 3A). We used Pearson's correlation method to compare two biological replicates. As shown in Fig 3B, Pearson's correlation was very high (r = 0.96) between the two biological replicates for both the vehicle and estradiol treated groups suggesting that our RNA-seq is robust. Next, we applied a bioinformatics approach, Ingenuity Pathways Analysis (IPA) (Ingenuity, CA), to identify a potential network based on the regulated genes to unveil the molecular mechanisms of estradiol action in the mammary gland. Using this approach, we identified βestradiol (p = 9.46E-21) and TGFβ (p = 5.75E-19) as the top upstream regulators (data not shown). As expected, these signaling cascades were predicted to be activated by acute treatment with estradiol treatment. The most enriched gene networks were centered on FOS, a classical estrogen target gene (Fig 3C). Studies have shown that FOS has been implicated to regulate cell proliferation, differentiation, and transformation [44,45]. Furthermore, DAVID functional annotation analysis was also used to identify the biological function of differentially regulated genes. Gene signatures-induced by estradiol were enriched for GO terms associated with negative regulation of protein kinase activity, epithelium development, regulation of cell proliferation and morphogenesis of a branching structure. However, genes repressed by estradiol were enriched for GO terms associated with muscle contraction, mammary gland development, epithelial cell differentiation, gland morphogenesis and epithelium development ( Fig  3D). The complete GO terms analysis is included in S6 Table.

Validation of RNA-seq by qPCR
To validate our RNA-seq data, we measured the expression of well-established estradiol responsive genes by qPCR. To examine this, 6-week old female mice were ovariectomized and then rested for 10 days. After 10 days, mice were treated with vehicle or estradiol for 2 hours

Fig 2. Representative screen shots of ChIP-seq data showing gene ERα recruitment at 2 hours after exposure to estradiol IGV screen shots showing ERα binding sites in relation to the TSS of (A) Greb1, (B) Pgr, (C) Fos (D) Gata3
and (E) Areg. Peak locations relative to the TSS are listed below each screen shot and the numbers indicate peak value of each gene. Red boxes represent peaks that were validated by ChIP-qPCR. Graphs representing validation of ERα occupancy using ChIP followed by qPCR for Greb1, Pgr, Fos, Gata3 and Areg. A total of six mice per replicate; Results are means ± SEM of three independent experimental replicates. �� , p < 0.01; ��� , p < 0.001.
https://doi.org/10.1371/journal.pone.0220311.g002 and #4 mammary glands were harvested. We used triplicate pools of mammary glands from vehicle or estradiol treated mice (4-5 mice per pool) for RNA isolation and subjected to qPCR. First, we tested estradiol-induced genes. Studies have shown that Greb1, Pgr, and Fos are canonical ER target genes [27][28][29]. Our RNA-seq analysis revealed that these genes were induced more than 2-fold by acute treatment with estradiol. Similar to what we observed in the RNA-seq data, acute estradiol treatment caused a significant increase in mRNA levels of Greb1, Pgr, and Fos expression in the mammary gland (Fig 4A). Our ChIP-seq analysis also identified ERα binding sites near these genes and confirmed direct hormone-dependent ER recruitment to these genes (Fig 2A-2C). Interestingly, our RNA-seq also identified several known genes (Krt4 and Csf1) as well as previously unidentified estrogen target genes such as Bmp8a, Heyl, Errfi1, Six1, and Gata6. These estrogen-induced genes were independently confirmed by qPCR (Fig 4A). Our RNA-seq analysis also revealed that acute in vivo treatment with estradiol significantly downregulated Gata3, Foxa1, Areg, and Ccnd1 expression. These downregulated genes were independently confirmed by qPCR. As shown in Fig 4B, acute estradiol treatment caused more than a 2-fold reduction of Gata3, Foxa1, Areg, and Ccnd1 mRNA levels when compared with vehicle treatment. The possible explanation for this observation is that a secondary transcriptional response may be involved in induction of these genes in response to estradiol signaling. Our acute treatment (2 hours) makes it unlikely that we will detect secondary changes in transcription. Further, it has been shown that estradiol treatment (24h) did not increase the Ccnd1 and Areg mRNA levels. However, in the presence of progesterone treatment, these genes were dramatically increased suggesting that progesterone may also be required for these gene expression in the normal mouse mammary gland [46]. Indeed, Nuclear Receptor Signaling Atlas (NURSA) transcriptomine database revealed that Foxa1 and Gata3 are downregulated in response to estradiol. Studies have shown that GATA3 is necessary for luminal epithelial cell differentiation and the gene is often mutated in human breast cancer [47][48][49][50][51][52][53][54]. Several direct downstream targets of GATA3 in the luminal epithelium have been identified including FOXA1, an important regulator of ERα expression. Notably, our ERα ChIP-seq identified a peak near the Gata3 gene and confirmed direct estrogen-dependent ERα recruitment to Gata3. In addition, we also observed significant downregulation of several luminal cell makers such as Cdh1, Muc1, Krt7, Krt8 and KRT18 gene sets in acute treatment with estradiol treatment suggesting that maturation of luminal cells may be reduced (Fig 4B). Collectively, our data suggests that acute treatment with estradiol represses luminal cell differentiation leading to expansion of a de-differentiated epithelial cell population.

Identification of ERα targetome in mammary epithelial cells
To identify subsets of estrogen-regulated genes that are direct targets of ER positive luminal epithelial cells, we integrated ERα ChIP-seq and RNA-seq data to uncover the subset of estrogen regulated genes that directly recruit ERα in ER positive luminal epithelial cells. In this analysis, we used the 220 estradiol-induced and 273 repressed genes identified after acute treatment with estradiol (Fig 3) to interrogate the 3686 unique ER binding genes. This integration showed that 36% (177 genes) of the 493 differentially regulated genes (220 induced and 273 repressed genes) directly recruit ERα (Fig 5A and 5D). These genes (177) are shown in S7 Table. Since ERα is expressed only in a subset of luminal epithelial cells, these genes could potentially represent a global estradiol-regulated luminal epithelial cell targetome of ERα. We next analyzed ERα binding site distribution in the ERα targetome in mammary epithelial cells. Interrogation of the ERα binding site distribution in upregulated genes showed that only 7% ERα binding sites were distributed in the proximal region, defined as -1 kb to +400 bp. However, the vast majority of ERα binding sites were located in the distal upstream (defined as -1kb to -100kb) and downstream regions (defined as +400 bp to +100kb), 34% and 59%, respectively (Fig 5B). We performed the same analysis on genes repressed by estradiol treatment. As shown in Fig 5E, only 4% of genes are associated with ERα binding in the proximal region, 38% in the distal upstream, and 58% in the distal downstream. Taken together, our ERα targetome revealed that the majority of ERα binding is located in distal regions of genes in the mammary gland, which is consistent with previous reports of ER binding sites in ERpositive breast cancer cells [21][22][23].
To identify the network of transcription factors linked with the ERα targetome, we employed the Cistrome tool SeqPos and identified enriched binding motifs within the estradiol-induced and repressed genes associated with ER binding regions. ER target genes induced by estradiol contained a canonical ERE motif as the most highly enriched cis-element. Additionally, ESRRB, ATF2, de novo, PAX2, and PAX6 motifs were enriched in estradiol-induced genes ( Fig 5C). However, ERα binding sites in genes repressed by estradiol showed a significant enrichment of ERE, ESRRA, and SF1 motifs (Fig 5F). The co-enrichment of PAX2 sites with ERα binding is consistent with a recognized role for PAX2 in human breast cancer cells [55]. PAX2 was shown to be expressed in a subset of breast cancers and was recruited to the ER binding site after both estrogen and tamoxifen treatment [55].
To identify the biological function of the ERα targetome, we again used DAVID functional annotation analysis. Interrogation of the 135 estrogen-induced ERα targetome identified a significant enrichment for GO terms associated with regulation of small GTPase mediated signal transduction, regulation of specific transcription from RNA polymerase II promoter, epithelial cell differentiation, growth and morphogenesis of a branching structure. On the other hand, analysis of the 42 genes comprising the estrogen-repressed ER targetome showed enrichment for GO terms associated with muscle contraction, response to inorganic substance, and response to inorganic substance (Fig 5G). The complete GO terms analysis is included in S8 Table.

Discussion
In this study, we report for the first time a comprehensive analysis of ERα binding events which elucidates the molecular mechanism of ERα action in the normal mouse mammary gland in vivo in response to acute estradiol treatment. We have identified 6237 high-confidence ERα binding sites in two biological replicates. These binding regions also corresponded with 3686 unique genes in the mouse genome that recruit ERα in response to acute estradiol treatment. Recent studies have shown that ER interacts with DNA sites in the absence of ligand in breast cancer cells and mouse uterus [56,57]. However, these ER chromatin binding sites are relatively lower than ligand-induced ERα binding sites. It has also been shown that ER binds the same chromatin binding sites regardless of the ligand; however, signal intensity was highest for estradiol treatment [29]. The distribution of ERα binding sites revealed that the majority of these sites (95%) are located in the distal enhancer, which is in agreement with other studies in ER-positive breast cancer cell lines [21,22]. Our data also suggest that ERα interactions occur at distal sites under normal physiological conditions. Interrogation of ERα-DNA interaction sites revealed ERE, PAX2, SF1, and AP1 motifs were highly enriched at sites where ERα binds, suggesting that these transcription factors interact with DNA and contribute to stabilizing the ER complex on chromatin. Furthermore, our ER binding sites also identified the FOXA1 motif, whose expression is required for luminal epithelial ERα expression and post-pubertal development of the normal mammary gland [58]. FOXA1 is an established pioneer factor for the majority of ERα binding to the genome in human breast cancer cells [21,29]. Whole transcriptome analysis by RNA-seq revealed that 493 genes are differentially regulated by acute treatment with estradiol in the mouse mammary gland in vivo. By integrating genome-wide ERα binding and global gene expression signatures, we uncovered a novel ER targetome of 177 genes that are associated with estrogen-dependent ductal elongation in the mouse mammary gland in vivo. This cistromic and transcriptome resource revealed the cisregulatory elements and cofactors that are involved in estrogen signaling in the normal mouse mammary gland.
It is well documented that the ovarian hormones estrogen and progesterone act as master regulators of mammary gland development. Specifically, estrogen triggers ductal elongation during puberty whereas progesterone is involved in mammary gland side-branching [46,[59][60][61]. Studies have shown that estradiol-mediated activation of ERα signaling is required for ductal mammary epithelial cells to proliferate, leading to ductal outgrowth. Deletion of ERα results in normal mammary glands before puberty. However, after puberty, terminal end buds remained absent, and the failure of ducts to invade into the fat pad, suggesting that ERα action is essential for mammary gland development [16,62]. Only subset of luminal epithelial cells are ERα positive, and release of paracrine signals from these ER-expressing cells permit other nearby epithelial cells, both luminal and myoepithelial, to participate directly in ductal outgrowth [16,18].
In our whole genome transcriptome profiling using RNA-seq, we identified that ER-regulated genes were affected by acute treatment with estradiol in the mouse mammary gland in vivo. We discovered that 493 genes are differentially regulated by estradiol. According to DAVID functional annotation analysis, these differentially regulated genes are primarily involved in molecular and cellular processes ranging from epithelium development to morphogenesis of a branching structure. Furthermore, we validated the RNA-seq results independently by qPCR analysis to confirm that these genes were indeed regulated in the mammary gland in response to acute treatment with estradiol. For validation, we selected canonical ER target genes as well as previously unidentified ER target genes. Studies have shown that Greb1, Pgr, and Fos are classical estradiol target genes and are strongly induced by estradiol. ERα directly controls GREB1 expression, and GREB1 is required for breast cancer cell growth. Clinically, the loss or reduced expression of GREB1 is predictive of poor outcome and decreased relapse free survival, similar to ER status [63][64][65][66][67]. It is well documented that PR is an ER target gene and induces proliferation of neighboring cells primarily by a paracrine mechanism [68]. It has been shown that progesterone induces receptor activator for nuclear factor κB ligand (RANKL) secretion from ERα+/PR+ mammary cells, and then this factor can directly binds to its receptor RANK to induce side-branching and alveolar development during ductal development and pregnancy [59,69].
Fos is an immediate early gene acutely induced by intracellular signaling cascades [70]. Fos, along with other members of the Fos family, dimerize with Jun to form the AP-1 transcription factor complex which regulates many genes involved in proliferation, differentiation, and survival. In breast cancer cells the Fos gene also plays a key role in tumorigenesis and invasive least one ER binding site. (B) ERα binding distribution at genes induced by estradiol at 2 h. (C) SeqPos Motif enrichment in direct ER target genes was induced at 2 h after estradiol exposure. (D) Proportional Venn diagram representing unique estradiol-repressed genes at 2 h as detected by RNA-seq (orange), unique genes with at least one ERα -binding region as detected by ChIP-seq (purple), and estradiol-repressed genes at 2 h with at least one ERα binding site (overlap). (E) ERα binding distribution at genes was repressed by estradiol treatment. (F) SeqPos Motif enrichment in the direct ER target genes repressed at 2 h after estradiol exposure. (G) Enriched GO terms of ER targetome by using the DAVID functional annotation tool.
https://doi.org/10.1371/journal.pone.0220311.g005 growth. Our data showed that estradiol induces expression of Greb1, Pgr, and Fos genes in the normal mouse mammary gland. Furthermore, our ERα ChIP-seq data confirmed that ER binds these genes at the distal region.
BMP8a is a member of transforming growth factor β (TGFβ) and plays a central role in cell proliferation, differentiation, and survival [71]. A recent study showed that BMP8a acts as paracrine factor for estrogen-dependent regulation of epithelial cell proliferation in the uterus [72]. In the present study, we observed acute exposure to estradiol induced Bmp8a mRNA levels in the mammary gland, suggesting that Bmp8a may also act as a paracrine effector of estrogeninduced mammary epithelial cell proliferation. Furthermore, our ERα ChIP-seq data also identified ERα binding near Bmp8a gene, which was independently confirmed by ChIP qPCR (data not shown), suggesting that Bmp8a is also a direct target of ERα in the normal mammary gland.
GATA3 has been implicated in breast tumorigenesis and its highest expression was observed in the luminal subtype of breast cancer [73,74]. Mutations in GATA3 have also been identified in a subset of breast cancers [75]. In the normal mouse mammary gland, Gata3 is the most highly expressed transcription factor in the mammary epithelium [76]. Using a mammary epithelium-specific knockout of Gata3, it has been shown that Gata3 is essential in maintaining luminal epithelial cell differentiation [48,49]. Several direct downstream targets of GATA3 in the luminal epithelium have been identified including FOXA1, a key regulator of ER binding in the breast cancer cells. Furthermore, ER binds near the GATA3 gene in breast cancer cells, as shown through ChIP assays [21,23,77].
Our ERα ChIP-seq data also revealed that ER binds near the Gata3 gene in a distal enhancer, also confirmed by ChIP qPCR, suggesting that Gata3 is a target of ER in the normal mammary gland. Furthermore, we showed that Gata3 is a direct target for ERα through down regulation after acute estradiol treatment. Our results suggest that acute treatment with estradiol represses Gata3 gene expression leading to expansion of a de-differentiated epithelial cell population in the mammary duct. Our study also identified that estradiol-induced down regulation of Gata3 was associated with repression of its target genes such as Foxa1 and Ccnd1. It has been also shown that GATA3 is required for CCND1expression [78]. In addition, our data suggest that Foxa1 is involved with ERα binding in the normal mouse mammary gland, since motif analysis of ERα binding regions in the normal mammary gland showed Foxa1 motif enrichment. Studies have shown that FOXA1 acts as a pioneer factor for ER binding in the distal enhancer sites in breast cancer cells [26,29]. Our findings also support the notion that FOXA1 is not only involved in ER binding of breast cancer cells but also involved in ER binding in the normal mouse mammary gland. Studies have shown that the majority of ERα binding sites are located in the distal enhancer sites in breast cancer cells, similar to our results from the normal mammary gland. Indeed, it has also been shown that distal enhancer sites are known to interact with transcriptional complexes located near gene through DNA looping [41,79]. Thus it is possible that ER binds distal enhancer regions and induces its target gene transcription through chromatin looping.
In conclusion, our study offers identification of ERα binding events across the genome under normal physiological conditions in the developing mouse mammary gland and provides a useful dataset to allow further study of estrogen signaling through its receptor. Furthermore, our findings identify cis-regulatory factors that cooperate in ERα-mediated control of gene expression in the normal mammary gland.