Rewired RNAi-mediated genome surveillance in house dust mites

House dust mites are common pests with an unusual evolutionary history, being descendants of a parasitic ancestor. Transition to parasitism is frequently accompanied by genome rearrangements, possibly to accommodate the genetic change needed to access new ecology. Transposable element (TE) activity is a source of genomic instability that can trigger large-scale genomic alterations. Eukaryotes have multiple transposon control mechanisms, one of which is RNA interference (RNAi). Investigation of the dust mite genome failed to identify a major RNAi pathway: the Piwi-associated RNA (piRNA) pathway, which has been replaced by a novel small-interfering RNA (siRNA)-like pathway. Co-opting of piRNA function by dust mite siRNAs is extensive, including establishment of TE control master loci that produce siRNAs. Interestingly, other members of the Acari have piRNAs indicating loss of this mechanism in dust mites is a recent event. Flux of RNAi-mediated control of TEs highlights the unusual arc of dust mite evolution.


Absence of Piwi proteins in dust mite genome
We obtained a genome sequence for the American house dust mite Dermatophagoides farinae using Illumina and PacBio platforms. The HGAP pipeline was used through PacBio SMRT analysis portal to filter and assemble PacBio reads, which resulted in 1,828 contigs producing a total length of 93,777,723 bp [32]. Then, Illumina reads were used to connect and extend the PacBio contigs using SSPACE scaffolding, which produced a 93,804,520 bp assembly in 1728 scaffolds [33]. After removal of bacterial contamination, the final contig number was reduced to 1706 with a N50 read length of 19,371. The assembled and filtered final genome was~92 Mb compared to a 53 Mb genome that was previously reported [34]. Using mRNA-seq datasets we annotated~18,500 transcripts through the Cufflinks program [35]. 47% of the genic transcripts exhibited similarity to S. scabiei and/or D. melanogaster protein coding genes or to the NCBI conserved domain collection (Materials and Methods) (S1 Table) [34,36]. Ago/Piwi proteins were identified in the D. farinae genome using RNA-seq annotations and amino acid sequences of seven Ago and six Piwi proteins from Tetranychus urticae-the closest relative of D. farinae with a high-quality genome and experimentally supported annotations [29]. Eight confident Ago homologs were found (Ago1-GenBank ID: KY794591, Ago2-GenBank ID: KY794592, Ago3-GenBank ID: KY794593, Ago4-GenBank ID: KY794594, Ago5-GenBank ID: KY794595, Ago6-GenBank ID: KY794596, Ago7-GenBank ID: KY794597, Ago8-GenBank ID: KY794598). Ago proteins from T. urticae, D. melanogaster, C. elegans, and Ascaris suum were compared to D. farinae Agos using amino acid sequences of Paz, Mid, and Piwi domains (Fig 1A). Our phylogenetic analysis recovered two Ago family members likely involved in miRNA (DfaAgo1) and siRNA (DfaAgo2) pathways [37]. The remainder belong to a divergent clade specific to dust mites (DfaAgo3- 8). Surprisingly, none of the Agos from D. farinae belong to the Piwi clade. We examined D. farinae Agos for the presence of slicer motifs. The DEDH slicer motif, which is common in metazoan Ago and Piwi proteins, was found in DfaAgo1 (miRNA) and DfaAgo2 (siRNA). The divergent Agos have an uncommon DEDD catalytic motif (S2 Fig). Orthologs containing a DEDD motif can be found in scabies (S. scabiei), social spiders (Stegodyphus mimosarum), and in C. elegans Ago family members of unknown function; which emphasizes the unusual nature of this Ago clade [36,38,39].

