Estrogen receptor α (ERα)-binding super-enhancers drive key mediators that control uterine estrogen responses in mice

Estrogen receptor α (ERα) modulates gene expression by interacting with chromatin regions that are frequently distal from the promoters of estrogen-regulated genes. Active chromatin-enriched “super-enhancer” (SE) regions, mainly observed in in vitro culture systems, often control production of key cell type–determining transcription factors. Here, we defined super-enhancers that bind to ERα in vivo within hormone-responsive uterine tissue in mice. We found that SEs are already formed prior to estrogen exposure at the onset of puberty. The genes at SEs encoded critical developmental factors, including retinoic acid receptor α (RARA) and homeobox D (HOXD). Using high-throughput chromosome conformation capture (Hi-C) along with DNA sequence analysis, we demonstrate that most SEs are located at a chromatin loop end and that most uterine genes in loop ends associated with these SEs are regulated by estrogen. Although the SEs were formed before puberty, SE-associated genes acquired optimal ERα-dependent expression after reproductive maturity, indicating that pubertal processes that occur after SE assembly and ERα binding are needed for gene responses. Genes associated with these SEs affected key estrogen-mediated uterine functions, including transforming growth factor β (TGFβ) and LIF interleukin-6 family cytokine (LIF) signaling pathways. To the best of our knowledge, this is the first identification of SE interactions that underlie hormonal regulation of genes in uterine tissue and optimal development of estrogen responses in this tissue.

Estrogen hormones are intricately involved in key molecular events that underlie development of female reproductive tract tissues, optimizing responses to ovarian hormones essential for successful reproduction later in life (1). Many studies of estrogen response mechanisms have utilized in vitro cell culture models, which have advanced understanding but are limited in their application to more biological contexts. In vivo models of hormone response include the rodent female reproductive tract, which is highly hormone-responsive, facilitating study of mechanisms of tissue development and hormone response. For example, female reproductive tissues of estrogen receptor ␣ (ER␣)-null and estrogen-deficient Cyp19 (aromatase) KO mouse models, both of which resemble their WT littermates at birth, lack pubertal uterine development, and thus exhibit a hypoplastic uterine phenotype (1)(2)(3)(4). These experimental characteristics duplicate those found clinically in patients with insensitivity syndromes due to mutations in their ESR1 or CYP19 genes (5)(6)(7)(8). Estrogen's activity involves interaction with its nuclear receptor, ER␣, which localizes to enhancer regions of chromatin, driven by high affinity for estrogen response element (ERE) DNA motifs. Enhancer regions that ER␣ interacts with are often distal from promoters of estrogenregulated genes (1,9), and current models incorporate a mechanism of chromatin looping to facilitate contacts between enhancer and promoter regions (10,11). Earlier studies have noted that some enhancers have particularly high enrichment of indicators of active chromatin; these are classified as "superenhancers" (12,13). Frequently, these super-enhancers control production of key transcription factors that characterize a particular cell type (12)(13)(14). In mice, prior to production of ovarian hormones, which begins with the onset of estrus cyclicity at puberty (postnatal day 28), the prepubertal uterus is hypoplastic, has ER␣ in all cells, and exhibits rapid transcriptional and growth responses of all uterine cell types to exogenous estradiol (E2) (15)(16)(17). The response, however, differs from what is observed in an ovariectomized adult mouse uterus, where the uterus response is restricted to only the epithelial cells (15,17). Through the process of pubertal maturation, the reproductive tract acquires optimal reproductive functionality (3,15,16,18,19). The molecular details of this important process have not been explored; therefore, we sought to define the estrogen-dependent enhancer landscape of the developing mouse uterus by analyzing uterine chromatin isolated from prepubertal (21day-old) or from adult (10-week-old) mice that were ovariectomized (ovexed) to remove the source of endogenous E2 and then treated for 1 h with vehicle (V) or with E2 to gain insights into molecular components and responses impacted during development of a hormonally regulated tissue.

