Cytoplasmic polyadenylation-mediated translational control of maternal mRNAs directs maternal-to-zygotic transition

A period of transcriptional silence follows fertilization in animal embryogenesis, during which development is mainly driven by maternally deposited mRNAs and proteins (Korzh, 2009). At a certain time point, the transition of developmental control from maternal to zygotic occurs (Tadros and Lipshitz, 2009). This event is referred to as the maternal-to-zygotic transition (MZT) and is represented by the occurrence of zygotic genome activation (ZGA; Newport and Kirschner, 1982), the mass degradation of maternal mRNAs (Hamatani et al., 2004; Piko and Clegg, 1982; Walser and Lipshitz, 2011) and major epigenetic reprogramming (Lindeman et al., 2011; Vastenhouw et al., 2010). Transcriptional silence necessitates gene regulation at a different level to ensure appropriate activity of maternal mRNAs. Translational regulation of maternal mRNAs by cytoplasmic polyadenylation (CPA) is crucial for the timely activation of maternal mRNAs during oocyte maturation and the early embryonic development of several model organisms (Aoki et al., 2003; Evsikov et al., 2004; Hwang et al., 1997; Lieberfarb et al., 1996; Oh et al., 2000; Potireddy et al., 2006; Salles et al., 1992). However, there is still a lack of understanding of its biological function and molecular targets during embryogenesis.

In fish and frogs, the mid-blastula transition (MBT) represents a key event within the MZT. At the cellular level, notable changes include the asynchronization of cell division, the acquisition of cell motility and gastrulation (Newport and Kirschner, 1982; Rott and Sheveleva, 1968). The MZT is thus a continuous process that encompasses MBT events, including zygotic genome activation (ZGA), which marks the initiation of gastrulation, together with maternal mRNA clearance. In zebrafish, MBT is initiated at the 10th cell cycle and represents a continuous process lasting approximately until the 13th cell cycle (Kane and Kimmel, 1993). We previously identified a large cohort of maternal mRNAs that undergo CPA during the course of embryogenesis, leading to the MBT in zebrafish (Aanes et al., 2011). Strikingly, CPA occurred almost exclusively prior to the MBT, which led us to question whether these RNAs are implicated in the regulation of this important developmental transition.

Several key elements regulate ZGA. First, maternal inhibiting factors need to be titrated by an increase in DNA content as a result of cell division that increases the nucleocytoplasmic ratio (Dekens et al., 2003; Evsikov et al., 1990; Kane and Kimmel, 1993; Newport and Kirschner, 1982). Concentrated extracts of Xenopus eggs contain active transcription machinery, but these are outcompeted by the abundant supply of maternal histones, resulting in the overall repression of transcription (Amodeo et al., 2015; Prioleau et al., 1994). Histones have also been implicated in regulating the accessibility of the genome to activating transcription factors that determine the onset of ZGA (Joseph et al., 2017; Pálfy et al., 2017). In addition to the large supply of histones for chromatin assembly, the bivalent modification of these histone proteins was observed both prior to the MBT and at the MBT in various organisms (Andersen et al., 2013). The nuclear-to-cytoplasmic volume ratio has also been implicated to play a role in determining the timing of ZGA independently of DNA levels (Jevtić and Levy, 2015, 2017). Second, the timely induction of activating factors serves to promote ZGA (Aoki et al., 2003; Oh et al., 2000). In zebrafish, several maternal factors have previously been implicated in this process, namely nanog, pou5f3 and sox19b (Lee et al., 2013; Onichtchouk et al., 2010; Xu et al., 2012).

Here, we investigate the association of maternal and early zygotic transcripts with polysomes during the first few hours of embryonic development over the course of the MZT. Our analysis revealed distinct polysome association dynamics of cytoplasmically polyadenylated maternal mRNAs. We found that developmental progression through the MZT depends on the precise translational activation of a subset of maternal mRNAs by CPA. Our data revealed the complexity of mRNA expression dynamics and showed that CPA is essential for two key developmental events that determine embryonic passage through the MZT: (1) ZGA and (2) the clearance of a large number of maternal transcripts.

To profile the translational regulatory landscape during early zebrafish development, we performed polysome profiling at several developmental stages: activated unfertilized egg (egg), pre-MBT (1-cell, 16-cell and 128-cell), MBT [high, 3.5 h post-fertilization (hpf)] and post-MBT (50% epiboly, 5.3 hpf). Two fractions that were distinguishable by their density due to the presence or absence of polysome association were isolated and sequenced (Fig. 1A, Fig. S1, Table S1). Unfractionated total RNA from corresponding stages were also sequenced to provide a measure of overall transcript level. Using the total unfractionated RNA-seq dataset, we classified the transcriptome into two broad categories based on their origin: maternal and zygotic. We defined RNAs with transcripts that were present at high (≥10TPM) levels at the egg stage and did not increase before 3.5 hpf as maternal, regardless of whether zygotic contributions exist later on. We defined RNAs that were expressed only from 3.5 hpf onwards as zygotic. We also assigned the pre-MBT cluster to accommodate transcripts whose expression increased prior to 3.5 hpf. Applying these criteria to the RNA-seq data of the total (unfractionated) mRNA, we obtained 9154 maternal genes, 831 pre-MBT genes and 667 zygotic genes (Fig. 1B, Table S2).

