In this study we set out to investigate the transcriptome dynamics of early zebrafish embryogenesis with a comprehensive experimental set-up that interrogated individual embryos about every minute during the circa 3 h period during epiboly from late blastula to the mid gastrula stage. We did so, because most studies into (zebrafish) embryogenesis are typically using pools of staged embryos plus quite dispersed time points, and we reasoned that a high-resolution time course could provide much insight into the true dynamics of the individual genes of the embryonic transcriptome. The fact that we were able to identify several oscillating genes shows the value of our approach. We set-out to prove that it should be feasible to pinpoint the starting and stopping points of individual genes. We were able to do so for the genes that started or stopped during the epibolic stage. Hence it is also possible to define consecutively starting, or stopping genes, which should be of great help reconstructing cellular processes and molecular pathways. In this study, we did not aim to unravel the role of individual genes nor that of the molecular pathways in zebrafish embryogenesis. This should be done by the appropriate experts on each molecular process. To facilitate them, we have not only uploaded our entire dataset to NCBI’s Gene Expression Omnibus with accession number GSE83395 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE83395), but also set-up a web resource in which experts can conveniently investigate the expression profiles of the genes of their choice: http://rnabiology.nl/Dr-Browser.html in combination with the expression in zebrafish oocytes from our previous study . To gain perspective, it is important that by interpreting individual genes, one is aware what the common transcriptome dynamics are for all genes during the same developmental stage.
To put our study in context, we compared our finding with three other studies: Junker et al.  on genome-wide RNA tomography, which has one sample point (shield stage) that overlaps with our time-course; Aanes et al.  of which the 5.3-hpf sample point overlaps with the start of our time course; and our own study Rauwerda et al.  concerning maternal RNA in zebrafish eggs. In general, the findings of the current study were in line with the previous studies (Additional file 15). The genome-wide RNA tomography by Junker et al. provides an exciting opportunity to link the temporal gene expression data to embryonic location. Given that genes primarily function in the context of a cellular network this would allow the decomposition of the temporal transcriptome data from multicellular embryos to clusters of genes that operate together in pathways or cellular mechanisms. It was possible to look into this, as the RNA tomography data overlapped with one embryonic developmental sample point (shield stage) with our data. A similar number of expressed genes (7,358), with 70% overlap in our study were found in the overlapping shield-stage samples from Junker et al. (2014) (Additional file 15-B). Hence, we clustered the tomography data in 16 spatial gene-expression patterns and analysed which pathways were overrepresented in these spatial expression clusters. Genes of the KEGG Pathway “Ribosome”, which also is overrepresented in our time course experiment, were indeed found to be overrepresented in the same spatial expression cluster (Additional file 16). However, for many other overrepresented pathways from our expressed genes, we failed to find them back using the tomography data. We feel that the “Ribosome” pathway probably succeeded as these genes are highly expressed. Unfortunately, although the overall data of the tomography study is excellent, the data of the overlapping shield sample point is rather noisy compared to the other sample points, with only a few distinct gene-expression patterns. So, we believe that new tomography data could help decompose the whole-genome temporal transcriptome data, which could support the discovery of new pathways or new genes into known pathways.
In another comparison, also a high-percentage overlap (77% relative to our study) was found with the Aanes et al. (2011) study with also a similar relative gene expression for the overlapping genes per Aanes-defined gene-expression clusters in both studies (Additional file 15-D).
A comparison of the current results with our previous study  on maternal RNA revealed that 87% of the expressed genes were also found to be present as maternal RNA in egg (Additional file 15-F). Given the assumptions that in epiboly almost all maternal RNA is cleared and maternal RNA primarily functions in the pre-gastrula development, this is in line with the finding of Lee et al. (2013) that 74% of the early-zygotically transcribed genes still had a maternal contribution. There is however, no clear relation between the expression level for each gene as observed in eggs and that at the late blastula stage (Additional file 15 E-F). There is also no correlation between the genes that showed a distinct mother-specific effect in egg expression and a spawn-specific effect in epiboly expression. Altogether, the clearance and subsequent transcription of these genes seems a complicated aspect of the MZT.
The most impressive result of the high-resolution embryonic time course is the extremely tight regulation of gene expression, even over different embryos from different spawns. A technical conclusion is that pooling of embryos for experimentation poses no problem, if the staging of the individual embryos is done properly. The observed spawn-specific differences will obviously be lost, if eggs from different spawns are combined. At the same time, if spawns are not evenly distributed between pools, spawn-specific gene-expression differences might lead to false positive differential gene-expression conclusions.
More importantly, the presence of such a tightly regulated embryonic program for transcription, which took us by surprise, signals a rather robust system that in essence is resilient to the differences at many cellular levels that must occur between individuals. Notably, even though we classified the expressed genes into ten different types, virtually each individual gene has its own unique expression profile with respect to intensity, change and variability.
Our unique experimental set-up also allowed us to identify the starting and stopping points of gene-expression at an omics scale. Only several hundreds of genes appeared to be switched on or off during this period of embryogenesis, which seems to be quite a quite limited number given the multitude of developmental processes that occur. We must however consider that even though we found many genes expressed, we might miss out on several important developmental regulating genes as they could be below the detection limit of our experimental technique due to low gene expression or only be expressed just in a small fraction of the embryo cells. Regardless, if we accept the premise that maternal RNA is cleared previous to the gastrula stage, this means that, except for the identified 157 starting genes (type 9), all zygotic gene expression has to start prior to this embryonic stage. Although for some genes the maternal and zygotic contributions can be deconvoluted on the basis of their 3’UTR [17, 18], for many genes, polyadenylation of maternal transcripts and zygotic gene expression are confounded. So the expression levels and differences for each individual gene are a net result of those poly-A+ generating mechanisms, as well as specific RNA-degrading mechanisms. Collectively, the gene-expression changing rates are quite modest, on average about 1.3 FC/h. They do not explain the observed gene-expression levels as extrapolation of present genes predicts for most of the expressed genes a starting point well before fertilization. As, besides the oscillating genes, no genes with peaking gene-expression were found, we conclude that gene expression in this embryonic stage is unexpectedly gradual, which might be the sign of a complex, strictly regulated system.
Our final observation with respect to transcriptome dynamics concerns the spawn-specific gene-expression effects. Much alike we observed in our oocyte study, we detected several genes with a substantial spawn-specific difference. Moreover, it appeared that this property might be present in the majority, if not all genes, albeit often with extremely small differences. As the eggs in our previous study were unfertilized, we could identify this effect to be mother-specific. In the present study it is undoubtedly a mixed effect of several maternal and paternal epigenetic factors that collectively produce these differences in gene expression. This would also explain why in eggs multilevel expressed genes are almost always expressed at one level in all eggs of one mother, whereas in this study embryos from one spawn often show multilevel gene expression. Although we were unable to model all the gene-expression differences between the embryos, it is clear that given the strict regulation of the transcriptome, they will be related. As such, the spawn-specific differences provide an additional layer of information that could support the understanding of transcriptome regulation in embryogenesis.