DNA Sequence Constraints Define Functionally Active Steroid Nuclear Receptor Binding Sites in Chromatin

Abstract

Gene regulatory programs are encoded in the sequence of the DNA. Since the completion of the Human Genome Project, millions of gene regulatory elements have been identified in the human genome. Understanding how each of those sites functionally contributes to gene regulation, however, remains a challenge for nearly every field of biology. Transcription factors influence cell function by interpreting information contained within cis-regulatory elements in chromatin. Whereas chromatin immunoprecipitation–sequencing has been used to identify and map transcription factor–DNA interactions, it has been difficult to assign functionality to the binding sites identified. Thus, in this study, we probed the transcriptional activity, DNA-binding competence, and functional activity of select nuclear receptor mutants in cellular and animal model systems and used this information to define the sequence constraints of functional steroid nuclear receptor cis-regulatory elements. Analysis of the architecture within sNR chromatin interacting sites revealed that only a small fraction of all sNR chromatin–interacting events is associated with transcriptional output and that this functionality is restricted to elements that vary from the consensus palindromic elements by one or two nucleotides. These findings define the transcriptional grammar necessary to predict functionality from regulatory sequences, with a multitude of future implications.

Steroid nuclear receptors (sNRs) are ligand-activated DNA-binding transcription factors (TFs) that induce or repress gene transcription following their interaction with specific DNA-binding sequences within enhancers associated with target genes (1). Specifically, sNRs mediate cellular response to their cognate hormones by interacting with DNA sequence-specific cis-regulatory elements: the estrogen response element (ERE) used by the estrogen receptors (ERs) and the hormone response element (HRE) used by the ketosteroid receptors (KRs): androgen receptor (AR), glucocorticoid receptor (GR), mineralocorticoid receptor (MR), and progesterone receptor (PR) (2, 3). These DNA elements not only serve to address sNRs to target genes but also contain information that influences the transcriptional activity of bound sNR-coregulator complexes. The functionality of these response elements is defined by the context in which they operate within chromatinized DNA. However, the structural and primary sequence features that distinguish functional binding elements from those that permit sNR binding but are not associated with transcription have not yet been defined, a situation that limits the interpretation of data from genome-wide binding studies of sNRs.

In contrast to the classical model of sNR function, the vast majority of sNR chromatin interacting events do not occur near promoters (4, 5), a characteristic of most TFs (6). Instead, sNRs interact predominantly in introns and distal intergenic regions (enhancers), often at distances exceeding 10-kilobase (kb) pairs from the nearest transcription start site (TSS) (7). Distal cis-acting enhancers are short [20 to 400 base pairs (bps)] noncoding DNA elements that direct cell type–specific patterns of gene expression, regulated by the binding of signal-dependent TFs and their associated coregulators (8). These remote elements can regulate transcription independently of their location, orientation, position, and distance from a promoter (9), making it difficult to link TF binding events to target gene regulation. Currently used bioinformatic approaches assume that binding events located up- or downstream and/or occurring most proximal to regulated genes are those most likely to be involved in sNR-mediated responses (10). As a result, bias is introduced into most analyses by selecting for sites that are closer to the TSS, an approach that does not adequately account for the important contributions of distal binding sites in specific gene regulation and prematurely assumes those sites and the DNA elements within those sites are involved in the functional outcome of the most proximal gene.

In addition to the difficulty in assigning functionality to TF binding events genome-wide, sNR binding analyses have also revealed that the number of binding events is at least an order of magnitude greater than the number of genes regulated by the hormone (11), raising the question of how many of these sNR binding events are functionally active. Furthermore, a substantial portion of the sNR binding sites identified in genome-wide binding studies do not contain a canonical consensus, palindromic element, challenging the notion that sNRs bind only to stringent paired half-sites with highly specific orientations (12). These observations have led to the suggestion that sNRs can also tether, in a nonclassic manner, to other DNA-bound TFs such as AP1, Sp1, and NF-κB. However, the functional relevance of these noncanonical chromatin interacting sites is unknown. More recent studies have revealed that most sNR binding events do not correspond to functional enhancers, suggesting that the majority of binding events may merely reflect DNA accessibility rather than functionality (11, 13). It is unclear why only some chromatin-interacting events are functional and actively involved in transcriptional regulation. These findings highlight the need for strategies to distill these protein-DNA interactions into functional and nonfunctional categories, the objective of this study.

For this study, sNR DNA-binding domain (DBD) mutant mouse models were used to dissect the relative biological roles of cis-regulatory DBD-dependent transcriptional regulation. Furthermore, we used our own data together with publicly available chromatin immunoprecipitation–sequencing (ChIP-seq) and global run-on–sequencing (GRO-seq) datasets of sNRs and their collaborating factors in mouse tissues and human cell lines. Using these data, we systematically analyzed the DNA sequence constraints of sequence-specific cis-regulatory elements and defined the elements within these sequences that track with and likely define functionality: sites we refer to as Nuclear Receptor Functional Enhancers (NRFEs).

Experimental Procedures

Quality control and peak selection criteria

In total, 1332 publicly available datasets (plus 4 new ChIP-seq datasets: EAAE-Input, EAAE-Vehicle, EAAE-E2, WT–Vehicle = 1336) were obtained from Gene Expression Omnibus and converted to fastq format using fastq-dump v2.4.5; original dataset names were retained. Reads were first selected using a cross-correlation analysis (trim_and_filter_SE.pl) (Supplemental Table 6) (14). The sequencing reads were then mapped uniquely, allowing for no more than two mismatches, to the reference genome [mouse genome build 38 (mm10) or human genome build 37 (hg19)] using Bowtie v1.1.2 (15). Mapped reads were deduplicated using MarkDuplicates.jar from the Picard tools package v1.96. Peak selection was performed using Hypergeometric Optimization of Motif EnRichment (HOMER) v4.7.2. Multiple peak selection criteria were used (L4, L8, L10, L15, L20), where Lx represents an x-fold greater tag density at peaks than in the surrounding 10-kb region. This performs a low to high stringency analysis of the data. To identify enriched regions of variable length for RNA polymerase II (RNAPII) studies, we required a 1000-bp gap (D1000) in the signal for adjacent peaks to be considered part of the same region. Similarly, peak selection criteria D2500, D5000, D10000, and D20000 were used. RNAPII peaks selected using the L parameter are enriched in promoter regions, whereas RNAPII peaks selected using the D parameter are enriched in introns (Supplemental Tables 4:2–3). Nascent RNA production at specific sites in the genome was selected when the body-fold change exceeded fourfold (B4). Peak selection criteria B8, B10, B15, and B20 were also used. By studying the data with respect to a multiple spectrum of peak selection criteria, we adjust for both the risk of excess background noise and the risk of filtering out any low-amplitude information. See Supplemental Table 6 for script codes and data analysis of the 1336 datasets.

ERE and HRE motifs contained within ChIP-seq peaks and GRO-seq reads

The location coordinates of the 3676 0- to 3-nucleotide (nt) variants of the ERE and HRE on the positive/sense strand were identified in the mouse (mm10) and human (hg19) genomes using OligoMatch (University of California, Santa Cruz). Here ERE is the estrogen response element (5′-GGTCAnnnTGACC-3′) and HRE is the hormone response element (5′-GAACAnnnTGTTC-3′). Also, the 0- to 3-nt variants include the 1 unique perfect consensus element, 30 possible 1-nt variants, 405 possible 2-nt variants, and 3240 possible 3-nt variants, for a total of 3676 sequences.

Four nucleotide possibilities in a motif with 10 primary positions:

• N = number of primary positions in the motif = 10

• L = number of nucleotide possibilities = 4

• K = number of variants

Binomial Equation:  

$${\left(x+y\right)}^n={\displaystyle \sum }_{k=0}^n\left(\begin{array}{c}n\\ k\end{array}\right)x^k{y}^{n-k}$$

 

Combinations: $$L^n={\left[\left(L-1\right)+1\right]}^n={\displaystyle \sum }_{k=0}^n\left(\begin{array}{c}n\\ k\end{array}\right){\left(L-1\right)}^k {\left(1\right)}^{n-k}$$

 

$$={\displaystyle \sum }_{k=0}^n\left(\begin{array}{c}n\\ k\end{array}\right){\left(L-1\right)}^k=2^{20}=4^{10}$$

A customized C++ program, cpp_custom_full.cpp, was used to overlap location coordinates of the ERE and HRE motifs and the ChIP-seq peaks or Gro-seq reads from each dataset to determine how many peaks contained a 0- to 3-nt ERE or HRE motif (the entire 13-nt motif was required to be within the peak boundaries) (Supplemental Table 6). Motif assignment was given to the perfect consensus motif or the motif with the least number of variants from the perfect consensus motif (0-nt variant > 1-nt variant > 2-nt variant > 3-nt variant) within each peak, thus defining each peak by a single (unique) motif. See Supplemental Table 5 for counts of the ERE and HRE motifs in the mouse (mm10) and human (hg19) genomes. See Supplemental Tables 1, 2, 3, and 4 for number of overlaps between motif location and ChIP-seq peaks or Gro-seq reads.

Results

Regulatory element motif discovery

TFs engage their cognate enhancers by interacting with specific, high-affinity cis-acting DNA sequences. Previous studies have evaluated the sequence requirements of cis-regulatory elements within specific enhancers, although the insights from these single-gene analyses have been difficult to apply to data collected from genome-wide surveys of cis-regulatory elements. Thus, in this study, we have evaluated gene regulation in response to steroid hormones as a model system for other nonsteroid transcription systems that can be exploited to define general rules that establish the functionality of sequences that confer steroid responsiveness upon target genes.

