A genome-wide screening for RNAi pathway genes in Acari


 BackgroundRNA interference (RNAi) is a highly conserved, sequence-specific gene silencing mechanism present in Eukaryotes. Three RNAi pathways critical for organismal development and survival are known, namely micro-RNA (miRNA), Piwi-interacting RNA (piRNA) and short interfering RNA (siRNA) pathways. Little knowledge exist about the genes involved in these pathways in Acari. Moreover, variable successes has been obtained in gene knockdown via siRNA pathway in functional genomics and management of Acari species. We hypothesized that the clue may be in the variability in the composition and the efficacy of siRNAi machinery among Acari. ResultsBoth comparative genomic analyses and domain annotation suggest that all the analyzed species have orthologs of genes that mediate gene silencing via the three RNAi pathways though gene duplication and/or loss have occurred in the different species. We also identified orthologs of Caenorhabditis elegans RNA-dependent RNA polymerase (RdRP) gene in all Acari species though no secondary Argonaute homologs that operate with this gene in siRNA amplification mechanism were found. This finding suggests that the siRNA amplification mechanism present in Acari may be distinct from that described in C. elegans. Moreover, the genomes of these Acari species encode no ortholog of C. elegans systemic RNAi defective 1 (Sid-1) that mediate systemic RNAi, suggesting that the phenomena of systemic and parental RNAi that has been reported in some Acari species probably occur through a different mechanism. Orthologs of RNAi spreading defective-3 (Rsd-3) gene and scavenger receptors namely Eater and SR-CI that mediate endocytosis cellular update of dsRNA in C. elegans and Drosophila melanogaster were found in Acari genomes. This result suggests that cellular dsRNA uptake in Acari is endocytosis-dependent. Detailed phylogenetic analyses of core RNAi pathway genes in the studied Acari species revealed that their evolution is compatible with the proposed monophyletic evolution of this group.ConclusionsTaken together, our in silico comparative analyses have revealed the potential activity of all three pathways in Acari. However, much experimental work remains to be done in elucidating the operating mechanisms behind these pathways in particular those that govern systemic/parental RNAi and siRNA amplification in Acari.

