{'fraction_non_alphanumeric': 0.0389215214251324, 'fraction_numerical': 0.025825710158883006, 'mean_word_length': 4.985016136468419, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 5, 'https://': 4, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Genome wide analysis of 3′ UTR sequence elements and proteins regulating mRNA stability during maternal-to-zygotic transition in zebrafish

Posttranscriptional regulation plays a crucial role in shaping gene expression. During the maternal-to-zygotic transition (MZT), thousands of maternal transcripts are regulated. However, how different cis-elements and trans-factors are integrated to determine mRNA stability remains poorly understood. Here, we show that most transcripts are under combinatorial regulation by multiple decay pathways during zebrafish MZT. By using a massively parallel reporter assay, we identified cis-regulatory sequences in the 3′ UTR, including U-rich motifs that are associated with increased mRNA stability. In contrast, miR-430 target sequences, UAUUUAUU AU-rich elements (ARE), CCUC, and CUGC elements emerged as destabilizing motifs, with miR-430 and AREs causing mRNA deadenylation upon genome activation. We identified trans-factors by profiling RNA–protein interactions and found that poly(U)-binding proteins are preferentially associated with 3′ UTR sequences and stabilizing motifs. We show that this activity is antagonized by C-rich motifs and correlated with protein binding. Finally, we integrated these regulatory motifs into a machine learning model that predicts reporter mRNA stability in vivo.

pinpoint distinct segments driving regulation. To address these issues, we recently developed a high-throughput RNA-element selection assay (RESA) to measure the regulatory activities of mRNA sequences in vivo (Yartseva et al. 2017). RESA uses endogenous RNA fragments in a parallel reporter assay that allows for high sequence complexity and high-density coverage of the transcriptome, or targeted regions of interest, providing near nucleotide resolution of the regulatory activity of RNA sequences.
RNA regulatory elements are recognized by trans-factors, including miRNAs and RBPs (Glisovic et al. 2008). Once bound, RBPs can regulate the processing, stability, and translation of their target mRNAs (Gerstberger et al. 2014). Interactome capture has revealed a set of proteins in intimate contact with mRNAs across different eukaryotic systems (Baltz et al. 2012;Castello et al. 2012;Kwon et al. 2013;Sysoev et al. 2016;Wessels et al. 2016;Despic et al. 2017). In vitro affinity selection methods, such as SELEX (Blackwell and Weintraub 1990;Ellington and Szostak 1990;Tuerk and Gold 1990), RNA affinity profiling (Tome et al. 2014), and RNAcompete (Ray et al. 2013), have been complemented with UV crosslinking and immunoprecipitation (CLIP) to provide the set of targets and the binding motifs in vivo for a number of RBPs (Ule et al. 2003;van der Brug et al. 2008;Chi et al. 2009;Hafner et al. 2010;König et al. 2010;Chan et al. 2014;Hansen et al. 2015;Murn et al. 2015;Sugimoto et al. 2015;Galloway et al. 2016;Scheckel et al. 2016;Rot et al. 2017). However, the presence of a specific sequence motif is not always indicative of regulation in vivo, suggesting that additional sequences, or combinatorial interactions between RBPs, influence the regulatory output on an mRNA. Current efforts have not yet linked the RNA regulatory maps with the RBP binding profiles to define the posttranscriptional regulatory network in vivo.
Identifying functional regulatory sequences together with RBP trans-factors is an essential step toward understanding mRNA posttranscriptional regulation. This regulation is of particular importance during the early stages of animal development, which are instructed by maternally provided mRNAs (Walser and Lipshitz 2011). During the maternal-to-zygotic transition (MZT), mRNAs deposited in the oocyte undergo coordinated remodeling. Individual pathways have been implicated in the regulation of maternal mRNAs (for review, see Lee et al. 2014;Yartseva and Giraldez 2015). For example, in Drosophila, the RBP SMAUG destabilizes maternal mRNAs (Dahanukar et al. 1999;Tadros et al. 2007). In Xenopus, AREs within the 3 ′ UTRs of maternal mRNAs trigger their deadenylation after egg activation and decay after the mid-blastula transition (MBT) (Audic et al. 1997;Voeltz and Steitz 1998). In zebrafish, zygotic transcription of microRNA miR-430 regulates ∼20% of destabilized maternal transcripts (Giraldez et al. 2006). Codon usage in the coding sequence of mRNAs influences differential mRNA stability during the MZT across several vertebrates (Bazzini et al. 2016;Mishima and Tomari 2016), and it influences mRNA half-life in yeast (Presnyak et al. 2015). mRNA methylation has been implicated in shaping mRNA stability during ES cell differentiation (Batista et al. 2014) and in MZT (Zhao et al. 2017). Despite these singular discoveries and the recent effort to develop models based on 3 ′ UTR sequence elements (Rabani et al. 2017), there is a paucity of quantitative models that integrate the various elements to predict mRNA deadenylation, degradation, and RBP binding. Thus, it is still poorly understood how sequence composition influences mRNA stability, how the sequence context affects the regulatory potential of each motif, what are the relative activities of different elements, and which RBPs mediate this regulation.
Here, we identified the RBPs bound to mRNAs during MZT, their preferential binding sequences, and the 3 ′ UTR cis-regulatory sequences shaping mRNA stability during zebrafish embryogenesis, which we have integrated into a quantitative model that explains the regulatory activity encoded within individual 3 ′ UTRs.

Results
Distinct pathways regulate maternal mRNA decay Fertilization triggers remodeling of the transcriptome required for the first steps of embryonic development. Maternal and zygotic posttranscriptional pathways regulate the dynamics of mRNA stability during MZT (Yartseva and Giraldez 2015). To identify which pathways regulate each mRNA, we analyzed transcript levels during the first 8 h postfertilization (hpf), at 30-to 60-min intervals in wild-type zebrafish embryos using mRNA-seq (Supplemental Table S1). We identified mRNAs whose decay was dependent on zygotic transcription or specifically on miR-430. These transcripts were significantly stabilized when zygotic transcription was inhibited with the RNA Pol II inhibitor, α-amanitin (Lindell et al. 1970;Kane et al. 1996), or when miR-430 was inhibited using an antisense tiny-LNA complementary to miR-430 (LNA 430 ) (Staton et al. 2013). The remaining transcripts that underwent decay when zygotic transcription was inhibited were classified as regulated by the maternal mode. To normalize RNA expression across stages, we used exogenous yeast spike-in RNA, allowing us to quantify global changes in mRNA levels ( Fig. 1A; Supplemental Table S2). Comparing early (2 hpf) to late (6 hpf) developmental stages, we defined the main regulatory mode for 5847 mRNAs undergoing decay: 3909 were regulated by the zygotic mode (67%), of which 616 were primarily dependent on miR-430 (11% of total unstable) and 1938 were primarily regulated by the maternal mode (33%) (Fig.  1B). In situ hybridization analysis of endogenous transcripts selected from each mode displays stabilization patterns in the absence of zygotic transcription (α-amanitin) and/or miR-430 function (LNA 430 ) that were consistent with behaviors observed in the global RNA-seq analyses (Fig. 1C).
Although we were able to identify the predominant mode of decay for each unstable mRNA, the stability of any individual mRNA is likely dictated by multiple regulatory mechanisms operating within the embryo. We identified multiple regulatory modes for numerous mRNAs. To dissect potential combinatorial regulation, we calculated the degree of stabilization conferred on individual mRNAs by the loss of miR-430 regulation, or zygotic transcription ( Fig. 1D; Supplemental Table S3). We find that the majority of unstable transcripts were degraded by the combinatorial regulation of multiple pathways (Fig. 1E). For example, we observe that up to 3688 mRNAs (63%) are partially stabilized after blocking miR-430 function. Together, these results define three regulatory modes of maternal mRNA turnover and their relative contributions to shaping posttranscriptional regulation across the zebrafish transcriptome after fertilization.

RESA identifies regulatory RNA elements
We hypothesized that the combinatorial regulation observed for maternal transcripts might be encoded in discrete elements within the mRNA. To identify these elements, we used RESA ( Fig. 2A; Yartseva et al. 2017), which assesses the ability of sequences to regulate mRNA stability in vivo when placed in the 3 ′ UTR of a reporter mRNA library. We generated two RNA reporter libraries (Supplemental Table S1), one composed of random fragments spanning the entire embryonic transcriptome (∼30 nt length; transcriptome library) and a high-density library enriched for 3 ′ UTRs from 434 genes regulated during the MZT (∼100-nt length; targeted library). To identify regions that mediate differential mRNA stability, we injected each RNA reporter library into one-cellstage embryos and quantified the abundance of each reporter within the RESA library using RNA-seq before and after zygotic transcription (2 and 6 hpf) in wild-type, α-amanitin, or LNA 430 -injected zebrafish embryos. Depletion or enrichment of sequences over developmental time revealed 1404 destabilizing and 295 stabilizing regions, respectively, across 3456 genes. Destabilizing regions were modulated by the maternal (593), zygotic (184), or miR-430 (627) modes of mRNA decay (Supplemental Fig. S1A). During early embryogenesis, deadenylation and decay are uncoupled (Voeltz and Steitz 1998), allowing us to distinguish between elements that predominantly cause deadenylation or decay. By inferring the relative poly(A) tail length from capture efficiency (i.e., depletion of reads from poly(A) selection over time relative to total reads indicates deadenylation) (Beilharz and Preiss 2007;Bazzini et al. 2012), we identified several hundred regions predominantly causing deadenylation (Supplemental Fig. S1A).
To test the regulatory activity of these sequences, we validated two reporter mRNAs containing 3 ′ UTR sequences identified by RESA using qRT-PCR (Fig. 2B,C). Each reporter was destabilized in wild-type embryos, and this effect was blocked when inhibiting zygotic transcription, consistent with the specific regulation of these regions by the zygotic mode. We observed that the mean destabilization across miR-430 target sites corresponded to the predicted microRNA target site strength (8-mer > 7-mer > 6-mer), whereas inserts containing reverse complement miR-430 target sites were not depleted (Fig. 2D), showing that RESA can accurately quantify regulatory strength across target sites. Together, these results indicate that RESA can identify several hundred regions that promote mRNA deadenylation and decay across the transcriptome.

Identifying destabilizing and stabilizing regulatory motifs
To identify short-linear motifs enriched in coregulated sequences, we used Find-A C D E B Figure 1. Distinct pathways regulate maternal mRNA decay. (A) Hierarchical protocol followed to separate genes into their predominant mode of decay. Among decaying genes (green) determined as stable early (2 hpf) or higher in absence of zygotic transcription (α-amanitin treatment), miR-430 decaying genes were first characterized using LNA 430 treatment (dark blue). Remaining genes, as indicated by the black tree, were split into maternal (purple) and zygotic (light blue) decaying genes based on their dependence on zygotic transcription. Transcripts with nonsignificant expression changes are colored in light gray (see DESeq2 analysis in Methods). (B) Biplots representing mRNA expression levels over developmental time separated by decay pathways. Time-course is shown in WT, α-amanitin, and LNA 430 conditions. Each line represents expression level of individual mRNA. (C) In situ hybridization illustrating three genes from the maternal, zygotic, and miR-430-dependent decay modes (n = 20/20). (D,E) Relative contribution of the three modes to RNA decay. Each gene total decay was assessed using fold-change between untreated early and late stages calculated in A. LNA 430 and α-amanitin conditions were used to attribute the contribution of the miR-430 and zygotic modes, respectively. Remaining decay was attributed to the maternal mode. This normalization was taking into account that α-amanitin treatment is also blocking miR-430 transcription. For example, if a transcript's decay profile was not affected by blocking miR-430 or zygotic transcription, then that transcript's degradation was wholly attributed to maternal decay mode. In contrast, decay that was prevented by blocking miR-430 or zygotic transcription was attributed to miR-430 or zygotic decay pathways, respectively. As multiple comparisons, used in A to determine each gene mode of decay, can be significant, the main (lowest P) and significant (all significant P) modes of decay are shown (E).
ing Informative Regulatory Elements (FIRE) (Elemento et al. 2007;Oikonomou et al. 2014). This method analyzes all possible 7-mers to then optimize these seeds into sequence logos by maximizing mutual information (Elemento et al. 2007;Oikonomou et al. 2014). We identified motifs associated with destabilization and stabilization in the three decay modes (confidence cutoff Z-score >20 and >10 for the transcriptome and targeted libraries respectively, see Methods) ( Fig. 3A; Supplemental Fig. S2). These motifs were consistent with those identified across mRNAs with different kinetics of decay (Rabani et al. 2017). We find that miR-430 seed target sequences were specifically enriched in unstable sequences that were stabilized by the loss of miR-430 function, confirming that our strategy identified known motifs in the expected mode of regulation. Independent of miR-430, the motifs most strongly associated with unstable sequences were CCUCCNC and CUGCNC (Z-score 144.1 and 77.8). To validate these findings independently from RESA, we analyzed the expression of a reporter mRNA containing multiple copies of the CCUC motif derived from the gene zc3h18 (that was in common between the transcriptomic and targeted libraries). We observed decreased stability for CCUC motifs compared with a reporter in which those motifs were mutated to CGUC, resulting in lower protein expression of a GFP reporter mRNA compared with a dsRed control ( Fig. 3B-D). These results suggest that CCUC motifs are responsible for sequence-specific decay elements in the early zebrafish embryo. FIRE also identified U-rich motifs such as UUUNUUU, UUUUANA, or AANUAUU overrepresented in stable regions with UUUNUUU displaying the strongest activity (Supplemental Fig. S2). Stabilization at multiple U-rich sites was consistently observed within individual transcripts (Fig. 3E). In contrast, RESA samples prepared with additional poly(A) selection revealed that deadenylated sequences were enriched in miR-430 sites and U[CA]UUUAUU, an ARE known to be involved in regulating poly(A) tail length ( Fig. 3F; Supplemental Fig. S2; Audic et al. 1997;Voeltz and Steitz 1998;Rabani et al. 2017;Yartseva et al. 2017). To validate this motif, we analyzed the poly(A) tail length of a reporter mRNA containing tandem copies of UAUUUAUU derived from the gene dennd6b (also known as fam116b). This analysis revealed that poly(A) shortening by AREs was dependent on activation of the zygotic genome, as it was abolished in α-amanitin-treated embryos (Fig. 3G).
To measure the regulatory activity of each identified motif, we analyzed the average decay across all the regions containing that motif in the RESA-targeted assay before (2 hpf) and after MZT (6-8 hpf). We performed two different assays: (1) We compared their effect in the total and the poly(A) selected library to distinguish effects on stability and deadenylation, and (2) we analyzed their dependency on the zygotic mode of decay by determining whether their regulatory activity was blocked in the absence of zygotic transcription (α-amanitin-treated samples). The U-rich motif showed a stabilizing effect (mean stabilization of 0.17 log 2 foldchange (log FC ) at 2 vs. 6 hpf) (Fig. 3H), whereas CCUCCNC had a destabilizing effect (mean destabilization −0.22 log FC 6/2 hpf) (Fig. 3I). Regulation by these elements was mainly controlled through the maternal mode on the total mRNA as it was still observed in α-amanitin-treated samples but had a weaker effect on the poly(A) libraries. Analysis of lower Z-score FIRE motifs CCUGC and UUAUU and Pumilio motifs (UGUA[AU]AUA) (Gamberi et al. 2002;Gerber et al. 2006) revealed weak regulation acting on total and poly(A) libraries (Supplemental Fig. S1B-D). As a control, meta-analysis over 123 miR-430 seeds (Fig. 3K) showed a strong zygotic-dependent regulatory effect on stability (−0.5 log FC 6/2 hpf) and deadenylation (−1.15 log FC 8/2 hpf), with an effect that was more than 2.5 times stronger than U-rich and CCUCCNC motifs. Consistent with the reporter analysis, meta-analysis of the U[CA] UUUAU ARE motif over 50 loci revealed that this motif causes zygotic-dependent deadenylation (−0.4 log FC 8/2 hpf) with a robust depletion of polyadenylated fragments, without a significant effect on the total mRNA abundance within the time frame analyzed ( Fig. 3J). Relative regulatory activity of these elements revealed that within this context, miR-430 provides stronger regulation than ARE motif UAUUUAUU, with a rapid coupling of deadenylation and decay that is not observed in AREs. Together, these results identify sequence motifs that regulate reporter mRNA stability and deadenylation during the MZT.

Mapping RNA-interacting proteins in the embryo
To identify the trans-factors binding mRNAs during the MZT, we adapted the interactome capture technique to zebrafish embryos (Baltz et al. 2012;Castello et al. 2012;Kwon et al. 2013;Sysoev et al. 2016;Wessels et al. 2016;Despic et al. 2017). This method uses UV to crosslink protein-mRNA interactions and then poly(A) purification followed by mass spectrometry to identify the proteins bound to mRNA ( Fig. 4A; Supplemental Fig. S3). We analyzed the interactome across three independent biological replicates with or without UV crosslinking. We queried the interactome at 4 hpf, a timepoint that precedes the widespread changes in maternal mRNA stability and deadenylation characterized here (Bazzini et al. 2012). By using label-free quantitative mass spectrometry, we identified 160 proteins with two or more peptides with at least one unique peptide. Of the 160 identified proteins, 112 were significantly enriched in the UV-crosslinked sample compared with controls (-UV) (Fig. 4B). From this core set of 112 proteins, 90 had also been identified in previous mRNA interactomes (Supplemental Table S4; Baltz et al. 2012;Castello et al. 2012;Kwon et al. 2013;Liao et al. 2016;Despic et al. 2017). We observed a significant enrichment of RBPs in the interactome, as 67.8% of the identified proteins had annotated RNA-binding domains compared with 8.2% of proteins detected in the input. Further, our analysis revealed an enrichment for Gene Ontology (GO) terms associated with RBPs (Fig. 4C) and a selective enrichment for bona fide RNA-protein interactions (Supplemental Fig. S4). These included RBPs involved in RNA processing and splicing such as Xrn2, Srsf2, or Celf1; mRNA translation such as Eif4a1a, Eif2a, and Eif4enif1; and RNA stability such as Pum1, Piwil1, Hnrnpd, Khsrp, and Stau2 among others. To validate their ability to bind RNA, we immunoprecipitated six of the identified proteins after UV crosslinking and radiolabeled the bound RNA for detection ( Fig. 4D,E). Although all the proteins interacted with RNA, a subset of these interactions was also dependent on the activation of zygotic transcription (Ythdf1, Khsrp, Khdrbs1), because the levels of RNA pulled down were reduced in α-amanitin-treated embryos, without affecting the RBP levels. Taken together, the zebrafish interactome identified a large set of RBPs that participate in a wide range of RNA processing pathways during MZT, providing an entry point to identify the effector proteins that regulate mRNA stability and decay.

Identifying putative regulatory factors driving mRNA stability
To identify potential effector proteins that mediate the regulatory activity observed in RESA, we performed iCLIP experiments on 24 RBPs identified in the zebrafish interactome and analyzed their target sequences (Supplemental Fig. S5; Supplemental Tables S1, S5). We reasoned that having a common tag would allow us to compare the signal between different proteins to identify specific binding events for each protein and distinguish them from background common to all samples. Thus, we analyzed FLAG-tagged versions of each protein expressed from an injected mRNA. To ascertain that this approach captures bone fide binding sites, we compared iCLIP signal of one candidate, Khsrp, from the pull down of the FLAG-tagged protein versus the endogenous protein using an endogenous antibody. We observed a similar enrichment of both the FLAG-tagged protein and the endogenous Khsrp in the 3 ′ UTR of endogenous mRNAs and similar motifs sharing the core UUUAU (Supplemental Fig. S6), suggesting that the FLAG-epitope does not significantly alter the binding motif of Khsrp. Next, we characterized the binding pattern and the sequence motif preferentially bound by each RBP across replicates (Supplemental Table  S6). Cumulative count of the iCLIP reads within the 5 ′ , 3 ′ UTRs and CDS revealed that the majority of RBPs displayed strong occupancy within the 3 ′ UTR (Fig. 4F). Within this class, we observed variable accumulation of reads in different regions of the 3 ′ UTR.
Among these, Celf1, Elavl1b (also known as Hug), and multiple hnRNPs displayed preferential occupancy toward the distal end of the 3 ′ UTR. In contrast, Pcbp2 binding was more frequently observed directly proximal to and downstream from the annotated stop codon. iCLIP reads from cnbpa (also known as zff9) and purbb were distributed throughout the CDS. In particular, khsrp was preferentially enriched across 3 ′ UTRs in the transcriptome and was preferentially excluded from coding sequences and 5 ′ UTRs. At the exon junctions, most RBP binding was observed within exons and close to acceptor and donor sites similar to control profiles (Supplemental Fig. S7A). Known splicing factors such as Srsf4 had similar binding profiles as observed by Änkö et al. (2012), whereas Khsrp (Supplemental Fig. S7B) and Hnrnpd displayed high intron binding.
To assess the RNA sequence specificity of each RBP, we searched for the top 10 enriched hexamers bound by each RBP compared with a negative control lacking a FLAG epitope (Fig.  4H). For most proteins, the identified motifs were similar and partially overlapping to those previously identified in vitro (Supplemental Fig. S8; Ray et al. 2013). The iCLIP binding pattern for each RBP resembled the distribution of the top identified motifs (hereafter in silico binding) (Fig. 4G), suggesting that for most proteins, the presence of the binding motifs explains the binding distributions observed in vivo. However, we observed (1) higher density of in silico binding in the 5 ′ UTR than observed in vivo for several RBPs and (2) higher density of in silico binding to the CDS for pcbp2 despite similar in vivo and in vitro profiles for pcbp2 within the 3 ′ UTR. These differences raise the possibility that additional factors, such as the ribosome, might contribute to the observed occupancy profiles in vivo (Supplemental Fig. S9), as it has been suggested for UPF1 (Zünd et al. 2013).
We found that several of the RBP binding motifs possess regulatory activity as measured by RESA ( Fig. 4H; Fig. 5D; right panels). For example, Elavl1b, Fubp1, Hnrnpc, and Tia1 were preferentially bound to poly(U) associated in RESA with stabilizing motifs. In contrast, RBPs such as Pcbp2, Ptbp2a, or Cnbpa were preferentially bound to destabilizing motifs. Taken together, these data correlate the regulatory activity of specific RNA motifs to the recognition by specific RBPs during zebrafish embryogenesis.

Modeling the effect of sequence on mRNA stability in vivo
Our analysis identifies several motifs that are enriched within regulated sequences in the mRNA. To determine whether sequence information can be used to model reporter mRNA regulation, we used machine learning and developed a random forest model (Breiman 2001). We reasoned that this model would allow us to capture the association between motif frequencies and their corresponding regulatory activities along RNA sequences (Fig. 5A). We analyzed the frequency of k-mers (1-8 nt) in 100-nt sliding windows across the RESA-targeted library. To build the model, we retained 387 k-mers out of 87,380 with an absolute Spearman's correlation coefficient above 0.1 between k-mer frequency and fold-change measured by RESA. Next, we used fivefold cross-validation to train and assess the performance of the random forest model. This model achieved a 0.68 Pearson correlation between the predicted change in mRNA abundance for each window and the corresponding change measured by RESA with asymptotic P = 0.0 (Fig. 5B,C). The random forest model selected 57 motifs with a dominant effect on RNA stability ( Fig. 5D; left; Supplemental Fig. S11A).
To further analyze the importance of the motifs selected by the random forest model, we compared the differential stability   between windows that contain or lack each motif. This analysis revealed that the top 25% stabilizing windows were significantly enriched in U-rich motifs, as well as UUAG and UGUA motifs (P < 1 × 10 −146 , 8 × 10 −21 , and 9 × 10 −20 , respectively, Mann-Whitney U test with Bonferroni multiple testing adjustment), confirming the trends observed above. On the other hand, the top 25% destabilizing windows were significantly enriched in miR-430 complementary sites and CCUC motifs (P < 2 × 10 −96 and 6 × 10 −95 respectively, Mann-Whitney U test with Bonferroni multiple testing adjustment) (Supplemental Fig. S10A,B). By intersecting these k-mers with the sequences identified in the iCLIP experiment, we identified potential RBPs that could act as trans-factors to regulate bound mRNAs (Fig. 5D, right) and revealed UUAG was significantly enriched among the top 50 k-mers bound by Hnrnpa1b (P < 2 × 10 −73 , χ 2 -test followed by Bonferroni multiple testing adjustment) and UUUU was significantly enriched in k-mers bound by Elavl1b and Hnrnpc (P < 1 × 10 −174 and 2 × 10 −92 , respectively). On the other hand, CCUC was significantly enriched in k-mers bound by Pcbp2 (P < 2 × 10 −23 ), and CUCU was significantly enriched in k-mers bound by Ptbp2a (P < 1 × 10 −30 ). Together these analyses provided a model that predicts the regulatory information encoded in the 3 ′ UTR sequence to regulate reporter mRNA stability in vivo. mRNAs can be controlled by maternally or zygotically encoded modes. To model the different modes of mRNA decay, we built additional random forest models in the absence of zygotic regulation (α-amanitin) or miR-430 regulation (LNA 430 ). We observed a strong correlation between the predicted and the observed regulation in RESA targeted for both models (α-amanitin: r = 0.70 with asymptotic P = 0.0; LNA 430 : r = 0.74 with asymptotic P = 0.0) (Supplemental Fig. S11B,C). Consistent with the FIRE analysis, CCUC and U-rich motifs were identified as regulatory sequences independent of zygotic transcription as part of the maternal model, whereas miR-430 was identified as the main element regulating mRNA stability of the zygotic mode (Fig. 5D). In contrast, models built using the RESA-targeted libraries using poly(A) selection were less accurate (wild type: r = 0.53 with asymptotic P = 0.0; αamanitin: r = 0.33 with asymptotic P = 1 × 10 −248 ; LNA 430 : r = 0.31 with asymptotic P = 6 × 10 −224 ) ( Fig. 5E; Supplemental Fig. S12) and revealed that UAUUUAUU (AREs) motifs were the strongest zygotic-dependent regulator of the poly(A) tail after miR-430. The lower accuracy of the model when the zygotic mode is blocked suggests that regulation of the polyadenylation status by 3 ′ UTR sequence is mainly zygotic dependent. Based on these results, we conclude that our current models capture the importance of CCUC and U-rich motifs to maternally regulate RNA stability and of miR-430 and UAUUUAUU (AREs) to zygotically regulate mRNA polyadenylation.

Predicting stability of endogenous mRNAs
The random forest model captures the association between motif frequencies and stability. Thus, we applied it to predict differential stability of endogenous mRNAs, determined as the fold-change in mRNA levels between 2 and 6 hpf. To this end, we predicted the stability of endogenous mRNAs using the random forest model trained on the RESA-targeted library by averaging the predicted stability of all 100-nt sliding windows of the 3 ′ UTR for each endogenous mRNA. Our model achieved a 0.29 Pearson correlation between the predicted and measured mRNA fold-changes with asymptotic P = 1 × 10 −158 (Supplemental Fig. S10D), indicating that a large fraction of mRNA regulation is not captured by the model de-veloped on the 3 ′ UTR reporter mRNAs, consistent with the model developed by Rabani et al. (2017). This suggests that although the primary 3 ′ UTR sequence is responsible for about a third of endogenous mRNA stability, other factors not captured by RESA such as codon bias (Presnyak et al. 2015;Bazzini et al. 2016;Mishima and Tomari 2016) or RNA structure and modifications (Batista et al. 2014;Ke et al. 2017;Zhao et al. 2017;Beaudoin et al. 2018) could drive the remaining two-thirds (see Discussion).
Among the RESA enriched motifs, several were associated with differential stability or deadenylation of endogenous mRNAs. Endogenous mRNAs containing U-rich motifs were specifically stabilized in total RNA libraries (Fig. 5F). In contrast, we find that UAUUUAUU sequences (AREs) are enriched within endogenous mRNAs that are deadenylated, consistent with a role for AREs in poly(A) tail shortening (Fig. 5G). To test the regulatory activity of endogenous sequences identified by RESA and the motifs identified by the random forest model, we used CRISPR-Cas9 editing to mutate a potential destabilizing sequence element identified in the trip10a gene. We observed that a 223-nt deletion in the endogenous gene, targeting the region regulated in RESA, caused stabilization of the mutant mRNA compared with the wild type, without any apparent developmental phenotype. This shows that RESA can identify de novo functional regulatory elements in vivo (Fig. 5H,I). Sequence analysis of the regulatory region in trip10a 3 ′ UTR revealed multiple AUUUA, AUUA, and AAAUAAA, which, when disrupted, stabilized the mRNA as measured using high-throughput sequencing (Fig. 5J) and increased protein output (Supplemental Fig. S13B-D). We predicted the reporter RNA levels using our model and obtained a 0.60 Spearman's correlation coefficient with RNA levels measured experimentally ( Fig. 5J; Supplemental Fig. S13A). We concluded that a significant part of the regulation is encoded in the 3 ′ UTR sequence and that it could be accurately predicted using machine learning.

Antagonistic effects of different RBP binding motifs
The regulatory effect of any particular motif can be expressed as the average activity across hundreds of loci with that motif. For example, loci containing the Pcbp2-binding motif displayed a mean destabilization of 0.985 across 10,305 loci (Fig. 6A). However, for most motifs, we observed a broad spectrum of regulation and RBP binding, suggesting that the mere presence of a given motif is not the unique determinant of regulatory activity or RPB binding. We hypothesized that the sequence context for each target site might explain differential regulation (RESA) and/or binding (iCLIP). To assess this, we ranked each locus possessing a Pcbp2, miR-430, or Elavl1b-binding motif according to its RESA activity. An analysis of flanking sequences revealed significant enrichment of specific 3-mers when comparing the most regulated versus least regulated loci. For Pcbp2, CCU, CUC, and CCC were significantly enriched within flanking sequences of the most destabilized sequences, which was associated with increased Pcbp2 binding as shown by cumulative iCLIP signal (Fig. 6B,C). In contrast U-rich sequences were significantly depleted within these sites. Conversely, the context of Elavl1b-binding sites was significantly enriched in UUU 3-mers within the most stabilized sequences, which were also among the most abundantly bound sequences by Elavl1b (Fig.  6D-F; Supplemental Fig. S14A-C). miR-430 binding sites did not appear to have a specific sequence bias associated with stronger regulation, consistent with the strength of the regulation being primarily modulated by the size of the seed (Supplemental Fig.  S14D,E). Favorable sequence contexts for Pcbp2 and Elavl1b resemble the original binding site and have opposing nucleotide preferences. These results suggest that U-rich and poly(C)-binding proteins might antagonize each other's activity influencing the stability of the mRNA.
To test the antagonistic effect of U-rich and C-rich motifs on mRNA stability, we first searched for reporters in the RESA library overlapping strictly one or both motifs. Because RESA is composed of a large library of individual reporters, each reporter can be directly assessed instead of averaging the signal across reporters. We analyzed 18 loci with more than 20 RESA reporters overlapping a CCUCC motif, a UUUUUU motif, or both (Fig. 6G, left). We found that the reporters with both motifs had an intermediate stability compared with reporters with each individual motif (P = 2.14 × 10 −3 and 4.29 × 10 −2 for CCUCC and UUUUU compared with both motifs, respectively, Wilcoxon signed-rank test) ( Fig.  6G; right). To validate the effect of these motifs independently from RESA, we analyzed the expression of reporters derived from the ccdc22 gene containing three copies of U-rich and one copy of C-rich motif and vice versa. U-rich reporters were more stable and resulted in higher protein expression compared with the C-rich reporter mRNA (P = 2.90 × 10 −3 , Mann-Whitney U test) (Fig. 6H). Together, these results show that RESA can measure the combined effects of multiple motifs and show that U-rich and C-rich motifs antagonize each other's activity during early embryogenesis.

Discussion
Posttranscriptional regulation plays a major role in shaping gene expression during cellular transitions in which mRNAs from the previous state undergo repression and decay. During the MZT, both maternal and zygotic modes regulate mRNA abundance. Here, we identified (1) cis-regulatory elements using RESA, (2) trans-factor RBPs using interactome capture, and (3) RBP target sequences using iCLIP. We integrated these regulatory sequences into a random forest prediction model, which largely recapitulated . RESA transcriptome coverage ratio and iCLIP signal were averaged over windows centered on RBP motif (RESA minimum coverage >0.05 CPM). The motif is represented with gray bar. SEM of RESA is shown by red shaded outline. (B,E) Sequence context of RBPs potential target sites. Volcano plot representing 3-mer enrichment 20 nt upstream of and downstream from the motif between the top 10% most destabilizing for Pcbp2 (B) and most stabilizing for Elavl1b (E) and the bottom 10%. P-values were calculated using a G-test. Red line indicates 1% significance cutoff after Bonferroni multiple test correction. (C,F) Motif-centered metaplots for Pcbp2 (C) and Elavl1b (F ) comparing potential target sites with favorable context (red) or not (blue). Favorable context defined as top 10% of sites most enriched with significantly enriched 3-mers and least enriched in significantly depleted 3-mers as shown on panels B and E. Scale on y-axis was adjusted, showing higher regulation and wider iCLIP binding signal of sites with favorable context. (G) Schematic describing how to select RESA inserts to measure antagonistic effects of two motifs. Inserts from RESA-targeted overlapping C-motif only, U-motif only, and both motifs are colored in blue, red, and green, respectively. Inserts with a partial overlap in gray are ignored. To estimate stability of inserts, the same procedure is applied for early and late timepoints, and the ratio of the number of inserts is calculated (left). Box-and-whisker plots showing stability of inserts overlapping C-rich and U-rich motifs separately and both together (right). ( * * * ) P < 0.001; ( * * ) P < 0.01; ( * ) P < 0.05, Wilcoxon signed-rank test. (H) Zebrafish 24-h embryos injected with [Crich (blue box) and 3× U-rich motifs (red box)] or [3× C-rich and U-rich motifs] reporters (n = 18 per inserts). Reporter sequence is derived from 3 ′ UTR of ccdc22 gene. In G and H, boxes span first to last quartiles, and whiskers represent 1.5× the interquartile range. ( * * ) P < 0.01, Mann-Whitney U test.
3 ′ UTR-mediated regulation as measured by RESA. These findings integrate important aspects of the posttranscriptional regulation shaping mRNA stability in vivo.
Our implementation of RESA allowed us to measure the relative strength of regulatory sequences in vivo. Among the regions in the 3 ′ UTR, miR-430 targets were the strongest deadenylating and destabilizing elements. UAUUUAUU elements provided a strong deadenylating activity that was dependent on the activation of the zygotic genome, with a weaker effect on mRNA decay. On the same mRNA, RESA identified multiple sequence elements with antagonistic effects. For example, the cpeb1b (also known as zorba) mRNA contains adjacent stabilizing U-rich motifs and destabilizing miR-430 target sites. Combining these elements on the same mRNAs may provide differential temporal or spatial regulation, creating specific patterns of gene expression. We propose that these transcripts are first stabilized by maternally provided poly(U)-binding proteins and are later deadenylated and degraded by the zygotic mode of decay triggered by miR-430 or UAUUUAUU elements (Audic et al. 1997;Voeltz and Steitz 1998;Giraldez et al. 2006;Wu et al. 2006;Bazzini et al. 2012). Initial stabilization is likely coupled with active cytoplasmic polyadenylation, which would increase mRNA translation efficiency early in embryonic development (Subtelny et al. 2014). We find a strong correspondence between the motifs identified by RESA and those recently reported in zebrafish by Rabani et al. (2017). We further characterized the sequence contexts that influence the strength of the regulation associated with these elements. We identified an antagonistic activity between U-rich and CCUC sequences. Indeed, transcriptome-wide analysis of 3 ′ UTR mediated stability reveal a global tendency of 3 ′ UTR sequences (rich in U) to promote stabilization of the mRNA, with destabilizing islands that contain three main elements: miRNAs, CCUC domains, and UAUUUAUU motifs. Furthermore, the fact that deadenylation and decay are initially uncoupled during embryogenesis (Voeltz and Steitz 1998) allows us to define elements that preferentially affect mRNA stability (U-rich and CCUC) (Stoeckius et al. 2014), elements that mainly affect mRNA deadenylation (UAUUUAUU-ARE), and elements that induce both, such as miR-430. ELAV-like proteins (HuR) have been shown to stall decay following deadenylation mediated by AREs (Fan and Steitz 1998;Peng et al. 1998) and miRNAs (Kundu et al. 2012). The favorable sequence contexts for Elavl1b and Pcbp2 are similar to their binding motif and show a broader accumulation of iCLIP signal across the binding sites, suggesting that favorable context sites are bound by multiple proteins.
By using iCLIP, we characterized the binding motifs of multiple RBPs identified in the interactome capture. We observed that similar motifs are recognized by multiple RBPs, which are associated with similar regulatory activities. Although a high level of functional redundancy could be part of a robust developmental system, these RBPs might also recruit other proteins to achieve a specific regulatory response, as proteins with families containing the same binding domain are combined with different functional domains. Functional redundancy impairs our ability to genetically dissect their activity, an issue further complicated by the existence of multiple homologs for each RBP. Nevertheless, the genetic elimination of the regulatory elements identified in RESA clearly reveals the importance of these sites in the regulation of endogenous mRNA levels (Fig. 5H,I). These regulatory elements are enriched at the primary sequence level and associated with a high density of binding. Furthermore, we observed reduced RBP binding in the absence of zygotic transcription (Fig. 4E). This might be caused by lower levels of the target mRNAs in α-amanitin-treat-ed embryos. Alternatively, RBP binding might be activated by the MZT. For instance, Khsrp shows differential binding to RNA in our interactome capture and in wild-type versus homozygous MK2 or MK3 kinase mutant cells (Boucas et al. 2015). The activation of MK2 is known to play a fundamental role during MZT (Holloway et al. 2009), which led us to hypothesize that the reduced binding of Khsrp to RNA might be caused by lack of posttranscriptional activation. This hypothesis opens the possibility that an additional layer of regulation is controlling the clearance of maternal RNAs, emphasizing the potential central role of the regulatory elements recognized by the RBPs in regulating RNA stability.
We integrated these regulatory activities into a random forest model. Our model was able to predict regulatory activity across the RESA-targeted library of reporters for maternal and zygotic signals that mediate mRNA stability. We also identified sequence motifs associated with lower polyadenylation status in the mRNA (i.e., miR-430 and UAUUUAUU) dependent on the zygotic mode. The stability of endogenous mRNAs was less accurately predicted (r = 0.29 with asymptotic P = 1 × 10 −158 ) than the RESA 3 ′ UTR reporter library (per transcript Pearson's r = 0.82 with asymptotic P = 0.0) (Supplemental Fig. S10C,D). Our model is solely based on the primary sequence of the 3 ′ UTR, and its lower predictive power in endogenous transcripts is consistent with recent studies (Rabani et al. 2017). Generally, differences in the RESA and mRNA-seq assays could explain this difference in accuracy. Although using the same reporter backbone for all fragments tested in RESA allowed us to eliminate any bias owing to different coding sequences, we have yet to probe other factors regulating mRNA fate not assayed by RESA. Notably, the codon bias of the coding sequence of mRNAs influences their stability (Hanson and Coller 2018), specifically during the MZT (Bazzini et al. 2016;Mishima and Tomari 2016). Also, our RESA libraries are synthesized in vitro and therefore lack any RNA modifications. Modifications such as N6-methyladenosine (m 6 A) have been proposed to help the clearance of maternal mRNAs (Zhao et al. 2017), yet it is unclear whether m 6 A-containing mRNAs are specifically stabilized upon loss of m 6 A reader proteins (Kontur and Giraldez 2017). In addition, RNA structure analysis has identified long-range RNA interactions that are currently not being integrated in the RESA reporters owing to the limited size of the fragments tested (Aw et al. 2016;Lu et al. 2016;Sharma et al. 2016;Beaudoin et al. 2018). Finally, our model is restricted to linearly adding the effect of multiple motifs. This design was motivated by the additive effects of the U-rich and C-rich motifs we observed (Fig. 6G), but it excludes the cooperative effects of multiple sites (Hon and Zhang 2007;Saetrom et al. 2007;Lai et al. 2012). Including cooperative effects will require more sophisticated models capable of identifying the combined effects of all possible pairs of selected elements along with distance between each pair. This exponentially larger feature space will require a significantly larger training data set to avoid underfitting. Explaining the decay dynamics of all mRNAs will require the integration of the regulatory activities found in the 3 ′ UTRs with additional elements in the mRNA, namely, codons, RNA structure, and modifications to achieve a global prediction of mRNA dynamics.

Early zebrafish embryo transcriptome
For the large developmental gene expression time-course, RNA-seq libraries were prepared after collecting developing embryos between 0 and 8 hpf. Yeast total RNA was spiked to allow for appropriate normalization of fold-changes. After extraction, RNA was subjected to poly(A)-selected RNA-seq library preparation and ribosomal RNA-depleted total RNA-seq library preparation. Sequenced reads were then mapped onto zebrafish genome Zv9 using STAR (Dobin et al. 2013) and gene annotation from Ensembl r78 (Aken et al. 2017). Significantly over-and underexpressed genes were determined using DESeq2 (Love et al. 2014). Decay mode of genes (Maternal, Zygotic: miR-430-dependent and Zygotic: miR-430-independent) and their contribution to gene decay were determined based on significant DESeq2 tests, which are described in the Supplemental Methods.
In situ hybridization for the org, trip10a, and dnajc5ga genes were performed as by Thisse and Thisse (2008).

RESA
The transcriptome-based reporter library was generated by overlapextension PCR with primers mapping to the SP6 promoter and downstream from the SV40 polyadenylation site (Yartseva et al. 2017). Libraries were injected into one-cell zebrafish embryos. After library preparation (Yartseva et al. 2017) and sequencing, RESA profiles were obtained for each transcript (transcriptomic library) or UTR (targeted library). They represented positional read coverage normalized to count per million (CPM) using total counts of all transcript/UTR profiles per sample. RESA peaks in UACA and tdrd7a loci were validated using qRT-PCR. For trip10a, CRISPRmediated mutagenesis (Moreno-Mateos et al. 2015) was used to obtain a 223-nt deletion. Stabilization of the trip10a transcript was assessed by in situ hybridization.
To evaluate the regulatory effect of short RNA elements, RESA profiles were split into fixed-length sequence segments (30 and 100 nt long for the transcriptomic and targeted libraries, respectively) with 33% overlaps (Oikonomou et al. 2014). For each segment, a score was computed and differences between scores (early and late and between conditions) were transformed to Zscores. The scores were used to categorize segments to destabilizing and stabilizing modes of regulation: maternal, zygotic/miR-430 dependent, and zygotic/miR-430 independent. Within each category, segments were grouped based on their category-specific P. Then, FIRE, a computational framework for the discovery of regulatory elements (Elemento et al. 2007), was used for de novo discovery of short motifs that are significantly informative of the different modes of regulation. To this end, FIRE was extended to include N-fold cross-validation, a restriction on the average degeneracy of the elements, and zebrafish-specific options.

Interactome capture
Wild-type zebrafish embryos were irradiated at 4 hpf with UV at 254 nm for 4 min and then collected. Interactome capture was conducted as described by Castello et al. (2012) with oligo(dT) magnetic beads. The -UV control group of embryos were not UV-irradiated before collection. After peptides were eluted, LC-MS/MS was performed to identify them. Detailed protocol is included in the Supplemental Methods. Finally, enrichment in the molecular function category of identified proteins was tested.

RBP binding: CLIP
FLAG-tagged RBPs genes were first cloned and in vitro transcribed to capped mRNA. UV crosslinking and pull-down was conducted on 4 hpf embryos. The iCLIP protocol described by Huppertz et al. (2014) adapted to zebrafish was performed. Khsrp was iCLIPed using endogenous antibody. After read mapping using STAR (Dobin et al. 2013), RBP binding profiles were computed. To summarize RBP binding within protein-coding transcripts, all UTR and CDS iCLIP binding profiles were transformed into metagene profiles using 50-nt binning. To define RBP motifs, enrichment of 6-mer sequences at RBP target sites was calculated while taking into account the iCLIP experimental background. Logo representation of the top-10 motifs was built by aligning them using MAFFT (Katoh and Standley 2013). Finally, to simultaneously analyze binding and regulatory activity, sequences matching RBP motifs were searched within the transcriptome. Within the 100nt window, iCLIP (5 ′ end) and RESA (coverage ratio) windows were averaged.
Modeling the effect of sequence on mRNA stability The frequency of all k-mers (1-8 nt) within the 100-nt sliding windows (10-nt step) of the RESA profiles was used to train a random forest model with 500 trees, aiming to build a set of decision trees able to capture the association between k-mer frequencies and their correlation with RESA stability values. To make the learning process more efficient, a preprocessing, unbiased, and fast k-mer filtering step was implemented. After filtering rare, not correlated, and not important (using feature importance function) k-mers, 57 k-mers with the most significant contribution were selected. Trained models were validated using fivefold cross-validation for each library (RESA targeted and poly(A) selected). For each of these two libraries, a separate model was trained and validated for each condition (WT, LNA-430, and α-amanitin). Finally, the stability score per each transcript in the RESA-targeted library was calculated as the average of all the sliding windows predicted scores across the 3 ′ UTR of that transcript.
Stability of endogenous transcripts was predicted using the random forest model trained on RESA-targeted library by averaging the predictions obtained using the same sliding window approach.
Independent validation using targeted mutations on the trip10a locus was performed comparing (1) GFP expression compared with the control dsRed expression using fluorescence microscopy on 24 hpf zebrafish embryos and (2) stability measured by high-throughput sequencing of wild-type versus mutated trip10a reporter RNAs.

Data access
The sequencing data generated in this study have been submitted to the NCBI Sequence Read Archive (SRA) database (https://www .ncbi.nlm.nih.gov/sra). The developmental time-course mRNAseq data are available in SRP189512, in addition to the previously published SRP072296 (Bazzini et al. 2016) and SRP149556 (Beaudoin et al. 2018). RESA sequencing is available in SRP189389, in addition to the previously published SRP090954 (Yartseva et al. 2017). Demultiplexed and raw reads for iCLIP data are available in SRP189499. SRP149368 contains demultiplexed and raw reads for endogenous Khsrp iCLIP (Beaudoin et al. 2018). For iCLIP multiplexed data, a column barcode in the SRA annotation indicates the internal barcode used to multiplex replicates (see Methods). To facilitate data download, internal to laboratory (AGR) and SRA (SRR), run IDs are listed in Supplemental Table S1 and at https://data.giraldezlab.org. The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium via the PRIDE partner repository with the data set identifier PXD009514. Source codes of the FIRE software (Supplemental Code S1) and the RESA random forest model (Supplemental Code S2) are also available at https://data.giraldezlab.org. Updated gene counts, data sets, and genome tracks are available at https:// data.giraldezlab.org.