As a starting place in these studies, we identified the location coordinates of each perfect consensus ERE motif (5′-GGTCAnnnTGACC-3′) and HRE motif (5′-GAACAnnnTGTTC-3′) in the mouse and human genomes using OligoMatch (University of California, Santa Cruz). In this manner, 2367 perfect consensus ERE motifs were identified in the mouse genome, and 2194 were identified in the human genome (Table 1). Similarly, 3444 perfect consensus HRE motifs were identified in the mouse genome and 3535 in the human genome (Table 1). The location coordinates for the three alternative nt possibilities (variants) at each of the 10 primary positions in the 13nt perfect consensus ERE and HRE motif were also obtained and referred to as 1-nt variants (30 sequences) (Table 1). This approach was also used to identify all 2-nt variants (405 sequences) and 3-nt variants (3240 sequences) of the perfect consensus motif in both the mouse and human genomes (Table 1). Considering the information content of a typical 10-nt sNR motif, the frequency that any 10-nt sequence that has a maximum possibility of 4 nts at each position will occur is once every 1,048,576 nts (4¹⁰) at random occurrence—specifically, 2525 and 2729 times in a genome made up of 2,647,531,771 nts (mouse) or 2,861,339,759 nts (human) (Supplemental Experimental Procedures-Genome Assemblies). Therefore, the occurrence of the perfect consensus ERE and HRE motif in the genomes is essentially that which would be expected to occur at random. This random distribution is also displayed by the 1- to 3-nt variant groups (Table 1). Furthermore, there is no genomic distribution preference of ERE or HRE motif sequences toward any specific genomic feature—3′UTR, 5′UTR, intergenic, TSS, exon, intron, noncoding, and promoter-TSS (Supplemental Fig. 1A)—compared to the whole genome (Supplemental Fig. 1B). Although many studies have analyzed the effects of multiple copies of sNR consensus elements, their spacing and location on the transactivation of a reporter gene (16), the near-random occurrence and distribution of these sequences within the genome, suggests that the information content of these sequences would not be significantly different from any random distribution in the genome. Moreover, the number of 1-nt variant EREs and HREs (Table 1) in the genomes far exceeds the number of sNR chromatin interaction sites detected in previously published ChIP-seq experiments. Based on information theory/combinatorics alone, sNRs bind far fewer regions than would be expected. That is, only a small fraction of the near-consensus binding sequences occurring in the genome is occupied by a TF in a given cell type at a given moment (17). This is most likely due in large part to the restrictive nature of chromatinized DNA (18). These observations emphasize the importance of chromatin accessibility and nucleosome occupancy, determined primarily by pioneering factors during cellular development and differentiation, as major determinants for cell type–specific sNR signaling.

Table 1.

Number of 0- to 3-nt Variant ERE and HRE Motifs in the Mouse and Human Genomes $\mathbf{ERE}\mathbf{=}\mathbf{}\underset{1 2 3 4 5 6 7 8 9 10}{\mathbf{G}\underset{¯}{\mathbf{GT}}\mathbf{CAnnnTG}\underset{¯}{\mathbf{AC}}\mathbf{C}}$$\mathbf{HRE}\mathbf{=}\mathbf{}\underset{1 2 3 4 5 6 7 8 9 10}{\mathbf{G}\underset{¯}{\mathbf{AA}}\mathbf{CAnnnTG}\underset{¯}{\mathbf{TT}}\mathbf{C}}$

Combinatorial Counts						Actual Motif Frequency in Genomes
A	B	C	D	E	F	Mouse Genome (mm10)		Human Genome (hg19)
	ERE	Comb. C = $\left(\begin{array}{c}n\\ k\end{array}\right)$	4 nts D = C × 4A	Total Unique E = ΣF	Unique F = C × 3A	Total Unique	Unique	Total Unique	Unique
0	Perfect ERE	1	1	1	1		2367		2194
1	1-nt Variant ERE	10	40	31	30	73,795	71,428	62,507	60,313
2	2-nt Variant ERE	45	720	436	405	971,950	898,155	977,233	914,726
3	3-nt Variant ERE	120	7680	3676	3240	7,722,557	6,750,607	8,493,417	7,516,184
A	B	C	D	E	F	Mouse Genome (mm10)		Human Genome (hg19)
	HRE	Comb. C =$\left(\begin{array}{c}n\\ k\end{array}\right)$	4 nts D = C × 4A	Total Unique E = ΣF	Unique F = C × 3A	Total Unique	Unique	Total Unique	Unique
0	Perfect HRE	1	1	1	1		3444		3535
1	1-nt Variant HRE	10	40	31	30	100,483	97,039	108,302	104,767
2	2-nt Variant HRE	45	720	436	405	1,437,999	1,337,516	1,447,845	1,339,543
3	3-nt Variant HRE	120	7680	3676	3240	11,899,196	10,461,197	12,219,004	10,771,159
A	B	C	D	E	F
		Comb. C =$\left(\begin{array}{c}n\\ k\end{array}\right)$	4 nts D = C × 4A	Total Unique E = ΣF	Unique F = C × 3A
0	Perfect	10!/(0! × 10!)	1 × 1	(3° × 1) + 0	(3° × 1)
1	1-nt Variant	10!/(1! × 9!)	10 × 4	(31 × 10) + 1	(31 × 10)
2	2-nt Variant	10!/(2! × 8!)	45 × 16	(32 × 45) + 31	(32 × 45)
3	3-nt Variant	10!/(3! × 7!)	120 × 64	(33 × 120) + 436	(33 × 120)

The number of 0- to 3-nt variants of the ERE and HRE on the positive/sense strand identified in the mouse (mm10) and human (hg19) genomes using OligoMatch (University of California, Santa Cruz). Here ERE is the estrogen response element (5′-GGTCAnnnTGACC-3′) and HRE is the hormone response element (5′-GAACAnnnTGTTC-3′). The 0- to 3-nt variants include the 1 unique perfect consensus element, 30 1-nt variants, 405 2-nt variants, and 3240 3-nt variants, for a total of 3676 sequences.

Use of sNR DBD mutants to distinguish functional from neutral DNA-binding events

Biochemical studies have in the past indicated that sNRs can interact with DNA in both a specific (high-affinity) and a nonspecific (low-affinity) manner. Indeed, often overlooked, these two properties of sNRs were exploited in DNA chromatography protocols to purify sNRs (19). These activities, specific vs nonspecific DNA binding, were shown to be functionally separable as mutations in the DBD attenuated the interaction of the receptor with specific DNA elements, but had minimal effect on general DNA-binding properties (19). Considering this established biochemical property, we used sNR-DBD mutants to assess the extent to which DNA-binding activities reported in ChIP-seq studies can be attributed to “specific” vs “nonspecific” binding events. For this analysis, we used two previously described ERα DBD mutant mouse models (20). In one model, two point mutations, E207A and G208A, were introduced into the P-box of the first zinc finger in the DBD [ERα DBD mutant (ERα-KIKO)] [Fig. 1(a)]. In the second mouse model [ERα DBD Mutant (ERα-EAAE)], four point mutations that flank the P-box of the first zinc finger were introduced in the DBD: Y201E, K210A, K214A, and R215E [Fig. 1(a)]. We chose to use mouse models for these studies as they have the advantage of allowing an examination of the impact of the selected mutations on physiological responses to 17-β estradiol (E2) in a natural in vivo biological setting. E2 does not induce mammary gland growth or increase uterine weight in female ERα-KIKO or ERα-EAAE mice (20). Although the regulation of classical ERα target genes is diminished in the ERα-KIKO model, this particular receptor mutant exhibits an aberrant DNA-binding specificity with reduced binding of ERE motifs but has acquired the ability to bind HRE motifs and thus inappropriately activates some non-E2-regulated genes (20), a utility that we have exploited in these studies. Of specific importance to this study was the observation that the phenotypes observed in the ERα-EAAE mice are indistinguishable from those observed in the loss-of-function ERα-null mouse model (20, 21). Furthermore, genome-wide gene expression profiles (validated by quantitative reverse transcription polymerase chain reaction) revealed that gene expression in the uterus and liver of ERα-EAAE mice was completely unresponsive to E2 (20, 21). These two experimental model systems confirm in vivo that a functional ERα-DBD is required for all E2-mediated transcriptional responses (in the uterus and liver), both transcriptional induction and repression, and that indirect DNA-binding activity such as “tethering” cannot independently initiate gene transcription in response to E2. Previously, others have reported the development of a mouse model in which the DBD of the GR was mutated with a single point mutation, A458T, in the second zinc finger [GR DBD mutant (GR-Dim)]. [Fig. 1(b)] (22). Analogous to what was observed in the ERα-EAAE model, tissues derived from the GR-Dim model are refractory to endogenous glucocorticoids (23–27). Thus, we conclude from this study and our own previously published data that a functional DNA-binding domain is required for GR- and ER-dependent transcriptional activity. It also suggests that the general low-affinity DNA-interacting activity of sNRs is not, in and of itself, sufficient for transcriptional activity and affords us the opportunity to evaluate which sNR DNA-binding activity is required for the recognition of each of the DNA-binding elements we have defined above.

Figure 1.

sNR DBD mutant mouse models. (a) Two ERα DBD mutant mouse models of ERα: two point mutations at amino acids E207A and G208A in the P-box of the first zinc finger (ERα-KIKO) and four point mutations at amino acids Y201E, K210A, K214A, and R215E that flank the P-box of the first zinc finger (ERα-EAAE). (b) GR DBD mutant mouse model of the GR: an A458T point mutation introduced into the second zinc finger (GR-Dim).

