Inducible super-enhancers are organized based on canonical signal-specific transcription factor binding elements

Published online 19 December 2016 Nucleic Acids Research, 2017, Vol. 45, No. 7 3693–3706 doi: 10.1093/nar/gkw1283 Inducible super-enhancers are organized based on canonical signal-specific transcription factor binding elements D´ ora Bojcsuk1, Gergely Nagy2and Balint L. Balint1,* 1Department of Biochemistry and Molecular Biology, Genomic Medicine and Bioinformatic Core Facility, University of Debrecen, Debrecen 4032, Hungary and 2Department of Biochemistry and Molecular Biology, MTA-DE ‘Lendulet’ Immunogenomics Research Group, University of Debrecen, Debrecen 4032, Hungary Received August 12, 2016; Revised December 07, 2016; Editorial Decision December 08, 2016; Accepted December 13, 2016 ABSTRACT Super-enhancers are established through the interactions of several enhancers and a large number of proteins, including transcription factors and coregulators; however, the formation of these interactions is poorly understood. By re-analysing previously published estrogen receptor alpha (ER␣) ChIP-seq data sets derived from the MCF-7 cell line, we observed that in the absence of stimulation, future super-enhancers are represented by one or a few transcription factor binding event(s) and these extraordinary enhancers possess a response element largely specific to the ER␣dimer. Upon hormonal stimulation, these primary binding sites are surrounded by a large amount of ER␣and the critical components of active enhancers, such as P300 and MED1, and together with neighbouring sites bound by newly recruited ER␣, they generate the functional super-enhancers. To further validate the role of canonical elements in super-enhancer formation, we investigated some additional signaldependent transcription factors, confirming that certain, distinguished binding elements have a general organizer function. These results suggest that certain signal-specific transcription factors guide superenhancer formation upon binding to strong response elements. These findings may reshape the current understanding of how these regulatory units assemble, highlighting the involvement of DNA elements instead of protein–protein interactions. INTRODUCTION Estrogen receptor alpha (ER␣) is a well-studied member of the nuclear receptor (NR) superfamily and is a key hormone-regulated transcription factor (TF) in threequarters of breast cancer cases (1). The activation of ER␣ through 17␤-estradiol (referred to as estradiol or E2) results in the translocation of the receptor from the cytoplasm to the nucleus and strengthens the binding of the receptor and its collaborating factors to regulatory regions of DNA. These events contribute to the altered, particularly increased expression of the specific target genes (2). ER␣ primarily binds the estrogen response element (ERE), initially described as a conserved regulatory element of the promoter of Xenopus and chicken vitellogenin genes in 1986 (3). ERE is typically an inverted repeat of the NR half site (AGGTCA) that is separated by a 3-bp-long spacer (IR3). There are several variations of this element, which also binds the ligand-induced ER dimer (4). The most widely used human cell line to study the behavior of ER␣is MCF-7. MCF-7 is a well-established in vitro model for the investigation of estrogen-dependent biological processes of breast cancer development, as this model is an ER+breast cancer-derived cell line isolated from the pleural effusion of a patient with metastatic breast cancer (5,6). There are several putative ER␣transcription factor binding sites (TFBSs) in the reference human genome, but only a portion of these sites is functionally relevant. Using general bioinformatic motif scans, more than 1 million of these putative sites were identified, while the number of optimal or near-optimal sequence motifs was predicted to be approximately 230 000 according to previous studies (7,8). Among these, only a small fraction, 5105 elements, was bound by ER␣in MCF-7 cells, while approximately 11 000 additional ER␣binding events were identified without the presence of any ERE in the DNA. The discrepancy between the putative ER␣binding sites and those identified in experimental models may reflect the accessibility of the sites. Chromatin accessibility has been suggested to drive, in general, the site selection of TFs and chromatin regions with increased DNase I hypersensitivity were indeed enriched in cell lineand tissue-specific TFBSs (9). The emergence of these accessible sites may be primed by cell line-specific factors. More- *To whom correspondence should be addressed. Tel: +36 52 411717 (Ext. 50015); Fax: +36 52 314 989; Email: [email protected].hu C The Author(s) 2016. Published by Oxford University Press on behalf of Nucleic Acids Research. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact [email protected] 3694 Nucleic Acids Research, 2017, Vol. 45, No. 7 over, priming events at the chromatin level further prepare the conformational context for subsequent protein–DNA and protein–protein binding events that ultimately lead to transcriptional activation. The major collaborative factors of ER␣are Forkhead box protein A1 (FoxA1) and activator protein 2 gamma (AP2␥)(10,11). FoxA1 plays a role as a pioneer factor, and it has been proposed that FoxA1 binds ∼50% of the regions occupied by ER␣, even in the absence of estradiol (12–14). Moreover, FoxA1 is indispensable for any ER␣recruitment (15). AP2␥also joins to ER binding sites and may be involved in stabilizing ER␣-chromatin interactions (10). Recent studies have suggested that enhancers may form a higher-order structure within regulatory units, also called super-enhancers (SEs) (16). Several reports have suggested that the expression of cell type-specific genes is controlled through these markedly active enhancer groups. Thus, relevant genes have been published with respect to different types of cells, with genes defining the nature of a given cell (17). In cancer cells, SEs are associated with key oncogenes, such as c-MYC (16,18) and TAL1 (19); therefore, the genome-wide determination of SE-associated genes would be useful for mapping additional determinative genes that may be involved in tumor pathogenesis (20). Thus, an understanding of how these SEs work would be useful for bridging the gaps in cell differentiation and tumorigenesis. Moreover, an understanding of the mechanisms underlying the manner in which SEs influence the genetic programmes of a cell could bridge the phenotype–genotype gap as the foundation for the missing heritability. According to a current approach of HOMER bioinformatic analysis suite (21), active SEs are defined based on the elevated local sequence read density of such genomic regions where the individual ChIP-seq (chromatin immunoprecipitation followed by high-throughput DNA sequencing) peaks of the key TFs are closer to each other than 12.5 kb. Each ChIP-seq peak or group of peaks is given a read density score, sorted based on score values and those that are located on the region of the graph with a slope >1are considered SEs. Although relevant data have accumulated concerning how these SEs function, little is known regarding the detailed mechanisms that drive the activities of these complexes. SEs are bound by master regulators (17,22), and their prediction can also be done based on the simultaneous presence of BRD4 and/or the active histone mark H3K27ac (23,24). The Mediator complex, which is a key component of transcription initiation, is involved in the looping of enhancers to transcription start sites (25,26). The largest member of the Mediator complex, MED1, binds to several nuclear receptors, such as thyroid hormone receptor (TR), vitamin D receptor (VDR), peroxisome proliferator-activated receptor gamma (PPAR␥), hepatocyte nuclear factor 4 alpha (HNF4␣), ER and glucocorticoid receptor (GR), and it is also a co-regulator of SEs (17,27). Another subunit of the complex, MED14, also interacts with GR, HNF4␣and PPAR␥(27). Previous studies have demonstrated that the reduced level of Mediator members also affects the expression of lineage-specific genes, reflecting the absence of major TFs (28). This effect is consistent with the role of MED1 in forming a bridge between collaborating promoters and enhancers through DNA loop formation (29). In the present study, using previously published and publicly available ChIP-seq data sets, we conducted a detailed examination of how ER␣-driven SEs are formed. We compared the sites with the highest level of ER␣binding that upon estradiol treatment, together with the newly activated regions, were organized into SEs. Unexpectedly, we identified a subset of the enhancers within SE clusters that plays distinct roles in the organization of these functional units, and only these enhancers possess canonical EREs. Moreover, these enhancers have a different profile of protein– protein interactions upon ligand treatment. Further examination of the interdependence of ER␣, FoxA1 and AP2␥ showed that different TFs have specific subsets of SEs, in which the cooperative binding events are incidental. To further validate the role of canonical elements in SE formation, we investigated several additional NRs and TFs, confirming the pivotal role of strong binding sites in SE formation. MATERIALS AND METHODS Data selection We investigated seven TFs in five different human and mouse cell types: ER␣, FoxA1 and AP2␥in the MCF-7 cell line (10,30,31), androgen receptor (AR) in the prostate cancer-derived LNCaP cell line (32), JUNB in primary bone marrow-derived macrophages (BMDM) of the C57BL/6 mouse strain (33), VDR in mouse intestinal epithelial cells and RAR in the F9 mouse embryonic testis carcinoma cell line (34). Control and ligand-treated ChIP-seq samples for the above-mentioned TFs were selected from the publicly available gene expression omnibus database (35) (Supplementary Table S1). We selected seven additional E2-treated ER␣ChIP-seq samples (8 in total, from which at least 10 000 ER␣TFBSs were predicted) for further comparison (Supplementary Table S2) (30,36–40). P300, BRD4, MED1 and H3K27ac ChIP-seq and DNase-seq experiments were also performed as follows (Supplementary Table S3) (41–43). For strengthening our model, several further sets samples of ER␣ChIPseq experiments were included in our analysis that cover: vehicle-treated, untreated samples, vehicleand E2-treated FoxA1 knock-down samples, tamoxifenand fulvestranttreated samples as well as a time-course experiment of E2treatment (Supplementary Table S3) (8,31,44–47). Data processing Raw sequence files of the selected ChIP-seq samples were downloaded and processed using a previously published computational pipeline with the hg19 and mm10 reference genomes (48). ChIP-seq peaks were predicted using MACS2 (49), and SEs were predicted from the ligandtreated samples using HOMER findPeaks (21). Artifacts, based on the blacklisted genomic regions of the Encyclopedia of DNA Elements (50), were removed from the peak sets using BEDTools (51). Reads per kilobase per million mapped reads (RPKM) values for both the control and ligand-treated samples were Nucleic Acids Research, 2017, Vol. 45, No. 7 3695 calculated on the summit ±50 bp region of the peaks determined from the corresponding ligand-treated samples. Peaks with the highest read density (referred to as ‘mother’ peaks) of each future SE were selected from the control samples. The emerging SE peaks (referred to as ‘daughter’ peaks) were identified in the ligand-treated samples in a subsequent analysis. The detailed information about the sites of predicted SEs, mother and daughter enhancers was presented in Additional material (http://dx.doi.org/10.13140/ RG.2.2.22940.49286). Motif enrichment analysis was carried out by findMotifsGenome.pl and was performed on the −50, 50 bp flanking region of the peak summits. The individual and average read density of the different peak sets was determined by annotatePeaks.pl (21). Read distribution and average density heat maps were displayed by Java TreeView (52). Histograms and box plots were performed using GraphPad Prism version 6.00 for Windows (GraphPad Software, La Jolla, CA, USA, www.graphpad.com). RESULTS Putative ER␣binding sites in MCF-7 cells To determine the ER␣cistrome, we collected public ChIPseq data sets that included ligand-treated samples with the highest IP efficiency available. By re-analysing the top eight estradiol-treated MCF-7 samples derived from six laboratories (Supplementary Table S2), we identified a surprisingly large number of ER␣binding sites (∼189 000 in total), although the consistency of these sites was rather low; specifically, in at least two samples, only 56 931 ER␣TFBSs could be predicted (Supplementary Figure S1A and B). Among these, only 4387 TFBSs were present in all of the investigated samples, and we named these sites the ‘consensus peak set’ (the consensus ER␣TFBSs were listed in Additional material). The observed heterogeneity may reflect genetic alterations, the number of cell passages, the estradiol concentration used in the specific experiments, the induction time of the estradiol response or other technical differences, e.g. the number of cells or the protocol used for ChIP. Measuring DNase I hypersensitivity to identify chromatin regions accessible for TFs revealed the consensus peak set with the highest chromatin accessibility (Supplementary Figure S1C). While the majority of the detected ER␣binding sites are variegated, ER␣SEs showed a high degree of consistency within the reported experiments, highlighting the importance of these regulatory regions within the genome. Notably, more than 90% of the SEs contained at least one peak from the consensus set (Supplementary Table S2), and, interestingly, these peaks were typically present also in the vehicle-treated samples (Figure 1A). Primary enhancers designate future super-enhancers A comparison of the vehicleand estradiol-treated samples (GSM614610 and GSM614611 ER␣ChIP-seq samples were used during the main analysis) revealed that in the absence of ligand, each SE was represented by one or a few peak(s) and ligand-induced SEs formed around these activated genomic regions. In the selected sample, 392 highly covered regulatory regions could be predicted by the SE finder of HOMER (Supplementary Figure S1D). Figure 1A concentrates on the new binding events but cannot emphasize the increase on pre-occupied regions; the exact coverage of the examples is shown on Supplementary Figure S1E. Notably, these individual peaks could not be predicted as SEs in the vehicle-treated samples, as the currently used prediction criteria rely mainly on the presence of clusters of enhancers. These peaks are individual or a few separate peaks, but upon ER␣activation, a cluster of peaks assembles and is clearly identified as SEs. Based on this observation, we considered that these preferential primary enhancers may play a key role in the formation of SEs, such as regulatory units; therefore, we further investigated these genomic regions. We examined the size distribution of ER␣peaks (4042 peaks were investigated), plotting each member of ER␣SEs based on their tag density (in deciles) (Figure 1B and C). As a control, we used other ER␣TFBSs that do not cluster in SEs (75 693 peaks) (Figure 1D and E). The top 10% of SE peaks had an extremely high tag density compared with the other sites in both vehicleand estradiol-treated samples (Figure 1B and C). Although the distribution of the tag density was similar under both conditions, a more than 5-fold increase in the median tag density was observed upon ligand treatment. To exclude technical issues as the cause of this observation, we performed the same analysis using another setofER␣ChIP-seq samples and observed the same results (Supplementary Figure S1F and G). Plotting the ER␣enrichment of all ER␣-bound single enhancers that are not members of a SE revealed one-tenth of the tag density of these sites compared with the peaks within SEs. These results suggest that the binding events outside of SEs are less pronounced (Figure 1D and E). According to the assumption that primary enhancers form the basis of SEs, we referred to these elements as ‘mother enhancers’, and the subsequently appearing secondary enhancers were referred to as ‘daughter enhancers’. The determination of mother and daughter enhancers is described in the data processing description in the Materials and Methods section. To globally determine whether the binding of ER␣to the so-called mother enhancers precedes binding to the subsequent so-called daughter enhancers within the SEs, we sorted all ER␣SEs based on the RPKM values of the mother enhancers (calculated from the vehicle-treated samples). In Figure 1F, mother enhancers were well separated, and a relatively low coverage of the daughter enhancers was observed. However, upon ligand treatment, the existence of enhancers was detected on the majority of ER␣SEs (Figure 1F and G, Supplementary Figure S2A and B). We also plotted the average tag density of the ER␣SE peaks, indicating that mother enhancers show even higher signals in the vehicle-treated samples than daughter enhancers upon estradiol treatment (Figure 1H and I, Supplementary Figure S3A and B). To validate the precedence of mother enhancer occupancy, we compared vehicle-treated and untreated ER␣ ChIP-seq samples (Supplementary Figure S4A). These data suggest that vehicle treatment of MCF-7 cells changing the ER␣binding neither on mother nor on daughter enhancer 3696 Nucleic Acids Research, 2017, Vol. 45, No. 7 Figure 1. MCF-7-specific estrogen receptor alpha (ER␣) super-enhancers are represented by one (or a few) enhancer(s) even in the absence of estradiol treatment. (A) Integrative Genomics Viewer (IGV) snapshot of ER␣ChIP-seq coverage representing six ER␣super-enhancers (SEs) upon vehicle and estradiol (E2) treatment. The interval scale is 50 in both cases. (Band C)ER␣tag density upon vehicle or E2 treatment in the deciles determined based on ER␣recruitment at peaks within SEs. (Dand E)ER␣tag density upon vehicle or E2 treatment in deciles determined based on ER␣recruitment at peaks not overlapping with SEs. The boxes represent the first and third quartiles, the horizontal lines indicate the median reads per kilobase per million mapped reads (RPKM) values, and the whiskers indicate the lower and upper extremes per decile. (Fand G)ER␣tag densities of the mother (M) and the top 6 daughter enhancers (1–6) within the same super-enhancer region in vehicleand E2-treated MCF-7 cells. Enhancers were vertically sorted based on the RPKM values of the mother enhancers (in the first column), and the individual enhancers within a SE region were subsequently horizontally aligned based on the read enrichment of the vehicle-treated samples. (Hand I) Histograms show the average tag density of mother (392) and daughter enhancers (3650) in the presence or absence of E2. (Jand K) Histograms show the average tag density of mother (392) and daughter enhancers (3650) in the presence vehicleor E2-treatment and upon tamoxifenand fulvestrant-treatment. Samples derived from the same experiment. Nucleic Acids Research, 2017, Vol. 45, No. 7 3697 regions. The processed time-course experiment also suggested that mother enhancers are indeed occupied prior to daughter enhancers, and show a prolonged binding at later time points, too (Supplementary Figure S4B–D). Thus, there is a difference between the dynamics of the occupancy of the distinct SE peak sets. Importantly, the presence of ER␣on mother enhancers is remarkable after ER␣antagonist treatments such as tamoxifen and fulvestrant (Figure 1J and K). Moreover, in FoxA1-depleted state, as shown on Supplementary Figure S4E–G, the binding of ER␣is reduced, but not abolished neither on mother nor on daughter enhancers. As MED1 was reported to be a key component of the Mediator complex bridging SEs with transcription start sites, we visualized the presence of this component within the ER␣SEs. MED1 is recruited to ER␣-bound sites with a high binding affinity to the buds or initiator(s) of SEs, namely, mother enhancers, upon estradiol treatment (Figure 2A). We also investigated the MED1 densities based on the previously defined and sorted deciles of ER␣peaks, and as shown in Figure 2B and C, Mediator is preferentially recruited to ER␣peaks with higher tag densities, particularly in estradiol-treated samples (Figure 2C). Collectively, we observed collaborating enhancers clustered in SEs, showing a much higher ER␣occupancy than single enhancers, and even more prominent enrichments accompanied by Mediator binding were observed within SEs. The read distribution plot around the individual enhancers demonstrated that mother enhancers are located in the most accessible chromatin regions with high levels of DNase I signal and show the highest MED1, P300, H3K27ac and BRD4 coverage upon induction, suggesting that the top ER␣enhancers (392 mother enhancers in total) represent the most active regulatory regions, while the daughter enhancers (3650 in total) did not show this pattern (Figure 2D, Supplementary Table S3). Based on the previous observation, namely, that each super-enhancer region is indicated by a markedly active region, we assessed whether mother and daughter enhancers were observed for other TFs and cell types. First, we examined the other two lineage-determining factors associated with MCF-7 cells, namely, FoxA1 and AP2␥(Supplementary Figure S2C–F). Second, we extended these investigations to other NRs in different cell types (AR in LNCaP cells, RAR in mouse F9 cells and VDR in mouse intestinal epithelial cells) and to JUNB in mouse BMDMs (Supplementary Figure S2G–N). In all experimental systems investigated (except for FoxA1, as described later in this manuscript), patterns similar to those for ER␣in MCF7 cells were observed, indicating that a more general phenomenon occurs in the presence of several TFs (Supplementary Figure S2A–N, Figure S3A–N). Canonical elements drive transcription factors To better understand the factors that discriminate ER␣ mother and daughter enhancers, we applied a motif enrichment analysis for these regions. While mother enhancers showed strong canonical ERE motif enrichment from a relatively low number of target sequences, daughter enhancers were not enriched for the ER␣dimer-specific elements but were enriched for NR half sites and FoxA1 motifs (Figure 3A and B and Supplementary Figure S5A and B). Interestingly, at the sites of mother enhancers, a unique NR direct repeat element with overlapping half sites was also identified. Considering the strong P-value (1e-200) and the high enrichment of ERE (62.2%) (compared with the background, which was 3.37%), despite the small number of target sequences (n =328), we concluded that the high level of ER␣recruitment at mother enhancers (even in the absence of estradiol) reflects strong canonical DNA elements. However, for the emergence of daughter peaks, other TFs, such as FoxA1 and activator protein 1 (AP-1), act in concert with the increased level of ER␣after estradiol treatment. To validate these observations, we examined the motif dependence of binding frequency. We plotted the motif scores of the remapped EREs as a function of the deciles derived from the vehicle-treated peak size distribution in the case of mother enhancers and from the estradiol-treated peak size distribution in the case of daughter enhancers and observed that in all cases, the extremely high mother peaks had significantly stronger binding elements than the daughter enhancers (deciles 8–10 have P<0.0001) (Figure 3C). This finding suggests that the strongest elements of accessible euchromatic regions irresistibly attract the ER␣dimers. To determine whether this effect is also observed for the binding of other TFs, we investigated these correlations for other types of SEs. Motif enrichment analysis at the sites of the mother and daughter enhancers predicted from other investigated samples showed the same phenomenon. For example, we detected better motif enrichment of mother enhancers in MCF-7 for FoxA1 than that of their daughters, and this effect was also observed for AP2␥(Supplementary Figure S5C–F). Interestingly, AP2␥daughter enhancers harboured the motifs of the main collaborating partners, FoxA1, ER␣and AP-1, together with the CCCTC-binding factor (CTCF). In the case of AR in the LNCaP cell line, FoxA1 is needed for the establishment of mother enhancers, as shown by the extremely high motif enrichment in the mother enhancers (Supplementary Figure S5G and H). Importantly, daughter enhancers showed lower IR3/ARE (ARE: androgen response element) and FoxA1 motif enrichment, but a composite element of FoxA1 and AR was also identified, indicating a tight collaboration between these TFs, which may be less affected by other regulators. In the case of JUNB, mother enhancers of mouse BMDM cells showed all the main specific motif enrichments, namely, AP-1, AP-1::IRF (IRF: interferon regulatory factor) composite element (AICE), cAMP response element and the motif of the lineage-determining CCAAT/enhancer-binding protein (C/EBP). These motifs showed weaker enrichment in daughter enhancers, and further motifs of collaborating TFs, namely, purine-rich box-1 (PU.1), nuclear factor kappa-light-chain-enhancer of activated B cells (NF-␬B) and signal transducer and activator of transcription (STAT) appeared (Supplementary Figure S5I and J). The RAR and VDR mother enhancers of mouse F9 and epithelial cells, respectively, also showed stronger specific motif enrichment than the corresponding daughter peaks, and GATA protein(s) may also occupy VDR mother enhancers. However, daughter enhancers also showed pu- 3698 Nucleic Acids Research, 2017, Vol. 45, No. 7 Figure 2. MED1 is recruited to ER␣-bound sites with an extremely high binding affinity. (A) IGV snapshot of MED1 and ER␣ChIP-seq coverage, representing six ER␣SEs upon vehicle and estradiol (E2) treatment. The interval scale is 50 in both cases. (Band C) MED1 tag density upon vehicle or E2 treatment in deciles determined based on ER␣recruitment at peaks within SEs. The boxes represent the first and third quartiles, horizontal lines indicate the median RPKM values, and whiskers indicate the lower and upper extremes per decile. Paired t-test, * significant at P<0.05; ** at P<0.01; *** at P <0.001; **** at P<0.0001. (D) Read distribution plot of ER␣, DNase I, MED1, P300, H3K27ac and BRD4 upon vehicle or E2 treatment, relative to SE peaks in 2-kb frames. Mother (392 in total) and daughter (3650 in total) peaks are sorted according to ER␣tag density. tative AP-1 and Kr¨uppel-like factor (KLF) binding sites (Supplementary Figure S5K–N). All examined TFs playing roles in SE formation showed a phenomenon similar to that of ER␣, namely, these primary regulatory regions possessed canonical elements specific to dominant TF(s), while the further occupied regions had fewer specific elements together with their collaborative factors. The binding of canonical elements precedes, and upon treatment, likely facilitates the occupation of the nearby regions by dominant TFs and interacting partners. In summary, the primary ChIP-seq peaks, which were present prior to estradiol stimulation, had significantly stronger binding elements than the activated peaks, suggesting that certain elements have high DNA–protein interaction affinities and that there is no need for cooperative binding with other factors, e.g. FoxA1 or AP2␥. Therefore, we examined the differentially enriched ER␣and FoxA1 binding on ER␣SEs. As shown in Figure 3D and E, the occurrence of ER␣on the mother enhancers within ER␣ SEs is a predominant event (P<0.0001), but interestingly, Nucleic Acids Research, 2017, Vol. 45, No. 7 3699 Figure 3. Canonical elements provide higher DNA-binding affinity than non-canonical elements. (Aand B) Motif enrichments of mother and daughter enhancers. The P-value and target and background (Bg) percentages are included for each motif. (C) Estrogen response element (ERE) motif scores in the deciles determined based on ER␣recruitment upon vehicle (in the case of mother enhancers) or E2 treatment (in the case of daughter enhancers) within SEs. (Dand E)RPKMvaluesofER␣and FoxA1 coverage at ER␣mother and daughter enhancers. In the case of panels C, D and E, the boxes represent the first and third quartiles, the horizontal lines indicate the median RPKM values, and the whiskers indicate the 10th to 90th percentile ranges per decile. Paired t-test, * significant at P<0.05, ** at P<0.01, *** at P<0.001, **** at P<0.0001. FoxA1 also shows a similar pattern (P<0.0001). The results suggested that the binding affinity of ER␣was definitely higher at the sites of mother enhancers, and FoxA1 showed a stronger binding at these sites despite the lack of its binding elements that might mean the importance of protein–protein interactions, namely the binding to ER␣. ER␣, FoxA1 and AP2␥form distinct SEs Previous studies have suggested that interactions with collaborating pioneer factors are necessary for ER␣function; therefore, we conducted a detailed investigation of how ER␣, FoxA1 and AP2␥recruitment correlates with each other at the SE regions and examined the extent of the overlap between the SEs of these MCF-7-specific TFs. For this purpose, we plotted the FoxA1 and AP2␥tag density of mother and daughter ER␣enhancers (Figure 4A), and the resulting read distribution heat map showed a positive correlation between the three TFs following estradiol treatment, but the recruitment of the collaborating TFs was much lower than that of ER␣. Examining the mother and daughter FoxA1 peaks revealed an increase in AP2␥binding but not in FoxA1 and ER␣binding (Figure 4B), and we consistently obtained a similar result for AP2␥SE peaks (Figure 4C). Measuring the tag density of the three TFs on their SEs showed that ER␣is typically much less frequently recruited at the mother (and daughter) peaks of the other two TFs (Figure 4D and E) compared with the recruitment of FoxA1 and AP2␥at the sites of ER␣SEs. This effect was also observed for the recruitment of FoxA1 and AP2␥at each other’s SEs (Figure 4F–I), suggesting that SEs are typically dominated by a single TF and that other TFs are likely collaborating partners in their activity. TF specificity was also confirmed by the overlap of SEs determined for ER␣, FoxA1 and AP2␥(Figure 4J). The Venn diagram shows that of these three proteins, ER␣and FoxA1 have a more distinct SE profile, while AP2␥has its own, less prominent, SE profile. This finding is consistent with the coverage results shown in Figure 4A–C, as AP2␥was induced at all SEs examined. The observed dominances may reflect the canonical elements because not only the ER␣-bound mother enhancers but also those of FoxA1 and AP2␥showed stronger specific motif enrichment compared with that in the corresponding secondary regions (Figure 4A and B, Supplementary Figure S5A–F). On FoxA1and AP2␥-driven SEs FoxA1 showed no notable change upon E2 treatment, but interestingly, it did show higher density on ER␣-dominated SEs (Figure 4A–C). The former is not surprising as no significant change could previously been observed in FoxA1 binding upon E2 stimulation (46), but the latter is striking as we found a set of genomic regions, where E2 could affect FoxA1 binding. Focusing on some ER␣-driven SEs, we observed that there is typically at least one FoxA1-bound enhancer per SE in the absence of estradiol; however, these enhancers rather 3700 Nucleic Acids Research, 2017, Vol. 45, No. 7 Figure 4. ER␣, FoxA1 and AP2␥form distinct super-enhancers. (A–C) Read distribution plot of ER␣, FoxA1 and AP2␥upon vehicle or E2 treatment, relative to each other’s SE peaks in 2-kb frames. The number of mother and daughter peaks, sorted based on tag density, is indicated. (D–I)ER␣, FoxA1 and AP2␥tag density of each other’s mother and daughter enhancers in E2-treated MCF-7 cells. (J) Area-proportional Venn diagram illustrates the overlaps between the ER␣, FoxA1 and AP2␥SEs. occasionally overlap with mother enhancers (Figure 5and Supplementary Figure S6A and B). In contrast, AP2␥did not show enrichment prior to treatment and emerged together with ER␣, typically at mother enhancers and specific daughter enhancers occupied by FoxA1 (Figure 5and Supplementary Figure S6C and D). There is a similar enrichment at ER␣mother enhancers as of AP2␥, which suggests an ERE/ER␣-dependent FoxA1 and AP2␥recruitment by protein–protein interactions. Importantly, MED1 is specific only to mother enhancers, while P300 may also be recruited to regions bound by FoxA1. H3K27ac and BRD4, together with the DNA accessibility, also followed these tendencies, thus showing elevated activity upon treatment (Figures 5, 6A, D, G, J, M). By examining the presence of DNase I, MED1, P300, H3K27ac, BRD4 and collaborating TFs on FoxA1and AP2␥-driven SEs, we observed a striking phenomenon: as ER␣recruitment is neither induced at these sites, P300, H3K27ac and BRD4 nor show enrichment, while DNase I signal and MED1 binding is induced upon treatment (Figure 6B and C, E and F, H and I, K and L, N and O, Supplementary Figure S7A and B). Notably, according to the enrichment of H3K27ac and BRD4, the slight decrease in P300 recruitment did not affect these regions, as they were already at least as active in the untreated state as the induced ER␣SEs, and this activation did not change (Figure 6G–I). These results strongly suggest ER␣-dependent MED1/P300 recruitment, followed by histone acetylation Nucleic Acids Research, 2017, Vol. 45, No. 7 3701 Figure 5. Recruitment of FoxA1 and AP2␥at ER␣super-enhancers. IGV snapshot of ER␣ChIP-seq coverage, representing eight ER␣SEs upon vehicle (veh) or E2 treatment and the simultaneous presence of FoxA1, AP2␥, MED1, P300, H3K27ac, BRD4 and DNase I upon vehicle and E2 treatment. The interval scales are indicated in the upper left corners. Peaks, highlighted in grey, represent the sites of ER␣mother enhancers, while dashed lines indicate the location of FoxA1 peaks. and BRD4 binding (Figure 6A, D, G, J, M), while at FoxA1and AP2␥-dominated SEs, only MED1 recruitment and elevated chromatin accessibility were observed upon estradiol treatment (Figure 6E and F). On FoxA1and AP2␥- dominated SEs estradiol treatment results in MED1 binding. However, because P300 is not further recruited, neither histone acetylation nor the reading module of the acetylated histones, namely BRD4 will be further recruited to these sites (Figure 6E, F, H, I, K, L, N, O and Supplementary Figure S7A and B). These events suggested that in the case of TF-specific SE recruitment and the increase in MED1 binding, an increase in DNase I hypersensitivity is not sufficient for the initiation of these events. However, it is important to note that AP2␥-dominated SEs showed the most enriched active enhancer marks (Figure 6C, F, I, L, O), while those of FoxA1 were the least active (Figure 6B, E, H, K, N). Moreover, ER␣bound in a scattered manner to FoxA1and AP2␥-specific SEs (Supplementary Figure S6E–L). DISCUSSION In this study we first describe ligand-inducible superenhancers. Due to the almost prompt biological response (Supplementary Figure S4A–D) and the large variety of