Re-Arrangements in the Cytoplasmic Distribution of Small RNAs Following the Maternal-to-Zygotic Transition in Drosophila Embryos

Small ribonucleic acids (RNAs) are known to regulate gene expression during early development. However, the dynamics of interaction between small RNAs and polysomes during this process is largely unknown. To investigate this phenomenon, 0–1 h and 7–8 h Drosophila melanogaster embryos were fractionated on sucrose density gradients into four fractions based on A254 reading (1) translationally inactive messenger ribonucleoprotein (mRNP), (2) 60S, (3) monosome, and (4) polysome. Comparative analysis of deep-sequencing reads from fractionated and un-fractionated 0–1 h and 7–8 h embryos revealed development-specific co-sedimentation pattern of small RNAs with the cellular translation machinery. Although most micro RNAs (miRNAs) did not have a specific preference for any state of the translational machinery, we detected fraction-specific enrichment of a few miRNAs such as dme-miR-1-3p, -184-3p, 5-5p and 263-5p. More interestingly, we observed changes in the subcellular location of a subset of miRNAs in fractionated embryos despite no measurable difference in their amount in unfractionated embryos. Transposon-derived endo small interfering RNAs (siRNAs) were over-expressed in 7–8 h embryos and associated mainly with the mRNP fraction. In contrast, transposon-derived PIWI-interacting RNAs (piRNA), which were more abundant in 0–1 h embryos, co-sedimented primarily with the polysome fractions. These results suggest that there appears to be a complex interplay among the small RNAs with respect to their polysome-cosedimentation pattern during early development in Drosophila melanogaster.


Introduction
In Drosophila melanogaster, early developmental processes until the maternal-to-zygotic transition (MZT) require the function of nearly 7000 maternally-loaded RNAs [1,2]. The MZT transfers the developmental control to the zygotic genome and induces the degradation of maternal messenger RNAs (mRNAs). Accumulating evidence suggests that small non-coding RNAs (ncRNAs) play an essential role in regulating the MZT both in insects and mammals [3][4][5][6][7].
There are at least three types of well-documented ncRNAs, miRNAs, endogenous silencing RNAs (endo-siRNAs) and piRNAs [8][9][10]. Micro RNAs are ubiquitously expressed small RNAs of~22 nucleotide in length that regulate gene expression post-transcriptionally by decreasing