To validate our gene clusters, we compared our data with that of Harvey et al. (2013), in which crosses between different zebrafish strains were performed to enable distinctions between maternal and paternal alleles by single-nucleotide polymorphisms (SNPs). To allow direct comparisons, we re-analyzed the raw sequencing reads from this dataset and identified 4393 genes that had both maternal and paternal alleles (maternal-zygotic and zygotic genes) and 487 genes that had only maternal alleles (maternal genes). We found that 337 genes (69%) and 21 genes (4%) out of 487 genes with maternal-only alleles were within our maternal and pre-MBT clusters, respectively, whereas none of these were in the zygotic cluster (Table S18). The lack of overlap between the maternal-only genes and our zygotic cluster confirmed the zygotic origin of these genes. A total of 3755 genes (85%), 127 genes (3%) and 41 genes (0.9%) out of 4393 genes that had both maternal and paternal alleles in Harvey et al. (2013) dataset were within the maternal, pre-MBT and zygotic clusters, respectively, reflecting the zygotic contribution for these genes. Altogether, the results of these comparisons suggested that the method used in our study was able to distinguish maternal from zygotic genes.

To assess the association between translational regulation and CPA, identity needs to be assigned to transcripts based on their polyadenylation state. We focused on the maternal group of transcripts to identify those that underwent CPA during pre-MBT stages. As we employed mRNA isolation methods based on rRNA depletion, the information on the lengthening of the poly(A) tail of mRNAs that underwent CPA could not be inferred from this dataset. Therefore, we revisited our previous RNA-seq dataset (Aanes et al., 2011), in which mRNA isolation was performed using oligo d(T) selection that is biased towards those with a longer poly(A) tail. Thus, in the case of transcripts that underwent CPA, this results in an artificial increase in transcript level over time prior to the MBT. This particular group was represented by the ‘pre-MBT’ cluster in this dataset. We re-analyzed this dataset using the same criteria as those used for the current dataset to determine maternal, pre-MBT and zygotic transcripts, and re-mapped the raw reads to the latest genome assembly (zv10). Re-analysis of this dataset resulted in 3959 pre-MBT, 4598 maternal and 484 zygotic genes (Fig. 1B, Table S2).