As the physiological and transcriptional activities of the ERα-EAAE mutant mice were similar to those of the ERα-null mouse (20), we evaluated its DNA-binding activity at a genome-wide level using chromatin isolated from the uteri of E2-treated ERα-EAAE (EAAE-E2) mice. The dataset generated in this study was compared with ChIP-seq data derived in the same manner from uteri of E2-treated mice expressing intact ERα (WT-E2-1hr) or ERα-KIKO (KIKO-E2) (28). For comparative and control purposes, we also included ChIP-seq data from a study used to identify PR binding sites in the uteri of ovariectomized mice treated with progesterone (P4) in this analysis (Uterus-PGR-P4) (29). Irrespective of the ChIP-seq peak selection criteria (L), where Lx represents an x-fold greater tag density at peaks than in the surrounding 10-kb region, it was determined that the extent (i.e., number of binding events in the presence of hormone relative to vehicle conditions) of ERα recruitment to DNA in response to E2 was indistinguishable from that observed in parallel studies which assessed the recruitment of the ERα-EAAE to DNA (Supplemental Table 1:1). This result indicated that the majority of the ERα-chromatin binding events recorded in these ChIP-seq experiments did not appear to require a functional ERα-DBD. Furthermore, because ERα-EAAE did not support any measurable E2-dependent transcription, it suggested that the majority of the chromatin interacting events identified in ERα ChIP studies are not linked directly to transcriptional events. A comparison of (1) ERα and ERα-KIKO and (2) GR and GR-Dim ChIP-seq and ChIP-exonuclease (ChIP-exo) data also revealed that although the disruption of the DBD within these receptors significantly attenuated their transcriptional activity, it had a minor effect on their ability to interact with chromatin (Supplemental Tables 1:2–7). Of note, it was observed that the total number of peaks between replicate samples is highly variable, indicating that a normalization methodology will be required for further comparative analysis. With these data in hand, we then examined whether there were differences in the sequence constraints in ERE/HRE sequences that supported transcription vs those that did not appear to do so (i.e., serve as nonfunctional binding sites).

Only a subset of DNA elements occupied by ERα can be attributed to functional binding events

The functional relevance of the different classes of ERα-binding elements described above was further explored by evaluating the interaction of ERα (or ERα-DBD mutants) with chromatin isolated from the uteri of ovariectomized mice following acute exposure to E2. In this evaluation, the absolute number of ChIP-seq peaks containing each specific type of motif was counted. Because ChIP-seq peaks can contain multiple sNR-binding motifs, assignment was given to the perfect consensus element or the motif with the least number of variants from the perfect consensus motif (0-nt variant > 1-nt variant > 2-nt variant > 3-nt variant) within each peak, thus defining each ChIP-seq peak by a single motif. Under the least stringent peak selection criteria (L4, default peak selection criteria for a TF), 1201 peaks out of a total of 76,163 peaks (2%) in the ERα dataset (WT-E2-1hr) contain a perfect ERE motif [Fig. 2(a); Supplemental Table 2:1]. The proportion of peaks that contain a perfect ERE motif increased as the peak selection criteria (L) were increased, confirming that the most robust ERα-associated ChIP-seq peaks are enriched in perfect ERE motifs [Fig. 2(a); Supplemental Table 2:1]. This profile was also seen in a previously published ERα mouse uterine ChIP-seq dataset (Uterus-ER) (30) [Fig. 2(a); Supplemental Table 2:1]. In contrast, when analyzed in the same manner, the ERα DNA-binding domain mutants, ERα-KIKO (KIKO-E2) or ERα-EAAE (EAAE-E2), had a significantly smaller proportion of peaks that contained a perfect consensus ERE motif across all peak selection criteria [Fig. 2(a); Supplemental Table 2:1]. As expected, very few of the ChIP-seq peaks identified in the ChIP-seq study that evaluated the interaction of PR with chromatin isolated from the uteri of P4-treated mice (Uterus-PGR-P4) contained a perfect consensus ERE motif [Fig. 2(a); Supplemental Table 2:1]. A similar analysis revealed that peaks containing a 1-nt variant ERE were only apparent in the datasets generated with wild-type ERα [Fig. 2(b); Supplemental Table 2:1]. Peaks that contained a 2-nt variant ERE were also found in the wild-type ERα-datasets [Fig. 2(c); Supplemental Table 2:1]. However, they were also observed, albeit at a reduced proportion, in datasets generated with the ERα DBD mutants [Fig. 2(c); Supplemental Table 2:1]. Of note, the proportion of peaks that contained a 3-nt variant ERE remained uniform as the peak selection criteria were increased in the ERα-EAAE, ERα-KIKO, and Uterus-PGR-P4 datasets, whereas those identified in the wild-type ERα-datasets decreased [Fig. 2(d); Supplemental Table 2:1]. A similar pattern was observed when all remaining peaks (i.e., did not contain a 0- to 3-nt variant ERE) were examined [Fig. 2(e); Supplemental Table 2:1]. Motif discovery analysis, performed on all five groups of the elements under study, suggested that ERα-EAAE did not gain any detectable cryptic DNA-binding function or lose any other enriched motif that could confound the analysis of the data (Supplemental Fig. 2). Considering that ERα-EAAE was incapable of supporting E2-dependent transcription yet was able to interact with 2- and 3-nt variant EREs suggests that ERα binding to perfect consensus EREs, 1-nt variant EREs, and a subset of 2-nt variant EREs are required for E2-regulated transcriptional response in the uterus. It is further concluded that motifs that vary from the perfect consensus element by 2 nts represent a transition point from specific (high-affinity) to nonspecific (low-affinity) binding elements, with those containing more than 2-nt variants being considerably less likely to represent functional elements.

Figure 2.

Defining functionally active binding sites using sNR DBD mutant mouse models. (a–e) The proportion of ChIP-seq peaks that contain a 0- to 3-nt variant ERE (and remaining peaks) from five mouse uterine ChIP-seq datasets at L4 to L20 peak selection criteria: WT-E2-1hr, Uterus-ER, KIKO-E2, EAAE-E2, and Uterus-PGR-P4. (f–j) The proportion of ChIP-seq peaks that contain a 0- to 3-nt variant HRE (and remaining peaks) from five mouse uterine ChIP-seq datasets at L4 to L20 peak selection criteria: WT-E2-1hr, Uterus-ER, KIKO-E2, EAAE-E2, and Uterus-PGR-P4. (k–o) The proportion of ChIP-seq peaks that contain a 0- to 3-nt variant HRE (and remaining peaks) from four mouse liver ChIP-seq datasets at L4 to L20 peak selection criteria: GR-WT-pred-1, GR-WT-pred-2, GR-Dim-pred-1, and GR-Dim-pred-2. See also Supplemental Tables 2:1, 2:2, and 2:5.

Functional KR binding elements exhibit sequence constraints similar to those of EREs

The KRs (AR, GR, MR, PR) are a closely related family of sNRs that are structurally similar to ER. Thus, it was of interest to define the extent to which the DNA sequence constraints for functionality of an ERE were also apparent for an HRE. For this analysis, two ChIP-seq datasets were used: (1) PR interaction with chromatin from the uteri of P4-treated mice (Uterus-PGR-P4) and (2) ERα-KIKO (KIKO-E2) interaction with chromatin from the uteri of E2-treated mice. As observed for ERα/EREs, the proportion of Uterus-PGR-P4 and ERα-KIKO peaks that contained a perfect HRE or 1-nt variant HRE increased as the peak selection criteria were increased [Fig. 2(f) and 2(g); Supplemental Table 2:2]. Peaks that contained a 2-nt variant HRE were also found in these datasets [Fig. 2(h); Supplemental Table 2:2]. However, they were also observed, at a reduced proportion, in datasets generated with ERα and ERα-EAAE (serving as a control in this case) [Fig. 2(h); Supplemental Table 2:2]. The proportion of peaks containing a 3-nt variant HRE in the Uterus-PGR-P4 and ERα-KIKO datasets was uniform with the ERα and ERα-EAAE datasets [Fig. 2(i); Supplemental Table 2:2]. All remaining peaks (i.e., did not contain a 0- to 3-nt variant HRE) in the Uterus-PGR-P4 and ERα-KIKO datasets decreased as the peak selection criteria were increased while remaining uniform with the ERα and ERα-EAAE datasets [Fig. 2(j); Supplemental Table 2:2]. The acquired ability of ERα-KIKO to bind HRE motifs but reduced ability to bind ERE motifs was further confirmed in studies performed in ERα-KIKO–expressing mouse mammary epithelial cells (Supplemental Tables 2:3 and 2:4).

We also evaluated the architecture of the DNA elements that enabled the transcriptionally inactive GR-Dim mutant to interact with chromatin. This analysis revealed that GR-Dim had a significant reduction in the proportion of binding peaks that contained a perfect HRE or 1-nt variant HRE compared with wild-type GR, demonstrated decreased binding at 2-nt variant HREs, but had no such deficit at 3-nt variant HREs and all remaining peaks compared with wild-type GR [Fig. 2(k)–2(o); Supplemental Tables 2:5–10]. This suggests that GR binding to perfect HREs, 1-nt variant HREs, or a subset of 2-nt variant HREs is required to support all glucocorticoid-regulated transcriptional responses. A substantially similar conclusion was made in a recent study that determined that only those GR binding events that contained a HRE cis-regulatory element were able to confer glucocorticoid responsiveness upon a heterologous promoter via a STARR-seq reporter assay (31). Likewise, two DBD mutants in the AR, C562S and G568W, also exhibited decreased binding at perfect, 1-nt variant, and 2-nt variant HREs compared with wild-type AR, with no such deficit at 3-nt variant HREs and all remaining peaks compared with wild-type AR (Supplemental Table 2:11). Furthermore, this binding profile was also observed in primary prostate tumors that had acquired resistance to therapy, compared with tumors that were still sensitive to androgen deprivation (Supplemental Table 2:12). Conversely, a constitutively active AR mutant, ARv567es, notably involved in cancer progression and used to model prostate cancer resistance, as well as the constitutively active AR-Q640X variant, showed increased recruitment of AR only to sites that possessed a perfect, 1-nt, or 2-nt HRE compared with wild-type AR (Supplemental Tables 2:13 and 2:14). This binding profile was also observed in the AR-K2R SUMOylation mutant, K386R/K520R, which proliferate faster than wild-type AR (Supplemental Table 2:15). The similar DNA sequence constraints, observed in our studies of EREs and HREs, indicating that functionality is only associated with perfect consensus motifs or those that vary by 1-nt and, in some cases, by 2-nt variants, identify previously unappreciated constraints in the architecture of functional cis-regulatory elements. To distinguish these functional elements from those that support binding alone, we refer to them as NRFEs.

