Ago proteins are the key component of the siRNA and miRNA pathway and are indispensable binding proteins for the function of many other small RNAs. For example, tRF must bind Ago2 to direct the cleavage of target RNAs, which contain a perfect complementary sequence to the tRF . Using BmAgo2 polyclonal antibody, endogenous BmAgo2 was identified in BmN cells (Additional file 1: Figure S1). The Ago2 protein usually displays a Mg2+-dependent endonuclease activity, and Mg2+ provides a favorable contribution to the binding of the DNA-RNA to the Ago2 protein . Our work demonstrates that recombinant BmAgo2 could exhibit an endonuclease activity and cleave the IP RNA products when EDTA is removed from the lysis buffer (Additional file 1: Figure S6), suggesting recombinant HIS-tagged BmAgo2 might exhibit the same physiological activity as endogenous BmAgo2.
Using a HIS monoclonal antibody, we successfully isolated BmAgo2-associated small RNAs by immunoprecipitation from the BmN cells. We used an RIP-seq method to obtain a total of 817,302 BmAgo2-associated unique small RNAs from the ovary-derived BmN cells, 91,092 of which we could not match to the reference genome of the male silkworm. The unmatched reads might be derived from the W chromosome or might map to potential gaps in the sequenced female genome. To date, many species of small RNAs have been identified in Bombyx mori, including miRNAs [24–27, 29, 30, 91], piRNAs [35, 37, 38], TE-derived small RNAs [18, 34] and snoRNA/snRNAs of between 50nt and 500nt . In the present work, we identified BmAgo2-associated small RNAs between 18nt and 50nt in length, including the first identified tRFs, rRNA-, snoRNA- and snRNA-derived small RNAs in this species, as well as previously identified miRNAs, piRNAs and TE-derived small RNAs. Importantly, the BmAgo2 protein was found to recruit tRFs at an extremely high abundance (Figure 2B), suggesting that tRFs could function by binding the BmAgo2 protein.
We further employed Northern blotting to identify tRFs, and the results showed two obvious bands responding to the tRNAs and tRFs, respectively, suggesting that tRFs are a stable and novel class of small RNAs (Figure 3). Previous studies have shown that the ANG nuclease can generate 30-50nt long sitRNAs under conditions of stress from the anti-codon loops, and it has been reported to be a tRNA-specific enzyme [92–94]. The Bombyx 5′ and 3′tRFs identified here by Northern blotting were approximately ~33nt and ~40nt long, respectively. Northern blotting with overlapping probes further confirmed that the ~40nt 3′tRFs were generated by cleaving the anti-codon loops of the 72nt tRNAs (Additional file 1: Figure S7). These 5′ and 3′tRFs might be sitRNAs that are produced by the ANG nuclease. Our deep sequencing results further suggested that the ~33nt 5′tRFs and the ~40nt 3′tRFs could be further cleaved on their D and TψC loops by a certain nuclease, thereby generating the ~18nt 5′ truncated tRFs and ~21nt 3′trailer tRFs (Additional file 1: Figure S3), respectively. However, we failed to detect these BmAgo2-associated tRFs by Northern blotting from among the total RNAs of the normal or BmNPV-infected BmN cells, possible due to the low abundance of these tRF species. ANG and other RNases may be involved in generating these small tRF species, which can cleave tRNAs to produce ~20nt terminal tRFs . Furthermore, a previous study also showed that the major determinant in the processing of terminal small tRFs by ANG and possibly by related RNases is the TψC loop within the tRNA. The ANG tRF cleavage site is known to be located between the first and second base after the GTTC motif on the TψC loop . Our deep sequencing data show that a conserved GTTC domain exists at the termini of most 21nt 3′trailer tRFs, suggesting there is a potential tRF cleavage site between the GTTC cytosine and the first base after the GTTC motif on TψC loop (Additional file 1: Figure S3). The 3′trailer tRFs could not be detected by Northern blotting as they had a low total RNA abundance; however, the 3′ trailer tRFs were significantly more abundant than the 3′tRFs in the BmAgo2-associated total RNA. In the sequenced BmAgo2-associated small RNAs, the read copies of the 3′trailer tRFs were far greater than those of the 3′tRFs (Additional file 2: Table S2), suggesting that the 3′trailer tRFs preferentially associated with BmAgo2 for their functions, which might be similar to those of the siRNA pathway. For the tRNAs with a variable arm, the D loop could be cleaved at the GGCCGAGCGG cleavage site, generating the ~18nt 5′truncated tRFs (Additional file 1: Figure S3D). This result implies a possible relationship between the D loop cleavage and the variable arm during the biogenesis of 5′truncated tRFs.
The tRF biogenesis that occurs under conditions of viral stress provides an important piece of evidence about the functional role of tRFs. The present work shows that when the BmN cells were infected with BmNPV, many tRFs were generated and recruited onto the BmAgo2 protein. Therefore, these tRFs were enriched among the BmAgo2-associated small RNAs. The majority of these 5′tRFs identified by Northern blotting were expressed only in BmNPV-infected cells (Figure 3A), whereas the other tRF species identified by Northern blotting were usually expressed in both the normal and the virus-infected cells (Figure 3B), suggesting the main tRF species generated by viral stress may be the 5′tRFs. The asymmetric expression of 5′tRFs implies a relationship between the 5′tRFs and BmNPV infection in the BmN cells. Furthermore, 5′tRFs accounted for an overwhelming proportion of the BmAgo2-associated small RNAs. This was especially true for Bm-AspGTC-5′a, which accounted for 59.58% of the total BmAgo2-associated small RNAs. These observations suggest that the 5′tRFs are likely to bind the BmAgo2 protein and thus influence the BmNPV infection.
The piRNA system is thought to provide a germline defense against TE activity. A similar function was also recently noted in the siRNA pathway [51, 68, 95–97]. Drosophila utilizes two small-RNA systems to restrict transposon activity in its germline (mostly via piRNAs) and soma (mostly via siRNAs); in this system, Ago2 deficiency results in increased levels of transposon transcripts . In mouse oocytes, the loss of Ago2 can result in decreased levels of siRNAs and increased levels of retrotransposons . In S. japonicum, sjAgo2 strongly associated with TE-derived siRNAs and is believed to regulate transposons . These observations indicate that the Ago2 protein might be involved in regulating transposons at the transcriptional level in germline cells. The BmN cell line was obtained from the Bombyx ovary. TE-derived small RNAs accounted for almost half of the total unique BmAgo2-associated small RNAs. The peak of their length distribution was the typical siRNA length. Therefore, we suggested that BmAgo2 could execute important functions in the maintenance of genome stability in germline cells through an RNAi-like pathway.
piRNApredictor is a software program for piRNA prediction that uses a novel 'dynamic’ algorithm with a precision of over 90% and a sensitivity of over 60% . Using piRNApredictor supplemented with a perl script for screening Ping-pong pairs, 54,033 >25nt piRNA candidates and 11,337 typical Ping-pong pairs were identified from among the TE-derived small RNAs associated with BmAgo2. In addition, we further identified 56,131 >25nt piRNA candidates and 11,898 Ping-pong pairs with 17,533 sense piRNAs and 18,239 antisense piRNAs from the remaining unannotated reads (Additional file 2: Table S10). These piRNA candidates could not be mapped to the identified TEs and could be derived from some heterochromatic regions, other repeat sequences and unknown transposons in the Bombyx genome. piRNA is thought to defend the host genome against transposons in Bombyx mori. Recently, an increasing number of piRNAs have been identified in Bombyx mori[18, 35, 37], and their biogenesis has been further studied, providing evidence for their "Ping-pong Circle" and biosynthesis mechanisms [37–40]. We further compared our data with the previously identified piRNAs [18, 35, 37]. The results showed that our data has 0.83% (909), 14.95% (16,474), 30.62% (33,733), 36.54% (40,254), 17.55% (19,335) and 17.72% (19,517) overlap (mismatch not allowed) with the TE-, BmN4-, Siwi-, BmAgo3-, ovary- and testis-derived piRNAs, respectively. This overlap serves as an indicator of the high reliability of our dataset. We first identified piRNAs associated with the BmAgo2 protein, which also contains PAZ and PIWI domains. BmAgo2 might bind to piRNA-like RNA molecules through the PAZ domain, the biogenesis of which could be similar with that of piRNAs. The observations above suggest that BmAgo2 can associate with the 20 ~ 22nt TE-derived siRNAs and the >25nt piRNA-like small RNAs during the regulation of transposon expression.
We must emphasize that many more reads did not align to previously known ncRNAs, such as miRNAs in miRBase, tRNAs, rRNAs, piRNAs, snoRNAs, snRNAs and other reported Bombyx mori small RNAs , indicating that many small RNAs remain to be discovered in Bombyx mori. The unannotated small RNAs associated with BmAgo2 usually exhibit a higher abundance than the miRNAs. Only 49 miRNAs had a read copy of at least 100 (Additional file 2: Table S6), but 848 unannotated small RNAs had a read copy of at least 100.