Defining uterine super-enhancers that bind ER␣
We defined super-enhancers (SEs) that bind ER␣-in prepubertal uterine samples from mice treated for 1 h with vehicle (V) or with E2, as well as in uterine samples from ovexed adult mice treated for 1 h with V or E2 as described under "Experimental procedures." Briefly, ER␣-binding enhancers were identified using ER␣ ChIP-Seq data by stitching together regions with more than one ER␣ peak within 12,500 bp in at least one of the samples. The 4634 enhancer stretches containing multiple ER␣-binding sites were then classified as "typical" or "super" enhancers (TEs or SEs, respectively) by plotting ranked H3K27Ac signal of each; enhancers beyond the slope ϭ 1 elbow of the curve (Fig. 1A) were considered SEs. Therefore, locations of ER␣ ChIP-Seq peaks from all four sample types were either classified as single peaks (Fig. 1B) or were within enhancer regions with multiple peaks (Fig. 1B). The multipeak enhancer regions were then classified as super-enhancers or typical enhancers (Fig. 1B). Several hundred SEs were identified in each sample type (Fig. 1A). Comparison of the enhancers called in each sample indicates that many SEs overlap across the sample types ( Fig. 1A and Fig. 2), and most that are SEs, but not in all samples, are ranked in the top 25% of enhancers (Fig. 1A, red lines in the heatmap below each curve). The SEs in the prepubertal V-treated mice are present before the cells have been exposed to ovarian E2 beginning at puberty, whereas the SEs in the ovexed adult samples represent super-enhancers after the uterus has undergone pubertal maturation. For each super-enhancer curve, genes located at the top 10 SEs are listed. We noted that three of the genes at the SEs (Zmiz1, Ncor2, and Rara) were among the 10 highest ranked in all samples (Fig.  1C). Comparison of the SEs using heatmaps of ER␣ and H3K27Ac ChIP-Seq signals centered on each of the 2431 ER␣binding sites that is within a SE in at least one of the four samples ( Fig. 2) indicates that the H3K27Ac signal intensities of the adult and the 21-day-old are similar in V-and E2-treated samples, respectively, showing that the SEs seen in prepubertal samples, prior to pubertal uterine development, are largely unchanged by pubertal development. E2 leads to increased ER␣ binding at each site, as well as increased H3K27Ac signal flanking SE ER␣-binding locations. Interestingly, although E2 also increases ER␣ binding at the 21,894 non-SE ER␣-binding sites (TE and single ER␣-binding sites), H3K27Ac signal is more static (Fig. 2), with only the highest ranked ER␣ binding sites exhibiting E2-dependent H3K27Ac signal increase. Some researchers use H3K4Me1 ChIP-Seq signal intensity, also associated with active enhancers, to classify SEs. Using H3K4Me1 signal to classify uterine enhancers revealed that 65-87% of SEs by H3K4Me1 were also classified as SEs using H3K27Ac (Table  S1). Additionally, Ͼ98% of SEs by H3K4Me1 signal were within the top quartile of the H3K27Ac signal. Locations of SEs called using H3K27Ac in the adult ovexed 1-h E2 uterus relative to genes and chromosomes are illustrated in Fig. S1.

Chromatin interactions in uterine tissue
Distal enhancers can exert transcriptional control by physically contacting gene promoters through "looping" mechanisms (10). It is not known whether contacts are dynamically formed in order to exert regulatory signals between distal enhancers and promoters, whether loops are preformed and the activity of the promoters is mediated via dynamic transcription factor interactions, or whether a combination of mechanisms occurs (20). Assigning ER␣-interacting regions to regulated genes is important for understanding how ER␣ binding to distal enhancers regulates gene promoters; therefore, we analyzed ovariectomized adult mouse uterus 1 h after treatment with V, E2, or progesterone (P4) using Hi-C. We used the Juicer tool to call loops from the Hi-C data of each treatment as well as differential loops in pairwise comparisons (Fig. 3A). We observed little variability between replicates or between different treatments, suggesting that loops are not altered globally by the hormone treatments we administered (Fig. 3A). Therefore, we combined the loops from all WT samples into an "atlas" of uterine loops. Then we used cohesin subunit SMC1a ChIP-Seq to assess interactions and potential impacts of hormone treatment. Because cohesin slides along chromatin until it contacts a pair of converging CTCF sites (10,21), we compared SMC1a peaks at loop anchors selected to have cohesin at both ends as a surrogate for the presence of a loop in the sample. We examined SMC1A signal intensity, both at selected loop anchors and at other locations, and compared V-and E2-treated samples (Fig. 3B). When we compared V versus E2 SMC1a signal intensity, we did not observe a notable impact of E2 treatment (Fig.  3B) at selected loop anchors or at other locations, suggesting once again that loops are preformed and not impacted by our hormone treatments. To assess whether specific transcription factors might be associated with V versus E2 loops, we analyzed the DNA motifs associated with SMC1a peaks that were located at both ends of loops. The CTCF motif is highly enriched in SMC1a peaks at loop anchors of both V and E2 samples ( Table  1). The ER␣ motif estrogen response element and ER␣ "tethering" factor motif SP1 are both selectively enriched at SMC1a peaks of the E2 sample. These observations suggest E2-induced mechanisms in which ER␣ and associated transcription factors such as SP1 interact with preformed chromatin loops. Consistent with this mechanism, we note that Hi-C analysis of uterine tissue from mice lacking ER␣ shows that loops are present and do not differ from WT samples any more than the variation observed between replicates (Fig. 3A).