KRs exhibit a lower specific DNA-binding capacity than ERs

Next, the extent to which sequence constraints in its cis-regulatory element, observed with ERα in the uterus, was conserved in different tissues and species was examined. Analysis of 157 publicly available ER ChIP-seq datasets revealed agreement with the observations made in the studies in the mouse uterus. The proportion of peaks that contained a perfect ERE or 1-nt variant ERE increased as the peak selection criteria were increased [Fig. 3(a) and 3(b); Supplemental Table 2:16]. There was a subtle increase in the proportion of peaks that contained a 2-nt variant ERE [Fig. 3(c); Supplemental Table 2:16] and a decrease in the proportion of peaks that contained a 3-nt variant ERE and all remaining peaks in each dataset as the peak selection criteria were increased [Fig. 3(d) and 3(e); Supplemental Table 2:16]. On average, 4% to 7% (L4 to L20) of the peaks in each dataset contained a perfect ERE, 12% to 19% contained a 1-nt variant ERE, 21% to 26% contained a 2-nt variant ERE, 35% to 30% contained a 3-nt variant ERE, and 28% to 18% of the peaks did not contain a perfect, 1-nt, 2-nt, or 3-nt variant ERE (Supplemental Table 2:16). Thus, the signal intensity of ER ChIP-seq peaks does not distinguish between peaks containing an NRFE ERE vs peaks that do not contain an NRFE ERE: at L20 peak selection criteria (i.e., peaks with the highest signal intensity), ∼52% of the peaks contained an NRFE ERE, whereas ∼48% did not contain an NRFE ERE. Therefore, although we show that peaks with highest signal intensity include NRFE EREs at a higher proportion, many individual examples (∼48% of the ChIP-seq peaks within a given ER dataset) exist where this is not the case. Consequently, motif status of a given peak cannot be inferred through visual inspection. Of note, the total number of peaks across these 157 ER datasets ranged from 216 to 84,494 peaks (L4) or 126 to 33,122 peaks (L20) (Supplemental Table 2:16). That is, the proportion of peaks containing these elements was scale invariant over two to three orders of magnitude (Supplemental Table 2:16). Thus, the relative peak numbers were stable, whereas the actual peak numbers were highly variable.

Figure 3.

KRs exhibit a lower specific DNA-binding capacity than ERs. (a–e) The proportion of ER ChIP-seq peaks that contain a 0- to 3-nt variant ERE (and remaining peaks) from 157 ER ChIP-seq datasets at L4 to L20 peak selection criteria. (f–j) The proportion of KR ChIP-seq peaks that contain a 0- to 3-nt variant HRE (and remaining peaks) from 194 ChIP-seq datasets at L4 to L20 peak selection criteria. See also Supplemental Tables 2:16 and 2:17.

A similar profile (i.e., increase in proportion of peaks containing NRFE HREs and decrease in proportion of peaks containing 3-nt variant HREs and all remaining peaks in each dataset as the peak selection criteria were increased) was also seen in 194 publicly available KR (AR, GR, MR, PR) ChIP-seq datasets [Fig. 3(f)–3(j); Supplemental Table 2:17]. In comparison with the ER, however, on average, only 1% to 2% (L4 to L20) of the peaks in each dataset contained a perfect HRE, 6% to 11% contained a 1-nt variant HRE, 20% to 27% contained a 2-nt variant HRE, 40% to 38% contained a 3-nt variant HRE, and 32% to 22% of the peaks did not contain a perfect, 1-nt, 2-nt, or 3-nt variant HRE. Thus, the signal intensity of KR ChIP-seq peaks does not distinguish between peaks containing an NRFE HRE vs peaks that do not contain an NRFE HRE: at L20 peak selection criteria (i.e., peaks with the highest signal intensity), ∼40% of the peaks contained an NRFE HRE, whereas ∼60% did not contain an NRFE HRE. Therefore, although we show that peaks with highest signal intensity include NRFE HREs at a higher proportion, many individual examples (∼60% of the ChIP-seq peaks within a given KR dataset) exist where this is not the case. Consequently, motif status of a given peak cannot be inferred through visual inspection. Of note, the total number of peaks across these 194 datasets ranged from 328 to 122,550 peaks (L4) or 20 to 72,726 peaks (L20) (Supplemental Table 2:17). That is, the proportion of peaks containing these elements was scale invariant over two to three orders of magnitude (Supplemental Table 2:17). Thus, the relative peak numbers were stable, whereas the actual peak numbers were highly variable.

Indicators of enhancer activation are specifically enriched at NRFEs in response to hormone

Although a large portion of the genome possess enhancer-related chromatin features, only a fraction of enhancers are functionally active (32). ChIP-seq has been used to not only understand the genome-wide binding profiles of TFs but also to study the behavior of RNAPII, coregulators, and chromatin remodeling proteins. Because it is difficult to link individual enhancers to specific genes, we took the alternate approach of analyzing publicly available datasets to examine the extent to which genomic features that have been correlated with enhancer activation are enriched at NRFEs. For this analysis, a variety of genomic features, including enhancer RNA (eRNA) transcription, RNAPII occupancy, and coregulator and TF recruitment, were examined (Table 2).

Table 2.

Indicators of Enhancer Activation Are Specifically Enriched at NRFEs in Response to Hormone

Figure/Table ID	Technology	Comparing	Motif	Peak Selection Criteria	Datasets
Hormone-regulated eRNA transcription
Fig. 4(a); Supplemental Table 3:1	Gro-Seq	Veh vs E2	ERE	B4	74
Fig. 4(b); Supplemental Table 3:3	Gro-Seq	Veh vs DHT or Dex	HRE	B4	10
Supplemental Table 3:2	Gro-Seq	Veh vs E2	ERE	B4–B20	74
Supplemental Table 3:4	Gro-Seq	Veh vs DHT or Dex	HRE	B4–B20	10
Hormone-regulated RNAPII occupancy
Fig. 4(c); Supplemental Table 3:5	RNAPII ChIP-Seq	Veh vs E2	ERE	D1000	8
Fig. 4(d); Supplemental Table 3:12	RNAPII ChIP-Seq	Veh vs Cort, Dex, or Pred	HRE	D1000	3
Supplemental Table 3:6	RNAPII ChIP-Seq	Veh vs E2	ERE	D1000–D20000	8
Supplemental Table 3:7	RNAPII ChIP-Seq	Veh vs E2	ERE	L4–L20	8
Supplemental Table 3:8	RNAPII ChIP-Seq: ERα-WT and ERα-KIKO	Veh vs E2	ERE	L4–L20	2
Supplemental Table 3:9	RNAPII ChIP-Seq: ERα-WT and ERα-KIKO	Veh vs E2	HRE	L4–L20	2
Supplemental Table 3:10	RNAPII ChIP-Seq	Veh vs E2, Tam, or Fulv	ERE	L4–L10	3
Supplemental Table 3:13	RNAPII ChIP-Seq	Veh vs Cort, Dex, or Pred	HRE	D1000–D20000	3
Hormone-regulated recruitment of coregulators and chromatin remodeling proteins
Fig. 4(e); Supplemental Table 3:14	p300 ChIP-Seq	Veh vs E2	ERE	L4	2
Fig. 4(f); Supplemental Table 3:16	p300 ChIP-Seq	Veh vs Dex, T, P4, or R5020	HRE	L4	32
Supplemental Table 3:15	p300 ChIP-Seq	Veh vs E2	ERE	L4–L20	2
Supplemental Table 3:17	p300 ChIP-Seq	Veh vs Dex, T, P4, or R5020	HRE	L4–L20	32
Fig. 4(g); Supplemental Table 3:18	FoxA1 ChIP-Seq	Veh vs E2	ERE	L4	20
Fig. 4(h); Supplemental Table 3:20	FoxA1 ChIP-Seq	Veh vs Dex, T, or R1881	HRE	L4	10
Supplemental Table 3:19	FoxA1 ChIP-Seq	Veh vs E2	ERE	L4–L20	20
Supplemental Table 3:21	FoxA1 ChIP-Seq	Veh vs Dex, T, or R1881	HRE	L4–L20	10
Supplemental Table 3:22	FoxA1 ChIP-Seq	AR-C562S vs AR-WT	HRE	L4–L20	1
Supplemental Table 3:23	FoxA1 ChIP-Seq	AR-WT vs AR-Q640X	HRE	L4–L20	1
Fig. 4(i); Supplemental Table 3:24	DNase HSS ChIP-Seq	Veh vs E2	ERE	L4	12
Fig. 4(j); Supplemental Table 3:26	DNase HSS ChIP-Seq	Veh vs Cort, Dex, or T	HRE	L4	20
Supplemental Table 3:25	DNase HSS ChIP-Seq	Veh vs E2	ERE	L4–L20	12
Supplemental Table 3:27	DNase HSS ChIP-Seq	Veh vs Cort, Dex, or T	HRE	L4–L20	20
Supplemental Table 3:28	Cohesin ChIP-Seq: STAG1 and RAD21	Veh vs E2	ERE	L4–L20	2
Supplemental Table 3:29	Cohesin ChIP-Seq: STAG1 and RAD21	CTCF-knockdown vs CTL	ERE	L4–L20	2
Supplemental Table 3:30	PIAS1 ChIP-Seq	Veh vs R1881	HRE	L4–L20	2
Supplemental Table 3:31	AR ChIP-Seq	siPIAS1 vs CTL	HRE	L4–L20	2
Supplemental Table 3:32	FoxA1 ChIP-Seq	siPIAS1 vs CTL	HRE	L4–L20	2
Supplemental Table 3:33	Sumo2/3 ChIP-Seq	Veh vs R1881	HRE	L4–L20	2
Supplemental Table 3:34	Sumo2/3 ChIP-Seq	GR-K721R vs GR-WT	HRE	L4–L20	2
Supplemental Table 3:35	GR ChIP-Seq	GR-K721R vs GR-WT	HRE	L4–L20	2
Supplemental Table 3:36	GR ChIP-Seq	siHic5 vs CTL	HRE	L4–L20	2
Supplemental Table 3:37	HDAC 2 and 3 ChIP-Seq	Veh vs DHT	HRE	L4–L15	2
Supplemental Table 3:38	EZH2 ChIP-Seq	Veh vs DHT	HRE	L4–L20	1
Supplemental Table 3:39	TOP1 ChIP-Seq	Veh vs DHT	HRE	L4–L10	1
Supplemental Table 3:40	MRE11 ChIP-Seq	Veh vs DHT	HRE	L4–L15	1
Supplemental Table 3:41	Gro-Seq: Veh and DHT	siTOP1 or siMRE11 vs CTL	HRE	B4–B20	8
Hormone-regulated transcription factor recruitment
Supplemental Table 3:42	GATA3 ChIP-Seq	Veh vs E2	ERE	L4–L20	5
Supplemental Table 3:43	GATA3 ChIP-Seq: Trans	GATA3-WT vs GATA3-C318R or GATA3-N320K	ERE	L4–L20	2
Supplemental Table 3:44	ERα ChIP-Seq	siGATA3 vs CTL	ERE	L4–L20	1
Supplemental Table 3:45	AR ChIP-Seq	siGATA2 vs CTL	HRE	L4–L20	2
Supplemental Table 3:46	RARα/γ ChIP-Seq	Veh vs E2	ERE	L4–L20	2
Supplemental Table 3:47	RARα/γ ChIP-Seq: Trans	RARα/γ-WT vs RARα/γ-pBox or RARα/γ-ZF	ERE	L4–L20	4
Supplemental Table 3:48	Gro-Seq: Veh and E2	shRARs vs CTL	ERE	B4–B20	4
Supplemental Table 3:49	p300 ChIP-Seq	shRARs vs CTL	ERE	L4–L20	1
Supplemental Table 3:50	ERα ChIP-Seq	shRARs vs CTL	ERE	L4–L20	1
Supplemental Table 3:51	pERα-S118 ChIP-Seq	shRARs vs CTL	ERE	L4–L20	1
Supplemental Table 3:52	Med1 ChIP-Seq	Veh vs E2	ERE	L4–L20	2
Supplemental Table 3:53	Med1 ChIP-Seq: shCTL and shRARs	Veh vs E2	ERE	L4–L20	2
Supplemental Table 3:54	Gro-Seq	siSMC3 vs CTL	ERE	B4–B20	2
Supplemental Table 3:55	AP2γ ChIP-Seq	Veh vs E2	ERE	L4–L20	1
Supplemental Table 3:56	Gro-Seq: Veh and E2	shAP2γ vs CTL	ERE	B4–B20	4
Supplemental Table 3:57	Stat3 ChIP-Seq	Veh vs Dex	HRE	L4–L20	3
Supplemental Table 3:58	ERG ChIP-Seq	Veh vs DHT	HRE	L4–L20	1
Supplemental Table 3:59	CEBP-β ChIP-Seq and ChIP-Exo	Veh vs Dex or Pred	HRE	L4–L20	4
Supplemental Table 3:60	CREB1 ChIP-Seq	Veh vs Dex	HRE	L4–L20	2
Supplemental Table 3:61	DNase HSS ChIP-Seq	dn-CEBP-β vs CTL	HRE	L4–L20	2
Supplemental Table 3:62	AP1 ChIP-Seq	Veh vs Dex	HRE	L4–L10	1
Supplemental Table 3:63	AP1 ChIP-Seq	Veh vs Compound A	HRE	L4–L15	1
Supplemental Table 3:64	Gro-Seq	Veh vs Dex or Compound A	HRE	B4–B15	2