Loss of piRNAs in dust mite
A dust mite small RNA library of nearly 400 million reads was generated to investigate whether piRNA-class small RNAs could be identified (S1 Text) (S1 Fig). To accommodate the repetitive nature of piRNA targets all mapping events were captured for reads that mapped fewer than 100 times. An overall rate of~80% mapping was observed with~0.69% discarded due to mapping >100 times (S1 Fig). Next an algorithm that determines small RNA read overlap probabilities in mapping data was used to characterize biogenesis of dust mite small RNAs [40]. When applied to either all mapping or mapping in discreet size ranges no clear bias for 10nt overlapping reads was uncovered, showing an absence of ping pong processing (Fig 1B). Instead a strong signal seen in a register 2nt shorter than the length of read sizes. This is congruent with 2nt overhangs left by Dicer cleavage. Overlaps seen in dust mites starkly contrast with those seen in spider mites and Drosophila. In spider mites, ping pong signatures could be seen in longer reads (23-28nt) and the dicer-associated 2nt register in shorter (20-22nt) reads ( Fig 1C). Likewise, in RNAs sequenced from Drosophila female bodies a prominent ping pong signature is evident (Fig 1D). siRNA processing was not evident when considering whole genome mapping, but could be seen in a group of Drosophila IDEFIX retroelements that had  Table). Drosophila endo-siRNAs are a relatively small proportion of total small RNAs, and are frequently produced from inverted repeat loci which are not captured by the overlap probability calculation used here [41]. Moreover, this Dust mite RNAi genome integrity highlights a correlation between the presence of Rdrp in spider mites and an expanded population of Dicer products. Together this shows a dramatic departure in the composition of dust mites small RNA populations relative to those in the Piwi protein possessing spider mite. The difference is even more stark when comparing dust mites to the more distantly related fruit fly. The configuration of RNAi in spider mites is likely ancestral due to the similarities to Drosophila, which is supported by clear orthology of spider mite Piwis to distinct Drosophila pingpong partners: dmeAgo3 and Piwi/Aub (Fig 1A). Thus, RNAi pathways have diverged in dust mites and appear to be dominated by siRNA-like Dicer products, and lack the signature of amplifying ping pong piRNAs [42].
To functionally characterize dust mite small RNAs, we sought to identify genomic loci that generate and/or are targeted by these transcripts (S1 Fig). To achieve this, annotation of the dust mite genome was extended to find ncRNAs and repetitive elements using Repeatmasker. Additionally, non-miRNA, small RNA producing loci were annotated that exhibited >1000 read density and were longer than 200nt. The identities of regions were determined using blas-t2go [43]. Nearly a third of the loci were rRNA or mRNA. The remainder showed homology to either TEs or lacked similarity to known sequences. Together this permitted segmentation of the dust mite genome into mRNA, TE, rRNA, tRNA, snRNA, and unknown small RNAmapping loci (S2 Table). Small RNA reads were then mapped to these regions using multiple mapping conditions described above, as well as unique mapping. To ensure multi-mapping events were specific to loci groups, datasets were cleaned before mapping by removing reads that mapped to non-target genomic features (S1 Fig).
Multi-mapping alignments showed considerable enrichment at TEs relative to other classes, consistent with their repetitive nature (Fig 2A). Both multi-and uniquely mapping TE reads also exhibited lower strand bias with only a single locus showing 100% bias after unique mapping (S4 Fig). This is consistent with processing from dsRNA. Higher bias was seen at other loci, suggesting that some mapping events may be due to capture of RNA degradation fragments and not functional small RNAs. This was supported when overlap probabilities were calculated; which, with exception of TEs and mRNAs, did not show consistent processing signatures (S5 Fig). This includes the unknown loci, suggesting these transcripts may be degradation products of uncharacterized ncRNAs and are not generally siRNA or piRNA class small RNAs. Closer inspection of per locus overlaps did show Dicer processing at a minority of loci (S6 Fig). There was no clear ping pong processing at unknown loci. Small RNA mapping coverage was calculated per locus to understand siRNA production from TEs and mRNAs ( Fig 2B). On average, small RNA coverage was even across TE loci, while mRNAs had greater coverage at transcript 3' ends. This pattern at mRNA loci is suggestive of cis-NAT siRNAs [44]. Depth of coverage at TE loci varies, showing that active targeting is occurring at a subset of loci.
The absence of purely single-stranded small RNA producing loci that have homology to TE sequences suggests that dust mites also lack a Zuc-dependent piRNA-like pathway that is involved in genome surveillance. This does not rule out the existence of dual strand piRNA clusters; however, piRNAs produced from these loci are found to participate in the ping pong cycle, which we did not observe [45]. These data suggest that the piRNA pathway has been lost in dust mites and control of TE's is likely under the purview of a siRNA-like pathway.