Super-enhancers are predominantly located in loop ends and are associated with uterine genes
Next, we evaluated the relationships between SEs and chromatin loops and whether SE and uterine genes are brought into proximity in 3D space. For this part of our analysis, we focused on the adult ovexed E2-treated samples. First, we examined average H3K27Ac signal per loop at all loop ends or "anchors" (Fig. 4A), which indicated that this histone modification was centered on the ends of loops. By further comparing loop anchors overlapping TEs versus SEs, we observed significantly more H3K27Ac signal at TE-associated loop anchors than at all loop anchors (p Ͻ 0.0001; Fig. 4A). Signal was further enriched at SE-associated loop anchors relative to all or to TE-associated loop anchors (p Ͻ 0.0001), reflecting the more robust enhancer activity and confirming these regions as SEs. Of the 281 SEs identified in adult E2-treated tissue, 94% (263 SEs) were at the end of a loop, suggesting that they were likely to contact genes in distal regions (Fig. 4B). To get an indication of the gene regulation that might occur by looping between genes and SEs, we selected genes that either were at one of the SEs in a loop end or were at the other end of a loop that formed from a SE. Altogether, we found this entailed 1600 genes (Fig. 4C). Our previous studies have indicated that ER␣ interaction with chromatin occurs within 1 h and that resulting gene regulation occurs

Super-enhancers drive uterine estrogen responses
subsequently (9,22). Therefore, we used RNA-Seq of uterine RNA from ovariectomized mice treated with V or treated with E2 for 2 or 6 h to determine whether these 1600 ER␣-binding SE-associated genes were transcribed and E2-regulated in the uterus. Indeed, 963 of these 1600 SE-associated genes are expressed (transcripts per million (TPM) Ն 1 in at least one treatment condition) in uterine cells (Fig. 4C). Overall, 9975 uterine genes are regulated by E2 (TPM Ն 1, adjusted p value Ͻ0.01) in at least one treatment condition (Table S2). 569 of the 963 uterine SE-associated genes are regulated by estrogen ( Fig.  4D and Table S3). This suggests that ER␣-binding SEs may be involved in regulation of 5-6% of E2-responsive uterine genes. The remaining E2-regulated uterine genes are likely to be mediated by ER␣ interactions at TE or at single ER␣-binding sites, as we have previously demonstrated that the ER␣ is required for uterine E2 genomic response (4,22). To assess how the 569 SE-associated estrogen-regulated uterine genes impact biological functions, we used INGENUITY PATHWAY ANALYSIS. Genes involved with the promotion of gene transcription, cell survival, and genitourinary tract development and suppression of mortality and apoptosis were apparent (Table 2). Notably, upstream regulators with known roles in uterine function, such as estrogen and cytokines, including LIF, a critical mediator of embryo implantation (23), as well as IGF1 and TGF␤ (24 -32), are identified as activating upstream regulators (Table 3).
Therefore, these ER␣-binding SE-associated genes are enriched for important uterine processes.

Super-enhancer-associated genes acquire estrogen responsiveness during postpubertal development
To assess whether postpubertal development impacted E2 regulation of the SE-associated genes identified above, we analyzed gene regulation in 21-day-old and ovexed adult samples from WT mice and mice with deletion of the ER␣ gene (ER␣KO). We treated the mice for 2 h with V or E2 and isolated uterine RNA for microarray. We observed that 245 SE-associated genes were differentially expressed after 2-h treatment with E2. Hierarchical clustering reveals that the responses are ER␣-dependent ( Fig. S2 and Table S4). What is surprising is that, despite the presence of the SEs prior to pubertal development (Figs. 1 and 2), estrogen regulation is much more robust in adult ovexed samples ( Fig. S2 and Table S4), indicating that these genes acquire the transcriptional response as a result of pubertal development.

Super-enhancers are associated with key uterine factors
Our analysis has indicated an association between ER␣-mediated responses and looping between SEs and estrogen-regulated genes. Because genes at the SE often encode factors critical for specific cell types (14), we evaluated examples of four

. Estrogen increases ER␣ and H3K27Ac enrichment at super-enhancers. Shown is a heatmap of ER␣ ChIP-Seq and H3K27Ac
ChIP-Seq signal centered on 2431 ER␣ peaks that are at super-enhancers (top) or 21,894 ER␣ peaks that are not at super-enhancers but are at typical enhancers or are single ER␣-binding sites (bottom) in uterine chromatin samples from adult ovexed or 21-day-old female mice treated with saline V or with E2 for 1 h. For each set, enhancers are in the order of ER␣ signal in the adult E2 sample. Each panel shows ChIP-Seq signal at Ϯ1 kb relative to the ER␣ peak midpoint, with signal depth-normalized to 20 million uniquely mapped nonduplicate reads/sample.