Summary of genomic features analyzed. See also Supplemental Tables 3:1–64.

Abbreviations: Cort, corticosterone; CTL, control; Dex, dexamethasone; DHT, 5α-dihydrotestosterone; dnCEBP-β, dominant negative-CEBP-β; Fulv, fulvestrant; HSS, hypersensitivity site; Pred, prednisolone; R1881, methyltrienolone; R5020, promegestone; shAP2γ, knockdown of AP2γ; shRAR, knockdown of RARα/γ; siGATA2, knockdown of GATA2; siGATA3, knockdown of GATA3; siHic5, knockdown of Hic5; siMRE11, knockdown of MRE11; siPIAS1, knockdown of PIAS1; siSMC3, knockdown of cohesin subunit SMC3; siTOP1, knockdown of TOP1; T, testosterone; Tam, tamoxifen; Veh, vehicle.

Hormone-regulated eRNA transcription

Enhancer activity can be assessed by measuring the production of eRNAs, bidirectional long noncoding RNAs that are transcribed from sites within the enhancer (33–35). Production of eRNAs has been shown to be an inherent functional attribute of enhancers that is required for the hormone-dependent induction of target genes (36). Using Gro-seq data generated from E2-treated cells, compared with vehicle, we determined that only those sites that we have defined as NRFE ERE sites supported E2-dependent eRNA production [Fig. 4(a); Supplemental Tables 3:1 and 3:2]. Analysis of transcriptional response (GRO-seq) in response to 5α-dihydrotestosterone (DHT) or dexamethasone (dex) treatment revealed that the only eRNA transcripts that were upregulated in response to hormone were initiated at sites that possessed an NRFE HRE [Fig. 4(b); Supplemental Tables 3:3 and 3:4].

Figure 4.

Indicators of enhancer activation are specifically enriched at NRFEs in response to hormone. To facilitate the presentation of data derived from multiple datasets in a single figure, each data point represents the log of the proportion of peaks that contain the genomic feature in the presence of hormone relative to vehicle conditions at B4 (Gro-seq), D1000 (RNAPII), or L4 (p300, FoxA1, DNase hypersensitivity) peak selection criteria. Thus, data points above zero represent individual datasets that show increased presence of the genomic feature at those sites in response to hormone. All null event counts were set to unity to avoid zero divisors. In these analyses, we are trying to determine if there is a general consensus among datasets from various studies on the “parity” of the ratio under consideration, that is, above/below 1 (i.e., ± logarithm). In particular, the magnitude that the logarithm differs from zero in an individual dataset is just “strength of opinion” with regard to our interest in the general trend of the overall parity of the family of results. (a) The proportion Gro-seq reads that contain a 0- to 3-nt variant ERE (and remaining peaks) in the presence of E2 relative to vehicle conditions at B4 peak selection criteria. (b) The proportion Gro-seq reads that contain a 0- to 3-nt variant HRE (and remaining peaks) in the presence of DHT or dex relative to vehicle conditions at B4 peak selection criteria. (c) The proportion of RNAPII ChIP-seq peaks that contain a 0- to 3-nt variant ERE (and remaining peaks) in the presence of E2 relative to vehicle conditions at D1000 peak selection criteria. (d) The proportion of RNAPII ChIP-seq peaks that contain a 0- to 3-nt variant HRE (and remaining peaks) in the presence of cort, dex, or pred relative to vehicle conditions at D1000 peak selection criteria. (e) The proportion of p300 ChIP-seq peaks that contain a 0- to 3-nt variant ERE (and remaining peaks) in the presence of E2 relative to vehicle conditions at L4 peak selection criteria. (f) The proportion of p300 ChIP-seq peaks that contain a 0- to 3-nt variant HRE (and remaining peaks) in the presence of dex, T, P4, or R5020 relative to vehicle conditions at L4 peak selection criteria. (g) The proportion of FoxA1 ChIP-seq peaks that contain a 0- to 3-nt variant ERE (and remaining peaks) in the presence of E2 relative to vehicle conditions at L4 peak selection criteria. (h) The proportion of FoxA1 ChIP-seq peaks that contain a 0- to 3-nt variant HRE (and remaining peaks) in the presence of dex, T, or R1881 relative to vehicle conditions at L4 peak selection criteria. (i) The proportion of DNase hypersensitive sites that contain a 0- to 3-nt variant ERE (and remaining peaks) in the presence of E2 relative to vehicle conditions at L4 peak selection criteria. (j) The proportion of DNase hypersensitive sites that contain a 0- to 3-nt variant HRE (and remaining peaks) in the presence of cort, dex, or T relative to vehicle conditions at L4 peak selection criteria. See also Supplemental Tables 3:1, 3:3, 3:5, 3:12, 3:14, 3:16, 3:18, 3:20, 3:24, and 3:26.

Hormone-regulated RNAPII occupancy

RNAPII occupancy is a feature that unambiguously distinguishes active from inactive TF-bound sites (37). It was of significance therefore that we observed increased RNAPII occupancy at those sites that possessed an NRFE ERE in response to E2 treatment, independent of the peak selection criteria (L or D) parameter [Fig. 4(c); Supplemental Tables 3:5–7]. This E2-mediated increase in RNAPII occupancy at sites that possessed an NRFE ERE did not occur in ERα-KIKO following E2 treatment, consistent with the loss of E2-mediated gene expression in this ERα DBD mutant mouse model (Supplemental Table 3:8). However, an E2-mediated increase in RNAPII occupancy did occur at sites that possessed an NRFE HRE in ERα-KIKO, consistent with its ability to stimulate transcription using HRE motifs (Supplemental Table 3:9). Furthermore, in MCF7 cells, the E2-mediated increase in RNAPII occupancy at sites that possessed an NRFE ERE did not occur following treatment with the antagonist fulvestrant and very slightly with the antagonist tamoxifen (consistent with the different mechanisms of inhibition of these two ERα modulators), whereas chromatin occupancy of ERα was indistinguishable between these ligands (Supplemental Tables 3:10 and 3:11) (38). Correspondingly, analysis of RNAPII occupancy following corticosterone (cort), dex, or prednisolone (pred) showed an increase in RNAPII occupancy only at sites that possessed an NRFE HRE [Fig. 4(d); Supplemental Tables 3:12 and 3:13].