We then compared clusters that were derived from the two datasets. The numbers of genes in the zygotic clusters were comparable, but a substantial majority of Aanes pre-MBT transcripts [3667 genes (92.6%)] shifted to the maternal cluster in the rRNA-depleted dataset (Fig. 1C). This observation was consistent with the situation in which these 3667 mRNAs were truly undergoing CPA. Thus, we defined three distinct expression clusters in our current data: a maternal cluster [which was subdivided into a non-cytoplasmically polyadenylated subcluster (non-CPA, 5487 genes, 60%; and a cytoplasmically polyadenylated subcluster (CPA, 3667 genes, 40%)]; a pre-MBT cluster (831 genes); and a zygotic cluster (667 genes).

Among non-CPA subcluster transcripts, GO terms that were enriched (P<0.01) included ‘translation’, ‘DNA repair’, ‘cell cycle’ and ‘chromatin modification’ (Fig. S3, Table S3). These are relevant to embryonic development prior to the MBT when transcription is absent and embryonic cells undergo synchronous divisions. Interestingly, the enriched categories also included ‘mRNA splicing’. Maternal mRNAs are known to be spliced into a mature form prior to embryogenesis, and the presence of transcripts that regulate this process possibly represent a trace of earlier events that preceded fertilization. The CPA subcluster was enriched in several GO terms, including ‘protein transport’, ‘ubiquitin-mediated protein catabolic process’ and ‘cell cycle’ (Fig. S3, Table S4). Interestingly, this cluster also contained notable genes that were implicated in early developmental defects (Abrams and Mullins, 2009). These include buc (Dosch et al., 2004), piwil2 (Houwing et al., 2008), mapkapk2a and chuk (Wagner et al., 2004), and dhx16 (Pelegri and Mullins, 2004). In contrast, the zygotic cluster was enriched in GO terms related to embryonic patterning, such as ‘cell migration involved in gastrulation’, ‘dorsal/ventral pattern formation’ and ‘regulation of endoderm and mesoderm fate specification’ (Fig. S3, Table S5). Notably, this cluster was also enriched in transcription factors, likely reflecting the earliest transcriptional activity following ZGA.

Surprisingly, 831 transcripts were identified as being expressed pre-MBT. Libraries were prepared using an rRNA depletion method, and the presence of such transcripts suggests true increase in abundance prior to ZGA. These possibly represent the earliest wave of zygotic transcription that precedes the main wave at MBT. This group was enriched in the GO terms ‘nucleosome assembly’ and ‘chromatin silencing’, and contained transcripts encoding nucleosome core components (Fig. S3, Table S6), including multiple histone-encoding transcripts, which are non-polyadenylated (Marzluff, 2005). Our rRNA depletion-based strategy thus enabled their detection as a proof of principle of the improved detection of unique transcripts using this approach. To validate pre-MBT genes as true early zygotic transcripts, we compared them with those that were previously identified as first-wave genes (Lee et al., 2013). The comparison returned a very small overlap (14 of 269 first-wave genes and 831 pre-MBT genes; Table S18). We speculated that this might be caused by the later stage observed in the Lee study (4 hpf, dome stage), such that their gene list would also include zygotic genes that were initiated at the typical timing of ZGA (from 3.5 hpf in the present study). Indeed, a larger overlap (76 genes out of 269 first-wave genes and 667 zygotic genes) was observed with our zygotic cluster genes. Additionally, because our maternal cluster also contained genes with a zygotic contribution (i.e. maternal-zygotic genes), we obtained an equally sizeable overlap (86 genes out of 269 first-wave genes and 9154 maternal genes). Although the differences in the observed stages did not allow the validation of pre-MBT transcripts, it reaffirmed the zygotic origin of those within our zygotic cluster. We then turned to another dataset where early zygotic genes were detected through 4-sUTP labeling of nascent transcripts (Heyn et al., 2014). We found 257 out of the 592 genes reported as early zygotic genes in common with our dataset. However, these genes were distributed between the maternal, pre-MBT and zygotic clusters (113, 53 and 91 genes, respectively). It is worth noting that 440 of the early zygotic genes in Heyn dataset were also maternal-zygotic, some of which would unsurprisingly distribute to our maternal cluster. Additional differences were also caused by their ability to detect many miRNA species (56 in total), which is not the case in our dataset due to different methods used to isolate mRNAs. Owing to the unavailability of other comparable transcriptome data, we could not validate the remainder of the pre-MBT cluster genes. Nevertheless, their robust increase in transcript levels pre-MBT suggests contribution from zygotic transcription to this increase.

Hierarchical clustering and principal component analysis based on transcriptome profiles revealed tight grouping according to polysome association for samples prior to MBT and according to developmental stage for the MBT and post-MBT samples (Fig. 1D). Although this reflects the relatively constant nature of maternal mRNA composition, the distinction between polysome-bound and -unbound fractions suggests the presence of selective polysome engagement of transcripts prior to MBT. However, this selectivity disappeared at the MBT and post-MBT stages, indicated by the comparatively similar transcriptome profiles of both fractions. This is consistent with the observation that the length of poly(A) tail determined the translation rate only at 2 hpf and not at 4 or 6 hpf (Subtelny et al., 2014). The differential clustering behavior of the samples from the early and late stages thus reaffirms the notion that pan-embryonic translational regulation operates only on maternal mRNAs prior to the MBT.

To determine whether the CPA of maternal transcripts corresponds to their translational regulation, we performed a stepwise analysis of the polysome-bound and -unbound datasets. First, we sought to determine whether CPA correlates with the recruitment of transcripts to polysomes. We calculated the polysome association of each transcript species as a ratio of bound fraction over the total (bound+unbound) fraction (Fig. 2A). Taking the average values of all of the transcripts in a cluster, a decreasing trend of polysome association was observed in the non-CPA subcluster, whereas an opposite trend indicating an increase in polysome association was observed for the CPA subcluster (Fig. 2A). Next, we defined translational regulation as a disproportionate change in translation compared with transcript level. We reasoned that if translational regulation is exclusive to maternal mRNAs, then an increase in polysome association should be observed only in this subset and not in the subset of zygotic transcripts. We first defined the transcript level (TS) by measure of total mRNA levels and translation (TL) by measure of bound fraction (Fig. 2B). Using the egg stage as a baseline, we calculated the ΔTS and ΔTL of each transcript for each developmental stage. We classified transcripts based on their ΔTS: transcript level up or down (>2-fold compared with egg stage) or constant (<2-fold compared with egg stage). We then defined translational regulation by comparing the ΔTL of each transcript to its ΔTS and classified them as translationally up if ΔTL>ΔTS and translationally down if ΔTL<ΔTS. Thus, we obtained six different combinations of transcript level changes and translational profiles that we applied to transcripts of the non-CPA and CPA subclusters at each developmental stage (Table S7).

At the one-cell stage, the transcriptome consisted mostly of transcripts with constant transcript levels and modest changes in translation (<2-fold). As development proceeded, a decreasing correlation between ΔTS and ΔTL with modest changes in the former and more prominent changes in the latter was observed (Pearson's R=0.63 at the one-cell stage and R=0.21 at the 128-cell stage). This contrasts with what was observed at 5.3 hpf, in which a much larger correlation was found between ΔTS and ΔTL (Pearson's R=0.86). In agreement with the dynamics of the global polysome association (Fig. 2A), the two maternal subclusters had opposing profiles. The CPA subcluster exhibited a shift towards an increase in translation (Pearson's R=−0.005 at the 128-cell stage) and the non-CPA subcluster shifted toward a decrease in translation (Pearson's R=0.036 at the 128-cell stage; Fig. 2C-I). The peak of translational regulation appeared to occur at the 128-cell stage (Fig. 2H,I), in which the correlation between ΔTS and ΔTL was the lowest for both clusters (Fig. 2E). At this stage, the translationally increasing group represented 2464 of 3667 (67%) CPA subcluster transcripts, whereas the translationally decreasing group represented 2972 of 5487 (54%) non-CPA subcluster transcripts (Table S7). Interestingly, a considerable number [2079 (38%)] of non-CPA transcripts also translationally increased (Table S7), suggesting that the translation of these transcripts continued during embryogenesis despite their activation in an earlier CPA event prior to fertilization (Oh et al., 2000; Potireddy et al., 2006; Salles et al., 1992). At 3.5 hpf, fewer CPA cluster transcripts [1995 (54%)] were translationally upregulated, possibly reflecting a decline of the previous trend towards an increase in translation, which could be attributable to the initiation of maternal mRNA degradation in this population. Altogether, these observations suggest that maternal mRNAs that undergo CPA also experience positive translational regulation. This was manifested by an increase in the polysome association without a corresponding increase in the level of transcript. Our findings indicated that CPA is a translational regulatory mechanism of maternal transcripts that occurs exclusively during pre-MBT development.

To elucidate the biological significance of CPA in embryonic development, we used Cordycepin (3′-dA) to block the polyadenylation of maternal mRNAs. In the embryo, 3′-dA is converted to 3′-dATP, which lacks a 3′-OH group. Thus, its incorporation into the 3′ end of an mRNA molecule inhibits further elongation of its poly(A) tail (Aoki et al., 2003). One of the hallmarks of the MBT is gastrulation, characterized by the initiation of epiboly (Kane and Kimmel, 1993). We previously reported that the inhibition of CPA during pre-MBT development resulted in epiboly arrest and cytolysis of the embryo shortly after MBT [n=40; Fig. 3A-G; (Aanes et al., 2011)]. This phenotype was reminiscent of the phenotype of maternal mutants that affected epiboly (Kane et al., 1996; Wagner et al., 2004). To expand our analysis on the role of CPA in the MBT, we performed a shorter pulse of 3′-dA treatment (128-cell to 3.5 hpf). This treatment did not cause developmental arrest, but epiboly was delayed and embryonic patterning was significantly affected (n=30; Fig. 3H-K). This suggests that immediately prior to the MBT, development still depends on the CPA of maternal transcripts. Treatment from the one-cell stage up to 16-cell (n=40) or 128-cell (n=40) stages resulted in epiboly arrest and early mortality that is indistinguishable from that of one-cell to 3.5 hpf treatment, implying that possible effects of the latter treatment on the polyadenylation of zygotic transcripts are minimal to none (Fig. S6D-G).

We hypothesize that the translation of maternal mRNAs that undergo CPA is crucial for MBT. We treated embryos with cycloheximide (CHX) to inhibit translational elongation (Obrig et al., 1971) starting at different stages and ending at 3.5 hpf. CHX treatment starting from the one-cell stage blocked development prior to epiboly. These embryos contained fewer cells of uneven sizes (n=20; Fig. 3L), suggesting a failure of cytokinesis. Furthermore, as development proceeded, an acellular space formed between the blastoderm and yolk (Fig. 3M). This phenotype resembled cytokinesis-deficient maternal mutants (Dosch et al., 2004; Pelegri and Mullins, 2004; Yabe et al., 2007). Similarly, treatment that began at the 16-cell stage resulted in embryos with larger blastomeres at 3.5 hpf (n=61; Fig. 3N, Fig. S6). At later stages, these embryos remained arrested at the high-oblong stage (Fig. 3O,P). Cycloheximide treatment initiated closer to the MBT (128-cell stage) resulted in epiboly arrest followed by mortality, which resembled 3′-dA-treated embryos (n=66; Fig. 3Q-S). All of the CHX treatments caused mortality before 24 hpf. This suggests that similar processes were affected, supporting the hypothesis that the translation of CPA transcripts is necessary for MBT.

To characterize the phenotype resulting from 3′-dA treatment (one-cell stage to 3.5 hpf) in more detail, we looked at structures that indicate the progression of gastrulation. The yolk syncytial layer (YSL) is formed immediately before gastrulation (Kimmel and Law, 1985). The membrane of the YSL nuclei at the outermost periphery of the embryo [external YSL (e-YSL)] is connected to the overlying blastomeres by tight junctions (Betchaku and Trinkaus, 1978; Solnica-Krezel and Driever, 1994) and provides the driving force for epiboly (Behrndt et al., 2012). The e-YSL nuclei were clearly formed in 3′-dA-treated embryos at 3.5 hpf (n=5; Fig. 4A,D), but they were not visible at 5.3 hpf (n=5; Fig. 4B,E). Based on the expression of the YSL marker camsap2a2 at 6 and 8 hpf, the e-YSL was absent, although a weak signal from the inner YSL was still present (n=8; Fig. 4G-J). This suggests a delay in the progression of the deep cells compared with the overlying blastomeres. Therefore, 3′-dA treatment inhibited the progression of epiboly.

By 8 hpf, 3′-dA-treated embryos exhibited developmental arrest (Fig. 4C,F) and underwent cytolysis between 9 and 10 hpf. To determine the extent to which gastrulation proceeded in these embryos, we searched for the presence of dorsal shield (i.e. one of the earliest structures formed at gastrulation). In 3′-dA-treated embryos, no defined structure resembling the dorsal shield was observed, and the expression of the dorsal mesoderm marker (prechordal plate) goosecoid (gsc) was abnormal (n=8; Fig. 4K,L). Furthermore, the expression of wnt8a, which defines the ventrolateral mesoderm (Baker et al., 2010), was absent in 3′-dA-treated embryos (n=8; Fig. 4M,N). Together with the absence of the dorsal organizer, this suggests the lack of dorsoventral patterning and ventral mesoderm specification. Collectively, our findings showed that inhibition of CPA by 3′-dA resulted in a failure of developmental progression through the MBT, indicated by the absence of early hallmarks of gastrulation.

To elucidate the molecular events that underlie the control of MZT by CPA, we performed gene expression analysis using microarrays. We inhibited CPA by 3′-dA beginning at the one-cell stage until 3.5 hpf (long exposure), from the 16-cell stage until 3.5 hpf (medium exposure) and from the 128-cell stage until 3.5 hpf (short exposure), and assayed changes in gene expression profiles at 3.5 hpf. We identified genes that were consistently regulated (either upregulated or downregulated) in all three treatments because these are considered the most sensitive to 3′-dA treatment. This analysis consisted of a total of 3082 genes (2106 upregulated genes and 976 downregulated genes; P<0.05, fc≥2 in at least one treatment; Table S8). The downregulated group contained a large number of genes encoding transcription factors and were significantly enriched in GO terms related to transcription regulation, dorsoventral pattern formation (e.g. nog1, chd, gsc, ndr1, ndr2 and bambia), cell migration in gastrulation [e.g. mxtx2, apela, aplnra, cxcr4a, bon (mixl1) and rock2a] and germ layer formation [e.g. gata5, bon, wnt11, sox32, tbx16, ntla (ta), bmp2b, vent and oep (tdgf1); Table S9]. The downregulation of these functions agrees with the failure of gastrulation that was observed as a result of 3′-dA treatment.

Genes that were upregulated by 3′-dA were enriched in GO terms for DNA repair and chromatin modification, and included several components of chromatin and suppressive chromatin modifiers, such as ehmt2, hdac9b, hdac3, hdac4, setdb1a and smarcb1a (Table S10). The upregulation of these suppressive factors suggests a failure of their clearance. In order to confirm this, we looked at two degradation mechanisms for maternal mRNA clearance that were identified in zebrafish (Giraldez et al., 2006; Zhao et al., 2017). To identify possible interactions between CPA and these two degradation mechanisms, we searched for targets of each of these two degradation pathways among the genes that were affected by CPA inhibition. The targets of zygotically expressed miR-430 were identified as transcripts that were upregulated in MZdicer mutants and had the recognition sequence for miR-430 (Giraldez et al., 2006). Ythdf2, which is maternally expressed, recognizes the N⁶-methyladenosine (m⁶A) modification of RNAs that mediates their decay. Targets of this pathway were identified as those mRNAs that increased their stability in response to Ythdf2 loss of function (Zhao et al., 2017). We identified 580 genes in common between 3′-dA-upregulated genes and Ythdf2 targets; 33 of these were also targets of miR-430. Additionally, 25 genes upregulated by 3′-dA were miR-430 targets, although they were not Ythdf2 targets (Table S11). This observation supports the role of CPA in maternal mRNA clearance by both maternal and zygotic modes of degradation. Nevertheless, other degradation mechanisms also exist in the zebrafish (Bazzini et al., 2016), which may account for the rest of the maternal transcripts.

We then determined the expression cluster to which the 3′-dA-responsive genes belonged by intersecting them with the total RNA-seq clusters (Fig. 5A). A total of 1457 genes were in common between the two datasets. Among these, a large number of genes that were upregulated by 3′-dA were maternal [1090 out of 1140 (96%)], whereas very few genes were of zygotic origin [50 out of 1140 (4%)]. A total of 174 of the downregulated genes (55%) were of zygotic origin (zygotic and pre-MBT clusters) and 143 (45%) were maternal. Among the downregulated maternal genes, 83 of 143 (58%) had a zygotic contribution, as evidenced by the increase in their transcript level at 3.5 or 5.3 hpf in RNA-seq data. Therefore, the downregulation of these maternal-group genes may reflect failure of the activation of their zygotic counterpart. These observations are consistent with the scenario in which there is a concurrent failure of maternal transcript degradation and zygotic transcription initiation. Downregulated zygotic genes were enriched in similar GO terms to those of the overall microarray set of downregulated genes (Table S12), reflecting the earliest biological functions activated at the MBT, which includes gastrulation, cell fate commitment, convergent extension and the development of different germ layers. These interconnected processes are depicted in a network that shows the interactions between different functions (Fig. 5B).

Interestingly, several members of the core histone family, especially from the H2A family, were found among the downregulated pre-MBT genes (Table S13). Upregulated non-CPA genes presented similar GO term enrichment as the overall set of upregulated genes, including suppressive chromatin modifiers (Table S14). The observation of upregulated genes of maternal origin likely resulted from a failure of the transcriptional activation of maternal mRNA degradation components. To a lesser extent, some genes of maternal origin were downregulated as a result of CPA inhibition. These included pou5f3, a maternal (non-CPA) transcription factor implicated in ZGA (Lee et al., 2014). Upregulated CPA subcluster genes did not exhibit enrichment of any functional categories. Nonetheless, this group included a notable number of genes involved in protein degradation, including several ubiquitin ligases and their associated components (Table S15).

Altogether, our analyses indicated that the timely activation of maternal mRNA translation by CPA underlies major events associated with the MZT, namely the activation of expression of key zygotic genes that are required for gastrulation and germ layer formation, the degradation of maternal mRNAs and epigenetic adjustments to allow transcriptional activation.

Although more than 95% of the zebrafish genome has been assembled, its annotation is far from complete. A large fraction of our RNA-seq reads mapped to regions outside of those annotated as coding or noncoding transcripts (Fig. S5). We identified a total of 2176 novel transcribed fragments (transfrags), representing 1080 unique genomic loci based on combined RNA-seq reads across all samples. These possessed identifiable ORF characteristics that were at least 100 amino acids long. A considerable fraction exhibited junction-spanning reads, forming multi-exonic transcripts. To validate and annotate the novel identified loci, peptide sequences originating from the longest transfrags of each locus were aligned to the UniRef90 protein database to search for sequence homology to existing annotated transcripts. We identified 1922 transfrags (88%) that originated from 896 genomic loci (83%) (E-value threshold 1e-5). Taking the highest scoring match of each alignment, this corresponded to 878 unique UniRef90 entries, among which 570 were described as ‘predicted’ and 260 were described as ‘uncharacterized’ transcripts. Taxonomically, these homologous entries originated from 95 different species, the most prominent of which was Danio rerio (67.5%; Fig. 6A). In a parallel analysis, 1566 transfrags (72%) that originated from 770 loci (71%) corresponded to 1977 unique Pfam entries. Among these, the C2H2-type zinc-finger domain was the most prevalent (Fig. S4). Of these novel loci, we further identified 1952 transfrags (89.7%) that had at least one EST support. Altogether, a total of 2021 novel transfrags, representing 971 genomic loci, were annotated. Among these, 1740 transfrags (86%) and 782 genomic loci (80.5%) had both Uniref and EST support (Tables S16 and S17). To assess whether these novel transcripts had an expression profile distribution that was similar to the set of annotated transcripts, we calculated the total normalized read counts for novel genomic loci in each sample and performed clustering analysis. To our surprise, the clustering performed using only the set of novel transcripts presented an almost identical pattern to that of the annotated genes (Fig. 6B). The novel transcripts could also be clustered based on their expression pattern into maternal, pre-MBT and zygotic clusters, thus supporting their authenticity and biological relevance (Fig. 6C).

The high number of novel transfrags that possessed considerable sequence homology to known proteins and recognizable protein domains suggested that they were likely to be authentic and did not result from sample contamination. This also revealed that a considerable number of transcripts were missing from the current annotation, despite having EST support. One likely possibility is that these novel transcripts contain repetitive sequences that caused them to be filtered out during the data analysis. Therefore, our findings reiterate the need for a refinement of analytical strategies to increase the sensitivity of transcript identification. The discovery of novel transcripts would contribute to improving zebrafish genome annotation. More importantly, our datasets were generated without poly(A) tail bias, which enables the accurate detection of transcript levels, particularly in the case of maternal transcripts that are subjected to CPA at pre-MBT stages of development.

The polysome association analysis of the early embryonic transcriptome showed that maternal mRNAs that underwent CPA were increasingly associated with polysomes during pre-MBT stages, and their translation increased despite the lack of transcriptional input. One caveat of this analysis was that detection of the polyadenylation state relied on oligo d(T) selectivity based on the increasing length of the poly(A) tail. To define the occurrence of CPA with high confidence, we employed two strategies. First, we imposed a strict differential expression cutoff (see supplementary Materials and Methods) to define the increase in detected transcripts between two consecutive stages. Second, we cross-validated the oligo d(T)-generated dataset with another dataset from identical developmental stages that was generated without poly(A) bias. These strategies allowed us to obtain a robust gene set that could be directly validated by means of several methods (see Aanes et al., 2011) (Fig. S2).

Although translational activation through CPA is known to play an important role in oogenesis (Oh et al., 2000; Potireddy et al., 2006; Salles et al., 1994), its role in embryogenesis has been less appreciated, especially within the context of the MZT. We found that the CPA of a large cohort of maternal mRNAs was essential for the MZT in zebrafish. Cytoplasmic polyadenylation may regulate the MZT in at least two possible ways. First, it facilitates ZGA. Second, it is required for the activation of maternal mRNA degradation. Supporting the first mechanism, many zygotic genes were downregulated in response to pan-embryonic CPA inhibition, implying a failure of their transcriptional activation upon ZGA. Concurrently, upon CPA inhibition, many maternal mRNAs were upregulated, including targets of the two known maternal mRNA degradation mechanisms (Giraldez et al., 2006; Zhao et al., 2017). Defects in maternal mRNA clearance have been shown to affect early developmental events, including those that are associated with the MZT (Yartseva and Giraldez, 2015). Interestingly, Ythdf2 (Zhao et al., 2017) is cytoplasmically polyadenylated during pre-MBT stages, suggesting that its proper function in maternal mRNA degradation relies on translational regulation by CPA. Additionally, miR-430 is expressed zygotically, the activation of which likely depends on CPA. Therefore, one possibility is that CPA regulates the MZT, at least partially by ensuring the clearance of maternal mRNAs by these two degradation pathways.

Among the maternal factors that are known to play a role in zebrafish ZGA (Lee et al., 2013; Onichtchouk et al., 2010; Xu et al., 2012), nanog and pou5f3 were not subjected to CPA during embryonic development. However, sox19b falls within the CPA subcluster. Moreover, it exhibited a high degree of translational activation by the 128-cell stage, strongly suggesting that sox19b is under positive translational regulation during pre-MBT development. Based on loss-of-function analysis of these three factors, a combination of factors rather than the individual action of their activities likely contributes to ZGA (Lee et al., 2013). The stoichiometry between Sox2 and Oct4 (i.e. orthologs of Sox19b and Pou5f3, respectively) is known to affect the reprogramming of somatic cells into induced pluripotent stem cells (Carey et al., 2011; Papapetrou et al., 2009). Moreover, nanog, pou5f3 and sox19b, which are well known to be regulators of pluripotency, physically interact with each other, and the dynamics of their interactions correspond to the progression of cell differentiation (Perez-Camps et al., 2016; Rizzino, 2013). Whether the increase in Sox19b translation during development represents a similar mechanism for combinatorial activation of these factors in ZGA remains to be experimentally proven.

Altogether, our findings complement the currently accepted models of ZGA mechanisms by demonstrating a temporal component, in which CPA activates a subset of maternal mRNAs that are required for the MZT (Fig. 7A). Cytoplasmic polyadenylation might be considered an essential part of the mechanism that drives the MZT, which would expand the definition of MZT to encompass all processes that are essential for the transition of developmental control from maternal to zygotic. These would include the CPA of maternal mRNAs and other associated functions that are activated by this mechanism as a part of the MZT in addition to ZGA and maternal mRNA clearance (Fig. 7B).

Proteomics of the zebrafish egg and early embryo still poses technical challenge due to the lack of sensitive techniques and the presence of high molecular weight proteins in the yolk. Thus, there is insufficient knowledge of the factors that are provided maternally as fully functional proteins that might additionally contribute to the regulation of early embryogenesis. Polysome profiles that were identified in the present study provide crucial information on proteins that are synthesized from maternal factors that are deposited as mRNAs. Additionally, the identified novel transcripts are a valuable source of candidate regulators, the detailed functional analysis of which is currently under way. Complementing these data with knowledge of the maternal proteome that shapes the egg and early embryo is a next step towards a more comprehensive understanding of the control of the MZT.

Wild-type (AB) zebrafish were maintained in the IMCB zebrafish facility, Singapore. All animal experiments were performed according to the regulations of the Institutional Animal Care and Use Committee (Biological Resource Center, Biopolis). Embryos were maintained and staged according to standard protocols (Kimmel et al., 1995; Westerfield, 2000).

Embryos were collected at one-cell (0 hpf), 16-cell (1.5 hpf), 128-cell (2.5 hpf), high (3.5 hpf) and 50% epiboly (5.3 hpf) stages. Unfertilized eggs were harvested by squeezing the abdomens of spawning female zebrafish. Three hundred eggs or embryos from each stage were immersed in embryo medium containing 200 µg/ml cycloheximide for 10 min, followed by lysis in 500 µl lysis buffer (Sampath et al., 2011). The first four polysome fractions that were pooled as ‘polysome-unbound’ and the subsequent four polysome fractions that were pooled as ‘polysome-bound’ were precipitated with 2.5 ml of 100% ethanol for RNA isolation. Validation of expression dynamics in fractionated samples was performed by quantitative reverse transcription polymerase chain reaction (qRT-PCR) of representative genes in polysome-bound and -unbound fractions (Fig. S2).

RNA was isolated using RNEasy mini kit (Qiagen) with two additional rounds of column-washing steps. Total RNA was rRNA-depleted using Ribo-Zero Magnetic Gold Kit (Human, Mouse, Rat; Epicenter) and quantified using Qubit (Life Technologies). Its quality was assayed using RNA 6000 Nano chip on Agilent 2100 Bioanalyzer. Stranded RNA-seq libraries were generated using ScriptSeq v2 RNA Library Preparation Kit (Illumina) and sequenced on Illumina HiSeq2500 for 75 cycles to a depth of 100 million reads per sample. For detailed bioinformatics analysis methods and discovery of novel transcripts, see supplementary Materials and Methods.

Gene expression profiling was performed on a custom zebrafish array (Agilent Technologies, GIS V2) as previously described (Winata et al., 2013). Data analysis was performed in Genespring (Agilent Technologies). The statistical significance of changes in gene expression was determined using one-way analysis of variance (ANOVA); threshold values were set as P<0.05 and fold-change ≥2. Genes were annotated using BioMart (Smedley et al., 2009).

3′-dA treatment was performed as previously described (Aanes et al., 2011). Briefly, embryos were placed in glass Petri dishes containing 4 mM 3′-dA (Sigma-Aldrich, C3394) dissolved in embryo medium. Cyclohexamide treatment was performed with 20 µg/ml cyclohexamide (ROTH, 8682.3) in dimethylsulfoxide (DMSO). The control for CHX treatment was the corresponding amount of DMSO (0.04%) in embryo medium. To terminate the treatment, embryos were rinsed four times with embryo medium and transferred to a new dish containing fresh egg water. Observations were made of at least three independent experiments.

Immunofluorescent staining was performed according to established protocols (Chu et al., 2012). Confocal imaging was performed using an Olympus FluoView FV1000 (Olympus, Japan). Immunofluorescently labeled embryos were mounted in 0.5% low-temperature melting agarose in embryo medium. Two- or three-dimensional reconstructions of imaging data were performed using the standard Olympus software package.

Whole-mount in situ hybridization was performed as described previously (Winata et al., 2013). Digoxigenin-labeled RNA antisense probes were synthesized from partial cDNA clones of gsc, wnt8a and camsap2a2 (details are in the supplementary Materials and Methods), using DIG RNA labeling mix (Roche, Switzerland) and SP6 or T7 polymerase enzyme (Fermentas, Lithuania). Images were captured with a Leica M205FA camera and Zeiss AxioImager.M2.

We thank D. Solter, B. Knowles, W. Rybski, L. T. Chu and S. Fong for insightful discussions and technical assistance.