Super-enhancers drive uterine estrogen responses
genes at the SE that are potentially important for uterine development and E2 transcriptional response and that are among the top-ranked SEs (Fig. 1A). Retinoic acid receptor a (Rara) is a ligand-dependent transcription factor in the nuclear receptor superfamily (33). Homeobox (Hox) transcription factors are key mediators of anterior-to-posterior developmental patterning of tissues (34). Ncor2 (SMRT; nuclear receptor corepressor 2) interacts with and activates histone deacetylases, thereby reducing transcriptional activity (35). Zinc finger MIZ domaincontaining protein 1 (Zmiz1) is a novel transcription factor. We first used RT-PCR of RNA samples from prepubertal (21-dayold) or adult ovariectomized uteri to examine the impact of pubertal development on expression of these four genes. All genes showed increases in E2 response in adult tissues versus prepubertal tissues (Fig. 5A).
Hi-C analysis indicates distal chromatin loops with SEs at Rara (Fig. 5Ba), the Hoxd cluster (Fig. S3), Ncor2 (Fig. S4A), and Zmiz1 (Fig. S4B) genes. We wondered whether the expression of genes that either directly coincide with the SEs or are within loops would be impacted by E2. Therefore, we looked at expression and estrogen regulation of genes within loops formed at these four SEs. We used our RNA-Seq data to evaluate the levels and E2 regulation of these genes (Table 4). We assessed the relative level of expression of each gene by computing the average TPM of nine uterine RNA samples (3 replicates each of uterine RNA from animals treated with V or with E2 for 2 or 6 h). We examined the E2 regulation of each gene based on fold change of E2 2-h RNA versus V RNA or E2 6-h RNA versus V RNA. The Rara SE forms loops that include several genes (Table 4 and Fig. 5Ba), but only Rara and Igfbp4 were regulated by E2 (TPM Ն 1, FDR p value Ͻ0.01; Table 4). The Hoxd cluster SE coincides with the coding genes for Hoxd8, -9, -10, and -11 and contacts two additional SEs (Fig. S3A). Some of the genes within loops that include the SE are regulated by E2 (Table 4). There is a loop between the Hoxd SE and a third SE nearly 30 megabases away, in the Wilms tumor 1 (Wt1) gene (Fig. S3B). Thus, the Hoxd cluster in the uterus appears to be at the core of a very large 3D chromatin structure of SEs. Two SEs are localized at the Ncor2 coding sequence (Fig. S4A). Ncor2 is one of the most highly expressed genes in the loops that include the SE within it (Fig. S4A and Table 4). The Zmiz1 SE loops with two distal SEs (Fig. S4B). Zmiz1 and several other genes within the loop are expressed (Fig. S4B and Table 4). loops called from each sample in parentheses and differential loops called between the indicated pairs. Minimal differences are seen between replicates (V-1, V-2, V-3; P4-1, P4 -2) as well as between V and hormone treatment (E2 or P4) or between mice lacking ER␣ (ER␣KO) and WT litter mates (V-2). B, an "atlas" of uterine loops was built from the six combined WT sample Hi-C data sets. This collapsed set of loops was then filtered to retain only those with SMC1a ChIP-Seq peaks at both loop ends using either V-or E2-treated SMC1a ChIP-Seq data. The heatmap shows SMC1a signal centered on V (left) or E2 (right) SMC1a peaks that are (top) or are not (bottom) at loop ends of the selected loops. Each panel shows ChIP-Seq signal at Ϯ1 kb relative to the SMC1a peak midpoint, with signal depth-normalized to 20 million uniquely mapped nonduplicate reads/sample

Super-enhancers drive uterine estrogen responses
To experimentally determine whether a SE ER␣-binding site can impact expression of genes within Hi-C defined chromatin loops, we examined expression levels of genes within loops from a previously described SE distal from the Igf1 gene (36). The Igf1 SE forms a loop that includes Tyms-ps, Dram1, Pmch, Parpbp, Nup37, and Ccdc53 within it (Fig. 5Bb). Pmch and Tym-ps are not expressed (average TPM Ͻ 1); Dram1, Igf1, and Nup37 are regulated by E2 (Table 4). In our previous study, we disrupted the Igf1 distal SE by deleting one of its ER␣-binding sites (36); using this IGF1enh4KO model, we examined whether the distal SE is important for E2 regulation of genes in the loops that emanate from it. RT-PCR analysis of uterine RNA samples from IGF1enh4KO and WT littermates shows that although E2 regulation of Igf1 itself was affected by disrupting the SE (Fig. 5C), no change in expression of Dram1, Nup37, or Parpbp occurred, indicating that the deleted site selectively impacts only Igf1 expression.