Hormone-regulated recruitment of coregulators and chromatin remodeling proteins

The interaction of the sNR with specific coregulators is the key step in the assembly of transcriptional complexes at enhancers that regulate target gene transcription (39, 40). Thus, we considered that assessment of coregulator recruitment to sNR-occupied sites on chromatin would be a useful way to ascribe activity and evaluated the recruitment of p300 in response to hormone treatment. It was noted that, following treatment with E2, recruitment of p300 only occurred at sites that possessed an NRFE ERE [Fig. 4(e); Supplemental Tables 3:14 and 3:15]. Similarly, treatment with dex, testosterone (T), P4 or promegestone (R5020) only enabled p300 recruitment at sites that possessed an NRFE HRE [Fig. 4(f); Supplemental Tables 3:16 and 3:17].

The TF FoxA1 has been shown to function both as a pioneering factor, by facilitating access of sNRs to functional enhancers by opening up local chromatin during cellular differentiation (i.e., passively enabling sNR access to chromatin), and as a coregulator (41, 42). Our analysis revealed that (1) E2-dependent and (2) dex-, T-, or methyltrienolone (R1881)-dependent recruitment of FoxA1 occurred only at NRFE EREs and NRFE HREs, respectively [Fig. 4(g) and 4(h); Supplemental Tables 3:18–21]. This hormone-mediated NRFE HRE recruitment of FoxA1 was reduced in the AR-C562S DBD mutant and was increased in the constitutively active AR-Q640X variant, consistent with the binding profiles of these AR mutants at NRFE HREs (Supplemental Tables 3:22 and 3:23). Although most sNR binding occurs at preaccessible regions of chromatin, a small number of sites are created de novo subsequent to hormone action (43). Accordingly, E2 treatment resulted in an increase in DNase hypersensitivity only at those sites that possessed an NRFE ERE, and cort, dex, or T treatment increased DNase hypersensitivity only at sites that possessed an NRFE HRE [Fig. 4(i) and 4(j); Supplemental Tables 3:24–27].

In addition to p300 and FoxA1, we used published ChIP-seq datasets to probe for the interaction of other coregulators to sites that possessed an NRFE ERE or NRFE HRE and found cohesin (STAG1 and RAD21) to be recruited to NRFE EREs in response to E2 (Supplemental Table 3:28). This E2-mediated recruitment of cohesin was independent of its relationship with CTCF as hormone-mediated recruitment of cohesin was unaffected by CTCF knockdown (Supplemental Table 3:29). PIAS1, a sumo ligase and known AR coregulator, was recruited in a hormone-dependent manner to NRFE HREs; PIAS1 depletion had a minimal impact on the chromatin occupancy of its interacting partners AR and FoxA1 (Supplemental Tables 3:30–32). Correspondingly, Sumo2/3 was recruited to NRFE HREs in response to hormone; this recruitment was diminished in a SUMOylation-defective GR, GR-K721R, whereas chromatin occupancy of GR-K721R was indistinguishable from wild-type GR (Supplemental Tables 3:33–35). In addition, knockdown of Hic5, a GR coregulator, restricted GR occupancy only at NRFE HREs, and coregulators HDAC2 and 3 and EZH2 were all recruited in a DHT-dependent manner to NRFE HREs (Supplemental Tables 3:36–38). Furthermore, hormone-dependent recruitment of topoisomerase-I (TOP1) has been shown to facilitate enhancer activation by creating single-strand DNA nicks that are postulated to relieve transcriptional-induced torsional stress that results from eRNA synthesis (34). DHT-mediated recruitment of TOP1 and MRE11 occurred only at sites that possessed an NRFE HRE; knockdown of either TOP1 or MRE11 showed selective reduction of eRNA transcription at sites that possessed an NRFE HRE (Supplemental Tables 3:39–41).

Hormone-regulated TF recruitment

Recruitment of a mega-TF complex has been shown to distinguish active, functional enhancers from nonfunctional binding of ER (44). A key component of this complex, GATA3, was only recruited to NRFE EREs in an E2-dependent manner (Supplemental Table 3:42). It has also been shown that GATA3 DBD mutants can be recruited, in trans, to ERα-bound active enhancers (44). Correspondingly, we further demonstrated that these GATA3 DBD mutants were only recruited to sites that possessed an NRFE ERE, confirming that GATA3 was being recruited as a coregulator to NRFE EREs, rather than acting as a TF, binding its cis-regulatory element near NRFE EREs (Supplemental Table 3:43). In agreement, knockdown of GATA3 had no impact on the chromatin occupancy of ERα, and knockdown of GATA2 had no impact on the chromatin occupancy of AR (Supplemental Tables 3:44 and 3:45).

Similar to GATA3, RARα and RARγ have also been shown to be recruited to active, but not inactive, ERα-bound enhancers following E2 treatment (44). Increased RARα and RARγ recruitment in response to E2 treatment was observed only at peaks that contained an NRFE ERE (Supplemental Table 3:46). Four RAR DBD mutants showed increased recruitment, in trans, to sites that possessed an NRFE ERE, confirming that RARα and RARγ were being recruited as coregulators to NRFE EREs, rather than acting as a TF, binding its cis-regulatory element near NRFE EREs (Supplemental Table 3:47). Furthermore, knockdown of RARα and RARγ also inhibited E2-induced eRNA transcription and p300 recruitment only from sites that possessed an NRFE ERE (Supplemental Tables 3:48 and 3:49). Although knockdown of RARα and RARγ had no impact on ERα binding across the genome, it drastically reduced ERα-S118 phosphorylation only at sites that possessed an NRFE ERE (Supplemental Tables 3:50 and 3:51). In addition, E2-induced Med1 recruitment, part of the mediator complex, only occurred at sites that possessed an NRFE ERE; knockdown of RARα and RARγ inhibited this recruitment (Supplemental Tables 3:52 and 3:53). Med1 has also previously been shown to colocalize with cohesin at enhancers and promoters; cohesin contributes to E2-dependent gene activation, at least in part by stabilizing E2/ERα/eRNA-induced enhancer-promoter looping (36). Correspondingly, knockdown of cohesin resulted in a reduction of eRNA transcripts only from sites that possessed an NRFE ERE (Supplemental Table 3:54). Furthermore, E2-mediated AP2γ recruitment only occurred at sites that possessed an NRFE ERE; however, unlike RARα and RARγ, knockdown of AP2γ did not inhibit E2-mediated eRNA transcription (Supplemental Tables 3:55 and 3:56).

Examination of TF recruitment to NRFE HRE motifs revealed that hormone-mediated recruitment of the TFs Stat3 and ERG occurred only at sites that possessed an NRFE HRE (Supplemental Tables 3:57 and 3:58). Dex-mediated recruitment of CEBP-β to GR-bound enhancers has been shown to occur via mechanisms involving chromatin remodeling (45). We found hormone-dependent CEBP-β and CREB1 recruitment only to sites that possessed an NRFE HRE, and introduction of a dominant-negative CEBP-β reduced DNase hypersensitivity only at NRFE HRE sites (Supplemental Tables 3:59–61). Furthermore, the GR-interacting TF, AP1, was only recruited to sites that contained an NRFE HRE in response to dex; this recruitment did not occur following treatment with compound A, a selective GR modulator (Supplemental Tables 3:62 and 3:63). Correspondingly, dex-mediated eRNA transcription was increased at NRFE HRE sites; this increase in eRNA transcription did not occur following treatment with compound A (Supplemental Table 3:64).

Taken together, it is clear that only EREs and HREs, which have the characteristics of the NRFE defined before, are capable of productive associations with the TF/coregulator complexes that are associated with active enhancers. Of note, the components of the TF/coregulator complex recruited to individual NRFE loci are not equivalent; further investigation into the different components at individual NRFE loci will likely reveal the details involved in transcription expression, robustness, and sustainability of hormone-stimulated gene expression between different genes.

Furthermore, the reproducibility between replicate sNR ChIP-seq samples is less than ideal, generally within the range of 50% to 70% (Supplemental Fig. 3). In support of our conclusion, the overall reproducibility is highest at NRFEs and is significantly reduced for the remaining peaks in the datasets. This reduction in the continuity of the remaining peaks between replicates suggests that the remaining peaks largely account for the noise in the information and, hence, the reduced reproducibility between replicates (Supplemental Fig. 3). Moreover, although this study defines the DNA sequence constraints that distinguish active from inactive sNR-DNA interactions, our studies reveal that ERα binding does not correlate with E2-mediated gene expression of the closest TSS. This is independent of (1) peak robustness, (2) what type of motif (NRFE or not) is within the ChIP-seq peak, and (3) the proximity of ERα binding to or within a gene (Supplemental Fig. 4; Supplemental Table 4:1). Therefore, assigning each NRFE to the expression of specific genes requires the identification of the enhancer-promoter interactions, established during cellular differentiation, setting the stage for stimulus-specific transcriptional responses. In support of our observations, analysis of published ER chromatin interaction analysis by paired-end tag sequencing data has previously revealed that the percentage of sites bound by ER with distal chromatin interactions was directly correlated with ERE strength (46).

sNR cross-talk occurs between NRFE binding sites