siRNAs facilitate genome surveillance in dust mite
To investigate the role of dust mite small RNAs in genome surveillance we compared the biogenesis of TE-associated small RNAs to those found in spider mites. The size distribution of genome-aligned dust mite small RNAs is unimodal with a peak at 24nt, versus a bimodal distribution in spider mites ( Fig 3A). When TE-mapping reads are examined, the 24nt sized RNAs in dust mite were enriched by 10%, while in spider mites only larger size range RNAs were found (Fig 3B). In other locus groups, less read size bias was observed, consistent with the heterogeneity seen in strand bias and read overlap probabilities, further reinforcing that generally non-TE loci do not produce small regulatory RNAs (S7 Fig). Next we looked at the 5' nucleotide bias and found that dust mites TE siRNA reads have an equal prevalence of T and A residues versus spider mites where there was striking over representation of T ( Fig 3C). Then we examined per locus read size distribution and overlap probabilities to assess whether Dicer processed~24 nt small RNAs are common across dust mite TE loci ( Fig 3D). All loci exhibited mapping of predominantly 24 nt reads, and in the most prevalent size ranges (23-26nt) a clear pattern of overlaps could be seen that is consistent with Dicer processing ( Fig  3D). This contrasts with a similar analysis in spider mite where a ping pong signature was seen across all TEs. Together this suggests siRNAs are the main RNAi-based mode of controlling TEs in dust mites, accommodating the apparent loss of piRNAs. This is a clear departure from spider mites where stereotypical piRNAs target TEs.
In the D. farinae genome we found three Dicers (DfaDcr1-GenBank ID: KY794588, DfaDcr2-GenBank ID: KY794589, DfaDcr3-GenBank ID: KY794590). DfaDcr1 is a close ortholog of Arthropod miRNA-producing dicer (S8 Fig). The other two Dicer proteins are related to family members in other mites and lophotrochozoans, and are unrelated to Arthropod Dicer2 or nematode Dicer (S8 Fig). Unexpectedly, DfaDcr1 possesses an ATP binding helicase domain, which is implicated for processing of long dsRNA (S9 Fig) [46]. The more divergent Dicers, DfaDcr2 and DfaDcr3, lack both DUF283 and dsRNA binding domains, and have divergent PAZ domains (S9 Fig) [46][47][48]. Together this suggests that mites, and possibly other chelicerates, possess ancient Dicer biology present in basal protostomes that was lost both in nematoda and pancrustacea (insects and crustaceans). To verify whether TEs are controlled by Dicer-produced siRNAs we sought to inhibit the activity of dust mite Dicer proteins. To generate loss of Dicer function we elicited RNAi against each Dicer by feeding mites cognate dsRNA (Fig 4). Dust mites tolerate being soaked for several hours in aqueous solution, which they can be observed to ingest after 30 mins ( Fig  4A). Small RNAs (20-27nt) derived from dsRNA can be recovered from soaked mites ( Fig 4B).  Knockdown of target genes can also be observed (Fig 4C-4K). Depletion by RNAi of each DfaDcr protein resulted in derepression of multiple TEs ( Fig 4L) (S10 Fig). A strong effect was seen with loss of DfaDcr1 and DfaDcr2 function. The presence of processive helicase activity in DfaDcr1 suggests that long dsRNAs could be substrates. This combined with the lack of dsRNA binding motifs in DfaDcr2/3 suggests DfaDcr1 has a unique capacity to process dsRNA (S9 Fig), and therefore it is unsurprising that it has a significant role in the control of TEs (Fig 4L). Loss of DfaDcr2 showed a greater effect on TE expression compared to DfaDcr3. How these atypical Dicer proteins function is unclear; however, residues in the DfaDcr3 PAZ differ significantly from those in DfaDcr2 PAZ suggesting non-overlapping roles in the metabolism of dust mite small RNAs (S9 Fig). These results are consistent with reports that psoroptid mites are sensitive to dsRNA soaking, resulting in gene knockdown [49,50].
Investigation of RNAi in dust mites revealed loss of the piRNA pathway and replacement by siRNAs. This is similar to observations in nematodes and flatworms [21,22]. The loss of piRNA activity in dust mites, nematodes, and possibly in flatworms may be tolerated due to compensation by amplifying siRNAs produced by Rdrp [21,51]. The collective function of dust mite Rdrps, however, appears to be distinct from nematodes, as only processive versions are present, suggesting the de novo siRNA pathway may not be present in mites (S11 Fig). Substantial Rdrp activity does appear to be present in dust mites; dsRNA soaking results in elevation of target mRNA when reverse transcription is carried out with random hexamers (Fig 4E, 4G, 4I and 4K) but not oligo dT (Fig 4D, 4F, 4H and 4J). Increase of transcript abundance was not due to the presence of ingested dsRNA as the region cloned to generate dsRNA was distinct from the qPCR amplicon (S10 Fig). Random priming will capture Rdrp products, while oligo dT will only hybridize to the initial transcript. For all the genes tested an elevation of cognate transcripts could be observed after random priming that were poorly recovered from Oligo dT primed cDNA.