Mouse uterus super-enhancers differ from those in breast cancer cells
To evaluate how the mouse uterine SEs might be related to SEs in other E2-responsive cells, we compared our findings with a similar analysis that defined SEs that bind ER␣ in the MCF-7 breast cancer cell line grown in hormone-deprived medium (37). We compared the genes that are closest to 227 MCF-7 SEs with a list of genes closest to 281 uterine tissue SEs (Fig. 6). After removing genes that do not have both mouse and human orthologs (Fig. 6A), about 10% of the SE genes were common between the mouse uterus and MFC-7 cell systems (Fig. 6B). The small number of SEs located at genes shared in MCF7 cells and uterine tissue is likely due to differences in the biological responses E2 elicits in a normal E2-responsive tissue (mouse uterus) versus an immortalized breast cancer-derived cell line.

Discussion
The importance of E2 and ER␣ in uterine maturation is emphasized by the observation of hypoplastic uterine tissue in adult-aged ER␣-null mice (4), which circulate high levels of E2 but have no receptor protein, and in mice lacking Cyp19, which are therefore unable to synthesize E2 (2, 38) but have normal levels of ER␣. We defined the estrogen-dependent enhancer landscape of the developing mouse uterus and observed that overall, most of the 21-day-old and adult SEs had similar characteristics; in general, both ER␣ and H3K27Ac signals are increased by E2 at the SE (Fig. 2). Although E2 increases ER␣ signal of TE and single ER␣-binding sites as well, H3K27Ac is less impacted by E2 (Fig. 2), with only the highest-ranked ER␣binding sites showing an increase. Overall, this suggests that the enhancer landscape necessary for uterine maturation is already in place by postnatal day 21. It was clear, however, from our gene regulation analyses, that pubertal development leads to more robust regulation of genes associated with these SEs (Fig. S2).
Comprehensive analysis of ER␣ interaction with chromatin using ChIP-Seq has revealed that ER␣-binding regions are often distal from coding genes, which has led to a need to understand the 3D structure of chromatin within ER␣-expressing cells. We noted that all ER␣-binding loop ends had significantly more H3K27Ac signal, centered on the loop anchor, than loop anchors in general (Fig. 4A), especially those at superenhancers. This highlights the importance of ER␣ in the mechanism of activation via these 3D structures. We did note that H3K27Ac signal is greater at loop ends with TE (Fig. 4A), indicating a significant role of these ER␣-binding regions in E2 response as well. How responses mediated by single ER␣-binding sites or by TE versus SE impact uterine development and biological processes will be an important focus of future study. Our previous work characterized one such gene mechanism by deleting an ER␣-binding site in a super-enhancer distal from the Igf1 gene, resulting in loss of E2 induction of the coding transcript (36). Deleting one of the five distal ER␣-binding sites did not impact E2-dependent ER␣ binding to the four remaining distal sites or to a site near the promoter; however, E2-dependent recruitment of the histone acetyltransferase p300 and the cohesin subunit SMC1a was decreased at neighboring sites, and E2-dependent induction of enhancer RNA and coding mRNA transcription was lost (36). Thus, E2-dependent recruitment of factors necessary to optimize enhancer-promotermediated Igf1 transcription is governed by a critical site in the distal enhancer. The Igf1 distal SE also forms chromatin loops that include other genes, some of which are also regulated by

analysis showing de novo motifs in SMC1a peaks that overlap with anchors of loops with SMC1a at both ends (Fig. 3B) in V-treated (top) or E2-treated (bottom) samples
Shown is the average TPM (average of V, E2 2-h, and E2 6-h transcripts per million) of transcription factors that bind to motifs. ‫,ء‬ ratio of the percentage of target with motif to the percentage of background with motif.

Super-enhancers drive uterine estrogen responses
E2; however, deletion of the ER␣-binding site did not change the E2 response of these other genes, indicating that the primary target of the SE in this case is the Igf1 gene. Whether this relationship between SEs and regulated distal genes is generally seen, or whether there are cases in which expression of multiple genes is impacted by a single SE remains to be studied.
In this study, we identified SEs either in or looping to other well-known highly E2-responsive genes, including Greb1, Lif, Fos, Cebpb, Cyr61, and Inhbb (Table S3). The Greb1-interacting super-enhancer identified here was the focus of a previous study indicating the potential role of distal ER␣-binding regions in Greb1 regulation (39). The same ER␣-binding region was utilized to model mechanisms of distal interaction in a recent study (40), highlighting a role for the steroid receptor coactivator, SRC-3, in optimizing contacts between enhancer and promoter regions via interaction with an intronic SRC-3-binding sequence (40). Prior to E2 treatment, the enhancer and promoter are held in proximity to each other through their contacts with SRC-3 bound to the intron. E2 treatment leads to SRC-3-dependent formation of the enhancer-promoter transcriptional complex (40) and, consequently, increased Greb1 transcription. The formation of structures, such as the SRC-3/ intron/enhancer/promoter complex, may be a mechanism utilized by enhancer-promoter interactions.
We note that most of the SEs identified in adult ovexed E2-treated samples were at a loop end (263 of 281; Fig. 4B) and

