RNA interference (RNAi) technology has been shown to be applicable as a low-toxicity biopesticide to control agricultural insect pests and vectors of plant pathogens through silencing essential, biologically-relevant genes (Zotti and Smagghe 2015). RNAi shows great potential to be highly species-specific and thereby spares beneficial organisms, and is non-toxic to humans and other animal consumers. The RNAi approach for insect pest/vector control relies on the ingestion of long double-stranded RNAs (dsRNAs) to trigger gene silencing via siRNA production following Dicer processing (Head, Carroll et al. 2017, Knorr, Fishilevich et al. 2018). While a number of products are available, some arthropod pests exhibit moderate or only minor sensitivity to dsRNA upon ingestion (Yu, Christiaens et al. 2013, Zhu andPalli 2020). This suggests that to fully realize this strategy across most or all arthropods, RNAi triggers may require unique engineering relevant to each target species (Shukla, Kalsi et al. 2016, Parsons, Mondal et al. 2018).
Here, RNAi pathways were investigated for the whitefly Bemisia tabaci (Genn.) (Aleyrodidae, Hemiptera) to characterize the fundamental features that might be exploited to improve RNAi approach(es). As a group, this whitefly is considered a cryptic or sibling species. Although most B. tabaci are relatively benign, at least two variants/cryptic species transmit plant viruses and are among the most invasive species causing damage to crops grown in subtropical, tropical, and mild temperate parts of the world (J K Brown, D R Frohlich et al. 1995, Brown 2010, Chen, Hasegawa et al. 2016, de Moya, Brown et al. 2019, Grover, Jindal et al. 2019). Chemical pesticides can be toxic to the environment and consumers of these products, and regularly have become ineffective when resistance develops (Chen, Hasegawa et al. 2016). B. tabaci is closely related to greenhouse and spiraling whiteflies, and several other related phloem-feeding pests/pathogen vectors, including aphids, mealybugs, and psyllids. Previous RNAi studies with B. tabaci have demonstrated gene silencing in response to long dsRNA feeding, however, processing modes of these molecules and those in other non-holometabolous insects have not been characterized at the level of small RNA effector populations (Jaubert-Possamai, Le Trionnaire et al. 2007, Zha, Peng et al. 2011, El-Shesheny, Hajeri et al. 2013, Thakur, Upadhyay et al. 2014, Li, Zhang et al. 2016, Wang, Wang et al. 2016, Vyas, Raza et al. 2017, Grover, Jindal et al. 2019, Kanakala, Kontsedalov et al. 2019). In this study the behavior of these molecules is investigated in the context of an extant RNAi mechanism in B. tabaci.
Multiple biogenesis modes are reported for animal small RNAs, generally though three main classes are recognized: microRNAs (miRNAs), smallinterfering RNAs (siRNAs), and Piwi-associated RNAs (piRNAs) (Carthew and Sontheimer 2009). miRNAs are deeply conserved, short-hairpin derived RNAs that are present in the cells of nearly all metazoans (Bartel 2018). In contrast, the biology of the other two small RNA classes is highly variable, likely due their role in defense against invasive nucleic acids like transposable elements and viruses (Okamura 2012). Indeed, their overall roles in animals has been observed to diverge even among family members (Ozata, Gainetdinov et al. 2019). Each small RNA variety is sorted to distinct Argonaute (Ago)/Piwi proteins with Ago's binding siRNAs and miRNAs, and Piwi's binding piRNAs. siRNAs are 20-23nt products of Dicer cleavage, usually from a long > 100nt dsRNA molecule. While important to antiviral response many endogenous siRNA (endo-siRNAs) species derived from hairpin RNAs (hpRNAs) or cis-NATs can be found in arthropod genomes (Fagegaltier, Bougé et al. 2009, Lau, Robine et al. 2009, Claycomb 2014).
In comparison, piRNAs, are Dicer independent and produced by several methods that include endonucleolytic "slicer" activity present in Piwi proteins (Yamaguchi, Oe et al. 2020). They and are typically 26-30nt long and have a "U" residue at their 5' end. In Drosophila two biogenesis modes of piRNAs have been observed that are mediated by three piwi proteins: PIWI, Aubergine (Aub), and Argonaute 3 (Ago3) (Ozata, Gainetdinov et al. 2019). The "ping-pong" mode involves alternating target RNA "slicing" by Piwi proteins. This is an amplifying mechanism where processed RNAs are recruited as new piRNAs. Aub cleaves RNAs that load into Ago3 as secondary piRNAs. Reciprocally, Ago3 substrates load into Aub. The other mechanism relies on piwi proteins cleaving designated transcripts that then become substrates for the RNase Zucchini (Zuc), which the subsequently load into PIWI/Aub in Drosophila (Gainetdinov, Colpan et al. 2018).
In Drosophila, piRNAs are associated with germline, however, in many other arthropods piRNAs are found in soma, including hemipterans (Lewis, Quarles et al. 2018).
Relative to Drosophila, whitefly has additional Ago/Piwi proteins and exhibits somatic expression of Piwis. Using small RNA sequencing datasets miRNAs, endo-siRNAs, and piRNAs were annotated. Many of the loci resembled those described in Drosophila, however, we found many instances where there was production of both piRNAs and siRNAs. Using characteristics of these endogenous loci we designed RNAi triggers that exploit piRNA pathways. Gene silencing triggered by the piRNA pathway were equally as efficient as the siRNA mediated silencing of endogenous genes. Somatic piRNAs are widespread among insects, however, their application as a pest control tool is yet to be developed. Other hemipteran insect pests and vectors of pathogens to plants and humans, such as pea aphid and the kissing bug, respectively, have been found to produce somatic piRNAs (Brito, Julio et al. 2018, Lewis, Quarles et al. 2018). Insights of this study are expected to apply directly to these pests as well as many others.
Due to the divergent nature of siRNA and piRNA biology, species specific design is necessary to fully exploit these pathways for effective gene silencing.
To characterize the RNAi pathways of whitefly we first sought to identify the collection of Ago/Piwi proteins encoded in the whitefly 'B biotype' (also known as MEAM1) genome (MEAM1v1.2) using existing annotations and BLAST to curate sequences (Fig 1A) (Chen, Hasegawa et al. 2016). These sequences were then compared to Ago/Piwi proteins from D. melanogaster, Tribolium castaneum (red flour beetle), and subset from C. elegans. We found seven members of this family encoded in the whitefly genome. Three of the genes belong to the Ago family with one clearly related to miRNA-loading factors. The two other Agos group with the siRNA-associated Agos and appear to be a clade specific duplication as they are not homologs of the two si-Agos seen in T. castaneum. This is consistent with the diverging nature of siRNA biology and opens up the possibility that si-Ago function in whitefly might be distinct from fruit flies and beetles. The other four members belong to the piwi clade. One of the whitefly Piwi proteins is a homolog of DmeAgo3 while the other three group with DmePIWI/DmeAub. This indicates that ping-pong biogenesis is likely present in these animals. The phasing piRNA pathway is also presumably operative with apparent homologs of Zuc (Bta02312) and the RNA helicase Armitage (Bta07189) (Ishizu, Kinoshita et al. 2019).
piRNAs are found in somatic tissues of many insect orders, including hemipterans (Huang, Fejes Toth et al. 2017, Lewis, Quarles et al. 2018). To verify if this is also pertinent for whitefly, we investigated RNAi factor expression in the whitefly guts, salivary glands, and whole body (Fig 1B) (Cicero and Brown 2011). PolyA sequencing libraries from extirpated whitefly guts, salivary glands, and whole body were mapped to the RNAi factor sequences from MEAM1v1.2, per above (Chen, Hasegawa et al. 2016). The Alignments were then used to calculate RPKM values for each transcript. Expression of BtaAgo1 and BtaAgo2 was found in all tissues, as well as two Piwi's (BtaPiwi3 and BtaAgo3), along with BtaZuc,and BtaArmi (Fig 1B). This contrasts with BtaPiwi2, which is enriched in whole body presumably due to the inclusion of RNAs originating from gonad tissues. This suggests that similar to other hemipterans, whiteflies have somatic piRNAs with both ping-pong and phasing piRNA biogenesis modes being present. Significantly, somatic piRNAs are likely present in whitefly gut, the tissue that would be the primary target of ingested RNAi trigger molecules.
To further investigate whitefly RNAi pathways endogenous small RNA populations from whole body mixed adults (male and female) were examined using small RNA sequencing libraries mapped against MEAM1v1.2. From this alignment we first annotated miRNAs using miRDeep2 (Friedlander, Mackowiak et al. 2012). Subtracting miRNA-derived reads from datasets would allow focus on non-miRNA small RNA loci such as endo-siRNAs and piRNAs, which unlike miRNAs might have whitefly specific biology. 202 miRNAs are identified with high confidence with 89 being conserved in Drosophila (Fig 1C and Table S1). We also identified 124 additional miRNAs which were classified as lower confidence miRNA candidates due to suboptimal features such as low expression or imprecise precursor cleavage patterns. The miRNA repertoire of the whitefly genome is similar in size to other insects (Kozomara, Birgaoanu et al. 2019).
These results expand on prior miRNA annotations in whitefly due to the increased depth of datasets featured in this study.
Next we examined the size distribution of reads and found a bi-modal read size distribution with peaks at 22nt representing Dicer products (siRNAs and miRNAs) and 29-30nt (piRNAs) (Fig 1D). Among the Dicer products, roughly 56% derive from miRNAs. This shows piRNAs are more abundant relative to siRNAs in whitefly. To examine the modes of piRNA production we analyzed the abundance of read overlap pairs and the distance to 1U trailing reads (Fig 1D).
During ping-pong biogenesis piRNA pairs are cleaved at the 10th base of the guiding RNA. Thus, when this mode is active piRNAs are found to overlap by 10 bases, which is clear in the dataset. Phasing piRNAs are biased to occur end to end and can be recognized by close proximity of trailing reads. Phasing is also evident in the alignment. The abundance of piRNAs is further reflected by the high proportion 1U reads in the size distribution. Simultaneously, a significant proportion of the reads also exhibit a "A" at the 10th base which would be found on ping pong pair reads due to pairing with 1U.
Using reads subtracted of miRNAs we annotated non-miRNA, small RNAproducing loci. 3873 regions were identified with a read depth greater than 40 and 500+ bp length (Fig 2A andTable S2). The ratio of the number of small (19-23nt) to long reads (25-30nt) was then calculated to distinguish whether the locus produced smaller siRNAs (19-23nt) or longer piRNAs (25-30nt). This showed the majority of loci appear to be piRNA generating. Only 50 loci had a ratio of small to long that was greater than "one". Interestingly, the piRNA loci spanned regions ranging in size from 500nt to 50kb, indicating diverse transcripts generate this small RNA class (Yamanaka, Siomi et al. 2014). Apparent siRNA producing loci tended to be shorter regions, of which the longest was about 4kb in size. We then examined the distribution of read sizes at each locus (Fig 2B). Accumulation was most clear in the piRNA range, with substantially less for siRNAs. However, a minor signature of siRNA-sized reads could be seen at many loci. To confirm that most loci were sites of piRNA production we examined read overlaps and trailing 1U read distance, which shows evidence of ping-pong due to 10nt overlap bias and phasing with juxtaposed trailing 1U reads, respectively. The exception was ~100 loci that showed a greater accumulation of siRNAs.
To verify if these loci are sources of Dicer produced siRNAs we sought the 2nt overhang signature of RNase III processing (Fig 2C). Overlapping read pairs between 15-31nt with this signature were quantified. Pairs were identified where one strand (query read) of a potential duplex overlapped by less than two of its entire length, which would occur with a 2nt overhang (Antoniewski 2014). All potential combination of query and complementary target reads were quantified, revealing that 22nt reads show the greatest evidence of Dicer processing, and that this is likely size of Dicer2 products. The abundance of apparent Dicer overlapping reads differed from the distribution of the reads in different size ranges, validating this method of characterizing biogenesis. Interestingly, some signal could also be seen in the 29-30nt sizes that likewise was not reflected in the all read size distribution. This suggests a potential interaction between siRNAs and piRNAs unlike what is reported in Drosophila, and is consistent with the frequent co-occurrence of siRNA sized and piRNA sized reads across all annotated loci (Fig 2B ).
Next, we focused on the filter loci by expression to focus on the top 50 long read biased or short read biased loci (Fig S1, Table S1). The size distribution of reads for each locus was determined, which showed 28-30nt reads at long read loci, consistent with production of piRNAs (Fig S1). This contrasts with the short read loci, which show signal at 22nts and 29nts and 30nts, consistent with co-occurence of piRNAs and siRNAs seen across all loci (Fig 2B). We also characterized the two groups of loci by length, expression, and 1U bias (Fig S1). This showed that long read loci are larger, have a greater bias towards 1U and greater expression, which are characteristics of piRNA clusters.
Examining strand mapping showed that high expressing long read loci exhibit bias towards small RNA production from one strand indicating likely singlestranded precursor transcripts converted into phasing piRNAs (Fig S1) (Gainetdinov, Colpan et al. 2018). The short read loci are predominantly dualstranded, which is suggestive of a double-stranded RNA precursor serving as a substrate for Dicer (Claycomb 2014). We also identified 22 hairpin RNA loci indicating that this variety of locus is present as a minority of the overall collection of whitefly siRNA-generating loci (Fig S2).
To predict the function of these 100 loci, reads aligning to these loci were mapped back to the whitefly genome permitting up to three mismatches. This alignment was then intersected to MEAM1v1.2 annotations (Fig 2D ,E). The number of intersections was determined for each locus keeping mRNAs and TEs separate. Both long and short reads target mRNAs and TEs indicating possible roles for piRNAs and siRNAs not only in genome surveillance but also in gene regulatory networks. This is consistent with a proposed role for piRNAs in regulation of protein coding gene expression (Shamimuzzaman, Hasegawa et al. 2019). Taken together, this suggests that whitefly siRNAs and piRNAs are gene regulatory factors alongside miRNAs. These observations further reinforce the potential for exploiting these pathways for genetic technology that silences genic transcripts.
Whitefly endo-siRNA loci are also sources of piRNAs.
Prior work in whitefly has shown effectiveness of long dsRNA in gene silencing (de Paula, de Faria et al. 2015, Malik, Raza et al. 2016, Luo, Chen et al. 2017, Vyas, Raza et al. 2017, Grover, Jindal et al. 2019). The presumption is that these molecules are processed by Dicer into siRNAs. To better understand small RNAs simulated by fed dsRNA we used the computational approach described S1). These groupings were determined by evaluating dominant small RNA size and the dominant processing signatures of read pairs-2nt overhangs for Dicer or 10nt overlaps for ping-pong piRNA. The siRNA group are located in intergenic regions and have a strong bias towards short reads that appear to be dicer processed.
The cisNATs were sites of siRNA production between opposing mRNAs as showing in Fig 3B . In addition to these predominantly siRNA producing regions that are similar to ones observed in Drosophila (Czech, Malone et al. 2008, Ghildiyal, Seitz et al. 2008), many loci produce both siRNAs and piRNAs. These dual-identity loci could be grouped into one of three categories. One where there was equal production of siRNAs and piRNAs (No bias), a second for which some siRNAs were present, but piRNAs are dominant (piRNA), and the third group that harbors large piRNA clusters with only a minor production of siRNAs, the latter being similar to the locus shown in Fig 3B . These observations show that despite an Ago repertoire similar to Drosophila small RNA biogenesis in whitefly is distinct. This provides an opportunity to exploit these divergent activities for gene silencing and pest management.
We then extended our evaluation of processing dsRNA transcripts to those introduced exogenously via feeding. Here we tested three off-target, synthetic dsRNAs dissolved in a sucrose solution fed through an artificial system.
The RNAs cloned from genes of the potato psyllid Bactericera cockerelli (Sulc.) were fed to adult whiteflies from which small RNA and messenger RNA sequencing libraries were generated. Significant accumulation of reads arose exclusively from dsRNAs and not from other sections of the psyllid gene from which they were cloned (Fig 4A). However, only a fraction of the reads show a signature of Dicer processing based on 2nt overhangs. This suggests that most of the synthetic RNA was likely degraded with only a minority entering the siRNA pathway, which is reflected in the distribution of reads produced from the dsRNA sequences (Fig 4B). Only a modest peak was seen at 22nt with many more at the 15nt size. The low efficiency is likely in part caused by dsRNA specific nucleases (dsRNases), which are common in hemipteran insects. Several dsRNAses from gut and other tissues have been identified in whiteflies (Luo, Chen et al. 2017, Singh, Singh et al. 2017). Further, the alkaline pH of sternorrhynchan midgut is also compromised RNA integrity (Cristofoletti, Ribeiro et al. 2003, Molki, Ha et al. 2019).
Using these datasets, we sought to identify similarities between small RNAs derived from fed dsRNA and endogenously expressed siRNAs.
Specifically reads were subsetted based on sequence content to find population where signatures of dicer processing were most evident (Fig 4C). This was guided by the seqlogo results of endogenous siRNAs that showed preference for 1U and depletion of 1G (Fig 3A). In unfiltered reads only a slight enrichment of 22nt RNAs was seen with no evidence of 2nt overhangs. Next, reads were extracted based on their 5' residue, which showed similar size distribution to the unfiltered library with the exception of 1G reads where there was no bias towards 22nt reads. However, for each subset no 2' overhang specific to 22nt reads was observed. The analysis was then extended to include not only the first base of the read but also the 20th base. When considering 1U/A/C(H) and 20A/U/G(D) greater abundance of 22nt reads was seen but still no substantial 2nt overhang Dicer signature. When 1U/A(W) and 20A/U(W) were examined an even greater enrichment of 22nt reads as well as 2nt overhang enrichment for this size was observed. For individual nucleotide pairs (1U-20A, 1A-20U, 1C-20G, 1G-20C) read size and overhang enrichment increased, particularly for 1U-20A in size distribution and 1A-20U for 2nt overhang. This analysis provides a framework for computationally isolating siRNA processing signatures from degradation products, which is essential when considering exogenous dsRNA processing due to cloning of digestive contaminants.
To understand the physiological consequences of ingesting dsRNA we examined the effect on expression of the small RNA loci annotated in this study and protein coding genes. Studies in other animals suggest there may be competition between exogenous and endogenous small RNAs for biogenesis pathways, and that dsRNA may be recognized as a viral motif resulting in activation of defense pathways (Jelinek, Leonard et al. 2011). After feeding dsRNA, small RNA sequencing showed no significant change in endogenous small RNA expression compared to control (Fig 4D) (Table S3). For protein coding genes we observed about 500 transcripts that were differentially expressed based on a p-value ≤ 0.001 (Fig 4E) (Table S3). Only 20 of these genes exhibited a log(fold2) value greater than 2 or less than -2. All genes in this group have very low expression with 14 having unknown function. The genes having a known identity appear to be involved in basic metabolism or development. The one exception is a RNase H containing gene (Bta15726) that could be involved in an antiviral response. However, it is down regulated, which is inconsistent with being deployed to combat perceived viral infection. Thus, it would seem that whitefly does not mount a antiviral-type response to dsRNA. We also observed no change in expression of RNAi factors such as the Ago and Dicer proteins, suggesting that whitefly can metabolize exogenous RNAi triggers without affecting its core RNAi processes. Taken together it appears that when ingested, the bulk of dsRNA is degraded with a small amount contributing to the siRNA pool, and that exposure to dsRNA has minimal impact on off-target gene expression in whitefly.
In this study we found a significant population of piRNAs, which are more abundant than the endogenous siRNAs-the species exploited by existing RNAi approaches. The piRNAs also appear to be expressed in soma and show potential widespread control of mRNAs and not just a role in genome surveillance. This suggests that the piRNA pathway might be exploited to silence endogenous gene expression in whiteflies as an alternative method to the classic dsRNA-based siRNA strategy.
To trigger ectopic production, we engineered recombinant nucleic acids that take advantage of the major principle of piRNA biogenesis-recruitment of Piwi cleaved fragments into the pathway (Fig S3). We fused sequences from two loci annotated in this study, a piRNA bias locus (piRB-6) and siRNA-piRNA no bias locus (No bias-14) to target gene sequences. Both loci were among those that showed evidence of Dicer processing as well as piRNA production (Fig 3) S1). We chose two different genes to target with these constructs: Aquaporin1 (AQP1) and Alpha glucosidase1 (AGLU1), which were used in a previous study that yielded high gene knockdown via dsRNA (Vyas, Raza et al. 2017). To explore design principles, the positive strand of the locus was fused to AQP1 and the negative strand to AGLU1.
Using these constructs, both synthetic dsRNAs and ssRNAs (singlestranded RNA) were generated and fed to whiteflies in the artificial system described above. The concentration of RNAs (30 ng/ul) used was similar to what was previously fed to whitefly in dsRNA experiments (Vyas, Raza et al. 2017).
Luciferase sequences fused to piRB-6/No_bias-14 were used as off-target controls. Following feeding access for 6 days, expression of target genes was assessed by qPCR (Fig 5A). As previously reported, dsRNAs elicited gene knockdown of 68-80%. Satisfyingly, the piRNA triggers showed a similar degree of gene silencing with reduction in target expression of 60-80%. This result was observed for both ssRNA and dsRNA triggers with both piRNA sequences (piRB-6, No_bias-14) and targets (AQP1, AGLU1). For AQP1, piRNA triggers were equal to dsRNA (conventional dsRNA), while AGLU1 was not as well downregulated by the piRNA triggers relative to dsRNA suggesting inclusion of positive strand sequence might lead to superior knockdown. However, by combining ssRNA and dsRNA piRNA triggers for either sequence, gene silencing became comparable to conventional dsRNA for AGLU1. These results provide robust evidence that piRNA triggers, even those comprised of ssRNA, are capable of gene silencing in organisms that share RNAi biology with whiteflies.
Small RNAs were then sequenced to characterize the processing of the piRNA triggers. Small RNAs were sequenced from animals fed piRB-6 dsRNAs and ssRNAs targeted to both AQP1 and AGLU1 (Fig Fig 6). Reads mapping to these triggers showed significant heterogeneity in read size with no accumulation of a specific size, indicating the bulk of fed RNAs were degraded.
To identify potential small RNAs among the detritus we determined the relative abundance of 1U reads, a characteristic of piRNAs as well as siRNAs (Fig 5B).
From this we found significantly more 26-30nt piRNA-sized reads. In the doublestranded treatment a small peak possibly corresponding to 22nt siRNAs could be observed, but not for the single-stranded piRNA triggers.
Next we focused on the identity of small RNAs produced against the target gene. piRNA biogenesis could be observed for both triggers but more so for the single-stranded versions (Fig 5C). Ping-pong processing was observable when comparing the number of overlaps for different nucleotide combinations that can form pairs: 1U/10A, 1A/10U, 1C/10G, and 1G/10C. Read pairs were determined for all ranges of reads, and for those in piRNA sizes (28-30nt). The greatest enrichment for 28-30nt reads was seen for those with the signature of ping-pong piRNAs: 1U/10A. Phasing was also assessed for each strand of the piRNA triggers (Fig 5 D&E). This biogenesis mechanism was evident for the transcribed strand of single-stranded triggers, which is complementary to the target gene (AQP1, AGLU1). For both strands of the double-stranded trigger and the potential target derived reads in the single-stranded fed condition, less phasing was evident though a noticeable trend toward close proximity of 1U reads was seen. Next, we investigated whether siRNAs were processed from the triggers by examining 2nt overhangs in read populations as in This same phenomenon is seen in AQP1 triggers which sport the sense strand of piRB-6 for the on-target strand. For single-stranded AQP1 nearly all the RNAs appear to be phasing piRNAs, and for the double-stranded version most of the phasing piRNAs are on target. Both of these trigger versions lead to robust gene silencing. These results indicate that a superior choice for piRNA trigger design is to select the phased strand of piRNA loci to fuse with gene targeting sequences. It also shows that the small population of endogenous antisense ping-pong piRNAs or possibly even the siRNAs have a heightened role in promoting phasing. This is an intriguing departure from Drosophila where trailing piRNAs are produced downstream of a site of piwi protein-initiated cleavage.
Here it seems phasing of piRB-6 can be initiated internally because the region cloned for these triggers only includes an interior section of the locus (Fig 6D). It is also clear from these results that double-stranded triggers, expectedly, lead to greater production of siRNAs.
This study provides an in-depth analysis of the RNAi pathways in B. tabaci, a hemipteran insect pest and plant virus vector, and offers a rationale design of piRNA-based gene silencing biotechnology. Most significantly, we show ingested RNAs can enter piRNA pathways, which opens up the possibility for an entirely new strategy for gene silencing and potentially commercial products. On a superficial level whitefly small RNAs seem similar to Drosophila.
There are three distinct types of small RNAs (miRNAs, siRNAs, and piRNAs), as occur in fruit flies. However, upon close inspection, the biogenesis and function of the endogenous small RNAs in whitefly are quite different. Our work reinforces the consistent observation that non-miRNA RNAi pathways are fluid; clade specific duplication of the RNAi factors is common, even loss of an entire class of small RNA has occurred in several metazoan clades (Sarkies, Selkirk et al. 2015, Calcino, Fernandez-Valverde et al. 2018, Mondal, Klimov et al. 2018). Further, these findings illustrate the benefits of in-depth dissection of the RNAi biology for an evolutionarily and biologically different organisms, beyond those examined in model study systems, for developing genetic technology.
Through this comprehensive annotation of whitefly small RNA loci over 200 novel microRNAs are described, as well as 3878 siRNA or piRNA loci.
Previously described configurations whitefly siRNA and piRNA loci were observed such as large single-stranded, phased piRNA loci and siRNA expressing cis-NAT and hpRNA loci (Fig 3) (Table S1,Fig S2). However, curation of loci found extensive evidence of siRNA and piRNA biogenesis occurring simultaneously at many loci. In fact, this appeared to be the rule for the majority of endogenous siRNA and piRNA genes, and seemingly, is related to a different biogenesis and function for whitefly siRNAs or piRNAs. In Drosophila and vertebrates, piRNAs mainly control TEs in germline, however, many of the piRNA pathway accessory proteins such as Rhino, Deadlock, Cutoff, Moonshiner from the Drosophilids are not conserved indicating that piRNAs are shaped to individual organism's biology in an evolutionary arm race between the piRNAs and their targets (Ozata, Gainetdinov et al. 2019). Indeed, it is predicted that abundant somatic piRNAs engage gene regulatory networks in many basal arthropods, such as hemipterans, suggesting that this is the ancestral piRNA biology (Lewis, Quarles et al. 2018). This combined with observations that whitefly piRNAs respond to viral infection suggest diverse roles for these small RNAs in this insect (Shamimuzzaman, Hasegawa et al. 2019). Moreover, we find that phasing biogenesis can be initiated in the interior of loci which is distinct from the trigger/responder/trailing piRNA arrangement seen for phased Drosophila piRNAs.
While RNAi has been successful for controlling some pests like coleopterans (beetles), many other pests such as some lepidopterans (moths and butterflies) are unresponsive to exogenous RNAi trigger (Shukla, Kalsi et al. 2016, Parsons, Mondal et al. 2018). Penetrance of RNAi in hemipteran insects is moderate and higher dosage of dsRNA is required (Joga, Zotti et al. 2016). pH in the gut of the hemipteran insects is basic and presence of the nucleases in gut has been reported in whiteflies, aphids, and other hemipteran insects (Luo, Chen et al. 2017, Singh, Singh et al. 2017). We have noticed in this study that only a minority of the reads produced from the dsRNA trigger are siRNAs (Fig 4A&B).
This could be attributed to low abundance of the intact dsRNA for uptake by the gut epithelium cells. We see a similar accumulation of degradation products with piRNA triggers. Interestingly, even with single-stranded RNA triggers we see significant accumulation of small RNA reads from synthetic RNA along with some antisense reads. The antisense reads we observe have dominant ping-pong and minor siRNA signature. How the single-stranded RNAs trigger production of these molecules is not clear but could involve the recruitment of target mRNA into small RNA biogenesis. While this is standard behavior for piRNAs, it is not typical for siRNAs in organisms that don't possess Rdrp activity (Sarkies, Selkirk et al. 2015, Almeida, Andrade-Navarro et al. 2019, Pinzon, Bertrand et al. 2019).
We view this result as a first glimpse at a heretofore unappreciated small RNA biogenesis mechanism that involves interaction between siRNA and piRNA biogenesis, consistent with the widespread co-occurrence at endogenous loci.
Encountering such an unknown interaction is not entirely surprising as this study represents the first effort to characterize small RNA biogenesis on a per locus level in a non-holometabolous insect. This further reinforces the value of knowledge-based RNAi design gleaned from investigating exogenous trigger processing. In this study we provide clear rules for maximizing piRNA production, which could be fundamental to potent gene silencing technology aimed at aphids, mealybugs, psyllids, whiteflies, and other hemipterans.
As hemipteran insects respond to exogenous long dsRNA-mediated RNAi trigger only moderately, utilizing the gene silencing function of the piRNA pathway is exciting. These results show that in whitefly while there is significant sensitivity to dsRNA, there is very little physiological response to dsRNA feeding.
Even the secreted gut dsRNases do not become transcriptionally activated by feeding. This will likely apply to other hemipteran herbivores with similar composition of RNAi pathways and dsRNases. We expect that piRNA triggers, single-stranded or double-stranded, will likewise be physiologically neutral. The most promising result we report is that exogenous piRNA triggers are as effective as the siRNA versions. This study provides the first report of the exploitation of piRNAs as a feeding-based insect pest control strategy. Thus, this approach could become key for designing effective RNAi approaches against many insect pests that are found to be resistant to dsRNA-mediated RNAi. Finally, dsRNAs are capable of activating interferon response in humans and other vertebrates through binding of TLR3 receptors (Zhang, Xiang et al. 2016). Deploying ssRNA piRNA triggers as a pest control approach would avoid activating this pathway.
As a result, beneficial, non-pest organisms in the field would also be spared from off-target effects of dsRNAs as piRNA triggers rely on the specific genomic sequence of the target species and would not be converted into siRNAs as would happened with dsRNA-based triggers. Taken together, these findings demonstrate the benefit of in-depth studies of non-model organismal RNAi biology, and demonstrate somatic piRNAs can be used for environmental RNAi.
Insects in this study came from the type B. tabaci Arizona B biotype (AZ-B) whitefly colony established in Brown lab in 1988 following its discovery on poinsettia plants in Tucson, Arizona (Vyas, Raza et al. 2017). For this study, AZ-B adult whiteflies were serially transferred to and reared on cotton Gossypium hirsutum L. cv Deltapine 5415) plants at the 8-10 leaf stage. 199 nt long region from No_bias-14 locus (Scaffold40734: 1537-1774, 1811-2009) were fused to left and right sites of the GOI respectively. From piRB-6 locus the left and right flanking sequences were 342 and 366 nt respectively PCR products with T7 promoter sites on both strands were used for dsRNA synthesis while for ssRNA, PCR was carried by allowing the T7 promoter site in one strand. PCR products were directly used to synthesize the synthetic RNAs using MEGAscript ™ T7 Transcription Kit (Thermo Fisher Scientific, Catalog number AM1334) following manufacturer's protocol.
Using a hand-held aspirator 100 adult whiteflies were collected for each biological replicate from the colony and transferred to a plastic feeding chamber.
200 μl, 30 ng/μl RNA in 20% sucrose solution was sandwiched between two sterile Parafilm M® layers and feeding access to the solution was given to the insects for six days. On day 6, insects were collected for RNA extraction.
Total RNA was extracted following the standard Trizol RNA extraction method. The extracted RNAs were DNase I-treated (Invitrogen, DNA-Free™ kit, Lot 00522653) and 2 μg RNA was used for cDNA synthesis using High Capacity cDNA Reverse Transcription Kit (Applied Biosystems, Lot 00692533). The TaqMan qPCR master mix (Applied Biosystems, Universal PCR Master Mix, Lot #1908161) was used for quantitative gene expression analysis using standard protocol. Whitefly 18S rRNA gene was used for normalizing the expression of the target genes. All qPCR primer sequences from a previous study were used in this study and can be found in S1 Text (Vyas, Raza et al. 2017). Each treatment and control groups of the synthetic RNA feeding was carried out using at least three independent biological replicates. ΔΔCt method was used for gene knockdown analysis. Student t-test and one-way ANOVA were used for statistical analysis in CFX Maestro software v1.1 mRNA library preparation, sequencing, and gene expression Total RNAs were extracted using conventional Trizol RNA extraction method from different manually dissected tissues of whiteflies (gut, salivary gland, and whole body) (Cicero and Brown 2011). RNA integrity was confirmed using an Agilent 2100 Bioanalyzer (Agilent Technologies Inc., Santa Clara, CA).
Sequencing libraries were constructed using Illumina's TruSeq RNA Sample Preparation Kit v2, Catalog # RS-122-2002 (Set B). Using magnetic oligo (dT) beads, only poly(A) tail containing RNAs were separated from total RNA. Next, the mRNAs were fragmented by zinc treatment, and first-strand cDNA was synthesized from the fragmented RNAs using SuperScript II reverse transcriptase and random primers from Invitrogen. Then second-strand cDNA was synthesized, and Illumina multiple indexing adapters were ligated to the fragments. The remaining library construction steps were carried out following manufacturer's protocol. Quality filter and processing of the sequenced reads were performed using Illumina CASAVA v1.7.0, FastQC, and Trimmomatic. For each of the RNAi factors analyzed (Fig 1B ), reads were mapped with bowtie2 to transcripts sequences from the whitefly genome database www.whiteflygenomics.org (Langdon 2015, Chen, Hasegawa et al. 2016).
Bedtools was used to count read alignments to each transcript (Quinlan and Hall 2010).
Total RNA was extracted from adult whiteflies using standard Trizol RNA extraction protocol. Following DNase treatment, small RNA-seq libraries were constructed using NEXTflex Small RNA-Seq Kit v3 (NOVA-5132-06).
Small RNA reads were quality checked using FastQC and the adapter sequences were cleaved and trimmed using fastx toolkit. Next, 15-35 nt size reads were mapped to whitefly genome (MEAM1 genome v1.2) using bowtie with default parameters (Chen, Hasegawa et al. 2016). The genome mapped reads were isolated for the downstream analysis. mirDeep2 was used to annotate the miRNAs (Friedlander, Chen et al. 2008). Initial calls by the algorithm were manually inspected for recognized features of miRNAs (Berezikov, Liu et al. 2010). Annotations that showed evidence of mature and star strands in the appropriate Dicer cleavage register as well as significant expression were placed in the confident category. Deviation from these characteristics resulted in placement of annotation in the candidate category.
For non-miRNA annotations, small RNA reads, either taking all reads, 19-23nt sized reads, and 25-30nt reads were aligned using bowtie multimapping (-a -m 100) options. Bowtie was also used to identify the targets by allowing three mismatches. Size distributions were calculated with basic unix commands: awk, sort, uniq etc. Using bowtie alignments ping-pong overlap, piRNA phasing, and Dicer siRNA overhangs signatures were calculated as previously reported (Antoniewski 2014, Han, Wang et al. 2015). Samtools and Bedtools were used to count read alignments and identify high expressing regions and bias towards short and long read loci, as well as determine potential targets (Quinlan and Hall 2010). The R packages Scatterplot3d, sushi, heatmap2, pheatmap, and ggplot2 were used to draw the read density graphs (Kolde 2012, Phanstiel, Boyle et al. 2014, Warnes, Bolker et al. 2016, Wickham 2016, Ligges, Maechler et al. 2018).
The seqlogo program was used to visualize nucleotide biases (Crooks, Hon et al. 2004). Read subsetting based on sequence content was carried out using standard linux tools (grep, awk, etc).