Cataloging restricted sequences in siRNA producing master loci
Dust mites differ from nematodes that lost piRNAs in the organization of siRNA producing loci. A key feature of piRNA biology is the cataloging of restricted sequences into master loci. In nematode lineages lacking piRNAs, master loci also appear to be absent [21]. This is not the case in dust mites (Fig 5A). Three loci were discovered that span 62 kb, contain sequences from multiple varieties of TEs, and exhibit homology to 70% of TE mapped small RNAs ( Fig  5B). Two of the loci, ML-283 and ML-95, appear to be generated by duplication; however, some sequence divergence indicates they are distinct loci. Similar regions could not be found in the S. scabiei genome [52]. Though, poor conservation is a characteristic of piRNA master loci [53]. The dust mite loci appear to be generated from a dsRNA precursor as both strands of the loci show similar rates of read mapping (Fig 5A). We found a tendency for 2nt overhangs along with little evidence for nucleotide bias (S12 Fig). The loci were inspected for common motifs using the meme suite [54]. Motifs recovered were primarily simple repeats with none being shared between loci suggesting dust mite master loci don't possess elements like the Ruby motif which is central to directing piRNA transcription in C. elegans [55]. Following knockdown of each of the individual dust mite Dicers significant (>80%) reduction in siRNAs with derf1 dsRNA (upper panel) and coomassie staining of the membrane (lower panel). Animals were soaked for 30 min and after 4 days lysates were prepared. D-K. qPCR for dust mite transcripts, all experiments were performed at least three times. Values represent four technical replicates. Reverse transcription was carried out with either oligo dT (D, F, H, J) or with random hexamers (E, G, I, K). Target transcripts were derf1 (D,E), dcr1 (F,G), dcr2 (H, I), and dcr3 (J, K). Cntrl represents no treatment, and KD soaking in the indicated dsRNA. L. Increased expression of numerous TE's (S2 Table) following RNAi against three dust mite Dicers relative to untreated control. Error bars represent SEM. https://doi.org/10.1371/journal.pgen.1007183.g004 Dust mite RNAi genome integrity exhibiting homology to these regions was observed, indicating a dependence on the activity of all dust mite Dicers for biogenesis (Fig 5C). Detection of the siRNAs was accomplished with a combination of oligonucleotide probes complementary to sites of highest small RNA density in the three master loci (S1 Text). They also have homology to other regions of the genome, specifically TEs. Thus, the Dicer sensitive siRNAs include master loci derived primary siRNAs Dust mite RNAi genome integrity and potentially secondary siRNAs generated from processed TE transcripts. This is consistent with loss of TE control after knockdown of each Dicer (Fig 4L). However, there is a clear difference in the magnitude of TE expression, which may point to roles for dust mite Dicer proteins outside the production of siRNAs and to involvement in targeting of TE transcripts. This could be similar to limiting of latent viral infection by Drosophila Dcr2 [56].
Next, we sought to characterize terminal moieties of master loci associated siRNAs through biochemical tests to gain greater insight into their biogenesis (Fig 5D and 5E). The primary goal was to determine if the siRNAs had characteristics of Dicer cleavage: 5'-monophosphates and 3'-OH groups. β-elimination showed a shift to a lower molecular weight indicating an unmodified 2'OH; therefore, unlike Drosophila Ago2 endo-siRNAs or C. elegans Prg-1 associated small RNAs, dust mite siRNAs are not 2'-OH methylated (2'OMe) (Fig 5D) [57,58]. Next, we identified groups on 5' ends of small RNAs using the 5' monophosphate specific terminator ribonuclease. After treatment, a 50% reduction in siRNAs could be observed (Fig 5E). Degradation by terminator could be abrogated by prior treatment with calf intestinal phosphatase (CIP). There is a noticeable lag in siRNA gel migration following CIP treatment, which is consistent with removal of 5' phosphate groups and loss of charge. These results also reinforce the absence of a de novo siRNA pathway. Small RNAs produced by non-processive Rdrps in C. elegans have 5' triphosphate groups. While treatment with terminator did not completely eliminate siRNAs there was no observable change in migration. If the remaining small RNAs were spared due to the presence of trisphosphate groups there would be shift towards a smaller molecular weight, relative to untreated. Together, dust mite master loci associated siRNAs appear to be Dicer products arising from a dsRNA precursor, possess the expected 5'-monophosphate, but differ from insect endo-siRNAs due to the absence of 2'-OMe groups. We were able to identify a dust mite gene with similarity to Hen1 methyltransferase proteins; however, inspection of potential open-reading frames revealed the absence of a common motif involved in recognition of 2 nt 3' overhangs characteristic of Dicer products (S12 Fig). This likely explains the lack of 2'-OMe groups on dust mite siRNAs.

DNA methylation is not involved in dust mite TE control
Extent of DNA methylation in CG widely varies across insect clades and can be as high as 40% in roaches, while other groups, like flies, show little evidence for this modification [59]. Here we investigated whether this epigenetic control mechanism is a component of TE control in dust mites, as the genomes of nematodes and platyhelminths that lack the piRNA pathway are frequently modified by cytosine methylation [21,60]. Dust mites differ from these organisms, as evidence for this modification seems minimal and it is not enriched at TE loci ( Fig 6A). Indeed, bisulfite sequencing showed potential CG and CHG methylation is underrepresented in TE sequences, despite these sites occurring at the same rate as other genomic loci. Furthermore, the overall rate of DNA methylation (0.5%) was very low, suggesting this base modification is not a major feature of dust mite chromatin regulation. Moreover, we found a single DNA methyltransferase in the D. farinae genome, a Dnmt1 homolog (Fig 6B and 6C). It is likely a pseudogene as it appears to be truncated and shows little evidence of expression. This further highlights the distinct, derived nature of small RNA-mediated genome surveillance in dust mites.

Discussion
This work provides insight into the elaborate nature of RNAi in chelicerates, many of which appear to have both Piwi proteins and Rdrps [29,30,39]. Loss of the piRNA pathway in dust mites probably occurred in the parasitic ancestor. Inspection of the scabies mite genome similarly failed to uncover Piwi proteins (S13 Fig) [36]. Members of the divergent dust mite Ago family; however, were found. Indeed, a deeper inspection of scabies mite RNAi factors uncovered further similarities to dust mites (Table 1). Thus, absence of the piRNA pathway in  Dust mite RNAi genome integrity dust mites is likely a consequence of descending from an ancestor that underwent dramatic genome changes, potentially during the acquisition of a parasitic life style. This highlights plasticity of RNAi pathways and how clade-specific biology might impact evolution of RNAi technologies. Dust mites exhibit a highly distinct RNAi biology, possessing both novel and ancient effectors that haven't been studied in popular ecdysozoan model organisms. Indeed, there seems to be wholesale changes to the small RNAome of these organisms. Dicer produced siRNAs are an unusually common feature of the dust mite small RNA populations, comprising approximately three-fourths of all small RNA species. This contrasts with many other organisms where microRNA-class small RNAs are the archetype. Dust mite siRNAs are, at least in part, involved in genome surveillance. They target TE's and depletion of Dicer proteins causes derepression of these elements. Control of TE's is typically carried out by piRNAs in flies, from which dust mite siRNAs are distinct. A common feature of nearly all piRNAs is a "U" residue at the first position. We do not observe this in any subset of dust mite siRNAs. Furthermore, welldescribed modes of piRNA biogenesis found in Drosophila and C. elegans are absent in dust mites. Loss of piRNAs seems specific to psoroptidian mites, as they are clearly present in other Acari, like spider mites. The divergent nature of dust mite siRNAs is particularly apparent in the absence of 2'-OMethylation of siRNAs-a common feature of siRNAs and piRNAs in other organisms. Interestingly, scabies mites also lack the requisite Hen-1 protein [36]. Inspection of syntenic regions of the dust mite and scabies mite genomes showed rearrangements at this locus, potentially linking the loss of this activity to the evolution of Psoroptidia-specific Ago proteins (S13 Fig) (Table 1). The highly divergent RNAi pathways of dust mites provide an evolutionary perspective not only on the utility of small RNAs to acquire roles in genome surveillance, but also that the precise mechanism may not be that important. This is supported by relatively similar composition of classes of TE's in spider mites, dust mites, and scabies mites (S15 Fig). While similar classes were observed their locations and specific identities are distinct. Furthermore, this indicates that the collection of dust mite TEs analyzed in this study accurately represent the overall TE population.
Flux of small RNA pathways correlates with evolutionary innovation; for example, higher arthropods lost Rdrp in favor of piRNA control of TE [61]. This also occurred when vertebrates diverged from basal chordates [62]. In both cases, loss of Rdrp accompanied innovation in body plan and sensory organs. In vertebrates, whole genome duplication occurred twice following descent from a Rdrp expressing chordate ancestor, affirming a period of genome instability [62]. TE mobilization may be fortuitous for adaptation, and dramatic evolutionary changes may require extreme events such as perturbation of surveillance mechanisms.

Genome assembly pipeline
The dust mite genome was assembled using reads produced by PacBio and Illumina platforms. The initial assembly was generated by PacBio HGAP. Illumina reads were preprocessed in three steps before using them for extending PacBio contigs: a) Using Trimmomatic [63], from both ends of reads, nucleotides with base quality lower than 15 were removed. b) Using FastUniq [64], duplicate pairs were removed from the PE library, and c) SOAPec [65] was used to correct read error [64,65]. Any initial genome sequence has bacterial contamination due, at least, to the presence of gut microbiota in DNA isolates. To remove bacterial DNA sequences from D. farinae genome sequence, 4,864,367 Bacterial genome sequences [66] were downloaded from RefSeq database at: ftp://ftp.ncbi.nih.gov/refseq/release/bacteria and a blast database was created using the sequences [66,67]. All the contigs were blasted against the created bacterial genome database to check bacterial contaminations in the sequenced contigs. Then the matched percentages were calculated for each of the contigs. If the matched percentages were higher than 10% of an individual contig length, the contig was considered as contaminated by bacterial DNA and was discarded. After this process, our final contig number was reduced to 1706, N50 Read Length of 19,371 with the total length of 91,947,272 bp. Finally, a published dust mite genome [34] was compared to our assembled contigs using QUAST [34,68]. 79.3% bases of the reference genome could be aligned in the new assembly.