Super-enhancers drive uterine estrogen responses
were likely to be at or form a loop to genes expressed in the uterus (Fig. 4C). Many of these SE-associated genes are involved in transcription and in pathways important for uterine function, such as p53 (41) and Tgf␤ (26) signaling; therefore, these SEs are associated with important mediators of uterine maturation and function. We did not observe a SE at the Esr1 gene, although ER␣ is detected in all uterine cells and SE-associated genes are often characteristic of a particular cell type. Perhaps we should have expected that ER␣ is not driven by a SE, as we know that the ER␣-null mouse develops all uterine cell types normally (4,42). We did see that Esr1 was found to be one of the genes closest to an ER␣-binding SE in both the MCF-7 cell and our uterus data sets (Fig. 6C). When we examined the impact of E2 on expression of the 963 SE-associated uterine genes, 569 were regulated by E2, a small portion of the 9975 uterine genes that are regulated by E2 (Fig. 4D and Tables S2 and S3). These SE-associated differentially expressed genes (DEG) are involved with activation of TGF␤ and estrogen responses and with cytokine pathways, including LIF, a cytokine that is critical for embryo implantation (23) ( Table 3).
Our finding of a SE overlapping the uterine Rara gene confirms the association of SE with key genes that determine cell fate (12,13), as critical roles for retinoic acid (RA) signaling in uterine development have been described. RAR␣ is activated by retinol derivatives (Vitamin A), and disruption of retinol signal-ing, either by disrupting the Rara gene or with a retinol-deficient diet, impacts female reproductive tract development (43,44). RA is required for Mullerian duct development and for maintenance of epithelia of adult female reproductive tracts (45). Mechanistically, retinoic acid signaling influences uterine versus vaginal epithelial differentiation (46); therefore, precise regulation of RA signaling mediators is important for normal development and function. This is partly achieved via regulation of enzymes involved in retinoic acid synthesis and degradation (47), but clearly RAR expression is needed to mediate the signal, suggesting that SE in the Rara gene plays an important role in female reproductive tract development and function. Additionally, inappropriate perinatal estrogen exposures are known to disrupt the normal patterning of female reproductive tract mesenchyme, but this effect can be prevented by co-treatment with RA (45), illustrating an interplay between estrogen and RA signaling during female reproductive tract development. Evaluation of the Rara SE showed that it is within chromatin loops that include two members of the DNA replication complex, Cdc6 and Top2a (Fig. 5Ba and Table 4), which may indicate a role for the Rara SE in expression of genes involved in growth responses to E2.
There are four Hox clusters in the mouse genome (Hoxa, Hoxb, Hoxc, and Hoxd), and Hox9 -13 are believed to have arisen from a single ancestral gene (48). The Hoxa cluster, in particular, has been shown to play a key role in reproductive tract development and function (34). We were interested, then, to note that one of the top ranked prepubertal V SEs was located in the Hoxd cluster. Hoxa factors are expressed in a progressive pattern in the developing female reproductive tract, with Hoxa9 in the oviduct, Hoxa10 in the uterus, and Hoxa11 in the cervix (34). Hoxd9, -10, and -11 are expressed similarly to their Hoxa counterparts in the female reproductive tract (49) (Fig.  S5A). A role for Hoxd in uterine function has been described using complementation studies with Hoxa-and Hoxd-deleted mice (49). In addition, a recent study evaluated the impact of deletion of one allele each of Hoxa9, 10,11,Hoxc9,10,11, and Hoxd9,10,11 on uterine development. This study indicated that the combined deletion of one allele each of these nine Hox genes greatly decreased uterine gland development (50), further showing a role for Hoxd genes in uterine function. Although the Hoxa cluster has been very thoroughly described and characterized for its role in uterine function, we did not find any SEs at this gene cluster (Fig. S5B). We do observe chromatin loops between the Hoxa cluster and a region more than 1 megabase from the Hoxa cluster that includes SEs (Fig. S5B).
Ncor2, also called SMRT (silencing mediator of retinoic acid and thyroid hormone receptor), was originally described as a corepressor for RAR and TR, mainly through its interactions with these receptors and its associated histone deacetylase activity (35). Further study has revealed its dual role as a geneselective co-activator and co-repressor for ER␣ (51,52). A role for NCOR2 in the uterus has not been investigated, but considering its potential to regulate responses of TR, RAR, and ER␣, it is noteworthy that we see two SEs in the Ncor2 gene (Fig. S4A). Recently, it was shown that removing Ncor2 disrupts formation of a retinoid gradient in mouse embryos and leads to perturbation of HOXC expression patterns, greatly impacting proper   Includes CSF1, CSF2, EDN1, IFNG, IL1, IL17A, IL1A, IL1B, IL2, Il3, IL4, IL6,  OSM, TNF, and TNFSF11. c Includes EGF, ERBB2, MAPK3, MAPK8, MAPK9, and p38 MAPK.

