Analysis of environmental RNAs competent to become trans-acting is complicated by fragmented transcripts in samples. Immunopurification of Ago complexes followed by isolation and sequencing of associated RNAs can be used to validate small RNA identity [19]. Unfortunately, validated antibodies necessary for this approach are unavailable for all but a handful of organisms. To address this problem, we used Hi Trap QFF chromatography, which can enrich for Ago/Piwi-associated small RNAs (Fig. 1) [20]. Unbound RNAs in animal lysates are retained on the resin while highly basic Ago/Piwi proteins pass through. When used on spider mite extracts, a single clear peak of small RNAs is observed (Fig. 1a). In comparison, total RNA samples show many additional RNAs (Fig. 1b). After library creation and sequencing, size distribution of column purified RNAs shows a lower proportion of reads under 20 nt and an enrichment of piRNA sized (25-28 nt) species (Fig. 1c). Degradation fragments are typically associated with smaller read sizes, thus greater representation of longer sized reads in the column purified libraries demonstrates the effectiveness of isolating functional small RNAs. We also observed three times the representation of T.urticae miRbase miRNAs (Fig. 1d) [21]. Using this purification strategy, bona fide trans-acting small RNAs processed from exogenous transcripts can distinguished from fragments.
Following previously reported methods, we sought to induce maximum gene silencing using an animal soaking protocol with dsRNA targeted to spider mite Actin and GFP (Fig. 2 & Supplemental Fig. 1). Following exposure to dsRNA, animals were processed into a lysate and small RNAs were isolated with the Hi Trap approach described above. Following sequencing, small RNA reads were aligned to actin and GFP mRNA sequences, which showed significant accumulation of 18–21 nt RNAs at the target region (Fig. 2a & Supplemental Fig. 1). For untreated conditions after Hi Trap enrichment the accumulation was absent with very few alignments. This contrasts with many apparent degradation fragments for total RNA (Fig. 2a). We also found actin and GFP siRNAs were not heritable as none were present in embryos sired by soaked mites (Supplemental Fig. 1b,c).
We then sought to examine Dicer processing of long dsRNA-derived siRNAs by documenting overlaps between complementary pairs of small RNA reads [22]. Specifically, the abundance of small RNA pairs with 2 nt overhangs, characteristic of RNase III cleavage, were summed and plotted in a matrix of pair lengths (Fig. 2b,c). This reveals the sizes of small RNAs found in duplexes that have a signature of Dicer cleavage. In total ~ 45% of Actin and ~ 20% for GFP mapping reads showed signs of Dicer activity consistent with the centrality of this enzyme in the processing of dsRNA, and further validates the purification strategy for examining functional siRNAs. Common pairs were the expected 18–21 nt length and tended to be offset by 1–2 nts (Fig. 2b,c). Actin siRNAs showed more diversity in pair lengths compared to GFP siRNAs, suggesting capture of target cleavage products associated with Ago complexes by Hi Trap purification. We also observed many antisense reads mapping outside the target region that were not present in the total RNA library. This potentially could be the result of Rdrp activity producing antisense transcripts in response to siRNA targeting (Fig. 2d). We then inspected 5′ end sequence biases of actin reads in the dsRNA flanking regions (Fig. 2e,f). While there are no clear sequence motifs in sense mapping reads, antisense have clear 5′ terminal “T” and “A” at the 10 position. This is suggestive of piRNAs, which exhibit these features as a result of piRNA “ping-pong” processing [23]. In spider mites siRNAs act upstream of piRNAs and stimulate their production [6]. These results suggest this may extend to siRNAs derived from exogenous sources. This is consistent with the presence of somatic piRNAs in mites and suggests that piRNA-based processing might participate in eliminating siRNA-targeted transcripts.
Next, we investigated another source of exogenous RNAs–those acquired from plant hosts to assess the universe of RNA species that infiltrate RNAi pathways. Spider mites infest many plants but are readily reared on bean plants (Phaseolus vulgaris). In sRNA sequencing libraries created from total RNA extracted from bean-raised animals, ~ 3.13% map perfectly to bean sequences. Of these reads, 85% align exclusively while 15% also map to the mite genome. When Hi Trap purified libraries are aligned to plant sequences there is a 10-fold decrease in mapping rates, and a substantial shift in the proportion of reads that align to both mites and the plant (67%). Hi Trap mite purified RNA was also extracted from animals raised on Arabidopsis thaliana, which is not a preferred host of spider mites [24]. We observe reduced plant aligning reads to a rate of 0.09%, consistent with reduced intake of plant materials. 57% of reads from this sample also mapped to the mite genome. Plant mapping RNAs derived from mites have features of endogenous mite small RNAs, such as a peak of read sizes expected for siRNAs (Supplemental Fig. 2). Comparison to public small RNA sequencing libraries from the plants themselves (P. vulgaris and A. thaliana) showed shifts in the dominant read sizes (Supplemental Fig. 2). Thus, the RNAs in the mite library are unlikely to be mature plant sRNA contaminants.
Simultaneously, a substantial portion of the RNAs do appear to be generated from plant transcripts. Dicer overlap pairs from mites derived RNAs that map to plant sequences exhibited a pattern similar to long dsRNA and the mite genome (Fig. 2b,c, Fig. 3a, Supplemental Fig. 3a). The same pattern was not seen with total RNA library mapping to plant (Fig. 3a). Endogenous plant sRNAs were quite different, showing the greatest number of pairs at 24 nt, which is consistent with the 24 nt Dicer products found in plants (Fig. 3a, Supplemental Fig. 3c) [25]. Strikingly, the plant mapping reads from both Hi Trap and total RNA samples do not recapitulate this pattern when aligned to the mite genome (Fig. 3a). Thus, it would seem the plant mapping RNAs derived from mite samples are for a large part produced from plant transcripts despite co-mapping to the T. urticae genome. Further, plant mapping reads that do map back to the spider mite genome overlap with longer-sized reads in mites, again potentially tying piRNA production to siRNA targeting in this animal.
To understand the origin of plant derived siRNA found in mite libraries we annotated loci that transmit RNAs to mites based on read alignment density to plant genomes (Fig. 3b, Supplemental Fig. 3d). Using the total RNA dataset, we uncovered 1335 loci and 199 with the Hi Trap sample. Nearly all the Hi Trap loci overlapped with the ones annotated with total RNA, however, this was not tied to relative expression level. Examining the size distribution of reads mapping for each Hi Trap annotated loci showed bias similar to synthetic dsRNAs (Fig. 3 c, Supplemental Fig. 1a, Supplemental Fig. 2c,d). The pattern was not seen for these loci when using total RNA library alignments or a plant endogenous small RNA library (Fig. 3 c). Next, we intersected endogenous plant small RNAs loci with those called from mite derived RNAs libraries (Fig. 3 d). Generally, there was significant overlap with Hi Trap libraries but not total RNA (Fig. 3 d). However, the Hi Trap RNAs are not generated from well-processed plant small RNAs like miRNAs, which were only found in the total RNA dataset. This suggests that while plant RNAs that have features of small RNAs precursors are more likely to enter mite RNAi pathways, it is probably a double-stranded form of the RNA that is taken up. We also observe enrichment of Dicer overlap reads as a portion of mapping reads in Hi Trap samples with an apparent further enrichment for animals exposed to the synthetic dsRNA (Fig. 3 e). Thus, it would appear super abundance of dsRNAs in diet could saturate turnover enzymes and leads to greater stability of the small RNAs derived from dietary sources.
To explore the greater accumulation of Dicer products after dsRNA feeding, we focused specifically on alignments of mite derived libraries to plastids as these organelles do not have Dicer activity. Potential dsRNA/shRNA substrates are not processed prior to ingestion [26]. Plastids also lack well defined transcriptional units, leading to potential widespread formation of dsRNA [27]. Consistent with precursors of small RNAs being transmitted we observe over-representation of chloroplast sequences in all mite-derived libraries (Fig. 4a). Greater accumulation was seen in bean fed mites versus A. thaliana where less plant material is ingested. Interestingly, we also observe dsRNA feeding correlated with the fold increase in plastid genome coverage. This appears to be in part caused by a shift to dsRNA substrates. The amount of genome coverage for single-stranded mapping is similar between dsRNA fed and non-fed mites while the portion covered by dual strand mapping nearly doubles in the fed animals (Fig. 4b).
Generally, small RNAs map to regions of high gene density and rRNA loci (Fig. 4c) (Supplemental Fig. 4). In total RNA libraries the majority of reads map to a single strand of the rRNA loci. In contrast, the Hi Trap dataset show more consistent dual strand mapping, even at rRNA, albeit at a greatly reduced rate (Fig. 4d). In animals with dsRNA feeding, dual-strand coverage is even more equal with a greater portion of siRNA-sized (19-23 nt) transcripts. This suggests that among plastid RNAs ingested by mites a small, yet clear subset enter RNAi pathway. Prior exposure to dsRNA appears to stabilize the on-pathway RNAs which we propose is the result of saturating processing machinery. Alternately, by targeting Actin the dsRNA treatment could affect gut mobility leading to retention of RNA molecules and potentially infiltration of dsRNA pathways. However, together it would appear that siRNAs generated from exogenous RNAs are very unstable being preferentially eliminated rather than incorporated into regulatory complexes.
Based on our finding that structured plant-derived RNAs enter spider mite RNAi pathways we investigated the activity of short hairpin RNAs that mimic distinct, known miRNA-type biogenesis (Fig. 5a). The ability of four different configurations of synthetic hairpin RNAs to trigger silencing of Actin was tested alongside long dsRNA also targeted to Actin by soaking in in vitro synthesized molecules. qPCR was used to quantify target knockdown (Fig. 5b). The short hairpins were designed to transit: canonical miRNA biogenesis (shRNA), Ago processed short loop RNAs that mimic miR-451 with and without a G-C clamp at the hairpin base (SL1 & SL2), and a G-C clamp stabilized RNA with a loop sequence complementary to actin (Loop) [28, 29]. GC clamps were added to the structure bases of the RNAs to increase stability as this high energy fold is a challenging substrate for RNases [30]. After soaking, relative abundance of actin transcripts was assessed showing the reported ~ 40% reduction for long dsRNA [13]. Three of the structured RNAs showed a greater degree of knockdown with a reduction of transcript accumulation reaching nearly 60%, suggesting that the siRNA pathway of T. urticae may not be the optimal mode of RNAi to exploit for gene silencing (Fig. 5b).
We also sequenced samples after feeding with the SL2 (Fig. 5c) and shRNA (Fig. 5d) RNAs. SL2 derived small RNAs were more abundant but showed smaller fragments than shRNA, consistent with being an atypical RNAi substrate (Fig. 5e,f). For both RNAs, however, cleavage sites were non-uniform, and for SL2 there was little evidence slicer-mediated precursor processing was occurring as seen with miR-451, but rather standard dicer processing, suggesting that additional engineering could yield an even more superior gene silencing molecule. Nevertheless, many reads were complementary to Actin that were able to trigger knockdown. This together with the appearance of long dsRNA-derived siRNAs shows that a variety of ingested RNAs enter spider mite RNAi pathways, and that the designated siRNA/long dsRNA pathway may not be the most potent option for eliciting knockdown.