In the physiologic context, cells are contained in an environment with complex mixtures of hormones, allowing for multiple sNRs to be activated concurrently. Previous studies indicate that P4 functions as a molecular rheostat to control ERα chromatin binding and transcriptional activity (47). In agreement, under full media conditions, P4 or R5020 treatment resulted in increased recruitment of ERα specifically to sites that possessed an NRFE HRE and loss of ERα at NRFE EREs (Fig. 5; Supplemental Tables 3:65 and 3:66). This progestin-mediated ERα redistribution was diminished in PR-deficient cells, demonstrating the critical role for activated PR in determining this progestin-mediated ERα binding (Supplemental Tables 3:67–71). Correspondingly, p300 was also redistributed from NRFE EREs to NRFE HREs following P4 or R5020 treatment (Supplemental Tables 3:72 and 3:73). This suggests PR is recruiting ERα, along with other coregulators, from sites that are undergoing active transcription to NRFE HREs following P4 or R5020 treatment. In agreement, it was shown that the P4/R5020 ERα binding regions correlate with P4/R5020 downregulated genes (47). Similarly, ERα was redistributed from NRFE EREs to NRFE HREs in response to dex treatment, and E2-mediated GR recruitment occurred at NRFE EREs (Supplemental Tables 3:74–77). Furthermore, expressing higher levels of FoxA1 in cells redistributed AR specifically from NRFE HREs (Supplemental Table 3:78) to FoxA1 binding sites (42), and NRFE HREs became more nucleosome rich (Supplemental Table 3:79). The redistribution of AR from NRFE HREs also occurred in an AR mutant that is unable to interact with FoxA1 (Supplemental Table 3:78). Taken together, these studies suggest that TFs can redistribute other TFs and coregulators from sites undergoing active transcription. This model is consistent with the approximate equivalency of gene induction vs gene repression on each chromosome in response to hormone (Supplemental Fig. 5A and 5B). This observation is also supportive of previously observed transcriptional antagonism between sNRs (48). Taken together, these findings suggest that hormone-mediated NRFE binding by sNRs results directly in gene induction; the coregulators that are recruited to those sites are being redistributed from sites that are undergoing active transcription. This redistribution process consequently results in an effective repression of genes driven by those sites. Furthermore, although STARR-seq can report transcriptional activation and repression, it has recently been shown that GR and other TFs work exclusively as transcription activators (i.e., hormone-mediated transcriptional “repression” is not an activation process but rather a deactivation process regressing to the background unactivated state) (31, 49).

Figure 5.

sNR cross-talk occurs between NRFE binding sites. (a–e) The proportion of ERα ChIP-seq peaks that contain a 0- to 3-nt variant ERE (and remaining peaks) in the presence of P4 or R5020 relative to full media conditions at L4 to L20 peak selection criteria. (f–j) The proportion of ERα ChIP-seq peaks that contain a 0- to 3-nt variant HRE (and remaining peaks) in the presence of P4 or R5020 relative to full media conditions at L4 to L20 peak selection criteria. See also Supplemental Tables 3:65 and 3:66.

Discussion

The question of how to distinguish functional and nonfunctional DNA in mammalian genomes has attracted considerable attention for a number of years. Until recently, the study of sNR action was limited to promoters of well-characterized target genes (4). This was due in part to the inherent difficulties in establishing a causal relationship between a sNR-bound enhancer and a specific gene. A common feature of sNRs is a dramatic increase of sNR chromatin interacting sites in response to hormone (Supplemental Table 1:8), revealing an excess in the number of hormone-mediated sNR chromatin interacting events relative to the number of regulated genes. Importantly, we have shown the extent to which sNRs interact with NRFEs and nonfunctional sites appears to be similar. Thus, previous causal relationships established by correlating the levels of a given transcript to the level of sNR recruitment were premature as they do not demonstrate that the binding event noted was causally linked to the transcriptional event. Thus, validation of the functionality of individual sNR binding sites in vivo is essential to understanding transcriptional regulation.

Most studies to date verify the functional activity of a single enhancer on a one-by-one basis. Although genome-wide chromatin immunoprecipitation (ChIP) studies provide a snapshot of sNR occupancy across the genome, these types of analyses are devoid of kinetic information and are dependent not only on the parameters of peak determination but also on the duration of chromatin exposure to the fixative agent (i.e., highly transient interactions will be captured with long exposure times). Thus, although functional binding events are certainly present within the thousands of genome-wide binding events, neutral binding is likely to be more common. Here, we present three sNR DBD mutant mouse models used to define enhancer functionality across the genome. These models revealed that DNA sequence constraints are a key component involved in distinguishing functional from nonfunctional binding events. Specifically, sNR-bound functional enhancers require the presence of a perfect consensus element or variation by only 1 nt or, in a limited number of cases, 2 nts from the consensus, palindromic cis-regulatory element; loss of sNR binding at these sites results in complete loss of all hormone-dependent transcriptional activity. Furthermore, only those sNR binding sites that contain a perfect consensus motif or an element that varies from the consensus by 1 or 2 nts exhibit the biochemical features associated with active enhancers in response to hormone. These findings confirm that it is the physical property of the DNA itself, dictated by sequence content, that guides specific recognition and functional activity of a binding site. Furthermore, these findings also provide the ability to assign functional probabilities to regions of DNA based on the number of variants within the consensus palindromic element, leading to a more accurate and mechanistic interpretation of genome sequences.

In addition to demonstrating that DNA sequence constraints define functionally active sNR binding sites in chromatin, this study also highlights the importance of analyzing the relative number of ChIP-seq peaks or Gro-seq reads that contain a certain motif (rather than the actual number of peaks or reads) when comparing datasets. The actual number of peaks between datasets can be highly variable due to any number of factors: immunoprecipitation efficiency; differences in the amount of starting material, affecting the duplication rate; differences in fragmentation, affecting the length distribution; differences in sequencing quality or read count; and so on. In contrast, by analyzing the relative number of peaks or reads, the impact of these variables is diminished by this normalization methodology.

Chromatin itself functions as an active player in transcriptional regulatory mechanisms, serving not only as a substrate but also as an inducer of specific recognition of cis-regulatory regions required for the regulation of gene expression. Our findings are consistent with previous suggestions that most sNR binding events are nonfunctional (43, 50). This observation raises the question of what is the purpose and/or cause of this large amount of independently nonfunctional binding. Or, assuming a purpose, what information is contained in this structure, and what is this information telling us. This topic is of much interest for future studies.

Conclusion

The accurate and comprehensive identification of functional regulatory sequences in mammalian genomes remains a major challenge for nearly every field of biology. ChIP-sequencing does not provide any information about regulatory function, only protein-DNA coassociation, highlighting the growing need for rational strategies to enable functional and nonfunctional DNA-protein interactions to be distinguished. Our study highlights the importance of DNA sequence constraint in hormone-mediated enhancer activation and begins to reveal the relationship between genome sequence and gene expression through sNR-mediated responses. These data confirm that sNR-binding sites within functional enhancers are significantly less variable than anticipated and call into question the role of the nonfunctional binding of both sNR half-sites and indirect tethering of sNR to enhancers in a DBD-independent manner. These findings are a first step in decoding the gene regulation instructions written in our genome sequence, a critical task performed by cells every day with astonishing reliability and precision (51).

Data Access

The new datasets described herein are available from the NCBI Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/) under accession numbers GSE94737 (ChIP-seq: EAAE-Input, EAAE-Vehicle, EAAE-E2, WT-Vehicle) and GSE94875 (Microarray: E2-mediated gene expression in mouse uteri, Project 584).

Acknowledgment of Data Resources (1332 Datasets)

Figures

• Figure 2(a)–2(e) ; Supplemental Table 2:1 (28) (20) (30) (29) (NEW: EAAE-E2)

• Figure 2(f)–2(j) ; Supplemental Table 2:2 (28) (20) (30) (29) (NEW: EAAE-E2)

• Figure 2(k)–2(o) ; Supplemental Table 2:5 (25)

• Figure 3(a)–3(e) ; Supplemental Table 2:16 (28) (30) (52) (53) (54) (55) (56) (38) (44) (47) (57) (58) (59) (60) (61) (62) (63) (64) (65) (66) (36)

• Figure 3(f)–3(j) ; Supplemental Table 2:17 (29) (67) (68) (69) (25) (53) (45) (70) (71) (72) (47) (73) (74) (75) (76) (77) (78) (79) (80) (34) (81) (82) (83) (42) (84) (85) (86) (87) (88) (89) (90) (30) (25) (27) (91) (61) (92) (93) (94) (95) (96) (97) (98)

• Figure 4(a) ; Supplemental Table 3:1 (99) (100) (101) (36) (35) (102) (44)

• Figure 4(b) ; Supplemental Table 3:3 (103) (34) (91)

• Figure 4(c) ; Supplemental Table 3:5 (28) (38) (54) (55)

• Figure 4(d) ; Supplemental Table 3:12 (67) (30) (25)

• Figure 4(e) ; Supplemental Table 3:14 (55) (44)

• Figure 4(f) ; Supplemental Table 3:16 (67) (104) (47)

• Figure 4(g) ; Supplemental Table 3:18 (54) (55) (44) (61) (57)

• Figure 4(h) ; Supplemental Table 3:20 (61) (105) (81)

• Figure 4(i) ; Supplemental Table 3:24 (106) (61) (53)

• Figure 4(j) ; Supplemental Table 3:26 (67) (45) (104) (53) (61) (107) (108)

• Figure 5(a)–5(e) ; Supplemental Table 3:65 (47)

• Figure 5(f)–5(j) ; Supplemental Table 3:66 (47)

Supplemental figures

• Supplemental Figure 2 (28) (NEW: EAAE-E2)

• Supplemental Figure 3 (56)

• Supplemental Figure 4 (28) (20) (109) (110) (111) (112) (113) (NEW: Project 584)

• Supplemental Figure 5 (28) (20) (109) (110) (111) (112) (113) (NEW: Project 584)

• See Supplemental Table 7 for overlap of eight microarrays (mouse uteri treated with E2 for 2 hours).

Supplemental tables (1:1–8, 2:1–17, 3:1–79, 4:1–3)

• Chromatin Occupancy of sNR and sNR DBD Mutant Mouse Models

• Supplemental Table 1:1 (28) (20) (29) (NEW: EAAE-Vehicle, EAAE-E2, WT-Vehicle)

• Supplemental Table 1:2 (53)

• Supplemental Table 1:3 (25)

• Supplemental Table 1:4 (25)

• Supplemental Table 1:5 (25)