Number of gene orthologs involved in the biosynthesis of miRNAs and their characteristic conserved domains. Asterick (*) following species name abbreviations indicates that its genome is fully sequenced whereas no asterisk indicates that its genome is partially sequenced. ( # ) Indicates that the gene contains no conserved domains. With the same search approaches mentioned above, we identi ed homologs of Dicer1 (Dcr1) and Dicer2 (Dcr2) proteins in the genome assemblies of all the studied Acari species. Phylogenetic analysis based on the alignment of the two catalytic conserved Ribonuclease domains of all Dicer orthologs identi ed in these species along with orthologs of this protein in C. elegans (Ce), D. melanogaster (Dm) and an Acari species Dermatophagoides farina (Df) was carried out. Additionally, we searched for the presence of known conserved domains of Dicer proteins in each identi ed ortholog as outlined in the method's section. The analysis revealed that orthologs of Dicer1 and Dicer2 in Acariformes species clustered separately from those found in Parasitiformes species (Fig. 1A). Moreover, we found that Dicer1 orthologs identi ed in V. destructor (Vd), V. jacobsoni (Vj), T. urticae (Tu), I. scapularis (Is) and D. pteronyssinus (Dp), clustered together with signi cant bootstrap support (Fig. 1A) and have similar domain architecture to Dicer1 of both D. melanogaster (Dm-Dcr1) and C. elegans (Ce-Dcr1) (Fig. 1B). In the same vein, Dicer2 orthologs in these species clustered together and have similar domain architecture to Dm-Dcr2 with the exception of the ortholog identi ed in T. urticae, which clustered with Ce-Dcr1 (Fig. 1A, B). A single Dicer1 ortholog was found in all these species while several copies of Dicer2 were detected in the genome assemblies of V. destructor (2), V. jacobsoni (2) and D. pteronyssinus (2). As shown in Fig. 1A, duplicates of Dicer2 in V. destructor and V. jacobsoni clustered together while those identi ed in D. pteronyssinus did not. It is worthwhile to note that duplicates of Dicer2 in V. destructor, V. jacobsoni and D. pteronyssinus have similar domain architecture to Dm-Dcr2 (Fig. 1B).
Dicer1 ortholog found in T. mercedesae (Tm) clustered with Vj-Dcr1 and Vd-Dcr1 with signi cant bootstrap values (Fig. 1A) though it did not have similar domain architecture with Dm-Dcr1 and Ce-Dcr1 ( Fig. 1B). Homolog of Dicer2 in this species has a similar domain structure with Dm-Dcr2 and clustered closely with Vj-Dcr2 and Vd-Dcr2 (Fig. 1A, B). On the other hand, orthologs of Dicer1 and Dicer2 identi ed in E. maynei and S. scabiei clustered closely with Dp-Dcr1, Df-Dcr1, Dp-Dcr2 and Df-Dcr2 (Fig. 1A) though they differed signi cantly in domain architecture from Dm-Dcr1/Ce-Dcr1 and Dm-Dcr2, respectively (Fig. 1B). As shown in Fig. 1A   (IPR000999/PF00636) of Dicer orthologs identi ed in the studied Acari species using MAFFT [44]. The studied species include: Varroa destructor (Vd), Varroa jacobsoni (Vj), Tetranychus urticae (Tu), Ixodes scapularis (Is), Dermatophagoides pteronyssinus (Dp), Tropilaelaps mercedesae (Tm), Euroglyphus maynei (Em), Metaseiulus occidentalis (Mo) and Sarcoptes scabiei (Ss). Orthologs of Dicer proteins in Caenorhabditis elegans (Ce), Drosophila melanogaster (Dm) and Dermatophagoides farinae (Dp) were also included in this tree. The species names were abbreviated for convenience and the letters (a, b and c) following the names of some species indicate the number of gene copies found in their genomes. Names in "Red" and "green" colors on the tree are Parasitiformes and Acariformes species, respectively while those in "Blue" and "Purple" are D. melanogaster and C. elegans, respectively. The tree was constructed using PhyML [45], with the model recommended by Lefort et al. [46] under the Akaike information criterion (LG+G) with 500 bootstrap replicates. The domain architecture of Dicer proteins in (B) was generated by searching for known conserved domains in the Pfam database using both InterProScan [42] and HmmScan [43] and the letters (a, b and c) indicate the number of gene copies found in their genomes.
We identi ed homologs of Loquacious in the genomes of all Acari species investigated herein and several copies of this gene were found only in V. jacobsoni (6) and T. urticae (2) genomes (see Additional le). All orthologs of this protein found in these Acari species have the characteristic conserved dsRNA binding domain (IPR014720/PF00035). However, we did not nd either R2D2 protein nor its ortholog C3PO protein previously identi ed in the genome of the red our beetle Tribolium castaneum [33], in the genomes of any of the studied Acari species (see Additional le).
Argonautes are the core effector proteins of the RNA-induced silencing complexes (RISCs) in the cytoplasm. They belong to the Argonaute superfamily, which is functionally divided into three families based on their interactions with speci c RNA substrates [19]. Members of the Ago family interact with miRNAs and siRNAs, those of the Piwi family interact with piRNAs and those of the WAGO family, which are restricted to nematodes, interact with secondary siRNAs produced by RdRPs [19,36]. Members of the Ago family are further divided into two distinct groups based on their involvement in either the miRNA-or siRNA-directed gene silencing. For instance, in D. melanogaster, Argonaute1 (Ago1) is required for miRNAdirected gene silencing while Argonaute2 (Ago2) is required for siRNA-directed gene silencing alone [47].
Meanwhile, in C. elegans, Alg1 and Alg2 are required for miRNA-directed gene silencing and Rde1 and Ergo1 are required for siRNA-directed gene silencing [21,36]. A typical Argonaute protein is comprised of four domains: The N-terminal, PAZ, Mid and PIWI domain and the latter domain is the slicing domain that degrades the target mRNA transcript [19].
In this study, we found homologs of Ago family Argonautes in Parasitiformes and Acariformes species using the same search approaches mentioned above. A maximum likelihood tree based on the alignment of the PIWI domain showed that orthologs of the Ago family identi ed in Acari belong to two distinct groups that appear to be specialized in either mi-or siRNA-directed gene silencing (Fig. 2). It was intriguing to nd that some of the Acari orthologs clustered together with Dm-Ago1, Ce-Alg1 and Ce-Alg2, whereas others clustered closely with Dm-Ago2 and Ce-Rde1. Therefore, we referred to these orthologs as Ago1 and Ago2, respectively. Except in E. maynei whose genome encodes a single copy of Ago1 and Ago2, we found more copies of Ago2 than Ago1 in the genomes of the remaining Acari species. The tree further showed that orthologs of Ago1 and Ago2 in Acariformes species clustered separately from those identi ed in Parasitiformes species, with signi cant bootstrap values (Fig. 2). Figure 2 Phylogenetic analysis of Agonaute proteins. The tree was constructed based on the alignments of the conserved PIWI domain of Argonaute orthologs identi ed in the studied Acari species: Varroa destructor (Vd), Varroa jacobsoni (Vj), Tetranychus urticae (Tu), Ixodes scapularis (Is), Dermatophagoides pteronyssinus (Dp), Tropilaelaps mercedesae (Tm), Euroglyphus maynei (Em), Metaseiulus occidentalis (Mo) and Sarcoptes scabiei (Ss) using MAFFT [44]. Orthologs of these proteins in Caenorhabditis elegans (Ce), Drosophila melanogaster (Dm), Tribolium castaneum (Tc), Psoroptes ovis (Po) and Dermatophagoides farinae (Dp) were also included in this tree. The species names were abbreviated for convenience and the letters (a, b, c, d, e, f, g and h) following the names of some species indicate the number of gene copies found in their genomes. Names with "Red" and "green" colors on the tree are Parasitiformes and Acariformes species, respectively while those in "Blue" and "Purple" are D. melanogaster/T. castaneum and C. elegans, respectively. The tree was constructed using PhyML [45], with the model recommended by Lefort et al. [46] under the Akaike information criterion (LG + G + I + F) with 500 bootstrap replicates.
The search for the presence of known conserved domains revealed that almost all the orthologs of Ago1 in Acari have the same domain architecture with Dm-Ago1 and Ce-Alg1 and 2, with the exception of the single ortholog identi ed in E. maynei, which lacks the N-terminal domain (Fig. 3). Contrary to Ago1 homologs identi ed in Acari, which mostly shared the same domain architecture with orthologs of this protein in C. elegans and D. melanogaster, most of the Ago2 homologs in Acari had a domain architecture that was quiet distinct from those of both C. elegans and D. melanogaster (Fig. 3). The majority of them lack the Mid domain or both the N-terminal and the Mid domains. Only all the homologs of Ago-2 identi ed in I. scapularis (Ago-2a, b, c and d) and one homolog of this protein identi ed in M. occidentalis (Ago-2c) and T. urticae (Ago-2 h) have similar domain structure with Dm-Ago2. Figure 3 Schematic domain architecture of Argonaute proteins of the Ago family in Acari. The domain architecture of the proteins was generated by searching for known conserved domains in the Pfam database using both InterProScan [42] and HmmScan [43]. The letters (a, b, c, d, e, f, g and h) following the names of the gene (Ago1 and Ago2) indicate the number of gene copies found in in the individual Acari species.
Argonautes of the Piwi family such as Argonaute-3 (Ago3), Aubergine (Aub) and Piwi were found in the genomes of I. scapularis, T. urticae, V. destructor, V. jacobsoni, M. occidentalis and T. mercedesae, but not in genomes of D. pteronyssinus, E. maynei and Sarcoptes scabiei. We identi ed homologs of Ago3 in the genomes of I. scapularis, T. urticae, V. destructor, V. jacobsoni, M. occidentalis and T. mercedesae (see Additional le). Notably, several copies of this protein were found only in the genomes of I. scapularis (2) and V. jacobsoni (2). We further identi ed putative homologs of D. melanogaster-Aub and Tribolium castaneum (Tc)-Aub/Piwi only in the genome assemblies of I. scapularis and T. urticae (see Additional le). It is important to note that one homolog of Drosophila Aub/Piwi and Ago-3 were previously detected in the genome of the red our beetle, Tribolium castaneum (Tc) [33]. In our study, we identi ed one and two copies of Dm-Aub and Tc-Aub/Piwi, respectively, in I. scapularis and one and four copies of Tc-Aub/Piwi and Dm-Aub, respectively in T. urticae (see Additional le). All homologs of the Piwi family Argonautes identi ed in T. urticae and I. scapularis had the PIWI and PAZ domains, while those identi ed in V. destructor, V. jacobsoni, M. occidentalis and T. mercedesae had only the PIWI domain. Phylogenetic analysis placed homologs of Ago3, Aub and Piwi in Acari with orthologs of these genes in D. melanogaster and T. castaneum (Fig. 2).
We found homologs of the WAGO family Argonautes (Ce-PPW1, -PPW2, -SAGO1 and -SAGO2) in some Acari species studied herein using Proteinortho5 [41] (see Additional le). One and two homologs of Sago2 and Sago1, respectively, were identi ed in I. scapularis's genome. Also, one homolog of PPW2 was found in T. mercedesae while one and three homologs of Sago1 were identi ed in V. jacobsoni and V. destructor, respectively. In T. urticae, homologs of Sago1 (1) and PPW2 (1) were also found. Interestingly, most of the orthologs of the WaGO family Argonautes were the same homologs, which were highly similar to Ago family Argonautes (Ago2) identi ed in these species (represented as Is-Ago2a and Ago2b, Tm-Ago2c, Vd-Ago2b, Vj-Ago2b and Tu-Ago2c) (see Additional le). Since all the orthologs of the WAGO family Argonautes clustered with orthologs of the Ago family Argonautes, with signi cant bootstrap support (Fig. 2), we tentatively referred to them as members of the Ago family Argonautes that is Ago2 homologs. It is worth mentioning that all C. elegans WAGO family Argonautes used as query sequences in our study have both the PAZ and the PIWI domains.
Orthologs of the Ago and Piwi families Argonautes identi ed in all these Acari species were further examined for the presence of the conserved residues, Aspartate-Aspartate-Histidine/Aspartate/Lysine (DDH/DDD/DDK). These residues that are present within the catalytic PIWI domain presumably enable some Argonaute proteins to cleave the target mRNAs (reviewed in [48,49]). As shown on Table 2, the DDH catalytic residue was found in all the Argonautes identi ed in I. scapularis. Similarly, Argonautes of the Ago family identi ed in E. maynei's genome contained this motif. The only member of the Piwi family that is Ago3 identi ed in the genomes of M. occidentalis, T. mercedesae and V. jacobsoni did not have any of these residues in the sequences of their PIWI domains, whereas all the Argonautes of the Ago family identi ed in these species have the DDH residue. Also, Ago3 and one homolog of Ago-2 (that is Ago2d) identi ed in V. destructor lacked these conserved residues, though the remaining proteins of the Ago family had the conserved DDH residue. Also, Ago1 of S. scabiei did not have these residues, whereas Ago2 homologs did have. In T. urticae, all the Argonautes of the Ago family had the DDH motif, while only four out of the six Argonautes of the Piwi family had the DDH motif. In D. pteronyssinus, Ago1 and Ago2a have the DDH motif, while the remaining Argonautes that were Ago2b to 2 g had the DDD motif. Table 2 Catalytic residues of Argonautes of the Ago, Piwi and Wago groups in Acari species, Caenorhabditis elegans, Drosophila melanogaster and Tribolium castaneum.                                       Other components of the RISC complex notably Vasa intronic gene (Vig-1) and Tudor-staphylococcal nuclease (Tsn-1), that also contribute to the degradation of the target mRNAs [50,51], were also identi ed in the genome of the Acari species in this study, with the exception of E. maynei, which had Vig-1 but not Tsn-1 (see Additional le).

DsRNA uptake and spread
We searched for Sid-1, Rsd-2, Rsd-3, Rsd-6, Eater and SR-CI genes that are involved in cellular dsRNA uptake and systemic spread of the silencing dsRNA molecules as mentioned above. We were unable to identify genes similar to Sid-1 nor its orthologs SilA, SilB and SilC in the studied Acari species, previously found in the T. castaneum genome [33] in the studied Acari species. Similarly, orthologs of Rsd-2 and Rsd-6 were not found in any of the studied Acari species, though orthologs of Rsd-3 protein were present in all of them (see Additional le). Duplicates of this protein were only found in the genomes of I. scapularis (2), V. jacobsoni (5) and V. destructor (9) respectively. Interestingly, all the identi ed orthologs had a single conserved domain, the epsin N-terminal homology (ENTH) domain that has been shown to be su cient to mediate the transport of dsRNA molecules into both somatic and germ cells of C. elegans [29]. We also identi ed orthologs of Dm-Eater in the genomes of some Acari species except in those of D. pteronyssinus, E. maynei and S. scabiei. Again, duplicates of this protein were found in T. urticae (3), V. destructor (5) and V. jacobsoni (2). Orthologs of this protein in the Acari species have several epidermal growth factor (EGF)-like modules. Furthermore, one homolog of Dm-SR-CI gene was identi ed only in the genome of T. urticae. The Drosophila SR-CI receptor and its homolog had a characteristic conserved Meprin, A5 protein, and protein tyrosine phosphatase Mu (MAM) domain.

siRNA secondary ampli cation
Systemic and trans-generational gene interference rely on the ampli cation of the initial siRNA trigger for e cient gene knockdown throughout the treated organism. This mechanism for enhancing RNAi potency necessitates the action of a cellular RNA-dependent RNA polymerase (RdRP) protein, which has been extensively studied in C. elegans [52,53]. Interestingly, we found that all Acari species from the two sister lineages, Acariformes and Parasitiformes studied here, poses such RdRP proteins. Duplicates of this protein were found in the genomes of all studied species. Moreover, all the Acari RdRP proteins identi ed in this study including those previously identi ed in C. elegans (Ego-1 and Rrf-1) [37] and other Acari species e. g. D. farinae [54] and Psoroptes ovis [55] had the characteristic conserved RdRP domain [34]. In order to infer the evolutionary history of Acari RdRP proteins, we conducted a phylogenetic analysis based on the amino acid alignment of RdRP domain using C. elegans as outgroup, as described in the methods section. The phylogeny revealed that proteins of Acariformes (indicated with "green" color on Fig. 4) clustered strongly together (494/500 bootstraps) and separately from those of Parasitiformes. The Parasitormes proteins (indicated with "red" color on Fig. 4) also clustered strongly together (408/500 bootstraps). However, three out of the total four proteins identi ed in I. scapularis were found to be evolutionary closer to Acariformes proteins than to those of Parasitiformes (479/500 bootstraps) (Fig. 4). Ixodes scapularis (Is), Dermatophagoides pteronyssinus (Dp), Euroglyphus maynei (Em), Tetranychus urticae (Tu), Sarcoptes scabiei (Sc) and Psoroptes ovis (Po). The nematode species, Caenorhabditis elegans (Ce) was used as outgroup. The species names were abbreviated for convenience and the numbers (1, 2, 3, 4 and 5) following their names indicate the gene copies found in each species. Names with "Red" and "green" colors on the tree are Parasitiformes and Acariformes species, respectively while those in "Purple" are C. elegans. The tree was constructed using PhyML [45], with the model recommended by Lefort et al. [46] under the Akaike information criterion (LG + G + I + F) with 500 bootstrap replicates and is based on the amino acid alignment of the conserved RdRP domain (IPR007855/PF05183) using MAFFT [44].

Discussion
In this study, we aimed to identify genes that are associated with the three RNAi pathways in the genome assemblies of nine Acari species that belong to diverse Acari orders of agricultural, veterinary and medical importance [2]. Our ndings demonstrate that the genomes of all these Acari species encode orthologs of the genes that mediate the silencing process in the three RNAi pathways. Regarding the siRNA pathway, we also found orthologs of the genes that mediate the cellular uptake of dsRNA molecules and the ampli cation of the initial RNA molecule trigger. However, in the genomes of these Acari species we could not detect orthologs of genes that are required for systemic RNAi. In addition, we demonstrated that duplication and/or loss of gene orthologs involved in all three RNAi pathways have occurred in the Acari species across different lineages (Table 3). Gene duplication could be informative in predicting specialization or new gene functions whereas gene loss could help predict the contrasting effects of RNAi gene silencing in the different organisms as was already found in insects by Dowling et al. [56]. However, as only the genomes of T. urticae, M. occidentalis and I. scapularis have been fully annotated to date, we cannot completely state with conviction that gene loss occurred in the other species, whose genomes have been only partially annotated so far. In particular, our phylogenetic analyses of the core RNAi genes: Dicer, Argonaute and RdRP revealed that orthologs of these genes grouped according to their respective orders and those identi ed in Acariformes species clustered separately from those found in Parasitiformes species. This nding thus provides additional support for this taxonomic division in the subclass Acari [3]. Table 3 Summary of numbers of gene orthologs associated with the three RNAi pathways identi ed in the nine investigated Acari genomes. Asterick (*) following species name abbreviation indicates that its genome is fully sequenced whereas no asterisk indicates that its genome is partially sequenced.  Genes associated with the three RNAi pathways in Acari In this study, we identi ed orthologs of Drosha, a member of the Ribonuclease III family, its dsRNAbinding partner, Pasha and Exportin1, 2 and/or 5 in Acari genomes. This nding indicates that these proteins; which play vital roles in the initial step of miRNA biogenesis that takes place inside the nucleus are conserved in Acari as in other animal species [37,39,57]. Our phylogenetic analysis of Dicer proteins showed two Dicer isoforms with distinct functions in Acari: Dicer1 and Dicer2 involved in miRNA and siRNA biosynthesis, respectively (Fig. 1A, B) as was previous shown in other two Acari species: P. ovis [55] and D. farina [54]. These ndings also align with those of all insect species analyzed so far [33,56,58,59], but differ from those in organisms such as nematodes and mammals. In the latter species, a single Dicer1 isoform is involved in both siRNA and miRNA pathways [19,59]. The degree of functional specialization of Dicer proteins was further supported by the presence of Loquacious orthologs in the Acari species. This protein is indispensable for pre-miRNA processing by Dicer1 as deletion of this protein by RNAi resulted in the accumulation of pre-miRNAs in Drosophila S2 cells [40]. In contrast, R2D2 which we did not nd in any of the Acari genomes analyzed is apparently not an essential partner of Dicer2 in dsRNA processing as demonstrated in Drosophila R2D2-defective S2 cells by RNAi that exhibited no processing defect [60]. Functional specialization of Argonaute proteins was also apparent in Acari [19]. Of the three known Argonaute protein families in eukaryotes [17,19], our phylogenetic analysis revealed representatives of two families (Ago and Piwi) in Acari (Fig. 2). Within the Ago family, we further detected representatives of two distinct functional groups: Ago1 and Ago2 specialized in miRNA-and siRNAdirected gene silencing, respectively as was previous shown in two other Acari species: P. ovis [55] and D. farina [54]. The fact that these two gene families clustered separately from WAGO family Argonaute speci c to C. elegans with signi cant bootstrap value indicates that Acari, as insects [33], lack clear orthologs of WAGO family Argonautes in their genomes. The match that occurred between C. elegans WAGO family Argonautes and some orthologs of Ago family Argonautes (Ago2) in our study might imply similar domain structures and/or a common ancestor for both protein families rather than similar functions. This assumption is supported by the presumption that Argonaute protein families were all present in the last eukaryotic common ancestor and conserved their distinct functions as well as their similar domain architecture in the course of evolution [61]. However, the number of isoforms for each Argonaute family varies across the different eukaryotic lineages. For instance, humans have eight Argonaute proteins (four for each Ago and Piwi family), D. melanogaster has four Argonaute proteins (two for each Ago and Piwi family) and C. elegans has at least twenty six Argonaute proteins ( ve are Ago-like, three are Piwi-like and eighteen are WAGO-like) (reviewed in [49]). Our results also showed clear differences in the number of orthologs for each Argonaute family among the studied species (Table 3).
Even though the miRNA, siRNA and piRNA pathways regulate distinct biological functions within eukaryotic organisms, they may compete for shared effector proteins, especially the rst two pathways, and regulate each other at certain levels [17]. It is therefore important to gain a better understanding of the interconnectivity that may exist among the different RNAi cascades. The nding that functional separation of core effector proteins involved in the exogenous siRNA and endogenous miRNA/piRNA machineries that have occurred in Acari as in insects [33,56], supports the previous hypothesis that duplication and specialization of Dicers and Ago family Argonautes have occurred in Arthropods [59].
Duplication and specialization of these genes likely occurred early before the split of extant Arthropods into two monophyletic groups: Chelicerata and Mandibulata dating back to the Precambrian era [62]. However, further analysis of homologs of these genes in the other Arthropod lineages is needed to ascertain whether these phenomena are widely conserved or merelylineage-speci c.
This functional specialization of the effector RNAi proteins may suggest that there is no competition for these proteins among the RNAi cascades in Acari and probably in Insects as well. On the other hand, the pathways may be interlinked. Is it possible that the exogenous activation of the siRNA pathway via dsRNA treatment will affect the endogenous regulatory RNAi machineries (miRNA/piRNA/siRNA), which in turn may have a positive or negative effect on the silencing e ciency of the target gene in Acari? A recent study conducted in the pea aphid, Aphis pisum (Hemiptera: Aphididae) showed that the activation of the exogenous siRNA pathway following dsRNA treatment led to a signi cant change in the expression levels of both the exogenous siRNA and endogenous miRNA/piRNA associated genes thereby suggesting that crosstalk among these pathways can occur [63]. In fact, crosstalk among these three RNAi pathways have already been reported in a few eukaryotic species [17]. In addition, the Loquacious gene identi ed in the Acari species is known to be required in both the endogenous siRNA and miRNA machineries in D. melanogaster [64], thus indicating that these pathways may compete for this protein [17]. Therefore, future studies to understand the interconnection that exists among these RNAi pathways in Acari, which will expedite the e cient application of RNAi to both functional genomics and pest management are warranted.
Phylogenetic analysis of Dicer and Argonaute proteins in Acari revealed a high divergence of exogenous siRNA-associated Dicer2 and Ago2 proteins, but not of endogenous miRNA-associated Dicer1 and Ago1 proteins and piRNA-associated Piwi Argonaute, thereby suggesting a strong selective pressure on the exogenous siRNA core genes. For example, the genomes of V. destructor, V. jacobsoni, M. occidentalis and D. pteronyssinus even encode for more than one copy of Dicer2 (Table 3). The genomes of all the studied Acari species, except for E. maynei, encode for more than one copy of Ago2 (Table 3). This high divergence suggests that the ancestral role of exogenous siRNA machinery in antiviral immunity is conserved in Acari [22]. In fact some Acari species including V. destructor, V. jacobsoni, I. scapularis, T. mercedesae, T. urticae are known to be intimately associated with several viruses [65,66,67]. The ancestral siRNA-based antiviral defense is also conserved in insects [18] and plants [68,69]. In contrast, siRNA-mediated anti-viral immunity may be lost in vertebrates as a recent study provided evidence that their genomes do not encode homologs of siRNA-class Argonautes [70]. This may have occurred because of evolution of an interferon system regarded as the primer route of defense against viruses in vertebrates rather than on RNAi-based antiviral defense system [71,72].
Slicing activities of Dicer and Argonaute proteins The slicing activity of Dicer and Argonaute proteins is determined by the two consecutive Ribonuclease domains and PIWI domains, respectively [19]. The phylogenetic analysis of Dicer proteins based on the alignment of these domains showed that orthologs of these proteins clustered according to their functional relatedness to either siRNA or miRNA pathway with the exception of Dicer2 ortholog identi ed in T. urticae, which clustered closely to Dicer1 of C. elegans (Fig. 1B), suggesting that this protein in T. urticae may be involved in both pathways. Since very little is known about slicing activity of Dicer proteins in Acari, we cannot conclude how unique is the slicing mechanism in T. urticae. Perhaps, the analysis of the ATPase hydrolysis activity of the conserved Helicase domain of Dicer1 and Dicer2 in T. urticae may provide an insight into their cleavage preference of either miRNA or siRNA substrates. It has been reported that the processing of RNA substrates by some Dicer homologs is in uenced differently by its Helicase domain across organisms [73,74,75].
For most of the miRNA-and siRNA-class Argonautes in Acari, we found the highly conserved DDH or DDD catalytic motif in their PIWI domains, which suggest their slicing abilities. However, the catalytic motif of some Argonaute forms has been demonstrated to be required but not su cient for slicing the target mRNA transcripts (reviewed in [49]). Moreover, the majority of the orthologs of the Piwi family Argonautes identi ed in V. destructor, V. jacobsoni, T. mercedesae, M. occidentalis, I. scapulari and T.urticae lack the conserved DDH/DDD/DDK motif in their PIWI domains. Serine-Arginine-Glycine (SRG), Serine-Asparagine-Serine (SRS) and/or Asparagine-Aspartate-Histidine/Tyrosine (NDH/NDY) have substituted this motif ( Table 2). With these substitutions, it remains unclear whether the catalytic activity of these proteins is still intact.
Endocytosis mechanism for dsRNA uptake and systemic RNAi in Acari Our ndings suggest that cellular uptake of dsRNA in Acari is strongly endocytosis-dependent. In fact, the genomes of all the Acari species contained a clear ortholog of Rsd-3 gene, which is an integral component of systemic RNAi in C. elegans. In this organism, it facilitates the importation of dsRNAs in both somatic and germ cells via endocytosis [29]. In some Acari like I. scapularis, V. jacobsoni and V. destructor, we even found more than one copy of this gene (Table 3). Marr and colleagues also found a copy of this gene in P. ovis genome [55], suggesting that this is a general situation in Acari species studied so far. Homologs of this gene are also present in insect genomes e.g. D. melanogaster and T. castaneum [33]. Additionally, genes similar to scavenger receptors Eater and SR-CI that mediate endocytosis-dependent dsRNA uptake in S2 cells of D. melanogaster [31,32] were found in some Acari species (Table 3). Although they appear to be speci c to the hemocyte-like S2 cells of D. melanogaster, which display high rates of endocytosis when compared to other cell types, it is possible that they function in a similar manner in Acari cells. It is worth mentioning that the EGF-like modules and MAM domain present in orthologs of Eater and SR-CI receptors, respectively in D. melanogaster and the Acari studied herein are evolutionary conserved structures that have been reported in several extracellular adhesive receptors that function in phagocytosis, a form of endocytosis in mammals, nematodes and invertebrates [76,77]. Our results corroborate with results of a previous study, which also demonstrated that the cellular uptake of dsRNA in the tick Haemaphysalis longicornis is also endocytosis-dependent [78].
Regarding the systemic transfer, our study revealed that the Sid-1 gene that mediates systemic RNAi in C. elegans somatic and germ cells [27,28] or its orthologs found in T. castaneum [33] are both absent from the genomes of all the studied Acari species, including that of P. ovis genome [55]. The apparent absence of Sid-1 gene in Acari may imply that systemic RNAi reported in M. occidentalis [79], T. urticae [80] and I. scapularis [81] occurred via a different mechanism. In fact, it has been found that in insects, systemic RNAi can occur independently of this gene [33].
siRNA ampli cation mechanism in Acari The potency of RNAi gene silencing throughout the organism depends on an ampli cation mechanism mediated by an RNA-dependent RNA polymerase gene as mentioned above. In this study, we con rmed the presence of homologs of RdRP gene along with their characteristic conserved domain in the genomes of all the Acari species (Fig. 4) as was commonly found in other Arthropods including insects, crustaceans, chelicerates and myriapods [35]. We also con rmed the fact that duplications of this gene has occurred in almost all the studied species, except in S. scabiei (see Additional le). Previous studies also reported the presence of several copies of this gene in the genomes of P. ovis [55] and D. farina [54].
Overall, this and previous genomic studies indicate that RdRP ampli cation activity seems to occur in Acari species studied to date, which may result in the production of secondary siRNAs for e cient gene silencing throughout the organisms as was suggested for Arthropod species in which this gene was found [82]. This is not surprising since systemic and trans-generational gene silencing have been reported in the few Acari species mentioned above. Interestingly, in the siRNA ampli cation mechanism described in the model organism C. elegans [36,52,53], WAGO family Argonautes associate with secondary siRNAs produced by RdRPs to increase the silencing of target mRNA transcripts [36,52]. Yet, our analysis of Argonaute proteins did not reveal any clear orthologs of WAGO Argonautes in any of the Acari species. This ndings supports the recent suggestion of Pinzon and colleagues that the role of RdRP genes in siRNA ampli cation is not conserved and probably varies among different animal species [35]. It remains to be determined, what is the siRNA ampli cation mechanism in Acari.

Conclusions
Using genomic data annotated at the protein level, we provide for the rst time information on the putative proteins in nine Acari species that maybe involved in the gene silencing processes of the three RNAi pathways in both Acariformes and Parasitiformes lineages. It is clear that Acari, like insects have evolved a set of specialized proteins for gene silencing in the three RNA pathways and an endocytosis mechanism for cellular dsRNA uptake. However, the mechanisms and molecular basis through which siRNA ampli cation and systemic/parental RNAi occur in Acari, still remain a puzzle and require follow up by structural and functional studies. A complete annotation of genome sequences of the Acari species studied herein as well of additional studies from other families might not only shed light into the underlying mechanism and molecular basis of siRNA ampli cation and systemic RNAi, but no doubt also Page 44/57 uncover additional RNAi pathway genes. So far, our results suggest that all the evaluated Acari species have a potential for active exogenous siRNA-directed gene silencing machinery, though difference in gene knockdown is likely to occur given that the number of copies of core RNAi genes vary among the species. Our ndings provide a basis for exploiting siRNA pathways in functional genomics and in the management of economically important species, especially species like T. mercedesae, D. pteronyssinus and E. maynei in which RNAi-gene silencing has not yet been exploited. Moreover, additional studies are needed to elucidate the regulatory mechanisms behind RNAi pathways as well as any possible cross talk mechanisms.
In the genome assembly of individual organism, we searched for RNAi pathway genes belonging to the following functional groups: small RNA biosynthesis genes, Argonautes and RISC components, siRNA ampli cation protein and dsRNA uptake and spread proteins as summarized in the Additional le using Proteinortho5 (https://www.bioinf.uni-leipzig.de/Software/proteinortho/) [41] and/or Blastp search against the NCBI database. As query sequences, we used orthologous proteins within the genomes of P. ovis, D. farinae, D. melanogaster, T. castaneum and C. elegans. The Amino acid sequences of the query genes from the rst four species were retrieved from the NCBI database, while those from the forth organism were retrieved from the WormBase database. The accession numbers of the query sequences are available in the Additional le. Only orthologs that yielded top blastp hits were retained. We further veri ed these hits by Blastp search in NCBI and annotated them as potential RNAi genes only when the Blastp search provided the query sequence as the best hit. Figure 2 Phylogenetic analysis of Agonaute proteins. The tree was constructed based on the alignments of the conserved PIWI domain of Argonaute orthologs identi ed in the studied Acari species: Varroa destructor (Vd), Varroa jacobsoni (Vj), Tetranychus urticae (Tu), Ixodes scapularis (Is), Dermatophagoides pteronyssinus (Dp), Tropilaelaps mercedesae (Tm), Euroglyphus maynei (Em), Metaseiulus occidentalis (Mo) and Sarcoptes scabiei (Ss) using MAFFT [44]. Orthologs of these proteins in Caenorhabditis elegans (Ce), Drosophila melanogaster (Dm), Tribolium castaneum (Tc), Psoroptes ovis (Po) and

Figures
Dermatophagoides farinae (Dp) were also included in this tree. The species names were abbreviated for convenience and the letters (a, b, c, d, e, f, g and h) following the names of some species indicate the number of gene copies found in their genomes. Names with "Red" and "green" colors on the tree are Parasitiformes and Acariformes species, respectively while those in "Blue" and "Purple" are D. melanogaster/T. castaneum and C. elegans, respectively. The tree was constructed using PhyML [45], with the model recommended by Lefort et al. [46] under the Akaike information criterion (LG+G+I+F) with 500 bootstrap replicates.

Figure 3
Schematic domain architecture of Argonaute proteins of the Ago family in Acari. The domain architecture of the proteins was generated by searching for known conserved domains in the Pfam database using both InterProScan [42] and HmmScan [43]. The letters (a, b, c, d, e, f, g and h) following the names of the gene (Ago1 and Ago2) indicate the number of gene copies found in in the individual Acari species. was used as outgroup. The species names were abbreviated for convenience and the numbers (1, 2, 3, 4 and 5) following their names indicate the gene copies found in each species. Names with "Red" and "green" colors on the tree are Parasitiformes and Acariformes species, respectively while those in "Purple" are C. elegans. The tree was constructed using PhyML [45], with the model recommended by Lefort et al. [46] under the Akaike information criterion (LG+G+I+F) with 500 bootstrap replicates and is based on the amino acid alignment of the conserved RdRP domain (IPR007855/PF05183) using MAFFT [44].

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download. Additional.xlsx