Super-enhancers drive uterine estrogen responses
development (53). The disruptions observed incorporate developmental pathways involving three SE-associated uterine genes (Ncor2, Rara, and Hoxd), suggesting that these SEs might underlie processes important for pubertal development of uterine cells.
Zmiz1 is a member of the PIAS (protein inhibitor of activated STAT) family and was originally identified due to its role in prostate cancer cells as a co-activator that facilitates androgen receptor sumoylation (54). ZMIZ1 regulates activities of certain transcription factors, several of which are important for uterine development and function, including p53 (41), Notch (55), SMAD (24), STAT (56), and androgen receptor (57). Its essential role is highlighted by the finding that its deletion in mice causes embryonic lethality (12). Embryonic fibroblasts isolated from Zmiz1 null mouse embryos fail to proliferate normally (58,59), emphasizing its critical cellular role. A missense mutation in ZMIZ1 has been identified from endometroid cancer samples (Catalogue of Somatic Mutations in Cancer; RRID:   SCR_002260 (60)), and additionally, decreased ZMIZ1 expression (61) has been reported in adenomyosis biopsies, indicating an essential role for ZMIZ1 in human uterine health as well. We could readily detect Zmiz1 in the mouse uterus ( Fig. 5A and Table 4). Its expression was recently noted in rat uterine stromal cells, with down-regulation observed in response to circadian synchronization using the glucocorticoid agonist dexamethasone (62). SE-dependent expression of Zmiz1 could ensure optimal uterine responses mediated by pathways it regulates; however, the embryonic lethality of the global deletion prevents evaluation of uterine function with current models. We plan to examine its role in uterine tissue via conditional deletion of its expression.

Super-enhancers drive uterine estrogen responses
Our approach of identifying super-enhancers in developing uterine tissue is validated by observations that some SEs are at genes encoding factors known to be critical for uterine development and function, such as RAR␣. The location of most super-enhancers at the ends of chromatin loops and the uterine expression and developmental acquisition of E2 regulation of genes either at SEs or connected by looping to a SE indicate an important role for these enhancers. Indeed, regulatory pathways known to be involved in uterine development and function are impacted by these SE-associated E2-responsive genes. These important regulatory regions likely serve to optimize appropriate, efficient, and timely responses to circulating E2 during each estrous cycle and in preparation for establishing pregnancy.

Animals
All mice were used in accordance with an NIEHS-approved animal study protocol and using the 2015 edition of the Public Health Service (PHS) Policy on Humane Care and Use of Laboratory Animals. For E2 response experiments, intact 21-dayold or ovexed adult (10ϩ-week-old) female C57bl6/J mice were purchased from Charles River Laboratories (Raleigh, NC). ER␣KO mice (4) (also called Ex3␣ERKO) and their WT and heterozygous (ER␣ ϩ/Ϫ) littermates were produced in our contract colony at Taconic Farms (Albany, NY). IGF1enh4 KO mice and their WT littermates were produced in our colony at NIEHS and ovariectomized at 10 weeks of age. Genotypes were determined from an ear biopsy by Transnetyx (Cordova, TN). 21-Day-old females were used the week they were received; ovariectomized females were housed for 10 -14 days before the experiments to allow endogenous ovarian hormones to diminish. There was no blinding, and mice were randomly assigned to treatment groups. Mice were given a single intraperitoneal injection of 250 ng of E2 (Research Plus Inc., Barnegat, NJ) dissolved in 0.1 ml of normal saline, which is a 10 g/kg dose. Some were injected subcutaneously with 1 mg of P4 (Sigma) dissolved in 100 l of sesame oil (Sigma). Control V animals were injected with 0.1 ml of normal saline. Uterine tissue was collected 1, 2, 6, or 24 h after the injections and was snap-frozen in liquid nitrogen. RNA was isolated and cDNA was synthesized for real-time RT-PCR as described previously (22,36). Primer sequences are listed in Table S5. For gene expression studies, at least three animals per group were used based on an at least 2-fold change, with a coefficient of variation of 0.2 and 90% power, tested by two-way analysis of variance.