• Supplemental Table 1:6 (25)

• Supplemental Table 1:7 (25)

• Defining Functionally Active ERE and HRE Sites

• Figure 2(a)–2(e) ; Supplemental Table 2:1 (28) (20) (30) (29)

• Figure 2(f)–2(j) ; Supplemental Table 2:2 (28) (20) (30) (29)

• Supplemental Table 2:3 (53)

• Supplemental Table 2:4 (53)

• Figure 2(k)–2(o) ; Supplemental Table 2:5 (25)

• Supplemental Table 2:6 (25)

• Supplemental Table 2:7 (25)

• Supplemental Table 2:8 (25)

• Supplemental Table 2:9 (25)

• Supplemental Table 2:10 (27)

• Supplemental Table 2:11 (42)

• Supplemental Table 2:12 (94) (95)

• Supplemental Table 2:13 (80)

• Supplemental Table 2:14 (42)

• Supplemental Table 2:15 (83)

• KRs Exhibit a Lower Specific DNA-Binding Capacity than ERs

• Figure 3(a)–3(e) ; Supplemental Table 2:16 (28) (30) (52) (53) (54) (55) (56) (38) (44) (47) (57) (58) (59) (60) (61) (62) (63) (64) (65) (66) (36)

• Figure 3(f)–3(j) ; Supplemental Table 2:17 (29) (67) (68) (69) (25) (53) (45) (70) (71) (72) (47) (73) (74) (75) (76) (77) (78) (79) (80) (34) (81) (82) (83) (42) (84) (85) (86) (87) (88) (89) (90) (30) (25) (27) (91) (61) (92) (93) (94) (95) (96) (97) (98)

• Hormone-Regulated eRNA Transcription

• Figure 4(a) ; Supplemental Table 3:1 (99) (100) (101) (36) (35) (102) (44)

• Supplemental Table 3:2 (99) (100) (101) (36) (35) (102) (44)

• Figure 4(b) ; Supplemental Table 3:3 (103) (34) (91)

• Supplemental Table 3:4 (103) (34) (91)

• Hormone-Regulated RNAPII Occupancy

• Figure 4(c) ; Supplemental Table 3:5 (28) (38) (54) (55)

• Supplemental Table 3:6 (28) (38) (54) (55)

• Supplemental Table 3:7 (28) (38) (54) (55)

• Supplemental Table 3:8 (28) (20)

• Supplemental Table 3:9 (28) (20)

• Supplemental Table 3:10 (38)

• Supplemental Table 3:11 (38)

• Figure 4(d) ; Supplemental Table 3:12 (67) (30) (25)

• Supplemental Table 3:13 (67) (30) (25)

• Hormone-Regulated Recruitment of Coregulators and Chromatin Remodeling Proteins

• Figure 4(e) ; Supplemental Table 3:14 (55) (44)

• Supplemental Table 3:15 (55) (44)

• Figure 4(f) ; Supplemental Table 3:16 (67) (104) (47)

• Supplemental Table 3:17 (67) (104) (47)

• Figure 4(g) ; Supplemental Table 3:18 (54) (55) (44) (61) (57)

• Supplemental Table 3:19 (54) (55) (44) (61) (57)

• Figure 4(h) ; Supplemental Table 3:20 (61) (105) (81)

• Supplemental Table 3:21 (61) (105) (81)

• Supplemental Table 3:22 (42)

• Supplemental Table 3:23 (42)

• Figure 4(i) ; Supplemental Table 3:24 (106) (61) (53)

• Supplemental Table 3:25 (106) (61) (53)

• Figure 4(j) ; Supplemental Table 3:26 (67) (45) (104) (53) (61) (107) (108)

• Supplemental Table 3:27 (67) (45) (104) (53) (61) (107) (108)

• Supplemental Table 3:38 (114)

• Supplemental Table 3:29 (114)

• Supplemental Table 3:30 (81)

• Supplemental Table 3:31 (81)

• Supplemental Table 3:32 (81)

• Supplemental Table 3:33 (81)

• Supplemental Table 3:34 (78)

• Supplemental Table 3:35 (78)

• Supplemental Table 3:36 (93)

• Supplemental Table 3:37 (86)

• Supplemental Table 3:38 (86)

• Supplemental Table 3:39 (34)

• Supplemental Table 3:40 (34)

• Supplemental Table 3:41 (34)

• Hormone-Regulated Transcription Factor Recruitment

• Supplemental Table 3:42 (55) (44)

• Supplemental Table 3:43 (44)

• Supplemental Table 3:44 (44)

• Supplemental Table 3:45 (88)

• Supplemental Table 3:46 (44)

• Supplemental Table 3:47 (44)

• Supplemental Table 3:48 (44)

• Supplemental Table 3:49 (44)

• Supplemental Table 3:50 (44)

• Supplemental Table 3:51 (44)

• Supplemental Table 3:52 (44)

• Supplemental Table 3:53 (44)

• Supplemental Table 3:54 (36)

• Supplemental Table 3:55 (44)

• Supplemental Table 3:56 (44)

• Supplemental Table 3:57 (70)

• Supplemental Table 3:58 (86)

• Supplemental Table 3:59 (45) (25)

• Supplemental Table 3:60 (30)

• Supplemental Table 3:61 (45)

• Supplemental Table 3:62 (91)

• Supplemental Table 3:63 (91)

• Supplemental Table 3:64 (91)

• Genomic Feature Association Analysis of DNase Hypersensitivity

• Supplemental Table 4:1 (67) (45) (106) (107) (108)

• sNR Cross-Talk Occurs Between NRFE Binding Sites

• Figure 5(a)–5(e) ; Supplemental Table 3:65 (47)

• Figure 5(f)–5(j) ; Supplemental Table 3:66 (47)

• Supplemental Table 3:67 (65)

• Supplemental Table 3:68 (65)

• Supplemental Table 3:69 (65)

• Supplemental Table 3:70 (65)

• Supplemental Table 3:71 (65)

• Supplemental Table 3:72 (47)

• Supplemental Table 3:73 (47)

• Supplemental Table 3:74 (53)

• Supplemental Table 3:75 (53)

• Supplemental Table 3:76 (53)

• Supplemental Table 3:77 (53)

• Supplemental Table 3:78 (42)

• Supplemental Table 3:79 (42)

• Hormone-Mediated sNR Recruitment to Chromatin

• Supplemental Table 1:8A (28) (53) (54) (55) (56) (38) (47) (115) (57) (116) (114) (30) (44) (60) (61) (62) (64) (65)

• Supplemental Table 1:8B (29) (67) (53) (45) (47) (117) (77) (80) (82) (83) (118) (42) (84) (85) (86) (87) (89) (30) (91) (61) (92) (94) (95) (96) (98)

• Genomic Feature Association Analysis of RNAPII

• Supplemental Table 4:2 (28) (20) (119) (67) (45) (71) (120) (38) (121) (54) (55) (73) (117) (81) (46) (30)

• Supplemental Table 4:3 (28) (20) (119) (67) (45) (71) (120) (38) (121) (54) (55) (73) (117) (81) (46) (30)

Appendix. Antibody Table

Peptide/Protein Target	Antigen Sequence (if Known)	Name of Antibody	Manufacturer, Catalog No	Species Raised in; Monoclonal or Polyclonal	Dilution Used	RRID
C-terminus of ERα of mouse origin		ERα (MC-20)	Santa Cruz Biotechnology, sc-542	Affinity purified rabbit polyclonal antibody	20 μL for 30 μg of chromatin	RRID:AB_631470

Abbreviations:

 
• AR

  androgen receptor

 
• bp

  base pair

 
• ChIP

  chromatin immunoprecipitation

 
• ChIP-exo

  chromatin immunoprecipitation-exonuclease

 
• ChIP-seq

  chromatin immunoprecipitation–sequencing

 
• Cort

  corticosterone

 
• DBD

  DNA-binding domain

 
• dex

  dexamethasone

 
• DHT

  5α-dihydrotestosterone

 
• E2

  17-β estradiol

 
• ER

  estrogen receptor

 
• ERE

  estrogen response element

 
• eRNA

  enhancer RNA

 
• ERα-EAAE

  ERα DBD Mutant

 
• ERα-KIKO

  ERα DBD Mutant

 
• GR

  glucocorticoid receptor

 
• GR-Dim

  GR DBD Mutant

 
• GRO-seq

  global run-on–sequencing

 
• HRE

  hormone response element

 
• kb

  kilobase

 
• KR

  ketosteroid receptor

 
• MR

  mineralocorticoid receptor

 
• NRFE

  Nuclear Receptor Functional Enhancer

 
• nt

  nucleotide

 
• P4

  progesterone

 
• PR

  progesterone receptor

 
• pred

  prednisolone

 
• R1881

  methyltrienolone

 
• R5020

  promegestone

 
• RNAPII

  RNA polymerase II

 
• sNR

  steroid nuclear receptor

 
• T

  tetosterone

 
• TF

  transcription factor

 
• TSS

  transcription start site.

Acknowledgments

The authors thank Drs. Tianyuan Wang, Frank Day, and David Fargo at the National Institute of Environmental Health Sciences (NIEHS) Integrative Bioinformatics and NIEHS computational biology facility for their help with writing script codes and server utilization. We also thank members of the NIEHS microarray core for their expertise and Drs. Franco Demayo, Leping Li, and Ching-yi Chang for critical reading of our study.

Financial Support: Studies in this paper were supported by the Division of Intramural Research of the National Institute of Environmental Health Sciences, National Institutes of Health [Projects 1ZIAES70065 (to K.S.K.)]. The funders had no role in study design.

Author Contributions: S.C.H. performed the ERα-EAAE and gene expression experiments. L.A.C. and A.B.B. designed and performed the data processing and motif analyses. L.A.C. wrote the manuscript. K.S.K. and D.P.D. supervised and funded the research. All authors reviewed the manuscript prior to submission.

Disclosure Summary: The authors have nothing to disclose.