Transcriptional and epigenetic signatures of zygotic genome activation during early drosophila embryogenesis
© Darbo et al.; licensee BioMed Central Ltd. 2013
Received: 22 June 2012
Accepted: 28 February 2013
Published: 5 April 2013
In all Metazoa, transcription is inactive during the first mitotic cycles after fertilisation. In Drosophila melanogaster, Zygotic Genome Activation (ZGA) occurs in two waves, starting respectively at mitotic cycles 8 (approximately 60 genes) and 14 (over a thousand genes). The regulatory mechanisms underlying these drastic transcriptional changes remain largely unknown.
We developed an original gene clustering method based on discretized transition profiles, and applied it to datasets from three landmark early embryonic transcriptome studies. We identified 417 genes significantly up-regulated during ZGA. De novo motif discovery returned nine motifs over-represented in their non-coding sequences (upstream, introns, UTR), three of which correspond to previously known transcription factors: Zelda, Tramtrack and Trithorax-like (Trl). The nine discovered motifs were combined to scan ZGA-associated regions and predict about 1300 putative cis-regulatory modules. The fact that Trl is known to act as chromatin remodelling factor suggests that epigenetic regulation might play an important role in zygotic genome activation. We thus systematically compared the locations of predicted CRMs with ChIP-seq profiles for various transcription factors, 38 epigenetic marks from ModENCODE, and DNAse1 accessibility profiles. This analysis highlighted a strong and specific enrichment of predicted ZGA-associated CRMs for Zelda, CBP, Trl binding sites, as well as for histone marks associated with active enhancers (H3K4me1) and for open chromatin regions.
Based on the results of our computational analyses, we suggest a temporal model explaining the onset of zygotic genome activation by the combined action of transcription factors and epigenetic signals. Although this study is mainly based on the analysis of publicly available transcriptome and ChiP-seq datasets, the resulting model suggests novel mechanisms that underly the coordinated activation of several hundreds genes at a precise time point during embryonic development.
KeywordsDrosophila Melanogaster Zygotic Genome Activation Transcriptional Regulation Epigenetic Regulation Transcriptome ChIP-seq
After fertilisation, Drosophila melanogaster embryos undergo a series of 13 fast mitotic divisions without cytokinesis (thus forming a syncytium, i.e. a single cell with multiple nuclei). The first seven mitotic cycles are fast (8 min/cycle) and synchronous, while the zygotic genome remains transcriptionally inactive. The 8 th cycle coincides with the migration of nuclei to the periphery of the embryo (forming the syncytial blastoderm). Concomitantly, a first wave of ZGA occurs, leading to the expression of about 60 genes , including most of the segmentation genes and the genes required for cellularisation at cycle 14. From then on, the duration of mitotic cycles progressively increases up to 20 minutes at cycle 13. The second wave of ZGA involves over a thousand genes [2, 3]. This massive transcriptional activation coincides with a long pause (about 1 h) during the interphase of the 14 th cycle. During the first thirteen cleavage divisions, the volume of the embryos remains stable while the amount of DNA increases exponentially.
Using haploid mutants (with a nucleo-cytoplasmic (NC) ratio amounting to the half of the wild type one), Edgar et al.  have shown that cellularisation was delayed by one mitotic cycle (cycle 15 instead of 14) and proposed that this phenomena was due to the titration of maternal repressors by the increasing amount of DNA. Pritchard et al.  highlighted that fushi-tarazu repression was dependent on maternal repressor Tramtrack, itself dependent on the NC ratio. More recently Lu et al.  have shown that a few zygotic genes are activated depending on the NC ratio. However, a large fraction of the ZGA wave appears to be independent from the NC ratio and rather depends on the maternal clock model, which assumes that the triggering of gene expression depends on the absolute time after fertilisation. The two afore mentioned mechanisms are not exclusive, and they may play complementary roles in ZGA.
Recently, a combination of genetic and functional genomic studies demonstrated a major implication of the factor Zelda between one and three hours after fertilisation . Zelda has been shown to play a role of general transcription amplifier collaborating with Dorsal , STAT92E , and some other maternal morphogens . This factor binds the TAGteam motif (CAGGTAG), which has been previously proposed to play a role in the activation of pre-cellular blastoderm genes [2, 11, 12]. The TAGteam motif is overrepresented in peaks obtained from ChIP-seq experiments targeting 21 transcription factors involved in embryonic segmentation . Apart from Zelda, which has been recently shown to be involved in the two waves of ZGA , all the other factors reported so far are related with the minor wave. Thus, other factors remain to be identified in order to understand the mechanisms underlying ZGA in Drosophila, including epigenetic regulation.
The goal of our study is to explore the regulatory mechanisms involved in the activation of zygotic genes during the MZT. For this, we started from three transcriptome studies in early Drosophila embryos [2, 3, 6], selected clusters of genes specifically activated during MZT, discovered over-represented motifs in their regulatory region and predicted cis-regulatory modules comprising combinations of these motifs. Interestingly, this “factor-centric” analysis suggests an important role for Trl, a chromatin-remodelling factor, which led us to further investigate the potential associations between ZGA-associated cis-regulatory modules and various epigenetic marks.
It has been recently established by numerous studies that various types of histone modifications affect transcriptional activation, including methylation and acetylation of histone tails to cite a few [14–17]. Using complementary computational tools, we therefore further investigated the relationship between the presence of binding sites for key transcriptional factors and the presence of different in-vivo histone modifications and DNA binding event, focusing on genomic loci associated with ZGA genes. Our computational results prompt a model that tentatively explains the onset of ZGA by a combination of genetic and epigenetic factors.
Results and discussion
Selection of ZGA-responding genes
Transcriptome studies used in this analysis
In order to identify novel factors involved in ZGA, we have used a series of computational analysis tools to revisit three transcriptomic studies: (1) The first study aimed at detecting genes involved in the process of cellularisation: Pilot et al. (2006)  extracted mRNAs at five time points corresponding to fertilisation (T0), slow (T1) and fast (T2) phases of cellularisation, early gastrulation (T3) and late gastrulation (T4), respectively (Figure 1B); (2) De Renzis et al. (2007)  compared the expression profiles of wild-type embryos to those of embryos deleted for half-chromosomes, in order to analyse the respective contributions of maternal and zygotic mRNA during early embryogenesis. They identified five main classes of early expressed genes: (i) maternal and zygotic; (ii) maternal degraded and zygotic; (iii) purely zygotic; (iv) early-activated zygotic; (v) secondary targets (Figure 1C); (3) Lu et al. (2009)  compared expression profiles in haploid mutants versus wild type embryos in order to distinguish genes regulated by the NC ratio from those controlled by the maternal clock (Figure 1D).
Although these studies addressed distinct questions, the three datasets can be re-analysed and combined to extract genes with marked transcription variations in order to identify specific ZGA regulatory features.
Discrete transition profiles as signatures of co-expressed gene clusters
Biological interpretation of the 34 clusters obtained from discrete transition profiles
Nucleo-cytoplasmic ratio NC
d d D d d H
Maternal mRNA degraded during cellularisation
d d D d s H
d s D d s H
Maternal mRNA degraded during slow phase of cellularisation
d s D s s H
Maternal mRNA degraded during slow phase of cellularisation
and zygotic mRNAs transcription during late phase of gastrulation
d s D s u H
Maternal mRNA degraded during slow phase of cellularisation and zygotic mRNAs transcription during early phase of gastrulation
Maternal mRNA degraded during slow phase of cellularisation and zygotic mRNAs transcription during fast phase of cellularisation
s d D s d H
s d D d d H
Maternal mRNA degraded from fast phase of cellularisation
s d D d s H
s s D s d H
Maternal mRNA degraded from early phase of gastrulation
Maternal mRNA degraded from late phase of gastrulation
Zygotic mRNAs transcription during late phase of gastrulation
Zygotic mRNAs transcription during early phase of gastrulation
s s D s u H
Zygotic mRNAs transcription during fast phase of cellularisation
s u D s u H
s u D u u H
Transient zygotic mRNAs transcription during cellularisation
u s D u s H *
Zygotic mRNAs transcription during slow phase of cellularisation
u s D s s H *
u u D u u H *
u u D u s H *
Regarding the analysis of the data of Lu et al. , the transitions between consecutive time points were named by appending the genetic background (denoted by D or H, for diploid or haploid) to the reached time point, with a suffix specifying an early or late stage (E or L respectively). As shown in Figure 2D, transition profiles obtained from Lu experiments in wild type and haploid embryos can be combined in order to distinguish genes responding to the nucleo-cytoplasmic (NC) ratio from those activated by a “maternal clock”. Indeed, genes that depend on NC ratio are expected to respond one mitotic cycle later in haploids than in diploids, since the former embryos contain half the amount of DNA. Thus, the profile “ u s D u s H ” (read “up, then stable diploids, up, then stable haploids”) (Figure 2E) regroups genes activated at transition to the early 14 th mitotic cycle in diploids (transition X D 14E between time points D13 and D14E), but one cycle later in haploids (transition X H 15E between time points H14 and H15E). In contrast, genes whose activation fit the maternal clock model vary at the same absolute time, irrespective of the DNA amount. For example, genes having the profile s u D u s H (Figure 2F) are activated at 165-190 minutes after egg laying in diploids (time point D14L), and at 165-185 minutes in haploid (time point H15E). In total, the 32 = 9 diploid profiles combined with the 32 haploid profiles can form 81 possible transition profiles. However, we obtained only 37 different transition profiles, 24 of which contained at least ten genes. Furthermore, only 16 of them were classified as NC ratio or maternal clock responding genes (Table 1). We left aside the nine remaining clusters because we were not able to interpret the discrete profiles, based on the rules defined in Figure 2E and F.
Grouping of co-expression clusters based on discovered motifs
In addition to the 34 clusters obtained from the discrete transition profiles described in previous section (Additional file 4: Table S1), we included six clusters resulting from the previous published studies: five clusters containing maternal and/or zygotic genes defined by De Renzis and co-workers , and one cluster containing genes activated dependently on the NC ratio, defined by Lu and co-workers .
In order to detect similarities between clusters containing the same type of genes (i.e. maternal, early or late zygotic genes, etc.) and to regroup the most relevant genes for ZGA regulation analysis, we performed a preliminary discovery of over-represented heptanucleotides  in the regulatory regions associated with each of the 40 clusters.
The resulting classification tree shows that the clusters containing only zygotically activated genes appear to have a coherent regulation since they cluster tightly, whereas maternal+zygotic clusters reveal a more complex pattern of regulation. Indeed, some motifs over-represented in first introns and 5’UTR of maternal+zygotic clusters are also over-represented in upstream sequences of zygotic clusters (yellow frame on Figure 4A), whereas the motifs discovered in upstream regions of the maternal+zygotic clusters are also found in upstream sequences of maternal clusters (blue frame on Figure 4A). Moreover, clusters containing genes activated during late cellularisation showed none or few motifs and are present at unresolved branches of the hierarchical tree.
We were mostly interested in zygotic activation; the coherent clustering shown by purely zygotic clusters, and the fact that we did not find any specific motif to NC dependent and independent genes led us to merge the ten u X X X and zygotic clusters (Pilot “usss”, “udss”, “uuss”,“ussd”, Lu “ u s D u s H ”, “ u s D s s H ”, “ u u D u s H ”, “ u u D u u H ”, De Renzis purely zygotic, early zygotic) into a single cluster totalizing 417 genes, hereafter denoted as “ZGA cluster”, for further analysis (Additional file 5: Figure S4B).
The 40 most significant associations between the GO terms and genes of the ZGA cluster
GO term definition
Nb genes in
GO class coverage
BP: anatomical structure morphogenesis
BP: multicellular organismal development
BP: organ development
BP: biological regulation
BP: anatomical structure development
BP: developmental process
BP: embryo development
BP: system development
BP: regulation of biological process
BP: cell fate commitment
BP: pattern specification process
BP: generation of neurons
BP: regulation of cellular process
MF: nucleic acid binding transcription factor activity
MF: sequence-specific DNA binding transcription factor activity
BP: embryonic morphogenesis
BP: cell fate determination
BP: regulation of transcription, DNA-dependent
BP: cellular developmental process
BP: cell differentiation
BP: regulation of RNA metabolic process
BP: organ morphogenesis
BP: tissue development
BP: regulation of cellular macromolecule biosynthetic process
BP: regulation of macromolecule biosynthetic process
BP: regulation of nucleobase, nucleoside, nucleotide and nucleic acid metabolic process
BP: regulation of nitrogen compound metabolic process
BP: multicellular organismal process
BP: regulation of gene expression
BP: tissue morphogenesis
BP: regulation of cellular biosynthetic process
BP: post-embryonic organ development
BP: regulation of biosynthetic process
BP: neuron differentiation
BP: regulation of macromolecule metabolic process
BP: epithelium development
BP: embryonic pattern specification
BP: regulation of cellular metabolic process
Zelda, Tramtrack and Trithorax-like binding motifs are over-represented in ZGA genes
Remarkably, a motif corresponding to the Tramtrack (TTK) binding motif was discovered with the de novo approach. TTK is a maternal repressor, which is progressively titrated as the NC ratio increases during early mitotic cycles, thereby releasing the expression of zygotic genes . Surprisingly, the TTK binding motif is found over-represented in the sequences of pre-cellular activated blastoderm genes and of the genes with the discrete signature “Lu u s D s s H ”, but not in the sequences of genes known to depend on the NC ratio, which might be explained by the intervention of some other factors in this mechanism .
The TTK protein has been reported to physically interact with TRL proteins and to repress TRL-mediated even-skipped activation . TTK could act either directly by binding DNA and repressing the transcription of specific target genes, or indirectly by repressing an activator such as Trl. Interestingly, the TTK motif is significantly under-represented (sig = 5) in upstream sequences of maternal+zygotic and maternal clusters. This is consistent with a repressing activity of TTK. Indeed, the presence of TTK binding sites would result in early inactivation in the presence of maternally expressed Ttk. A motif matching the binding motif of Caudal (a maternal factor involved in segmentation) was further detected as over-represented in purely zygotic genes, but not in the ZGA cluster.
Two motifs were discovered in zygotic clusters, as well as in the ZGA cluster, which do not match any annotated transcription factor binding motif (“AGATACA” and “AaAAGGATCG”). However “AGATACA” was previously reported to be involved in chromosome pairing between regulatory regions associated with the mechanism of transvection . It thus seems particularly relevant that the strongest over-representation of this motif was found in 5’UTRs, as well as in upstream sequences. Finally, the analysis of over-represented motifs in the ZGA cluster revealed four more unknown motifs (not discovered in separated zygotic clusters). Logos and significance of all these motifs are displayed in Figure 5. As a control, we performed motif discovery analyses on 410 randomly selected gene clusters (with 41 different sizes) which did not return any of these motifs (cf. Additional file material available on RSAT website at the address (http://rsat.bigre.ulb.ac.be/rsat/data/published_data/Darbo_2013). This confirms the biological relevance of the discovered motifs.
Based on these results, and in order to predict putative cis-regulatory modules (CRMs), we scanned each type of ZGA non-coding sequences with the nine discovered motifs and predicted cis-regulatory modules (CRMs) by detecting cis-regulatory elements enriched regions (CRERs) using matrix-scan around ZGA defined genes. We detected 528 CRERs in upstream sequences, 313 in the 5’UTR, and 553 in first introns. Because we retrieved non-coding sequences associated with all alternative transcripts, upstream sequences of the smaller transcripts may overlap first introns or 5’UTR sequences. Moreover, in some genes, the first intron is embedded in 5’UTR. About 70% of the upstream sequences, 50% of the first introns and 40% of the 5’UTR contain at least one CRER (Additional file 7: Table S3). Thus, after having merged the CRERs detected in the different types of regulatory regions, we obtained a final set of 1394 non-overlapping CRERs, hereafter denoted as “predicted CRMs”.
In addition to de novo motif discovery, we analysed the enrichment of the ZGA cluster for known motifs, using the program cisTargetX . This tool reveals enriched regulatory features (e.g. motifs or in-vivo datasets) in a set of regions, and ranks these features using a Z-score like enrichements score. Consistently, the results reveal a high enrichment for Zelda (score 12.8) and TRL (score 3.5) binding motifs (Additional file 8: Figure S5). In this analysis, binding motifs for Dorsal (DL), Krüpple (KR) and Bicoid (BCD) were also reported as significantly enriched, which is not unexpected, given the high level of correlation between the binding of these TFs and Zelda. Indeed, a first study of in-vivo binding of BCD, CAD and KR showed that by far the most over-represented motif under the binding peaks was CAGGTAG, hence the binding motif for Zelda . However, many of the CRER do not contain any binding sites for these factors: of the 1394 predicted CRMs, only 765 (54%) contain a predicted binding site for BCD, KR, DL or CAD, using a threshold of p = 0.0001 on the binding site. Hence, the observed enrichment for these transcription factors is restricted to a subset of predicted CRMs. This can be confirmed using in-vivo datasets for these factors including Zelda; we restricted our analysis to the CRER containing a predicted Zelda binding site (780 CRER), of which 599 overlap with an in-vivo Zelda binding event, using the ZLD-ChIP datasets published in . Of these Zelda-bound CRER, 423 (70%) do not have any overlap with any of the 21 TFs published in . Hence, Zelda bound in these CRER is not acting as a precursor for segmentation TFs, which are thus likely to be specifically involved in ZGA.
CRER composition gives insight into ZGA mechanisms
We analysed the motif composition of the defined CRERs, to get insight into the respective contribution of each motif to the ZGA mechanism. A first observation is that roughly 75% of individual binding sites are contained in CRERs, a percentage that is constant across regions (upstream, intron and 5’ UTR), with the exception of 3’ UTR where this proportion is around 60%. This percentage varies depending on which motif is considered (Additional file 9: Figure S6A). Given that the CRER regions span between 15% and 30% of the regions analysed for motifs, this proportion of motifs in CRER represent a significant enrichment over random expectation, and supports the fact that most of the discovered motif instances are indeed bona fide binding sites.
In order to unveil specific organisation patterns, we used a randomization procedure which shuffles the motif instances across CRERs, maintaining the total number of instances of each motif across CRERs, and the number of binding sites in each CRER. A first striking observation is the strong over-representation of homotypic CRER configurations, which is particularly strong in upstream regions (Additional file 9: Figure S6B). As homotypic clusters are known to play an important role in the response to morphogens during early embryogenesis, this is not unexpected, but further supports the validity of the CRERs. The motifs showing the highest over-representations of homotypic clusters are Zelda and the unknown motif AGATACA. For Zelda, the prevalence of homotypic clusters might be a way to respond to “temporal morphogens”, as has been suggested previously , while clusters of AGATAC motifs had been identified previously and hypothesized to play a role in chromosome pairing and DNA looping . While heterotypic configurations are globally under-represented, specific combinations are nevertheless found more often than expected, and might point at particular cooperative mechanisms (Additional file 9: Figure S6C). For example, Zelda is found in heterotypic clusters together with either TTK-like motifs or the previously mentioned AGATACA motif, suggesting a mechanism by which distant enhancers bound by Zelda might be brought into contact with promoter regions with the help of mediator proteins binding AGATACA- motifs .
Using in-vivo datasets to investigate epigenetic mechanisms
The motif discovery analysis described in the previous section, and in particular the presence of Trl-related motifs, suggests a possible involvement of epigenetic factors in the activation of zygotic gene expression. In order to complement the previous motif analysis, we thus decided to make use of recently published in-vivo datasets and to investigate epigenetic regulation by analyzing ChIP-seq and DNAse1 accessibility data, using read densities as well as peaks locations. We focused on the factors CBP (0-4 h), Trl (0-8 h) and Zelda (extracted at 1 h, 2 h and 3 h after fertilisation), modified histones (H3K4me1 0-4 h, H3K4me3, H3K9Ac, H3K27Ac, H3K27me3) and open chromatin (DNAse1 accessibility, stage 5).
Differential motif analysis reveals ZGA specific associations
After having shown that ZGA CRERs reveal specific associations between motifs, we wanted to investigate whether these associations were general, or ZGA specific. In order to do so, we systematically performed a differential motif analysis with the program peak-motifs (RSAT) [27, 28] between the ChIP peaks located in non-coding regions associated (hereafter denoted “ZGA-peaks”) and not associated (“non-ZGA peaks”) with genes of the ZGA cluster (Additional file 10: Figure S7 and Additional file 11: Figure S8).
The Zelda binding motif is over-represented in CBP, TRL and DNAse1 ZGA-peaks vs. non-ZGA peaks, confirming the importance of this factor for the control of zygotic genome activation (Additional file 10: Figure S7). The unknown motif AGATACA appears also to be systematically enriched in ZGA vs. non-ZGA datasets, confirming its relevance to ZGA specific processes.
The differential analysis of Zelda-bound regions at different time points shows that the TRL-related motif and AGATACA are highly differentially enriched, underlying the ZGA-specific association between these three motifs. As expected, and as a control of the differential analysis, the Zelda motif does not appear, being present in ZGA as well as non-ZGA peaks.
CBP does not establish direct interactions with DNA, but interacts with a large diversity of DNA-binding transcription factors. In a recent study, the importance of Dorsal for the recruitment of CBP has been shown . Interestingly, in this study, a strong correlation between CBP and TRL binding had also been shown. Here, we do not find Dorsal binding sites over-represented in ZGA vs. non-ZGA CBP peaks. However, a strong enrichment in Zelda binding motif might suggest that Zelda might take over the role of Dorsal for CBP recruitment in the case of ZGA. The TRL motif is found when motif discovery is performed independently on ZGA and non-ZGA CBP peaks (data not shown), showing that CBP and TRL are indeed associated, as noted previously . The fact that the TRL motif does not appear in the differential analysis is likely due to the fact that CBP and TRL co-localize also outside ZGA-specific regions. However, a much stronger overlap between CBP and TRL peaks appears around ZGA-genes: while 18% of TRL-peaks overlapp a CBP-peak between 0-4 h, the proportion reaches 46% when restricting the analysis to peaks located around ZGA-genes.
High enrichment of CRMs for marks of transcriptional and epigenetic regulation
The previous analysis indicates the prominent role played by CBP, TRL and Zelda around ZGA-specific genes. We then wanted to investigate in more details the importance of these factors at the precise locations of our predicted CRMs.
In order to detect specific associations, we analysed the densities of reads from ChIP-seq experiments under the 1394 predicted CRMs regions. To evaluate the level of enrichment, we ran the same analysis on a positive control set (114 curated blastoderm-specific CRMs from RedFly) and three types of negative sets: (i) regulatory regions of the 417 ZGA genes scanned with randomized (column-permuted) motifs, (ii) regulatory regions of 417 randomly selected genes, and (iii) 317 CRMs not supposed to be active in blastoderm, according to RedFly annotations.
For each of these datasets, we computed the density of reads under CRMs for various marks of transcriptional and epigenetic regulation: Zelda (global transcription factor), CBP (non-DNA binding cofactor) and TRL (chromatin remodelling factor), histone marks, and DNA accessibility profiles, and compared it with the density of reads under randomly selected regions of similar sizes and types (upstream, intron, …). We also computed a p-value using the Wilcoxon rank test in order to evaluate the difference of enrichment between ZGA CRMs and controls (Additional file 12: Table S4).
The ROC curves (Additional file 14: Figure S10) highlight a strong enrichment of ZGA predicted CRMs for Zelda (1 h, 2 h, 3 h), TRL (0-8 h), CBP (0-4 h) and H3K4me1 (0-4 h) as well as DNAse1 hypersensitive sites (stage 5) that together correspond to signatures of active enhancer. This alone confirms the biological relevance of our CRMs defined purely from sequence motifs around ZGA specific genes. Similar levels of association were found in blastoderm-specific CRMs for marks of active enhancers. However, TRL was found enriched for ZGA CRMs but not for blastoderm-specific CRMs (wilcoxon p-value 8e-3). Blastoderm-specific CRMs were also enriched for two repressive marks (H3K9me3 and H3K27me3). This might reflect the tight regulation of the genes controlled by these CRMs, which are active in few spatially located nuclei, but highly repressed by Polycomb-group proteins in the major part of the embryo, as indicated by a recent study by Negre and co-workers . Moreover these repressive marks remain associated with blastoderm CRMs at later stages (Additional file 15: Figure S11).
In contrast, during the time window corresponding to zygotic genome activation (0-4 h), the predicted CRMs of ZGA genes (red curves on Additional file 14: Figure S10) show a significant enrichment for some marks of transcriptional activity (H3K4me1, CBP) but not for repressive marks (H3K27me3, H3K9me3), where the red curve is intermingled with the negative controls (green, purple and blue curves). This seems consistent with a general activation of many genes in the whole embryo.
Previous studies have shown that some of these marks are correlated  and do not act independently from each other. Using a computational strategy developed previously , we used a ranking approach to compute the correlation between these marks for (i) random non-coding regions of the genome matching positional biases of ZGA CRMs, (ii) specifically for the ZGA predicted CRMs as well as for (iii) blastoderm CRMs from Redfly and (iv) central nervous system CRMs from Redfly (Figure 6B; Methods).
Most combinations show a global positive correlation, even in randomly selected regions. Since random regions have been sampled from locations characteristic of ZGA CRER, this reflects a positional effect specific to upstream or intronic regions. The combination CBP/H3K4me1 shows a higher correlation for all three classes of functional elements compared to random regions, as expected from previous studies .
However, some combinations show a much higher degree of correlation for ZGA CRERs compared to random regions or other CRMs, notably CBP/Trl and H3K4me1/Trl. The fact that Trl is involved in these ZGA-specific combinations is interesting, as Trl alone is not the best discriminant between ZGA CRERs and other regions (Figure 6A). While Trl and CBP are known to interact [29, 30], our results suggest that the synergy between them is even higher on ZGA-specific CRMs and might contribute to the activation of the zygotic genome.
From transcriptome data to CRMs prediction and epigenetic context characterisation
The goal of our study was to investigate the mechanism of zygotic genome activation. In order to do so, we (i) re-analysed published datasets to carefully define a list of ZGA-related genes, (ii) applied motif discovery approaches to uncover potential regulators of this process, and (iii) combined in-vivo datasets for various epigenetics factors to understand the interplay between the different regulators of the ZGA.
In particular, using published transcriptome data, we proposed a novel method to cluster gene expression profiles in time-course experiments, which does not require any parameter in order to define co-expression clusters. Functional analysis (expression profiles, non-coding sequence analyses, functional classes enrichment) of the different clusters allowed us to delineate a comprehensive and coherent cluster of genes activated during ZGA. The motifs discovered in the corresponding genes led us to propose several factors and co-factors potentially acting in trans, along with putative cis-regulatory modules.
Analyses of specific associations of predicted CRMs and epigenetic marks led us to propose a model combining different factors (Zelda, TRL, CBP and other unknown factors), which presumably bind accessible and active chromatin regions. In particular, we highlighted to importance of a DNA-motif, AGATACA, which is not yet characterized, but might correspond to a structurally important element or a DNA-binding motif.
From our results, we ranked the predicted CRMs combining TRL, CBP, DNAse1 accessibility, and H3K4me1 data to select the most relevant ones, which can be visualized in their genomic context using the UCSC genome browser. For example, Figure 6C shows the region upstream the TSS of the gene crocodile, a purely zygotic gene, whose activation is dependent on the NC ratio and which is involved in the specification of the most anterior head segment.
Tentative regulatory model and prediction of novel CRMs potentially involved in ZGA control
As TTK becomes titrated by the increasing NC ratio, TRL could be released and become active. Moreover, according to the recent RNA-seq data from Gelbart and Emmert, the amount of Trl mRNA increases from 2h to 4h after egg fertilisation . TTK could thus repress TRL while its abundance is still low, suggesting a mutually enforced effect of TTK titration (NC ratio-dependent) and TRL increase (NC ratio-independent). Binding of TRL could in turn trigger the recruitment of chromatin remodelling complexes. Consistently, we found a high association between predicted CRMs in ZGA-associated regulatory regions, ChIP-seq profiles of TRL binding, and H3K4me1 occupancy.
TRL is not a ZGA-specific factor. What is its exact role during ZGA? While the answer to this question would require experimental validation, our study suggests a mechanism analogous to what has been recently described for dorso-ventral patterning , namely that the specificity of TRL action during ZGA might be conveyed by Zelda. This transcription factor has been shown to be primarily involved in the very early stages of embryogenesis, and we find ZGA-specific over-representation of Zelda binding motifs in CBP bound regions around ZGA-genes (Additional file 10: Figure S7). Our model thus involves Zelda, TRL and possibility another factor binding AGATACA. Specific enrichment of the combination of TRL, CBP, H3K4me1 and open chromatin suggests that the global cofactor CBP could be recruited by these factors at the location of the ZGA CRMs. Consistently, both TRL and CBP contain a Q-rich protein-protein interaction domain [34–36], suggesting potential interactions between these two proteins. The high over-representation of TRL binding motif in CBP peaks reinforces this hypothesis. In contrast, the absence of association with acetylated histone 3 suggests that CBP might not act as an acetyl-transferase here, but instead could act as a bridge between the transcription factors and the basal machinery.
Transcriptome and ChIP-seq data were retrieved from GEO database (http://www.ncbi.nlm.nih.gov/geo/). Transcriptome data: GEO Reference Series IDs GSE3955 , GSE14287 . ChIP-seq data: Zelda dataset (GSE30757), other datasets belong to the GSE23537 super series from ModENCODE project .
DNAse1 accessibility data were retrieved from Berkeley Drosophila transcription Network Project (http://bdtnp.lbl.gov/Fly-Net/browseAccess.jsp). As the reads from these experiments were mapped on the dm2 assembly, genomic coordinates were converted from dm2 to dm3 assembly with liftOver (http://genome.ucsc.edu/cgi-bin/hgLiftOver).
Discrete transition profiles
Each gene is thus characterized by a discrete transition profile described as a vector of letters u,d,s (Figure 2C).
where n is the total number of comparisons between a GO class and a gene cluster. To avoid under-estimating the significance, only genes with at least one annotation in GO were considered for this analysis (the “population size” parameter of compare-classes was set to the number of Drosophila melanogaster genes annotated in GO, while the non-annotated genes were discarded from the clusters for this enrichment analysis step).
Analysis of regulatory sequences
We used the tool retrieve-ensembl-seq  to retrieve non-coding sequences associated to each Drosophila melanogaster gene (upstream, 5’UTR, 3’UTR, first intron). Upstream non-coding sequences were extracted up to the closest neighbor gene, with a maximal length of 5 kb. We activated the options to mask coding sequences and repeats, as well as options to retrieve non-coding sequences for all alternative transcripts and to merge overlapping ones.
To automatize motif discovery on the various non-coding sequence types for the different clusters defined during this study, we used the script gene-cluster-motifs, a task manager available in the stand-alone version of RSAT. Among the different motif discovery algorithms supported by this task manager, we ran oligo-analysis  and dyad-analysis .
These algorithms are based on words and dyads counting respectively. The number of occurrences of each word (dyad) is compared to the expected frequencies observed in a reference sequence set. Specific background models were built for each sequence type (upstream, first intron, 5’UTR, 3’UTR) by computing oligonucleotide and dyad frequencies in the whole set of genomics sequences of the same type. Significance of over-representation is estimated using binomial distribution by computing a nominal p-value.
Over-represented words (oligos) and spaced word pairs (dyads) were assembled and converted to position-specific scoring matrices with the tool matrix-from-patterns (RSAT).
An important advantage of word-based approaches is their scalability: the computing time increases linearly with sequence size, in contast with machine-learning approaches such as MEME or Gibbs motif sampler, whose complexity is quadratic or worse (see  for a quantitative evaluation of time efficiency).
We ran all motif discovery algorithms available in the web site (oligo-, position-, local-word- and dyadanalysis). We searched for over-represented 6- and 7-mers(oligo-, position-, local-word-analysis) and for pairs or trinucleotides spaced by 0 to 20 nucleotides (dyad-analysis). Background was computed from input sequences using a markov model of k - 2 with k representing the oligomer length (oligo-, dyad-analysis). We selected JASPAR Core Insects, DMMPMM and iDMMPMM motif databases for comparison of discovered motifs with known binding motifs.
CisTargetX was used with default parameters, excepting the parameter “Z-score threshold”, for which we selected the option “Determine threshold automatically” instead of the 2.5 default value.
Cis-Regulatory element Enriched Regions (CRERs)
CRERs were predicted with matrix-scan (RSAT) . To compute CRERs significance, we kept sites with a maximal p-value of 10-4, and imposed a distance of at least six nucleotides between consecutive sites to discard overlapping sites that would bias the computed significance. The CRERs length was allowed to vary from 30 and 800 bp. Only CRERs with at least a significance of 2 were further analysed. The background models were computed from input sequences using a Markov model of order 2.
Enrichment of CRMs in ChIP-seq reads
To compare predicted CRM and ChIP-seq profiles, we defined a method to integrate the density of reads over a given region. As a negative control, we measured the read density under random selections of genomic regions of the same sizes as the CRMs. The distributions of densities were compared with ROC (receiver-operating characteristic) curves. The random regions were generated from upstream, first intron, and 5’UTR locations in accordance to the analysed set of predicted CRMs.
Computation of the intensity under a region ( I r )
where x s is the starting position of the region of interest.
where x n and xn+1 are the discrete read positions just before and after x e , respectively, and d n and dn+1 the corresponding read densities.
Generation of random regions
For each CRM type (predicted or curated), we generated ten replicates of random regions of the same lengths as the original CRMs. For each sequence type (upstream, first intron, 5’UTR, 3’UTR), the random regions were retrieved from the whole set of sequences of the same type found in the Drosophila genome. For curated CRMs, random regions were retrieved from upstream sequences since they are almost all present in upstream sequences.
where k c is the resulting rank of a given region of a combination c of the experiments in set ω and is the median of the ranks assigned to the region for all experiments in ω.
Correlation between marks
Following a previous publication , we used a complete partition of the Drosophila non-coding genome representing about 136 K regions, and scored these regions with the marks of interest (CBP, Trl, H3K4me1 and DNAse1 HS sites). All 136 K regions were ranked according to these four features. Next, we extracted the subset of regions overlapping the ZGA CRERs, and computed the Spearman correlation between the ranks of these regions for all pairs of features over 1000 subsamples of 80% of the regions. For sake of comparison, we have extracted the CRMs annotated with the terms Blastoderm (226 regions) or Central nervous system (397 regions) from the Redfly database and performed the same analysis. The barplots on Figure 6B show the mean correlations over the 1000 subsamplings and the error bars indicate the standard deviations. As a negative control, one thousand random regions were sampled from the set of 136 K regions, such that the proportion of upstream and intronic regions matches those of the ZGA CRERs. For each pair of features, the mean and standard deviation of the correlation were computed and plotted.
Additional file information
All data used in this study and the results are available on the supporting Web site (http://rsat.bigre.ulb.ac.be/rsat/data/published_data/Darbo_2013).
This research was the central part of PhD thesis of ED, under the co-direction of DT and JvH in the TAGC laboratory. CH is Maître de Conférences at Aix-Marseille Université. His research activities consist in conceiving, developing, evaluating and applying bioinformatics approaches to analyse regulatory sequences. TL is a group leader in the institute of developmental biology of Luminy Marseille (France). His research activities consist in the understanding of the cell biological basis of cellular organisation and polarity using a combination of genetic, genomic, bio-physical and cell biological techniques. DT was Professor of Bioinformatics at the Université de la Méditerranée (Marseille, France) until January 2010, and is currently Professor of Systems Biology and group leader at the Institute of Biology of the Ecole Normale Supérieure (Paris, France). JvH was Professor at the Université Libre de Bruxelles (Brussels, Belgium) until Oct 2011, and is now Professor at Aix-Marseille Université (Marseille, France). His research activities consist in conceiving, developing, evaluating and applying bioinformatics approaches to analyse regulatory sequences and biomolecular interaction networks.
Area under curve
Chomatin Immuno-Precipitation and sequencing
Cis-regulatory element enriched region
Transcription start site
Transcription factor binding site
ED salary was funded by the Agence Nationale de la Recherche (ANR), through the project NeMo (ANR-07-BLAN-061). CH is supported by the ANR Young Researchers Grant “CardiHox”. The TAGC laboratory acknowledge support from the EU ERA-NET Plus scheme in FP7, through the ERASysBio+ project ModHeart. The BiGRe and TAGC laboratories acknowledge support from the EU-funded COST action [BM1006 “Next Generation Sequencing Data Analysis Network”]. Collaboration between DT and JvH was initially supported by the Belgian Federal Science Policy Office (Interuniversity Attraction Poles, project P6/25 - BioMaGNet) and further stimulated by a sabbatical year spent by JvH at the TAGC lab, and by a two-months Invited Professorship of JvH at ENS. We warmly thank Nathalie Dostatni and Morgane Thomas-Chollier for critical reading and thoughtful suggestions.
- Tadros W, Lipshitz HD: The maternal-to-zygotic transition: a play in two acts. Development. 2009, 136 (18): 3033-3042.View ArticlePubMedGoogle Scholar
- De Renzis S, Elemento O, Tavazoie S, Wieschaus EF: Unmasking activation of the zygotic genome using chromosomal deletions in the Drosophila embryo. PLoS Biol. 2007, 5 (5): e117-PubMed CentralView ArticlePubMedGoogle Scholar
- Pilot F, Philippe JM, Lemmers C, Chauvin JP, Lecuit T: Developmental control of nuclear morphogenesis and anchoring by charleston, identified in a functional genomic screen of Drosophila cellularisation. Development. 2006, 133 (4): 711-723.View ArticlePubMedGoogle Scholar
- Edgar BA, O’Farrell PH: Genetic control of cell division patterns in the Drosophila embryo. Cell. 1989, 57: 177-187.PubMed CentralView ArticlePubMedGoogle Scholar
- Pritchard DK, Schubiger G: Activation of transcription in Drosophila embryos is a gradual process mediated by the nucleocytoplasmic ratio. Genes Dev. 1996, 10 (9): 1131-1142.View ArticlePubMedGoogle Scholar
- Lu X, Li JM, Elemento O, Tavazoie S, Wieschaus EF: Coupling of zygotic transcription to mitotic control at the Drosophila mid-blastula transition. Development. 2009, 136 (12): 2101-10.PubMed CentralView ArticlePubMedGoogle Scholar
- Harrison MM, Li XY, Kaplan T, Botchan MR, Eisen MB: Zelda binding in the early Drosophila melanogaster embryo marks regions subsequently activated at the maternal-to-zygotic transition. PLoS Genet. 2011, 7 (10): e1.002266-View ArticleGoogle Scholar
- Reeves GT, Stathopoulos A: Graded dorsal and differential gene regulation in the Drosophila embryo. Cold Spring Harb Perspect Biol. 2009, 1 (4): a.000836-View ArticleGoogle Scholar
- Tsurumi A, Xia F, Li J, Larson K, Lafrance R, Li WX: STAT is an essential activator of the zygotic genome in the early Drosophila Embryo. PLoS Genet. 2011, 7 (5): e1.002086-View ArticleGoogle Scholar
- Kanodia JS, Liang HL, Kim Y, Lim B, Zhan M, Lu H, Rushlow CA, Shvartsman SY: Pattern formation by graded and uniform signals in the early Drosophila embryo. Biophys J. 2012, 102 (3): 427-33.PubMed CentralView ArticlePubMedGoogle Scholar
- ten Bosch JR, Benavides JA, Cline TW: The TAGteam DNA motif controls the timing of Drosophila pre-blastoderm transcription. Development. 2006, 133 (10): 1967-77.View ArticlePubMedGoogle Scholar
- Liang HL, Nien CY, Liu HY, Metzstein MM, Kirov N, Rushlow C: The zinc-finger protein Zelda is a key activator of the early zygotic genome in Drosophila. Nature. 2008, 456 (7220): 400-403.PubMed CentralView ArticlePubMedGoogle Scholar
- Satija R, Bradley RK: The TAGteam motif facilitates binding of 21 sequence-specific transcription factors in the Drosophila embryo. Genome Res. 2012Google Scholar
- Li L, Zhu Q, He X, Sinha S, Halfon MS: Large-scale analysis of transcriptional cis-regulatory modules reveals both common features and distinct subclasses. Genome Biol. 2007, 8 (6): R101-PubMed CentralView ArticlePubMedGoogle Scholar
- Suganuma T, Workman JL: Crosstalk among histone modifications. Cell. 2008, 135 (4): 604-607.View ArticlePubMedGoogle Scholar
- Nakayama T, Nishioka K, Dong YX, Shimojima T, Hirose S: Drosophila GAGA factor directs histone H3.3 replacement that prevents the heterochromatin spreading. Genes Dev. 2007, 21 (5): 552-561.PubMed CentralView ArticlePubMedGoogle Scholar
- Shimojima T, Okada M, Nakayama T, Ueda H, Okawa K, Iwamatsu A, Handa H, Hirose S: Drosophila FACT contributes to Hox gene expression through physical and functional interactions with GAGA factor. Genes Dev. 2003, 17 (13): 1605-1616.PubMed CentralView ArticlePubMedGoogle Scholar
- van Helden J, André B, Collado-Vides J: Extracting regulatory sites from the upstream region of yeast genes by computational analysis of oligonucleotide frequencies. J Mol Biol. 1998, 281 (5): 827-842.View ArticlePubMedGoogle Scholar
- Thomas-Chollier M, Sand O, Turatsinze JV, Janky R, Defrance M, Vervisch E, Brohée S, van Helden J: RSAT: regulatory sequence analysis tools. Nucleic Acids Res. 2008, 36: W119-W127.PubMed CentralView ArticlePubMedGoogle Scholar
- Pagans S, Ortiz-Lombardía M, Azorín F, Espinás ML: The Drosophila transcription factor tramtrack (TTK) interacts with Trithorax-like (GAGA) and represses GAGA-mediated activation. Nucleic Acids Res. 2002, 30 (20): 4406-4413.PubMed CentralView ArticlePubMedGoogle Scholar
- Lewis EB, Knafels JD, Mathog DR, Celniker SE: Sequence analysis of the cis-regulatory regions of the bithorax complex of Drosophila. Proc Natl Acad Sci U S A. 1995, 92 (18): 8403-8407.PubMed CentralView ArticlePubMedGoogle Scholar
- Turatsinze JV, Thomas-Chollier M, Defrance M, van Helden J: Using RSAT to scan genome sequences for transcription factor binding sites and cis-regulatory modules. Nat Protoc. 2008, 3 (10): 1578-1588.View ArticlePubMedGoogle Scholar
- Aerts S, Quan XJ, Claeys A, Naval Sanchez M, Tate P, Yan J, Hassan BA: Robust target gene discovery through transcriptome perturbations and genome-wide enhancer predictions in Drosophila uncovers a regulatory basis for sensory specification. PLoS Biol. 2010, 8 (7): e1.000435-View ArticleGoogle Scholar
- Li Xy, MacArthur S, Bourgon R, Nix D, Pollard DA, Iyer VN, Hechmer A, Simirenko L, Stapleton M, Luengo Hendriks CL, et al: Transcription factors bind thousands of active and inactive regions in the Drosophila blastoderm. PLoS Biol. 2008, 6 (2): e27-PubMed CentralView ArticlePubMedGoogle Scholar
- MacArthur S, Li XY, Li J, Brown JB, Chu HC, Zeng L, Grondona BP, Hechmer A, Simirenko L, Keränen SVE, et al: Developmental roles of 21 Drosophila transcription factors are determined by quantitative differences in binding to an overlapping set of thousands of genomic regions. Genome Biol. 2009, 10 (7): R80-PubMed CentralView ArticlePubMedGoogle Scholar
- Nien CY, Liang HL, Butcher S, Sun Y, Fu S, Gocha T, Kirov N, Manak JR, Rushlow C: Temporal Coordination of Gene Networks by Zelda in the Early Drosophila Embryo. PLoS Genet. 2011, 7 (10): e1.002339-10.1371\%2Fjournal.pgen.1002339.View ArticleGoogle Scholar
- Thomas-Chollier M, Herrmann C, Defrance M, Sand O, Thieffry D, Van Helden J: RSAT peak-motifs: motif analysis in full-size ChIP-seq datasets. Nucleic Acids Res. 2012, 40 (4): e31-PubMed CentralView ArticlePubMedGoogle Scholar
- Thomas-Chollier M, Darbo E, Herrmann C, Defrance M, Thieffry D, van Helden J: A complete workflow for the analysis of full-size ChIP-seq (and similar) data sets using peak-motifs. Nat Protoc. 2012, 7 (8): 1551-1568.View ArticlePubMedGoogle Scholar
- Holmqvist Ph, Boija A, Philip P, Crona F, Stenberg P, Mannervik M: Preferential genome targeting of the CBP co-activator by Rel and Smad proteins in early Drosophila melanogaster embryos. PLoS Genet. 2012, 8 (6): e1.002769-View ArticleGoogle Scholar
- Nègre N, Brown CD, Ma L, Bristow CA, Miller SW, Wagner U, Kheradpour P, Eaton ML, Loriaux P, Sealfon R, et al: A cis-regulatory map of the Drosophila genome. Nature. 2011, 471 (7339): 527-531.PubMed CentralView ArticlePubMedGoogle Scholar
- Kharchenko PV, Alekseyenko AA, Schwartz YB, Minoda A, Riddle NC, Ernst J, Sabo PJ, Larschan E, Gorchakov AA, Gu T, et al: Comprehensive analysis of the chromatin landscape in Drosophila melanogaster. Nature. 2011, 471 (7339): 480-485.PubMed CentralView ArticlePubMedGoogle Scholar
- Herrmann C, Van de Sande B, Potier D, Aerts S: i-cisTarget: an integrative genomics method for the prediction of regulatory features and cis-regulatory modules. Nucleic Acids Res. 2012, 40 (15): e114-PubMed CentralView ArticlePubMedGoogle Scholar
- Marygold SJ, Leyland PC, Seal RL, Goodman JL, Thurmond JR, Strelets VB: Wilson RJ and the FlyBase Consortium: FlyBase: improvements to the bibliography. Nucleic Acids Res. 2013, 41 (D1): D751-D757.PubMed CentralView ArticlePubMedGoogle Scholar
- Chopra VS, Srinivasan A, Kumar RP, Mishra K, Basquin D, Docquier M, Seum C, Pauli D, Mishra RK: Transcriptional activation by GAGA factor is through its direct interaction with dmTAF3. Dev Biol. 2008, 317 (2): 660-670.View ArticlePubMedGoogle Scholar
- Wilkins RC, Lis JT: DNA distortion and multimerization: novel functions of the glutamine-rich domain of GAGA factor. J Mol Biol. 1999, 285 (2): 515-525.View ArticlePubMedGoogle Scholar
- Suhara W, Yoneyama M, Kitabayashi I, Fujita T: Direct involvement of CREB-binding protein/p300 in sequence-specific DNA binding of virus-activated interferon regulatory factor-3 holocomplex. J Biol Chem. 2002, 277 (25): 22304-22313.View ArticlePubMedGoogle Scholar
- R Development Core Team: R: A Language and Environment for Statistical Computing. 2012, Vienna: R Foundation for Statistical ComputingGoogle Scholar
- Gentleman RC, Carey VJ, Bates DM, Bolstad B, Dettling M, Dudoit S, Ellis B, Gautier L, Ge Y, Gentry J, et al: Bioconductor: open software development for computational biology and bioinformatics. Genome Biol. 2004, 5 (10): R80-PubMed CentralView ArticlePubMedGoogle Scholar
- Irizarry RA, Bolstad BM, Collin F, Cope LM, Hobbs B, Speed TP: Summaries of Affymetrix GeneChip probe level data. Nucleic Acids Res. 2003, 31 (4): e15-PubMed CentralView ArticlePubMedGoogle Scholar
- Grumbling G: FlyBase: anatomical data, images and queries. Nucleic Acids Res. 2006, 34: D484-D488.PubMed CentralView ArticlePubMedGoogle Scholar
- van Helden J: Regulatory sequence analysis tools. Nucleic Acids Res. 2003, 31 (13): 3593-3596.PubMed CentralView ArticlePubMedGoogle Scholar
- Sand O, Thomas-Chollier M, van Helden J: Retrieve-ensembl-seq: user-friendly and large-scale retrieval of single or multi-genome sequences from Ensembl. Bioinformatics (Oxford, England). 2009, 25 (20): 2739-2740.View ArticleGoogle Scholar
- Defrance M, Janky R, Sand O, van Helden J: Using RSAT oligo-analysis and dyad-analysis tools to discover regulatory signals in nucleic sequences. Nat Protoc. 2008, 3 (10): 1589-1603.View ArticlePubMedGoogle Scholar
- Portales-Casamar E, Thongjuea S, Kwon AT, Arenillas D, Zhao X, Valen E, Yusuf D, Lenhard B, Wasserman WW, Sandelin A: JASPAR 2010: the greatly expanded open-access database of transcription factor binding profiles. Nucleic Acids Res. 2010, 38 (Database issue): D105-D110.PubMed CentralView ArticlePubMedGoogle Scholar
- Zhu LJ, Christensen RG, Kazemian M, Hull CJ, Enuameh MS, Basciotta MD, Brasefield JA, Zhu C, Asriyan Y, et al: FlyFactorSurvey: a database of Drosophila transcription factor binding specificities determined using the bacterial one-hybrid system. Nucleic Acids Res. 2011, 39 (Database issue): D111-D117.PubMed CentralView ArticlePubMedGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License(http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.