Identification of super-enhancers
Super-enhancers were identified based on the method described by Bojcsuk et al. (63). First, 1-h V and 1-h E2 ER␣ peak calls with high stringency were called by HOMER (parameters: -fdr 0.00001 -F 12 -style factor) for the ovexed adult and 21-day-old samples and then combined via BEDtools mergeBed (version 2.24.0), entailing more than 24,000 ER␣binding locations altogether. BEDtools mergeBed (version 2.24.0) was then rerun, this time to merge all peak calls within 12.5 kb. This set was subsequently filtered to retain only those regions with more than one contributing called peak, resulting in 4634 ER␣-binding enhancer regions, each with multiple ER␣-binding locations. Each enhancer region was scored by counting the number of overlapping uniquely mapped nonduplicate reads with BEDtools multiBamCov (version 2.24.0); reads in this calculation were H3K27Ac or H3K4Me1 ChIP-Seq data after extension of mapped read length to 200 nucleotides. After normalizing to 20 million reads, the input-subtracted signal was determined per region. The signal curve was plotted as H3K27Ac or H3K4Me1 signal versus region rank and rescaled to span 0 -1 on both axes. Definition of each enhancer region containing multiple ER␣-binding sites as a "typical enhancer" or "super-enhancer" was determined according to the elbow of the signal curve where slope ϭ 1.

SMC1a ChIP-Seq analysis
Processing of the SMC1a ChIP-Seq data was performed as described previously (36). Peak calls were made by HOMER version 4.10.3 via findPeaks with parameter "-style factor". The de novo motif analysis was done with HOMER version 4.10.3 findMotifsGenome.pl with parameter "-size given"; for consistency, query regions (SMC1a peaks at loop ends) were resized to 300 bp centered on the midpoint of the called peaks.

RNA-Seq analysis
Processing of the RNA-Seq data was performed as described previously (28). Read counts per gene were determined by Subread featureCounts version 1.5.0-p1 (64) with parameters "-s2 -Sfr -p" for gene models based on RefSeq transcripts as downloaded from the UCSC Table Browser on February 9, 2015. Differential gene analysis was performed with DESeq2 version 1.15.1 (65). Genes were considered expressed in the uterus if TPM Ն 1 in at least one condition. Genes were considered differentially expressed by applying a FDR cutoff of Յ0.01. SEassociated genes that were expressed or differentially expressed (FDR Ͻ 0.01) in the uterus were identified using the Partek List Manager Tool to find RefSeq gene symbols common to both sets.

Microarray
Two sets of uterine RNA samples were isolated. The first set was obtained from ovariectomized adult ER␣KO females and their WT littermates treated for 2 h with V or E2 as above. The second set was obtained from 21-day-old ER␣KO females and their WT littermates (produced by timed matings of heterozy-

Super-enhancers drive uterine estrogen responses
gous (ER␣ ϩ/Ϫ) males and females) treated for 2 h with V or E2. RNA was assessed by microarray as described previously (66) and analyzed in Partek using analysis of variance to find DEG and then filtering for signal intensity Ͼ100 in at least one sample and then filtering for super-enhancer-associated genes. Data from ovariectomized adults (68) (GSE100131) and data from 21-day-old mice (GSE148006) have been deposited in GEO.

Hi-C and data analysis
The Hi-C analysis was based on samples described previously (36) (V-1, E2, and P4-1 in Fig. 3A) together with additional samples. Samples V-2 and ER␣KO V were processed by Arima Genomics (San Diego, CA), and libraries were sequenced by the NIEHS Molecular Genomics Core, as described for the previous study (36). Samples V-3 and P4-2 were processed and sequenced by Active Motif, Inc. (Carlsbad, CA) using the Arima Genomics Hi-C kit. The Juicer (version 1.5.6) platform was used for processing the Hi-C samples as described previously (36,67). The quality control metrics from the Juicer output are shown in Table S6. Chromatin loops were identified with the Juicer hiccup utility, an algorithm for finding chromatin loops, by searching for clusters of contact matrix entries with enriched contact frequency relative to local background at default parameters (67). Simple overlap assessment via BEDtools inter-sectBed was used to determine localization of super-enhancers and SMC1a peaks at loop ends. The "atlas" of mouse uterine chromatin loops was generated by combining Juicer loop calls from the six individual WT samples, whereby all loops with both ends overlapping were collapsed.

Data availability
Hi-C and microarray data from samples are deposited in GEO under GSE147843 and GSE148006.
Acknowledgments-We are grateful to Drs. Alicia Chi and Joseph Rodriguez for critical reading of the manuscript. We thank the NIEHS Molecular Genomics Core for sequencing and the NIEHS Microarray Group for microarray. We greatly appreciate the animal care and surgical expertise provided by individuals in the NIEHS Comparative Medicine Branch. We thank Active Motif Inc. for providing complementary processing and sequencing of Hi-C samples V-3 and P4-2. Conflict of interest-The authors declare that they have no conflicts of interest with the contents of this article.