Enhancers reside in a unique epigenetic environment during early zebrafish development
- Lucas J. T. Kaaij1,
- Michal Mokry†2,
- Meng Zhou†3,
- Michael Musheev1,
- Geert Geeven4,
- Adrien S. J. Melquiond4,
- António M. de Jesus Domingues1,
- Wouter de Laat4,
- Christof Niehrs1, 5,
- Andrew D. Smith3 and
- René F. Ketting1Email author
© The Author(s). 2016
Received: 25 May 2016
Accepted: 20 June 2016
Published: 5 July 2016
Enhancers, not promoters, are the most dynamic in their DNA methylation status throughout development and differentiation. Generally speaking, enhancers that are primed to or actually drive gene expression are characterized by relatively low levels of DNA methylation (hypo-methylation), while inactive enhancers display hyper-methylation of the underlying DNA. The direct functional significance of the DNA methylation state of enhancers is, however, unclear for most loci.
In contrast to conventional epigenetic interactions at enhancers, we find that DNA methylation status and enhancer activity during early zebrafish development display very unusual correlation characteristics: hypo-methylation is a unique feature of primed enhancers whereas active enhancers are generally hyper-methylated. The hypo-methylated enhancers that we identify (hypo-enhancers) are enriched close to important transcription factors that act later in development. Interestingly, hypo-enhancers are de-methylated shortly before the midblastula transition and reside in a unique epigenetic environment. Finally, we demonstrate that hypo-enhancers do become active at later developmental stages and that they are physically associated with the transcriptional start site of target genes, irrespective of target gene activity.
We demonstrate that early development in zebrafish embodies a time window characterized by non-canonical DNA methylation–enhancer relationships, including global DNA hypo-methylation of inactive enhancers and DNA hyper-methylation of active enhancers.
KeywordsZebrafish development Enhancers DNA methylation Priming 4C
In recent years, the epigenetic environment in which enhancers reside has been extensively studied. This has been done primarily in cell culture but also in model organisms, including zebrafish, Xenopus, mice, and Drosophila [1–9]. The results of these studies have led to the identification of two main enhancer types: active and primed. Generally, active enhancers are marked by H3K4me1/2/3 and H3K27ac, have low levels of DNA methylation (hypo-methylated), are bound by transcription factors (TFs), are accessible, and are nucleosome depleted. On the other hand, primed enhancers are enriched for H3K4me1 (and sometimes for H3K27me3), contain more DNA methylation, are less accessible, and are often in the vicinity of key developmental genes (reviewed in [10, 11]).
It is largely unclear what the functional importance is of low DNA methylation levels at active enhancers. It has been documented that a subset of TFs preferentially/exclusively bind to hypo-methylated DNA [12–14]. Interestingly though, a recent study showed that the vast majority of protein DNA interactions are DNA methylation independent . Furthermore, enhancers are enriched for 5-hydroxy-methylcytosine (5hmC), suggesting that the hypo-methylation of enhancers is at least partly an actively regulated process [16–18]. A recent study investigating the role of TET proteins in enhancer biology suggested that the hypo-methylated state of enhancers facilitates fast induction of gene expression upon differentiation [17, 18], implying that the low level of DNA methylation has a direct influence on enhancer functionality. Still, if and how DNA methylation directly influences enhancer activity is not clear and, hence, the functional importance of low DNA methylation levels at active enhancers remains to be clarified.
The mechanisms by which enhancers increase the transcriptional output of genes are also not fully understood. It has been shown that enhancer–promoter contact by looping is essential, but not sufficient, for enhancer function since primed enhancers can be physically associated with target genes prior to their expression [19–24]. Intriguingly, the investigation of individual enhancers revealed that, in some cases, the binding of TFs to the enhancer is essential for looping, whereas in other cases the looping was independent of the presence of the enhancer sequence itself [21, 22].
Recently, the epigenetic landscape during zebrafish development has been studied with a focus on the midblastula transition (MBT), which happens approximately 3.5 h post-fertilization (hpf), and enhancers have also been investigated during later development [6, 25–27]. These studies revealed that H3K4me1 together with H3K27ac mark functionally active enhancers . Interestingly, whereas the mouse and human genomes undergo massive DNA de-methylation after fertilization, the DNA methylation landscape in zebrafish is largely stable post-fertilization [28, 29]. Only very few differentially methylated regions (DMRs) develop between MBT and 8 hpf, whereas prior to MBT and during later development more DMRs were detected . As expected, a subset of the identified DMRs overlapped with enhancers, which is similar to what had been previously described in other organisms [5, 27, 28, 30–32].
In this study we aimed to characterize enhancers during zebrafish early development in more detail by integrating a wide set of different genomic data sets. Surprisingly, we observed that active enhancers are generally hyper-methylated, whereas primed enhancers are mostly hypo-methylated in early zebrafish embryos. This is opposite to what has been described in many other model systems. We found that the inactive, hypo-methylated enhancers (from now on referred to as hypo-enhancers) reside in a unique epigenetic environment, characterized by the co-occurrence of H3K4me1 with H3K4me2/3, and sometimes H3K27me3. Furthermore, they are equally accessible compared with the active enhancers on hyper-methylated DNA (referred to as hyper-enhancers). Our investigation of hypo-enhancers throughout zebrafish development revealed that a subset is actually hyper-methylated in gametes and these specific loci show mild 5hmC enrichment prior to MBT. We also present evidence that hypo-enhancers are already marked with low levels of H3K4me1 prior to MBT, that the H3K4me1 mark is present throughout development, and that a subset is active in a tissue-specific manner later in development. Finally, 4C-Seq analyses of five hypo-enhancers reveals that these are physically associated with the transcriptional start site (TSS) of TF genes and remain stably associated with these genes throughout development, irrespective of the transcriptional status of the gene. In conclusion, we demonstrate that early zebrafish development comprises a time-window in which the “classic” correlations between DNA methylation and enhancer activity do not apply, with DNA hypo-methylation being specific for primed enhancers and DNA hyper-methylation found at active enhancers. These results question a direct relationship between the degree of DNA methylation of these loci and their activity in zebrafish.
DNA hyper-methylation associates with enhancer activity during early zebrafish development
To further study global DNA methylation dynamics at active enhancer elements, we identified stage-specific enhancers based on the presence of H3K27ac at specific developmental time points (4 hpf, 8 hpf, and 24 hpf) and probed the DNA methylation status of these loci at these three time points (Additional file 1: Figure S1e). Between 4, 8, and 24 hpf, the average DNA methylation level at these stage-specific active enhancers was largely stable at ~85–90 %, suggesting a minimal effect of enhancer activity on DNA methylation during early zebrafish development.
The global nature of the analyses presented thus far could miss DNA methylation changes at a limited set of specific enhancers. To probe the role of enhancers in shaping the DNA methylation landscape locus-specifically, we checked the amount of DMRs that overlap with an enhancer at 4 or 8 hpf. Again in agreement with recently published data , we only found a small number of DMRs (~0.5 % overlap of a DMR at the two time points) that overlap with an enhancer. Importantly, this stability of DNA methylation at enhancers is specific to early development as enhancers that are hyper-methylated during early development can lose DNA methylation in adult muscle (Additional file 1: Figure S1d). In summary, we detected only limited DNA methylation dynamics at active and inactive enhancer elements during early zebrafish development and found that, surprisingly, active enhancers tend to be hyper-methylated during this developmental phase.
Hypo-enhancers define inactive enhancers close to TFs
The vast majority of enhancers at hypo-methylated loci are not active in any of the embryonic data sets analyzed if we use H3K27ac as a proxy of enhancer activity. To follow up on this observation, we defined two sets of enhancers based on their DNA methylation: hypo-enhancers (<25 % DNA methylation and >5 CpGs, n = 774 at 4 hpf) and hyper-enhancers (>75 % DNA methylation and >5 CpGs, n = 7722 at 4 hpf). Examples of a hyper- and hypo-enhancers and their epigenetic environments are shown in Fig. 1b and c, respectively. To quantify the fraction of active enhancers among the hypo- and hyper-enhancers, we computed the fraction that overlaps with an H3K27ac peak at 4 and 8 hpf. As shown in Fig. 1a, hypo-enhancers show limited overlap with an H3K27ac peak (p < 2.2 × 10−16, hypergeometric distribution); only ~24 % of the hypo-enhancers do so (Fig. 1d). We sought to further support these findings by checking TF occupancy of hyper- and hypo-enhancers. For this, we analyzed publicly available ChIP-Seq data sets of three core pluripotency TFs (Nanog, Oct4 and Sox2) [41, 42] and computed the overlap with the two-enhancer types. Consistent with our H3K27ac analysis, hypo-enhancers also have significantly limited overlap with TF binding sites (p < 2.2 × 10−16, hypergeometric distribution; Fig. 1e).
Finally, we used the GREAT tool  to perform Gene Ontology (GO) term analysis to see if hypo- or hyper-enhancers are associated with specific classes of genes. Intriguingly, whereas hyper-enhancers are not associated with genes enriched for a molecular function, hypo-enhancers are only enriched for molecular function GO terms that point towards TFs (Fig. 1f). This shows that the set we marked as hypo-enhancers defines loci that mostly resemble inactive enhancers close to TFs.
Hypo-enhancers drive reproducible expression in enhancer assays
Previously, Bogdanovic et al.  performed a series of enhancer assays revealing that a subset of putative enhancers identified in their study are able to drive expression in a transgenic setting. Overlapping their validated enhancers with our two enhancer classes revealed that all (12/12) enhancers tested by Bogdanovic et al.  are hyper-enhancers. This overlap confirms that at least some hyper-enhancers identified in our study drive expression in a transgenic setting. To validate the functionality of the identified hypo-enhancer loci as enhancers, we performed transgenic enhancer assays using the ZED system . We cloned five hypo-enhancers close to genes with a tissue-specific expression pattern, of which four gave reproducible results (Additional file 1: Figure S2a–g). Of these, three resembled the expression pattern of the closest gene (dacha, tbx2a, and gsc). These data show that the identified loci, referred to as hypo-enhancers, are also capable of driving tissue-specific gene expression in transgenic settings. We conclude that at least some hypo-enhancers are enhancers.
A subset of hypo-enhancers show 5hmC enrichment prior to MBT
Hypo-enhancers and hyper-enhancers are equally accessible
Hypo-enhancers have a unique epigenetic make up
To further probe epigenetic differences between the two enhancer types, we generated ChIP-Seq profiles for multiple histone modifications at 4 hpf. We started out by computing the enrichment of both H3K4me3 and H3K4me2 over both the hypo- and hyper-enhancers. Unexpectedly, whereas the hyper-enhancers reveal minimal enrichment compared with the background control, the hypo-enhancers reveal significantly stronger enrichment with average enrichment of ~4.5- and ~3-fold for H3K4me2 and H3K4me3, respectively (p < 2.2 × 10−16, non-paired Wilcoxon test; Fig. 3c, d; Additional file 1: Figure S3b). Furthermore, we showed that the enrichment of H3K4me2/3 over hypo-enhancers is a feature of the majority of hypo-enhancers by plotting all read densities in box plots (Additional file 1: Figure S3c–f). We note that, in both cases, the enrichments are still considerably lower than those observed at TSSs (typically more than tenfold; data not shown). Previous work in ESCs showed that a subset of primed enhancers is enriched for the repressive H3K27me3 mark and this type of enhancer is considered to be poised for activation . We found that only the hypo-enhancers show a small enrichment (Fig. 3e), suggesting that only a few hypo-enhancers are poised for activation. In contrast to other work, we do not find any enrichment of Pol II at either the hypo- or the hyper-enhancers (Fig. 3f) .
Finally, we analyzed the epigenetic status of TSSs within 50 kb of either hypo- or hyper-enhancers. This analysis revealed that hyper-enhancers tend to be close to active genes as their proximal TSSs are more enriched for both Pol II and H3K4me3. Furthermore, the hyper-enhancer-associated TSSs are slightly more accessible and nucleosome-poor compared with those associated with hypo-enhancers (Additional file 1: Figure S3g–j). Altogether, these analyses show that hypo- and hyper-enhancers are loci with epigenetically distinct properties in early zebrafish embryos.
Hypo-enhancers are stable throughout zebrafish development
Hypo-enhancers are marked with H3K4me1 prior to MBT
Although the majority of the histone marks are deposited around MBT, low levels of H3K4me3 at TSSs have been observed prior to MBT [25, 26]. This, in combination with the stable nature of H3K4me1 at hypo-enhancers post MBT, prompted us to investigate the H3K4me1 status of the hypo- and hyper-enhancers prior to MBT by performing H3K4me1 ChIP-Seq at 2.5 hpf (256-cell stage). Initial peak calling revealed a low amount of high-confidence peaks but manual investigation of some of the peaks that were called suggested that these are bona fide H3K4me1 peaks (data not shown). Moreover, an H3K4me1 profile and a heat map generated over all refseq TSSs revealed mild enrichment and the well-known bimodality of this histone modification at TSSs (Additional file 1: Figure S5a, b). These results suggest that the chromatin immunoprecipitation (ChIP) procedure worked appropriately but that the amount of H3K4me1 present at this time point is simply very low, as shown for multiple other histone modifications previously [25, 26]. Indeed, conventional immunofluorescence coupled to confocal microscopy for H3K4me1 on 2-hpf embryos revealed very weak staining, if any, whereas at 4 hpf the presence of H3K4me1 is obvious (Additional file 1: Figure S5c).
We also generated H3K4me1 profiles over the identified hypo- and hyper-enhancers (called from 4-hpf embryos). In line with the stability of this mark at hypo-enhancers later in development, the hypo-enhancers display a clear, ~2.5-fold enrichment of H3K4me1 at 2.5 hpf (Additional file 1: Figure S5d). This enrichment is significantly higher than that observed over the hyper-enhancers (p < 2.2 × 10−16, non-paired Wilcoxon test). Furthermore, representation of these data in box plots reveals this is a global trend that includes a large fraction of the hypo-enhancers (Additional file 1: Figure S5e, f). Together, these data show that some hypo-enhancers are pre-marked by H3K4me1 at 2.5 hpf.
Hypo-enhancers overlap with H3K27ac peaks in adult zebrafish
Hypo-enhancers associate with TSSs of target genes independent of their transcriptional activity
Next, we analyzed hypo-enhancer–promoter interactions specifically. As expected for enhancers, for all five tested hypo-enhancers we found one or more 4C interactions close to a promoter. Importantly, we found that the hypo-enhancer interacts with the promoter of a TF in all five cases (Fig. 6c, d; Additional file 1: Figure S6), consistent with the observed GO term enrichment of TFs close to hypo-enhancers (Fig. 1f).
Because we performed the 4C-Seq experiments at two developmental time points and in adult brain, we could also assess the dynamics of the enhancer–promoter interactions throughout development and in adult brain. Focusing only on the hypo-enhancer interactions with a TSS of a TF gene, we detected promoter–enhancer interaction at all three time points in five out of five cases (Fig. 6c, d; Additional file 1: Figure S6).
As described previously [22, 23, 51], we could confirm that promoter–enhancer interactions are already present prior to gene expression (defined as RPKM values <1 in RNA-Seq experiments) in three out of five hypo-enhancer interactions with a TF promoter at 8 hpf (Fig. 6c; Additional file 1: Figure S6b, c). Surprisingly, sox5 and meis2a are not expressed significantly in the brain (RPKM <1) but are still physically associated with a hypo-enhancer in this tissue (Fig. 6d; Additional file 1: Figure S6c). In summary, we show that all tested hypo-enhancers physically associate with the promoter of a TF. These enhancer–promoter interactions are already observed prior to transcriptional activation and are maintained in adult tissue, irrespective of the transcriptional activity of the target gene.
The epigenetic landscape in which enhancers reside has received a lot of attention in recent years. In this study we integrate a wide range of epigenetic data sets and show that active distal regulatory elements are generally hyper-methylated (hyper-enhancers) early in zebrafish development, whereas inactive enhancers are enriched at hypo-methylated DNA (hypo-enhancers). These results are highly unexpected since mammalian systems display opposite global enhancer–DNA methylation correlations [17, 31, 32]. Furthermore, we show that the epigenetic make-up is strikingly different for the hyper- and hypo-enhancers . Even though primed enhancers can be hypo-methylated in human ESCs, primed states generally correlate with hyper-methylation, while active states generally correlate with hypo-methylation [17, 31, 32, 52]. Hence, the opposite correlations we find during early zebrafish development are surprising and represent a unique architecture of enhancers.
Enhancer activity does not influence DNA methylation
We describe that enhancer activity during zebrafish embryogenesis has minimal effect on the underlying DNA methylation state and vice versa. These findings are striking since there is compelling evidence that low levels of DNA methylation at enhancers do have functional importance in mammalian systems [17, 18]. Interestingly, the inverse association between DNA methylation and enhancer activity seems similar to promoter methylation during early Xenopus embryogenesis and mouse spermatogenesis, where some genes that are highly expressed are nevertheless hyper-methylated at the TSSs [53, 54]. These observations show that a temporal uncoupling of DNA methylation and transcriptional repression can take place. The question remains why active enhancers remain hyper-methylated during early zebrafish development. We currently cannot provide a good answer to this question, but it is of great importance to know that DNA methylation dynamics at enhancers in zebrafish can follow different rules to those described thus far in mammalian systems. Mechanistically, insufficient TET activity could be the basis for the lack of de-methylation of active enhancers, as studies on mESCs have shown a direct role for TET enzymes in enhancer demethylation [17, 18]. TET activity is more pronounced in adult zebrafish tissues and, indeed, active enhancers are hypo-methylated in adult tissues . Still, the functional relevance of these differences in DNA methylation behavior between developmental stages remains unresolved.
Hypo-methylation as a potential driver for enhancer priming
The activation of the zygotic genome and the priming of inactive enhancers are crucial for early development. In Drosophila, a single key factor, Zelda, has been identified to bind both active and inactive enhancers and the association of Zelda with inactive genes is crucial for their future activity [55, 56]. However, Zelda is not conserved in vertebrates. We hypothesize that the hypo-enhancers we have defined are in fact primed due to their hypo-methylated state. Interestingly, a similar mechanism has been described for non-methylated CpG islands that gain H3K4me3 due to their low levels of DNA methylation irrespective of the genomic context . Moreover, a study in mESCs lacking Lsh, a factor described to be required for DNA methylation at specific loci, has shown that locus-specific loss of DNA methylation coincides with a gain in H3K4me1 specifically at those regions , showing that low DNA methylation can trigger de novo acquisition of H3K4me1. The identification and characterization of the factors that set the H3K4me1 mark at hypo-enhancers in zebrafish will be of pivotal importance in addressing the hypothesis that enhancer priming is initiated by DNA hypo-methylation.
H3K4me1 is present at hypo-enhancers prior to MBT
In zebrafish, the majority of histone modifications are largely absent prior to MBT, although low levels of H3K4me3,for instance, have been described at these early stages [25, 26]. In this study we present evidence that hypo-enhancers bear low levels of H3K4me1 prior to MBT (at ~2.5 hpf). The significance of these low levels of H3K4me1 at primed hypo-methylated enhancers is currently unclear but they suggest that the marking of primed enhancers is not coupled to transcription since zygotic transcription of the nuclear genome is minimal at this stage . Furthermore, this limited zygotic transcription also implies that parentally inherited factors are responsible for the earliest priming of enhancers. Additional studies are needed to identify the factor(s) responsible for the positioning of H3K4me1 at hypo-methylated enhancers prior to zygotic gene activation.
Hypo-enhancers are stably associated with TF promoters
The mechanism(s) behind enhancer–promoter looping are poorly understood but multiple factors have been shown to play a role in the formation of enhancer loops, including cohesin, CTCF, and also TFs themselves [19, 20]. Our results demonstrate that hypo-enhancers can associate with promoters long before they become active, when they have epigenetic signatures that mark them as primed enhancers. Perhaps more surprisingly, we also found that these enhancers can stay associated with a target promoter long after the gene has been shut down. This situation has also been observed for the HoxA and HoxD loci previously [60, 61]. The mechanism behind this observation remains to be investigated but, in this light, it is interesting that the enhancer region of the SHH (sonic hedgehog) gene still contacts the SHH gene in the physical absence of the enhancer itself . This raises the interesting possibility that hypo-enhancers may interact with target promoters through enhancer-independent mechanisms.
During early zebrafish development, enhancer activity generally has an inverse association with DNA methylation. Furthermore, hypo-methylated enhancers identified early in zebrafish development are active in adult zebrafish and physically associated with target promoters irrespective of their expression.
Bisulfite sequencing library preparations
Total genomic DNA was extracted from carefully timed embryos at indicated developmental stages (2, 4, 8, and 24 hpf). Bisulfite sequencing (BS-Seq) libraries were prepared as described before .
Reads were mapped to the zebrafish genome assembly (Zv9). DMR and HMR calling was done as described elsewhere . DMRs were filtered by asking for at least five significant different CpGs (p < 0.05, Fisher exact test). Additional BS-Seq data sets were downloaded from Potok et al. , Jiang et al. , Stadler et al. , and Wang et al. . In the case where the 4-hpf hypo-enhancers were followed throughout development, only those hypo-enhancers that had more than an average fourfold coverage in all methylomes are displayed and further analyzed. Hypo-enhancers with 0.25 fractional DNA methylation more in the gametes than at 4 hpf were considered methylated in the gametes.
ChIP-Seq and Atac-Seq analysis
ChIP-Seq was performed essentially as described by Bogdanovic et al.  and Atac-Seq was performed as described by Buenrostro et al. . The following antibodies were used in this study: H3K4me1 ChIP grade (abcam), H3K27ac ChIP grade (abcam), and H3K4me2 ChIP grade (abcam). For the 2.5-hpf time point ~8000 embryos were used. The ChIP-Seq libraries were prepared using the Ovation Ultralow Library System V2 Workflow (NuGEN). In brief, the DNA ends were repaired, adapter ligated, and PCR amplified. Size selection of libraries between 250 and 500 bp was done on a LabChip XT instrument (PerkinElmer). These size-selected fractions were pooled and sequenced on a HiSeq 2500 (Illumina). Additional data sets were downloaded from Bogdanovic et al.  and Zhang et al. . Peaks of Oct4, Sox2, and Nanog were downloaded from Leichsenring et al.  and Xu et al. . For the initial enhancers, identification peaks called by Bogdanovic et al. were used. Enhancers were defined as a H3K4me1 peaks not overlapping a H3K4me3 peak or a RefSeq (zv9) or Ensemble (zv9) annotated TSS (±4 kb). For all other data sets used reads were mapped to the zebrafish genome assembly (Zv9) using Bowtie (v0.12.08) using the following settings: -m10 -n2 -l28 --best --strata --chunkmbs 256. In the subsequent analysis multiple reads mapping to the same location and strand have been collapsed to single reads and only uniquely placed reads were used for peak-calling. The Cisgenome v.2 software package was used for peak-calling (false discovery rate <0.05). A combination of custom PERL, Python, R scripts, BedTools, and Cisgenome functions were used for computational data analysis. In brief, all heat maps and profiles presented represent normalized (RPM) read densities and the background was set to 1 in the case of the profiles. HMRs not overlapping a TSS were defined as being more the 4 kb away from a TSS, similar to enhancers. To compute the stability of H3K4me1 throughout development over hypo- and hyper-enhancers, we computed the average read density ±200 bp from the mid-point of the enhancer. We only considered enhancers with >20 RPM combined over the three time points, and in the case where a time point had an average of 0 RPM, a 0.1 pseudo count was added to facilitate LOG transformation. In cases where the epigenetic status of TSSs close to hypo- and hyper-enhancers was evaluated, we only considered TSSs within 50 kb of a hypo- or hyper-enhancer. The Pol II density over TSSs was calculated as average RPM values in 100-bp bins within ±500 bp of the TSS.
4C-Seq was essentially performed as described in van de Werken et al. . In brief, single cell suspensions were made from the zebrafish embryos and brain using TrypLE (Life Technologies). Subsequently, the usual 4C protocol was followed using DPNII as the primary cutter and Csp6I as the secondary cutter. 4C template DNA from multiple (at least three) independent 4C experiments was pooled and subsequently 0.8–2 μg of template DNA was used in the 4C PCR reaction. 4C libraries were sequenced on the illumina HiSeq platform. The sequence of the primers used in the 4C-seq experiments can be found in Additional file 2.
4C-Seq data analysis
For peak-calling in a single 4C experiment, we performed explicit background modeling of the up- and downstream genomic regions. We assumed that, in a completely unstructured chromatin fiber, the contact probability monotonically decreases as a function of the distance to the viewpoint. We modeled this by performing monotonic regression of the 4C signal as a function of the distance to the viewpoint. For this we used the R package isotone, which implements the monotonic regression. We then compared the observed 4C signal to the predicted value from the background model and called the extremes that reach a significance threshold as peaks. For a given threshold q and a distribution F of residuals from the background model, every observation greater than Q3(F) + q × IQR(F), where Q3 is the third quartile of F and IQR(F) the inter-quartile range, is considered significant.
Parameters were as follows for calling a 4C interaction: q between 1.5 and 4 and minimal distance to viewpoint >20 kb.
For computing the overlap between 4C interactions and ChIP-Seq peaks, a non-redundant list of 4C interactions (while merging 4C interactions within 1 kb of each other) was intersected with non-redundant ChIP-Seq lists as indicated. The background overlap was computed by generating semi-random 4C peaks within 500 kb of the originally called 4C peak. These semi-random 4C peaks were thereafter intersected with the non-redundant ChIP-Seq list. The standard deviation for the observed and expected overlap was computed by bootstrapping 10,000 times.
Data sets were downloaded from Collins et al.  and Pauli et al. . Reads were mapped to zebrafish genome assembly version danRer7 (Zv9) by tools in the RMAP package . Gene expression was estimated by mapping reads to the transcriptome and computing RPKM for each refseq gene. Here we mapped the two ends of pair-end data separately, but if both ends of one pair were mapped to the same exon, only one was counted in the RPKM. The mappable length of one gene was adjusted based on mapping deadzones.
Enhancer assays were essentially performed as described by Bessa et al. . In the case of the Gsc and Tbx2a enhancer, the activity was scored in primary injected embryos. In the case of Gsc the specific expression pattern was observed in four out of seven transgenic embryos and the Tbx2a enhancers were observed in three out of five embryos. Dacha and Unxc4.1 enhancer activity was assayed in a transgenic assay. In total two out of two (Dacha), and two out of two (Unxc4.1) transgenic founders showed the presented expression pattern. The sequence of the primers used to clone the enhancer elements can be found in Additional file 3.
DNA was purified from staged embryos by phenol chloroform extraction and subsequently digested to nucleosides using nuclease P1 (Roche), snake venom phosphodiesterase (Worthington), and alkaline phosphatase (Fermentas) and subjected to stable isotope dilution LC-MS/MS analysis as described in detail elsewhere .
We would like to thank N. Kreim, from the IMB Bioinformatics Core facility, and S Boymans for bioinformatic support, C.-T. Han and H. Lukas from the IMB Genomics Core facility for Library preparation and sequencing, and M. Vermunt for providing the ZED vector.
This work was supported by a Marie Curie fellowship 623119 (LJTK) and ERC Grants from the Ideas Program of the European Union Seventh Framework Program (RFK and CN).
Availability of data and materials
The Gene Expression Omnibus accession numbers of the raw data reported in the paper are GSE74617 and GSE74789.
LJTK conceived the project, performed experiments, bioinformatics analysis, and data interpretation, and wrote the manuscript. MMokry, MZ, and AMJD performed bioinformatics analysis. MMusheev performed MS analysis. GG and AM helped with 4C analysis and interpretation. WdL supervised 4C analysis. CN contributed essential reagents and supervised MS analysis. ADS supervised bioinformatics analysis. RFK supervised the project and wrote the manuscript. All authors read and approved the final manuscript.
The authors declare that they have no competing interests.
Ethics approval and consent to participate
All experiments with animals were conducted according to the local regulations and with permission of local animal welfare officers.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Bonn S, Zinzen RP, Girardot C, Gustafson EH, Perez-Gonzalez A, Delhomme N, Ghavi-Helm Y, Wilczynski B, Riddell A, Furlong EE. Tissue-specific analysis of chromatin state identifies temporal signatures of enhancer activity during embryonic development. Nat Genet. 2012;44:148–56.Google Scholar
- Creyghton MP, Cheng AW, Welstead GG, Kooistra T, Carey BW, Steine EJ, Hanna J, Lodato MA, Frampton GM, Sharp PA, et al. Histone H3K27ac separates active from poised enhancers and predicts developmental state. Proc Natl Acad Sci U S A. 2010;107:21931–6.Google Scholar
- Gupta R, Wills A, Ucar D, Baker J. Developmental enhancers are marked independently of zygotic Nodal signals in Xenopus. Dev Biol. 2014;395:38–49.View ArticlePubMedPubMed CentralGoogle Scholar
- Heintzman ND, Hon GC, Hawkins RD, Kheradpour P, Stark A, Harp LF, Ye Z, Lee LK, Stuart RK, Ching CW, et al. Histone modifications at human enhancers reflect global cell-type-specific gene expression. Nature. 2009;459:108–12.Google Scholar
- Hon GC, Rajagopal N, Shen Y, McCleary DF, Yue F, Dang MD, Ren B. Epigenetic memory at embryonic enhancers identified in DNA methylation maps from adult mouse tissues. Nat Genet. 2013;45:1198–206.Google Scholar
- Bogdanovic O, Fernandez-Minan A, Tena JJ, de la Calle-Mustienes E, Hidalgo C, van Kruysbergen I, van Heeringen SJ, Veenstra GJ, Gomez-Skarmeta JL. Dynamics of enhancer chromatin signatures mark the transition from pluripotency to cell specification during embryogenesis. Genome Res. 2012;22:2043–53.Google Scholar
- Xie W, Schultz MD, Lister R, Hou Z, Rajagopal N, Ray P, Whitaker JW, Tian S, Hawkins RD, Leung D, et al. Epigenomic analysis of multilineage differentiation of human embryonic stem cells. Cell. 2013;153:1134–48.Google Scholar
- Shen Y, Yue F, McCleary DF, Ye Z, Edsall L, Kuan S, Wagner U, Dixon J, Lee L, Lobanenkov VV, Ren B. A map of the cis-regulatory sequences in the mouse genome. Nature. 2012;488:116–20.Google Scholar
- Pekowska A, Benoukraf T, Zacarias-Cabeza J, Belhocine M, Koch F, Holota H, Imbert J, Andrau JC, Ferrier P, Spicuglia S. H3K4 tri-methylation provides an epigenetic signature of active enhancers. EMBO J. 2011;30:4198–210.Google Scholar
- Calo E, Wysocka J. Modification of enhancer chromatin: what, how, and why? Mol Cell. 2013;49:825–37.View ArticlePubMedGoogle Scholar
- Spitz F, Furlong EE. Transcription factors: from enhancer binding to developmental control. Nat Rev Genet. 2012;13:613–26.View ArticlePubMedGoogle Scholar
- Wiench M, John S, Baek S, Johnson TA, Sung MH, Escobar T, Simmons CA, Pearce KH, Biddie SC, Sabo PJ, et al. DNA methylation status predicts cell type-specific enhancer activity. EMBO J. 2011;30:3028–39.Google Scholar
- Kim J, Kollhoff A, Bergmann A, Stubbs L. Methylation-sensitive binding of transcription factor YY1 to an insulator sequence within the paternally expressed imprinted gene, Peg3. Hum Mol Genet. 2003;12:233–45.View ArticlePubMedGoogle Scholar
- Perini G, Diolaiti D, Porro A, Della Valle G. In vivo transcriptional regulation of N-Myc target genes is controlled by E-box methylation. Proc Natl Acad Sci U S A. 2005;102:12117–22.View ArticlePubMedPubMed CentralGoogle Scholar
- Domcke S, Bardet AF, Adrian Ginno P, Hartl D, Burger L, Schubeler D. Competition between DNA methylation and transcription factors determines binding of NRF1. Nature. 2015;528:575–9.View ArticlePubMedGoogle Scholar
- Stroud H, Feng S, Morey Kinney S, Pradhan S, Jacobsen SE. 5-Hydroxymethylcytosine is associated with enhancers and gene bodies in human embryonic stem cells. Genome Biol. 2011;12:R54.View ArticlePubMedPubMed CentralGoogle Scholar
- Hon GC, Song CX, Du T, Jin F, Selvaraj S, Lee AY, Yen CA, Ye Z, Mao SQ, Wang BA, et al. 5mC oxidation by Tet2 modulates enhancer activity and timing of transcriptome reprogramming during differentiation. Mol Cell. 2014;56:286–97.Google Scholar
- Lu F, Liu Y, Jiang L, Yamaguchi S, Zhang Y. Role of Tet proteins in enhancer activity and telomere elongation. Genes Dev. 2014;28:2103–19.View ArticlePubMedPubMed CentralGoogle Scholar
- de Laat W, Duboule D. Topology of mammalian developmental enhancers and their regulatory landscapes. Nature. 2013;502:499–506.View ArticlePubMedGoogle Scholar
- Levine M, Cattoglio C, Tjian R. Looping back to leap forward: transcription enters a new era. Cell. 2014;157:13–25.View ArticlePubMedPubMed CentralGoogle Scholar
- Deng W, Lee J, Wang H, Miller J, Reik A, Gregory PD, Dean A, Blobel GA. Controlling long-range genomic interactions at a native locus by targeted tethering of a looping factor. Cell. 2012;149:1233–44.Google Scholar
- Amano T, Sagai T, Tanabe H, Mizushina Y, Nakazawa H, Shiroishi T. Chromosomal dynamics at the Shh locus: limb bud-specific differential regulation of competence and active transcription. Dev Cell. 2009;16:47–57.View ArticlePubMedGoogle Scholar
- Ghavi-Helm Y, Klein FA, Pakozdi T, Ciglar L, Noordermeer D, Huber W, Furlong EE. Enhancer loops appear stable during development and are associated with paused polymerase. Nature. 2014;512:96–100.Google Scholar
- Schoenfelder S, Sugar R, Dimond A, Javierre BM, Armstrong H, Mifsud B, Dimitrova E, Matheson L, Tavares-Cadete F, Furlan-Magaril M, et al. Polycomb repressive complex PRC1 spatially constrains the mouse embryonic stem cell genome. Nat Genet. 2015;47:1179–86.Google Scholar
- Lindeman LC, Andersen IS, Reiner AH, Li N, Aanes H, Ostrup O, Winata C, Mathavan S, Muller F, Alestrom P, Collas P. Prepatterning of developmental gene expression by modified histones before zygotic genome activation. Dev Cell. 2011;21:993–1004.Google Scholar
- Vastenhouw NL, Zhang Y, Woods IG, Imam F, Regev A, Liu XS, Rinn J, Schier AF. Chromatin signature of embryonic pluripotency is established during genome activation. Nature. 2010;464:922–6.Google Scholar
- Lee HJ, Lowdon RF, Maricque B, Zhang B, Stevens M, Li D, Johnson SL, Wang T. Developmental enhancers revealed by extensive DNA methylome maps of zebrafish early embryos. Nat Commun. 2015;6:6315.Google Scholar
- Potok ME, Nix DA, Parnell TJ, Cairns BR. Reprogramming the maternal zebrafish genome after fertilization to match the paternal methylation pattern. Cell. 2013;153:759–72.Google Scholar
- Jiang L, Zhang J, Wang JJ, Wang L, Zhang L, Li G, Yang X, Ma X, Sun X, Cai J, et al. Sperm, but not oocyte, DNA methylome is inherited by zebrafish early embryos. Cell. 2013;153:773–84.Google Scholar
- Kaaij LT, van de Wetering M, Fang F, Decato B, Molaro A, van de Werken HJ, van Es JH, Schuijers J, de Wit E, de Laat W, et al. DNA methylation dynamics during intestinal stem cell differentiation reveals enhancers driving gene expression in the villus. Genome Biol. 2013;14:R50.Google Scholar
- Stadler MB, Murr R, Burger L, Ivanek R, Lienert F, Scholer A, van Nimwegen E, Wirbelauer C, Oakeley EJ, Gaidatzis D, et al. DNA-binding factors shape the mouse methylome at distal regulatory regions. Nature. 2011;480:490–5.Google Scholar
- Ziller MJ, Gu H, Muller F, Donaghey J, Tsai LT, Kohlbacher O, De Jager PL, Rosen ED, Bennett DA, Bernstein BE, et al. Charting a dynamic DNA methylation landscape of the human genome. Nature. 2013;500:477–81.Google Scholar
- Meissner A, Mikkelsen TS, Gu H, Wernig M, Hanna J, Sivachenko A, Zhang X, Bernstein BE, Nusbaum C, Jaffe DB, et al. Genome-scale DNA methylation maps of pluripotent and differentiated cells. Nature. 2008;454:766–70.Google Scholar
- Lister R, Pelizzola M, Dowen RH, Hawkins RD, Hon G, Tonti-Filippini J, Nery JR, Lee L, Ye Z, Ngo QM, et al. Human DNA methylomes at base resolution show widespread epigenomic differences. Nature. 2009;462:315–22.Google Scholar
- Roadmap Epigenomics C, Kundaje A, Meuleman W, Ernst J, Bilenky M, Yen A, Heravi-Moussavi A, Kheradpour P, Zhang Z, Wang J, et al. Integrative analysis of 111 reference human epigenomes. Nature. 2015;518:317-330.Google Scholar
- Heintzman ND, Stuart RK, Hon G, Fu Y, Ching CW, Hawkins RD, Barrera LO, Van Calcar S, Qu C, Ching KA, et al. Distinct and predictive chromatin signatures of transcriptional promoters and enhancers in the human genome. Nat Genet. 2007;39:311–8.Google Scholar
- Rada-Iglesias A, Bajpai R, Swigut T, Brugmann SA, Flynn RA, Wysocka J. A unique chromatin signature uncovers early developmental enhancers in humans. Nature. 2011;470:279–83.View ArticlePubMedGoogle Scholar
- Andersen IS, Reiner AH, Aanes H, Alestrom P, Collas P. Developmental features of DNA methylation during activation of the embryonic zebrafish genome. Genome Biol. 2012;13:R65.View ArticlePubMedPubMed CentralGoogle Scholar
- Wang L, Zhang J, Duan J, Gao X, Zhu W, Lu X, Yang L, Zhang J, Li G, Ci W, et al. Programming and inheritance of parental DNA methylomes in mammals. Cell. 2014;157:979–91.Google Scholar
- Clouaire T, Webb S, Skene P, Illingworth R, Kerr A, Andrews R, Lee JH, Skalnik D, Bird A. Cfp1 integrates both CpG content and gene activity for accurate H3K4me3 deposition in embryonic stem cells. Genes Dev. 2012;26:1714–28.Google Scholar
- Leichsenring M, Maes J, Mossner R, Driever W, Onichtchouk D. Pou5f1 transcription factor controls zygotic gene activation in vertebrates. Science. 2013;341:1005–9.View ArticlePubMedGoogle Scholar
- Xu C, Fan ZP, Muller P, Fogley R, DiBiase A, Trompouki E, Unternaehrer J, Xiong F, Torregroza I, Evans T, et al. Nanog-like regulates endoderm formation through the Mxtx2-Nodal pathway. Dev Cell. 2012;22:625–38.Google Scholar
- McLean CY, Bristor D, Hiller M, Clarke SL, Schaar BT, Lowe CB, Wenger AM, Bejerano G. GREAT improves functional interpretation of cis-regulatory regions. Nat Biotechnol. 2010;28:495–501.Google Scholar
- Bessa J, Tena JJ, de la Calle-Mustienes E, Fernandez-Minan A, Naranjo S, Fernandez A, Naranjo S, Fernandez A, Montoliu L, Akalin A, Lenhard B, Casares F, Gomez-Skarmeta JL. Zebrafish enhancer detection (ZED) vector: a new tool to facilitate transgenesis and the functional analysis of cis-regulatory regions in zebrafish. Dev Dyn. 2009;238:2409–17.Google Scholar
- Tahiliani M, Koh KP, Shen Y, Pastor WA, Bandukwala H, Brudno Y, Agarwal S, Iyer LM, Liu DR, Aravind L, Rao A. Conversion of 5-methylcytosine to 5-hydroxymethylcytosine in mammalian DNA by MLL partner TET1. Science. 2009;324:930–5.Google Scholar
- Almeida RD, Loose M, Sottile V, Matsa E, Denning C, Young L, Johnson AD, Gering M, Ruzov A. 5-hydroxymethyl-cytosine enrichment of non-committed cells is not a universal feature of vertebrate development. Epigenetics. 2012;7:383–9.Google Scholar
- Kriaucionis S, Heintz N. The nuclear DNA base 5-hydroxymethylcytosine is present in Purkinje neurons and the brain. Science. 2009;324:929–30.View ArticlePubMedPubMed CentralGoogle Scholar
- Ito S, D'Alessio AC, Taranova OV, Hong K, Sowers LC, Zhang Y. Role of Tet proteins in 5mC to 5hmC conversion, ES-cell self-renewal and inner cell mass specification. Nature. 2010;466:1129–33.View ArticlePubMedPubMed CentralGoogle Scholar
- Zhang Y, Vastenhouw NL, Feng J, Fu K, Wang C, Ge Y, Pauli A, van Hummelen P, Schier AF, Liu XS. Canonical nucleosome organization at promoters forms during genome activation. Genome Res. 2014;24:260–6.Google Scholar
- Natoli G, Andrau JC. Noncoding transcription at enhancers: general principles and functional models. Annu Rev Genet. 2012;46:1–19.View ArticlePubMedGoogle Scholar
- Eijkelenboom A, Mokry M, de Wit E, Smits LM, Polderman PE, van Triest MH, van Boxtel R, Schulze A, de Laat W, Cuppen E, Burgering BM. Genome-wide analysis of FOXO3 mediated transcription regulation through RNA polymerase II profiling. Mol Syst Biol. 2013;9:638.Google Scholar
- Gifford CA, Ziller MJ, Gu H, Trapnell C, Donaghey J, Tsankov A, Shalek AK, Kelley DR, Shishkin AA, Issner R, et al. Transcriptional and epigenetic dynamics during specification of human embryonic stem cells. Cell. 2013;153:1149–63.Google Scholar
- Bogdanovic O, Long SW, van Heeringen SJ, Brinkman AB, Gomez-Skarmeta JL, Stunnenberg HG, Jones PL, Veenstra GJ. Temporal uncoupling of the DNA methylome and transcriptional repression during embryogenesis. Genome Res. 2011;21:1313–27.Google Scholar
- Hammoud SS, Low DH, Yi C, Carrell DT, Guccione E, Cairns BR. Chromatin and transcription transitions of mammalian adult germline stem cells and spermatogenesis. Cell Stem Cell. 2014;15:239–53.View ArticlePubMedGoogle Scholar
- Foo SM, Sun Y, Lim B, Ziukaite R, O'Brien K, Nien CY, Kirov N, Shvartsman SY, Rushlow CA. Zelda potentiates morphogen activity by increasing chromatin accessibility. Curr Biol. 2014;24:1341–6.Google Scholar
- Lin C, Garruss AS, Luo Z, Guo F, Shilatifard A. The RNA Pol II elongation factor Ell3 marks enhancers in ES cells and primes future gene activation. Cell. 2013;152:144–56.Google Scholar
- Thomson JP, Skene PJ, Selfridge J, Clouaire T, Guy J, Webb S, Kerr AR, Deaton A, Andrews R, James KD, et al. CpG islands influence chromatin structure via the CpG-binding protein Cfp1. Nature. 2010;464:1082–6.Google Scholar
- Yu W, Briones V, Lister R, McIntosh C, Han Y, Lee EY, Ren J, Terashima M, Leighty RM, Ecker JR, Muegge K. CG hypomethylation in Lsh-/- mouse embryonic fibroblasts is associated with de novo H3K4me1 formation and altered cellular plasticity. Proc Natl Acad Sci U S A. 2014;111:5890–5.Google Scholar
- Heyn P, Kircher M, Dahl A, Kelso J, Tomancak P, Kalinka AT, Neugebauer KM. The earliest transcribed zygotic genes are short, newly evolved, and different across species. Cell Rep. 2014;6:285–92.Google Scholar
- Montavon T, Soshnikova N, Mascrez B, Joye E, Thevenet L, Splinter E, de Laat W, Spitz F, Duboule D. A regulatory archipelago controls Hox genes transcription in digits. Cell. 2011;147:1132–45.Google Scholar
- Berlivet S, Paquette D, Dumouchel A, Langlais D, Dostie J, Kmita M. Clustering of tissue-specific sub-TADs accompanies the regulation of HoxA genes in developing limbs. PLoS Genet. 2013;9:e1004018.View ArticlePubMedPubMed CentralGoogle Scholar
- Song Q, Decato B, Hong EE, Zhou M, Fang F, Qu J, Garvin T, Kessler M, Zhou J, Smith AD. A reference methylome database and analysis pipeline to facilitate integrative and comparative epigenomics. PLoS One. 2013;8:e81148.Google Scholar
- Buenrostro JD, Giresi PG, Zaba LC, Chang HY, Greenleaf WJ. Transposition of native chromatin for fast and sensitive epigenomic profiling of open chromatin, DNA-binding proteins and nucleosome position. Nat Methods. 2013;10:1213–8.View ArticlePubMedPubMed CentralGoogle Scholar
- van de Werken HJ, Landan G, Holwerda SJ, Hoichman M, Klous P, Chachik R, Splinter E, Valdes-Quezada C, Oz Y, Bouwman BA, et al. Robust 4C-seq data analysis to screen for regulatory DNA interactions. Nat Methods. 2012;9:969–72.Google Scholar
- Collins JE, White S, Searle SM, Stemple DL. Incorporating RNA-seq data into the zebrafish Ensembl genebuild. Genome Res. 2012;22:2067–78.View ArticlePubMedPubMed CentralGoogle Scholar
- Pauli A, Valen E, Lin MF, Garber M, Vastenhouw NL, Levin JZ, Fan L, Sandelin A, Rinn JL, Regev A, Schier AF. Systematic identification of long noncoding RNAs expressed during zebrafish embryogenesis. Genome Res. 2012;22:577–91.Google Scholar
- Smith AD, Chung WY, Hodges E, Kendall J, Hannon G, Hicks J, Xuan Z, Zhang MQ. Updates to the RMAP short-read mapping software. Bioinformatics. 2009;25:2841–2.Google Scholar
- Welling M, Chen HH, Munoz J, Musheev MU, Kester L, Junker JP, Mischerikow N, Arbab M, Kuijk E, Silberstein L, et al. DAZL regulates Tet1 translation in murine embryonic stem cells. EMBO Rep. 2015;16:791–802.Google Scholar