Transcript annotations
Using available mRNA-seq datasets [34], transcripts were identified by the Tuxedo suite. Initial mapping with Tophat was followed by transcript annotation with cufflinks [69]. Transcript similarity was estimated using Blast2Go.

Small RNA analysis
Total RNA isolated via the trizol method from bulk collected dust mites in order to capture life stages of D. farinae. Small RNAs were cloned from total RNA with an Illumina small RNA truseq kit, and sequenced on the Illumina NextSeq platform. The dataset was comprised of nearly 400 million reads. Quality of the sequenced library was assessed by FastQC tool and the small RNA reads were analyzed using a custom pipeline (S1 Fig) [70].

dsRNA soaking of mites and northern blotting
Mites collected with the salt bath method were suspended in a solution of dsRNA dissolved in nuclease free water (S1 Text). After 6 hours, animals were washed in water and dried on filter paper. After that the animals were kept in 23˚C with relative humidity of 80%. After two days, total RNA was extracted using trizol method and resolved in a 12.5% urea-polyacrylamide gel. When animals were fed unlabeled dsRNAs, RNAs were transferred to nylon membranes and subject to northern blotting as previously described (S1 Text) [56]. If radiolabeled RNAs were fed, gels were directly exposed to phosphoimager screens.

β-elimination
20 μg of total RNA was oxidized at room temperature in borax/boric-acid buffer (60 mM borax and 60 mM boric acid-pH 8.6) containing 80 mM NaIO4 for 30 min. β-elimination reaction was carried out for 90 min using 200 mM NaOH at 45˚C. Following precipitation, RNA was resolved on a 12.5% urea-polyacrylamide gel, and subject to northern blotting as previously described [56].