Small RNA Deep-Sequencing Data Analysis
Embryo collection, polysome profiling and small RNA deep-sequencing data have been previously described ( [35], GEO accession number GSE35443). We used the small RNA-seq data to investigate the small RNA dynamics during early development in Drosophila melanogaster. To this extent, the 36-bp multiplex sequences were split into their corresponding samples according to their barcodes. Only sequences without any ambiguity were used for further analyses.
After barcode selection, the 3 adapter sequence was trimmed from the raw reads in 4 steps by using the 3 adapter sequence 5 -ATCTCGTATGCCGTCTTCTGCTTGT-3 : (1) Trimming full adapter sequences identified inserts with a minimum length of 7 nt, (2) If no adapter sequences were found, the last base of the 3 adapter was removed in a stepwise manner down to a minimum adapter size of 4 nt, which permit identifying inserts up to 8-28 nt, (3) Finally, the remaining reads were searched for adapter sequences with mismatches, and (4) Inserts with poly(A) and/or multiple unidentifiable nucleotides (Ns) were removed from the data. The resulting 15-29-bp inserts were collapsed into a fasta format and used for subsequent analyses.
For quality control analyses, nexalign program [36] was used to align all sequences to the genome (dmel-r6. 18, flybase.org) [37] and CCA-appended mature transfer RNAs (tRNAs) (flybase.org). The sequences were categorized in an order as those with exact match, one (insertion or deletion or mismatch), two or three mismatches. The remaining sequences were grouped as unmapped.
To calculate miRNA expression levels, we used the exact-matched mature miRNAs sequences and hairpin sequences in miRBase (Release 21). The fold-change for each miRNA was calculated using the formula (fold change = ((7-8h_RPM + 10)/(0-1h_RPM + 10)) in which 10 reads were added to each read to eliminate overestimations in lower reads. The resulting miRNAs were clustered using Cluster 3.0 [42] and visualized by java Treeview [43]. The transposon sequences from Flybase were used to identify the transposon-derived piRNA expression levels. We aligned all transposon-derived 23-29-nt reads to the Repbase collection to obtain the log 2 ratio of 0-1 h and 7-8 h RPM values. In order to validate that the 23-29-nt reads are indeed piRNAs, we searched for two features in these reads (1) a 10-nt 5 -5 complementarity and (2) a preference for U and A at the nucleotide positions 1 and 10, respectively. To analyse the piRNA clusters, the previously published piRNA cluster coordinates were used [41]. The 21-nt reads were not used for piRNA analysis, as they represent siRNA population. All scripts for small RNA sequence data analyses were written in Python and available upon request.
We downloaded 10 fastq files for 0-2 h and 11 fastq files for 6-8 h embryonic transcripts from modENCODE [44]. All fastq files were aligned to Drosophila melanogaster genome (dmel-r6.18) using GSNAP [45]. featureCounts [46,47] was used to assign reads to each gene using GTF file (dmel-all-r6. 18.gtf) downloaded from flybase.org. We used R package DESeq2 [48] to normalize reads, calculate fold changes and p-adjusted values. Differentially expressed genes were selected using 2-fold change and p adjusted value <0.1. The pipeline for RNA-seq of mRNA data analyses were written in R and available upon request.

Quantitative PCR Analysis
RNAs smaller than 200 nt were isolated from unfractionated 0-1 h and 7-8 h embryos using the miRVana miRNA isolation kit according to the manufacturer's instruction (Thermo Fisher Scientific, Waltham, MA, USA). Complementary DNA was prepared from 200 ng small RNAs using RT 2 miRNA cDNA kit (QIAGEN, Hilden, Germany). qPCR was performed in duplicates of three biological replicates (Lightcycler 480, Roche, Basel, Switzerland) using the 3 standard reverse primer (5 TAGGT 22 -3 ) and the 5 forward primer (mature miRNA sequence). Relative miRNA expression was calculated after normalizing the data with U6 ncRNA.

Results
Blastoderm cellularization, which takes place between 2.5 and 3 h after fertilization in Drosophila melanogaster, is the first developmental process that requires post-MZT zygotic transcription [2]. Small RNAs play an important role in the regulation of the MZT to set the stage for post-MZT embryonic events [3,4,10]. To gain insight into the expression levels of small RNAs during MZT, we analyzed the small RNA-seq data that was previously deposited to GEO ( [35], GSE35443). This data set includes small RNA-seq of 0-1 h and 7-8 h embryos. Half of these embryos were used to isolate total RNAs from unfractionated embryos ( Figure 1A, un-fractionated (UF)) as a reference. The other halves were then fractionated, on sucrose density gradients (SDGs), into 4 major sub-fractions based on A 254 absorbance: (1) mRNPs devoid of rRNAs, (2) 60S, containing 28S rRNA, (3) monosome, and (4) polysome. Here we report the comparative analysis of dynamic changes that we have observed in the intracytoplasmic localization of all small RNAs during early development in Drosophila melanogaster.
One major concern in polysome profiling is the potential for a greater RNA degradation during sample processing prior to downstream events. Thus, we first checked the RNA quality of total RNAs and small RNAs using Agilent RNA 6000 Nano kit and Agilent Small RNA kit, respectively (Agilent, Santa Clara, CA, USA) ( Figure S1). Total RNA capillary electrophoresis verified high quality of RNAs phenol-extracted from fractions [35]. Interestingly, each fraction appears to possess a unique small RNA profile especially based on small RNAs that are smaller than tRNAs. When we calculated the frequency of deep-sequencing reads based on the read size, a great majority (88.8-99.5%) of adaptor-containing reads contained inserts of 15-29 bp in size, further verifying a high quality of input RNA. The low percentage of 0-14-bp inserts in fractionated RNAs, compared to that in unfractionated total RNAs, was consistent with the high quality of fractionated RNAs (Table S1). Interestingly, 1 h total RNA contained the greatest percentage of 0-14 bp inserts, congruous with the destabilization of maternal mRNAs [2]. Of the 12,553,921 total reads, 71.95% of the fragments matched perfectly to the Drosophila melanogaster genome with a number of unique sequences in each sample ranging from 71,830 to 757,803. Based on the relative length distribution, we observed enrichment in 22-28-bp RNAs ( Figure S2).

Small RNA Populations Localize in Distinct Cytoplasmic Reservoirs
We used the nexalign program to align sense or antisense sequences to the known RNAs [36]. A pre-alignment to all known RNAs prompted us to align sequences in an order as described in the Methods section. We obtained a perfect match ratio of 71-82% whereas more than 90% of sequences mapped with exact or a single mismatch. Quite interestingly, each small RNA population appeared to co-sediment with polysomal fractions to a different extent ( Figure 1B). In the non-polysomal mRNP

Small RNA Populations Localize in Distinct Cytoplasmic Reservoirs
We used the nexalign program to align sense or antisense sequences to the known RNAs [36]. A pre-alignment to all known RNAs prompted us to align sequences in an order as described in the Methods section. We obtained a perfect match ratio of 71-82% whereas more than 90% of sequences mapped with exact or a single mismatch. Quite interestingly, each small RNA population appeared to co-sediment with polysomal fractions to a different extent ( Figure 1B). In the non-polysomal mRNP fraction, the majority of small RNAs contained tRNA-derived fragments (tRFs), with 31% and 37% in 0-1 h and 7-8 h, respectively. Based on northern blot results, we detected differentially expressed tRFs in accordance with the deep-sequencing data, which was published previously [35]. In the polysomal fraction, the majority of small RNAs was derived from rRNAs and transposons. The pattern of enrichment with respect to the developmental stage (1 h versus 8 h) was similar in fractionated and unfractionated samples, indicating the potential biological significance rather than random degradation introduced during fractionation.
A relatively high percentage of rRNA-derived fragments in 80S and polysomal fractions of both 0-1 h and 7-8 h embryos is conceivable considering the abundance of ribosomes in these fractions ( Figure 1B). We noticed a much higher rRNA-derived fragment in the total RNAs of 7-8 h embryos compared to 0-1 h embryos ( Figure 1B) despite quite similar RNA quality based on Bioanalyzer capillary electrophoresis ( Figure S1A). Interestingly, while the percentage of rRNA-derived fragments are similar in mRNP, 60S and 80S fractions, it is higher in polysomal fractions in 7-8 h embryos. Normally, degraded fragments would be expected to enrich in mRNP fractions, however these small RNAs co-sediment mainly with the polysomal fraction, suggesting a biological degradation or processing from rRNAs rather than random degradation. Additionally, the transcript-derived fragments are more abundant in 0-1 h embryos, which probably indicates either clearance of maternal mRNAs or biogenesis of novel zygotic small RNAs. We then focused on miRNAs and transposon-derived small RNAs as they constitute the bulk of small RNAs in our samples.

Micro RNAs Interact with Cellular Translational Machinery at All States
Although the temporal expression of miRNAs during early development is well documented [15], the extent to which each dysregulated miRNA is associated with polysomes is unknown. To this end, we firstly aligned the RNA-seq reads to the miRNA hairpin with perfect matches. A total of 263 mature miRNAs were detected with at least 1 read in at least one sample. Setting the minimal threshold at 50 RPM filtered out 169 miRNAs, yielding a total of 94 mature miRNAs for further analyses.
To identify fraction-specific miRNAs, we selected ten most abundant hairpin sequences from each sample (10 hairpins from each sample of total, mRNP, 60S, monosome and polysome, yielding a total of 50 hairpins) and identified 18 common hairpin sequences that constituted 82.5-93.6% of all hairpin sequences with perfect matches in each sample (Table 1). Thus, these 18 miRNAs represent the most abundant miRNAs in each sample. We observed dynamic changes in the expression of many miRNAs in unfractionated total RNAs in agreement with the published results (Table 2 [15]). Analysis of miRNA expression in fractionated samples yielded several interesting points (Table 1). Firstly, for the majority of miRNAs, the miRNA levels in fractionated samples were comparable to those in the unfractionated samples, probably suggesting non-selective distribution of these miRNAs throughout the cytoplasm. Secondly, although most miRNAs are distributed nearly evenly throughout the fractions, a few miRNAs were enriched in particular fractions. For instance, dme-miR-5-5p is enriched in the mRNP fraction both in 0-1 h and 7-8 h embryos whereas dme-miR-1-3p is enriched in the 80S fraction. Thirdly and more importantly, the extent to which a miRNA is associated with a particular fraction appears to be relatively similar for most miRNAs. However, for a few miRNAs, the degree of association is developmentally regulated. For instance, dme-miR-9c-5p makes up of 11.5 and 4.1% of miRNAs in 0-1 h and 7-8 h embryos, respectively.
It is well known that maternal miRNAs are degraded after the MZT is completed and zygotic miRNA transcription modulates zygotic gene regulatory networks [28]. It is unknown however whether the polysome status of zygotic miRNAs follow that of maternal miRNAs. There are potentially two possibilities (1) if the target sequence determines the extent of association to different components present in each fraction, we would expect differential polysome association pre-and post-MZT, assuming that the targets of miRNAs pre-and post-MZT stage are different, (2) if the association of miRNAs with different components is independent from the target sequences, then we would expect to see a similar polysome-association pattern. To differentiate between these two hypotheses, we first calculated the miRNA expression ratios (0-1 h/7-8 h) in the unfractionated and fractionated samples. We then checked how the ratio in the unfractionated samples is reflected upon that in the fractionated samples. We assumed that if the cytoplasmic fate of a miRNA does not change pre-and post-MZT, the ratios obtained from the unfractionated and fractionated samples should be similar between 0-1 h and 7-8 h embryos. In contrast, if there are any changes in the localization of miRNAs, the difference in the expression level obtained from the unfractionated samples should manifest in a particular sub-cellular fraction. This approach resulted in identification of 4 different miRNA groups with each having a unique expression behaviour.
The first group includes 41 miRNAs that are over-expressed in 7-8 h embryos ( Table 2, G1; Figure 2). While 17 miRNAs co-sediment with the complexes in the mRNP fraction, some co-sediment with the 60S fraction. Interestingly, for 9 miRNAs, we did not see a proportional increase in their expression in fractionated RNAs although their expression was increased in un-fractionated embryos, suggesting the involvement of nuclear retention or other unknown mechanisms. The relative number of miRNAs classified in the second group does not change in 0-1 h and 7-8 h un-fractionated embryos, indicating a similar transcriptional activity and/or miRNA stability. However, we detected dynamic changes in their subcellular locations following the MZT (Table 2, G2; Figure 2). For instance, despite no difference in the total dme-miR-9a-5p amount in un-fractionated embryos, this specific miRNA becomes more polysome-associated in 7-8 h embryos (Log 2 fold = 3.3). The third group includes 29 miRNAs, whose expression decreases in 7-8 h embryos (Table 2, G3; Figure 2). Interestingly, the decrease in the miRNA expression is not distributed evenly throughout the 4 fractions, suggesting a preference for a specific fraction. For instance, 20 out of these 29 miRNAs sediment specifically with non-polysomal fractions, particularly in the mRNP fraction. This indicates that the majority of small RNAs acts at mRNP complexes in early embryo compared to the 7-8 h ones. For some miRNAs, we detected equal distribution in the cytoplasm despite a decrease in the unfractionated embryos, suggesting perhaps more efficient nucleo-cytoplasmic transport. The fourth group includes the miRNAs whose relative abundance in the un-fractionated and fractionated embryos are different. These miRNAs behave similarly to those in the second group in that the transcriptional output from these miRNAs does not lead to proportional changes in their cytoplasmic location, suggesting differential localization pathways pre-and post-MZT (Table 2, G4; Figure 2). We then selected six miRNAs, representing at least one miRNA per group, to validate RNA-seq data by quantitative PCR (qPCR) analyses of total RNAs from unfractionated embryos (Table 3). Results from both methods were in parallel to each other.

Micro RNA Cluster Behaviours and miRNA Editing
Some miRNAs are known to be transcribed as members of gene clusters [50]. Expression studies revealed that miRNA clusters are co-expressed [51] and cluster members are coordinated during target regulation [52]. Assuming such coordination during target RNA regulation, we hypothesized that the cytoplasmic fate of the members should be similar. Thus, we compared the extent to which each cluster member is associated with four different fractions in our experimental design, each representing a different translational state of the cell. Our cluster analysis showed that miRNA cluster members behave similarly with respect to their cytoplasmic localization (Table S2).
We also checked the frequency of post-transcriptional miRNA editing events likely to occur during the early development in Drosophila melanogaster as miRNA editing is commonly used in eukaryotes to modulate the targets of miRNAs [53]. We first aligned our sequences to the known miRNA sequences and looked for the sequences that align to the known sequences with a single mismatch. Based on this approach, we identified one candidate editing event (dme-miR-986, C→T at 11th position) ( Figure S3). The PCR-amplification and sequencing of dme-miR-986 from P2 strain embryos and S2 cells showed that this particular difference in the sequence stems from a single nucleotide polymorphism (SNP) not an editing process.

Cytoplasmic Micro RNA Re-Arrangement Does Not Correlate with Target Abundance
Our observation in the cytoplasmic distribution of miRNAs is interesting in the sense that localized miRISC complexes may differentially regulate their target mRNAs. To test this hypothesis, we examined a potential correlation between the four miRNA groups that we identified (Table 2) and the target mRNA abundance in 0-2 h and 6-8 h embryos [44]. Thus, we first checked the quality of 0-2 h and 6-8 h embryo RNA-seq data ( Figure S4A,B [44]) and then identified the differentially expressed transcripts ( Figure S4C). Of these mRNAs differentially expressed in 0-2 h and 6-8 h embryos, we selected the experimentally validated targets of miRNAs included in our four miRNA groups. When we categorized the target transcripts based on their abundance (i.e., up-regulated, down-regulated, no-change), we did not detect a correlation between the cytoplasmic distribution of miRNAs and their effect on target mRNA abundance ( Figure S5), suggesting perhaps the presence of some other mechanisms.

Transposon-Derived Small Interfering RNAs and PIWI-Interacting RNAs Interact with Different Complexes
Some piRNAs are maternally deposited into the oocyte and they might be involved in the degradation of maternal mRNAs from the embryo [7]. Also, MILI/MIWI is associated with polysomal fractions [54]. However intra-cytoplasmic distribution and re-arrangements, if any, of piRNAs during early development are unknown. We used the Repbase [40] collection to calculate the small RNAs generated from transposons. The expression of transposon-derived small RNAs decreased towards 7-8 h, suggesting the significance of these small RNAs probably before the MZT. To ensure that we selected the transposon-derived piRNAs, we looked for two main features associated with piRNAs. The first feature is the 10-nucleotide complementarity at their 5 ends as generated by the ping pong cycle ( Figure 4A). The second feature in Drosophila melanogaster is the presence of a U nucleotide at the 1st position and an A nucleotide at the 10th position ( Figure 4B; reference [55]). Furthermore, we only used the reads in length of 23-29 bp to select the sequences that match the aforementioned two criteria.
We divided transposon-derived small RNA transcripts into two groups based on their sizes: transposon-derived piRNAs of 23-29 nt and transposon-derived siRNAs of 21 nt [55]. To globally compare the polysome association of transposon-derived piRNAs to that of transposon-derived siRNAs in embryos, we calculated the read frequency of transposon-derived small RNAs. We noticed a much higher peak at the 21 bp reads in 7-8 h embryo total RNAs ( Figure 4D) compared with 0-1 h ( Figure 4C). Moreover, the 21-bp read frequency is higher in mRNP-associated complexes decreasing towards polysomal fractions. This data suggest that transposon-derived siRNA expression is more abundant in later developmental stages although piRNA expression appears to be slightly higher in 0-1 h embryos. Additionally, the intracellular localization of transposon-derived siRNAs appears to be different from that of transposon-derived piRNAs ( Figure 4C,D). Although transposon-derived piRNAs are distributed throughout four fractions nearly equally, transposon-derived siRNAs are predominantly associated with the mRNP fraction. We then checked the frequency of reads from the 42AB cluster as it is one of the best characterized piRNA clusters in Drosophila melanogaster [56]. This analysis showed, in parallel to previous findings, production of transcripts from both strands: Interestingly, there appears to be a strand-bias in polysome-associated transcripts especially in 0-1 h embryos ( Figure 3B).  The reads were first aligned to the Repbase collection [40], all reads and uniquely mapped reads were used to calculate the transcript expression levels. Then the log2 fold change was calculated and clustered as in Figure 2A; (B) The reads that are mapped to one of the most abundant transposon cluster in Drosophila melanogaster, 42AB. The reads were first aligned to the Repbase collection [40], all reads and uniquely mapped reads were used to calculate the transcript expression levels. Then the log 2 fold change was calculated and clustered as in Figure 2A; (B) The reads that are mapped to one of the most abundant transposon cluster in Drosophila melanogaster, 42AB.
To substantiate our observation that transposon-derived piRNAs are highly expressed earlier compared to the siRNAs, we collected previously published deep-sequencing data from different developmental stages of Drosophila melanogaster [15,41,57]. We used the same strategy as in Figure 3 to trace the temporal expression of siRNAs and piRNAs. Our analysis showed that the siRNA expression is relatively low in early developmental stages but increases significantly during later developmental stages ( Figure 5A). Using the same data sets, we also checked the relative abundance of Drosophila melanogaster small RNAs in various developmental stages. This analysis validated the notion that the transposon-derived small RNA abundance drops in later developmental stages while miRNA expression level increases ( Figure 5B). Moreover, the 21-nt siRNAs derived from transposon increase while the piRNA levels (23-29 nt) decrease.
predominantly associated with the mRNP fraction. We then checked the frequency of reads from the 42AB cluster as it is one of the best characterized piRNA clusters in Drosophila melanogaster [56]. This analysis showed, in parallel to previous findings, production of transcripts from both strands: Interestingly, there appears to be a strand-bias in polysome-associated transcripts especially in 0-1 h embryos ( Figure 4B).
To substantiate our observation that transposon-derived piRNAs are highly expressed earlier compared to the siRNAs, we collected previously published deep-sequencing data from different developmental stages of Drosophila melanogaster [15,41,57]. We used the same strategy as in Figure 4 to trace the temporal expression of siRNAs and piRNAs. Our analysis showed that the siRNA expression is relatively low in early developmental stages but increases significantly during later developmental stages ( Figure 5A). Using the same data sets, we also checked the relative abundance of Drosophila melanogaster small RNAs in various developmental stages. This analysis validated the notion that the transposon-derived small RNA abundance drops in later developmental stages while miRNA expression level increases ( Figure 5B). Moreover, the 21-nt siRNAs derived from transposon increase while the piRNA levels (23-29 nt) decrease.    The reads were first aligned to the Repbase collection [40], all reads and uniquely mapped reads were used to calculate the transcript expression levels. Then the log2 fold change was calculated and clustered as in Figure 2A;

Discussion
We present here an in-depth analysis of cytoplasmic distribution of small RNAs in early development with respect to their association with the cellular translational machinery. The

Discussion
We present here an in-depth analysis of cytoplasmic distribution of small RNAs in early development with respect to their association with the cellular translational machinery. The comparative analysis of the most abundant small RNAs reveals that the spatial location of each type of small RNA is quite different (Figure 1 and Table S3). Interestingly, miRNAs appear to utilize numerous molecular mechanisms as they interact with the translational machinery at all states. piRNAs, which are expressed more abundantly at the pre-MZT stage, are associated with polysomal complexes. In contrast to piRNAs, transposon-derived siRNAs are more abundantly expressed at the post-MZT stage and are primarily found in the mRNP fraction.
During the maternal-to-zygotic transition two crucial events take place: (1) clearance of maternal mRNAs and (2) transcriptional activation of the zygotic genome. The clearance of maternal mRNAs is orchestrated by numerous RNA-Binding proteins, including Smaug [1,2] and small RNAs such as miRNAs, siRNAs and piRNAs [4,7,[59][60][61]. Differential expression and specific functions of miRNAs in early development were previously reported [15,22] and our findings are in agreement with these results. For instance, dme-miR-8 and 10 are up-regulated in 7-8 h embryos in total RNAs while dme-miR-310 and 311 are down-regulated.
Although the differential expression of miRNAs during the MZT stage is well documented, the polysome profile of miRNAs during this switch is unknown. Our findings suggest that the majority of miRNAs appear to be associated with distinct fractions and probably with distinct complexes as a result (e.g., mRNP or polysomal complexes) ( Table 2). For instance, dme-miR-263-5p, 5-5p and 9c-5p are primarily found in the mRNP fraction whereas dme-miR-184-3p is mainly part of polysomal fractions. Bantam-3p and dme-miR-1-3p are primarily enriched in 60S and 80S fractions, respectively ( Table 2). Although we cannot attribute any mechanistic insight into the fraction-specific functions of miRNAs yet, it might be noteworthy to speculate that this specific localization might be associated with association of miRNAs with different components in each fraction. For instance, since mRNP complexes are known to include translationally inactive mRNAs [62], the miRNAs in this fraction might be involved in sequestration of mRNAs (or storage) away from the translational machinery. The miRNAs in the 60S and 80S fractions are likely to interfere with the translation initiation while the polysome-associated miRNAs probably interfere with translation elongation.
One other interesting observation from our findings is that the change in the number of miRNAs in total RNAs (transcriptional or post-transcriptional increases/decreases) is not distributed equally throughout each fraction in the embryos. Some of differentially expressed miRNAs are directed towards the translationally inactive mRNP complexes whereas some others are destined for the translationally active polysomal complexes (Table 2). Although the polysome association of miRNAs has been documented before [63], the switch in their state suggests an interesting layer of regulation. We interpret all these observations to mean that miRNAs have a pre-defined cytoplasmic fate following their nucleo-cytoplasmic transport. When we looked at the validated target abundance of miRNAs categorized under four different groups, we did not see any correlation between the cytoplasmic localization pattern of miRNAs and their target abundance. We might not have detected any correlation for at least three reasons; (1) other regulatory mechanisms, e.g., transcriptional regulation, could be in effect, (2) target RNA stability may not change upon miRNA binding (i.e., translational regulation), (3) effect on target mRNAs could be residual, making it difficult to detect. Thus, further studies with reporter target sequences would be required to precisely trace the actual effect of cytoplasmic miRNA distribution on its target.
For a few miRNAs ( Table 2, dme-miR-1-3p, dme-miR-2b-3p, dme-miR-92a-3p, dme-miR-92b-39 and dme-miR-1012-3p) miRISC complexes appear to switch from one translational state to another. For example, mRNP-fraction-associated miR-1-3p switches to polysomal fractions following the MZT. We cannot, however, conclusively claim that miRNAs switch their intracellular location as some of these post-MZT miRNAs could be zygotically transcribed miRNAs. In this situation, there should be additional factor(s), e.g., miRNA-binding proteins that specify the polysome status of zygotically transcribed miRNAs. In either case, our data suggests that the activity of miRNP complexes may be spatially modulated through the differential cytoplasmic localization of these complexes pre-and post-MZT. More direct evidence is required, however, to demonstrate whether the eukaryotic cells utilize intracellular re-localization of miRNPs as a means to modulate miRNA function and target mRNA(s) as a result.
We also detected cytoplasmically localized transposon-derived small RNAs, which are classified as siRNAs (21 nt) or piRNAs (23-29 nt). PIWI-interacting RNAs protect the zygotic genome against infection by retroviruses that can come from the surrounding follicle cells or endogenous transposons within the female germline [4]. PIWI-interacting RNAs function together with the well-characterized Smaug to direct the CCR4 deadenylase to specific mRNAs, thus facilitating maternal mRNA deadenylation and decay [7]. The cytoplasmic localization of piRNAs in early development is consistent with their role in deadenylation of maternal mRNAs. The exact role of endogenous siRNAs in early development is not well-defined. Based on the observation that the mRNA profiles of wild type and Dgcr8 null mouse oocytes were identical, it was proposed that endo-siRNAs, rather than miRNAs, are responsible for the Dicer knockout phenotype observed in mice [64]. Thus, endo-siRNAs and piRNAs are primarily involved in gametogenesis and very early development whereas miRNAs appear to get involved in later developmental stages [65,66]. Our data is consistent with this view in that transposon-derived transcripts are more abundant in 0-1 h embryos (both fractionated and unfractionated, Figure 1B) whereas miRNA expression is up-regulated in 7-8 h embryos. Endogenous siRNAs are well-known for their role in heterochromatin formation in the nucleus [67]. The identification of endo-siRNAs in the cytoplasm points to two possibilities (1) the process of maturation as they are produced from double-stranded RNAs (dsRNAs) in cytoplasm, (2) a potential role for these small RNAs in post-transcriptional gene regulation.
Transposon-derived siRNAs and piRNAs differ from each other in two ways: (1) 21-nt transposonderived siRNAs are highly expressed in 7-8 h embryos whereas piRNAs abound in 0-1 h embryos (Figures 2 and 4); and (2) siRNAs are mainly associated with the complexes in the mRNP fraction while piRNAs are largely associated with the polysomal complexes. It is quite interesting that both siRNAs and piRNAs are produced from the same transposons, yet they interact with the cellular translational machinery quite differently. Another interesting point is that we detected a strand-bias in the reads obtained from the dual-strand transposon 42AB in 0-1 h embryos especially in the polysome fraction ( Figure 5B). It remains to be investigated whether there is a functional relationship between the biased-production of piRNAs and development.
We previously reported that tRNAs also serve as templates for small RNAs, i.e.,~28-nt tRFs [35]. These small RNAs are expressed at both stages ( Figure 1B) and in fact throughout the embryonic development and in mature flies as well. They are mainly associated with non-polysomal fractions, resembling transposon-derived siRNAs. Studies on tRF function mainly focused on miRNA-like regulatory functions. In Drosophila melanogaster, tRFs are immunoprecipitated with anti-AGO1 antibody [68]. Interestingly, they inhibit cap-dependent translation initiation [69]. It is unknown however whether tRFs are specifically involved in modulation of translation at the MZT stage.
Supplementary Materials: The following are available online at www.mdpi.com2073-4425/9/2/82/s1. Figure S1: Total RNA quality control analysis, Figure S2: The insert length distribution of small RNAs, Figure S3: SNP on dme-miR-986-5p, Figure S4: Analysis of 0-2 and 6-8 hour old embryonic RNA-Seq data, Figure S5: Correlation between miRNA and their targets [67], Table S1: The percentage of hairpin and mature miRNA mapped reads in all reads mapped to miRNAs, Table S2: miRNA cluster expression analysis, Table S3: Alignment of the reads to known RNAs.