CIP and terminal exonuclease treatment
20 μg of total RNA was used for each of reaction. Terminator exonuclease (epicenter) was added to one tube and the tube was incubated at 30˚C for 60 minutes. After that the reaction RNA was purified by organic extraction protocol [71]. In the second condition, 1 μl CIP (Calf intestinal phosphatase, NEB) was added and incubated at 37˚C for 30 min. Terminator exonuclease was added followed by a second incubation at 30˚C for 60 minutes. Precipitated RNAs were resuspended in loading buffer and resolved on a 12.5% urea-polyacrylamide gel, and subjected to northern blotting as previously described [56].

Methylation analysis
A Methyl DNA seq library was created with Illumina Methyl-seq TruSeq Kit from dust mite DNA recovered by organic extraction followed by precipitation. Using the Bismark algorithm [72] base converted dust mite genome indexes were used to determine the rate of cytosine methylation. Using coordinates from cufflinks (mRNA), and RepeatMasker (TE) annotations, rates of methylation were determined for different genomic features. Reads were mapped uniquely and duplicated reads were discarded that resulting in an average 6X coverage depth [72]. Using bedtools, genomic regions that had >4 reads mapping were determined and the base conversion rate measured.  (18, 19, etc.) and together (18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29)(30). R heatmap2 package was used to draw the heatmap. 2nt dicer processing register is shown by blue arrow "D". Red arrow labeled "pp" shows 10nt ping-pong overlap signature. Blank areas in the heatmap are due to the absence of overlapping pairs. (TIF)