Small RNA Profiling in Dengue Virus 2-Infected Aedes Mosquito Cells Reveals Viral piRNAs and Novel Host miRNAs

In Aedes mosquitoes, infections with arthropod-borne viruses (arboviruses) trigger or modulate the expression of various classes of viral and host-derived small RNAs, including small interfering RNAs (siRNAs), PIWI interacting RNAs (piRNAs), and microRNAs (miRNAs). Viral siRNAs are at the core of the antiviral RNA interference machinery, one of the key pathways that limit virus replication in invertebrates. Besides siRNAs, Aedes mosquitoes and cells derived from these insects produce arbovirus-derived piRNAs, the best studied examples being viruses from the Togaviridae or Bunyaviridae families. Host miRNAs modulate the expression of a large number of genes and their levels may change in response to viral infections. In addition, some viruses, mostly with a DNA genome, express their own miRNAs to regulate host and viral gene expression. Here, we perform a comprehensive analysis of both viral and host-derived small RNAs in Aedes aegypti Aag2 cells infected with dengue virus 2 (DENV), a member of the Flaviviridae family. Aag2 cells are competent in producing all three types of small RNAs and provide a powerful tool to explore the crosstalk between arboviral infection and the distinct RNA silencing pathways. Interestingly, besides the well-characterized DENV-derived siRNAs, a specific population of viral piRNAs was identified in infected Aag2 cells. Knockdown of Piwi5, Ago3 and, to a lesser extent, Piwi6 results in reduction of vpiRNA levels, providing the first genetic evidence that Aedes PIWI proteins produce DENV-derived small RNAs. In contrast, we do not find convincing evidence for the production of virus-derived miRNAs. Neither do we find that host miRNA expression is strongly changed upon DENV2 infection. Finally, our deep-sequencing analyses detect 30 novel Aedes miRNAs, complementing the repertoire of regulatory small RNAs in this important vector species.


Introduction
Aedes mosquitoes are essential vectors for the transmission of important arthropod-borne viruses (arboviruses), including dengue virus (DENV), yellow fever virus, and chikungunya virus [1]. While several of these arboviral infections cause disease in humans, virus replication generally does not lead to severe pathology in vector mosquitoes. Infected mosquitoes thus serve as a persistent reservoir for arboviruses in the wild and they may transmit these viruses to vertebrate hosts throughout their entire lives [2].
After ingestion in a mosquito's blood meal, arboviruses need to overcome a number of anatomical and immunological barriers to reach sufficiently high titres in the saliva. Only then can transmission to a naive vertebrate host efficiently occur. One of the most important immune responses to arboviral infection is antiviral RNA interference (RNAi) [3][4][5]. This pathway is triggered by the presence of double stranded RNA (dsRNA), which is produced during the replication of RNA and DNA viruses [6,7]. The dsRNA is recognized and cleaved by the RNase-III enzyme Dicer-2 (Dcr2) into 21 nucleotide (nt) small interfering RNA duplexes (viral siRNA; vsiRNA) [8,9]. One of the siRNA strands is incorporated in Argonaute-2 (Ago2), the core protein of the RNA induced silencing complex (RISC) [10]. The siRNA-loaded RISC complex is guided to complementary viral RNA molecules and cleaves these target RNAs using the endonuclease (slicer) activity of Ago2 [11].
MicroRNAs (miRNAs) are a distinct class of small RNAs that are produced from genomeencoded stem loop-containing transcripts known as primary miRNA (pri-miRNAs). During the canonical miRNA biogenesis pathway, the stem loop structures, known as precursor miRNA (pre-miRNA), are released from the pri-miRNA by the microprocessor complex with at its core the RNase-III enzyme Drosha. After translocation into the cytoplasm, pre-miRNAs are cleaved by Dicer-1 (Dcr1) to produce a small RNA duplex comprised of the two mature miRNA strands. Usually, one of these strands is then preferentially incorporated into the Argonaute-1 (Ago1) containing miRNA-induced silencing complex (miRISC), whereas the other strand (the passenger or miRNA Ã strand) is usually discarded [12]. Loaded miRISC complexes are able to bind to specific target sites within mRNAs. This miRNA-mRNA interaction is initiated by nucleotide two to seven of the miRNA, the so-called seed sequence [13]. Stable binding of miRISC to an mRNA target, generally causes down-regulation of gene expression via translational inhibition and mRNA destabilization [14]. Importantly, infecting viruses can directly or indirectly, as a consequence of the immune response, reshape the host miRNA expression landscape. While quite a number of studies have reported on this matter in mammalian systems [15], little is known about virus-induced changes in miRNA levels in mosquito vectors. In Aedes mosquitoes, miRNA levels or modifications have been reported to be changed upon infections with DENV, West Nile virus, and chikungunya virus [16][17][18][19][20]. For most of these differentially expressed miRNAs, the biological relevance as well as the targeted mRNAs still await experimental validation.
Besides modulation of host miRNAs, some DNA and retroviruses encode their own miR-NAs to regulate viral and host mRNAs [21]. The expression of miRNAs from cytoplasmic RNA viruses has been controversial. However, functional introduction of artificial miRNAs into the genomes of Sindbis virus (SINV) and tick-borne encephalitis virus provides evidence that miRNA production from cytoplasmic RNA viruses may in principle be possible [22,23]. Yet, the presence and biological relevance of miRNAs encoded in the genomes of flaviviruses such as DENV is still an issue of debate [24][25][26].
The third, most enigmatic class of small RNAs are PIWI interacting RNAs (piRNAs). These are processed from long RNA precursors that are transcribed from genomic loci known as piRNA clusters. In sharp contrast to siRNAs and miRNAs, their biogenesis into mature piR-NAs is Dicer-independent. In Drosophila, piRNA maturation involves endonucleolytic cleavage of precursor transcripts by the Zucchini nuclease and the three PIWI proteins Piwi, Aubergine (Aub) and Argonaute-3 (Ago3) [27][28][29][30]. The primary function of the piRNA pathway in this model organism is the defence against transposable elements, mainly in germ-line tissues. Interestingly, piRNAs of viral origin (vpiRNA) have been found in somatic tissue of Aedes mosquitoes, suggesting that they contribute to the regulation of virus replication [31]. At present, vpiRNAs have been discovered upon infection with a number of Alphaviruses, Bunyaviruses and Flaviviruses, including DENV [31][32][33][34][35][36][37][38][39]. However, with the exception of SINV (Alphavirus), their molecular biogenesis has not been investigated [37].
Here, we make use of small RNA deep-sequencing in the siRNA, miRNA, and piRNA competent Aedes aegypti Aag2 cell line to investigate the production of small RNAs during DENV infection. We find that in addition to the well-characterized vsiRNAs, specific vpiRNAs are produced from DENV, which for their biogenesis in Aag2 cells rely on Piwi5 and Ago3 and, to a lesser extent, on Piwi6. We do not detect DENV-derived miRNAs, or prominent changes in host miRNA levels upon infection. Finally, we identify novel host miRNAs in our small RNA deep-sequencing libraries, complementing the currently annotated miRNA repertoire in Aedes aegypti vector mosquitoes.

Infection of Aag2 cells with DENV2
Aag2 cells were seeded one day prior to infection and infected with DENV2 at a multiplicity of infection (MOI) of 0.5 by directly adding the virus to the culture medium. Three days post infection the culture medium was removed and cells were harvested for RNA and protein isolation as detailed below.

Western blot
For the detection of the DENV NS1 protein in samples used for small RNA deep-sequencing, 5% of the cells were harvested in 50 μl lysis buffer (50 mM Tris-HCl pH 7.8; 150 mM NaCl; 1 mM EDTA; 0.5% NP-40; 1x Protease inhibitor cocktail (Roche); 1 mM DTT). 12.5 μl of 5x Laemmli buffer was added to each sample, incubated at 95°C for 5 min, and 30μl of each sample was loaded on a 12.5% polyacrylamide gel. After gel electrophoresis, proteins were transferred to a nitrocellulose membrane (Bio-Rad) using a semi-dry blotting system (Bio-Rad). The membrane was blocked in 5% non-fat dry milk (Bio-Rad) in 0.1% Tween20 in PBS (PBS-T) for 30 min at room temperature. Mouse anti DENV NS1 antibody was kindly provided by Dr. Peter Mason [41]. The antibody was added to the membrane in a 1:1,000 dilution in 5% blocking buffer. After an incubation for 1.5 hours at room temperature, the membrane was washed three times in PBS-T. IRdye680 conjugated goat anti mouse antibody (1:15,000 dilution in PBS-T; Licor) was then added to the membrane and incubated at room temperature for 1.5 hours. After three washing steps, the membrane was imaged on an Odyssey infrared image system (Licor).
dsRNA production and transfection of Aag2 cells dsRNAs targeting PIWI/AGO transcripts or luciferase as a negative control were produced by in vitro transcription from T7-promoter flanked PCR products as detailed in [37]. Primers to produce T7-flanked PCR products are indicated in S1 Table. For dsRNA transfection, 7.5x10 5 Aag2 cells were seeded in one well of a 24-well plate. For each condition, three wells were plated. The following day, transfection mixes containing 300 μl non-supplemented L-15 medium, 450 ng dsRNA and 1.8 μl X-tremeGENE HP (Roche) were prepared according to the manufacturer's recommendations. 100 μl of the mix was added dropwise to one well. After 2-3 hours the medium was replaced with fully supplemented L-15 medium. 48 hours later, the transfection was repeated to enhance knockdown efficiencies. Where indicated, the cells were infected with DENV2 which was added to the L-15 medium used to replace the transfection medium as specified above.

RNA isolation
Aag2 cells were lysed in Isol-RNA Lysis reagent (5 PRIME) as described in the manufacturer's instructions. Briefly, 200 μl of chloroform was added to 1 ml of Lysis reagent and mixed well. After centrifugation, the aqueous phase was collected and total RNA was purified using isopropanol precipitation. RNA was quantified on a Nanodrop spectrophotometer and RNA integrity was checked by ethidium bromide staining of ribosomal RNA bands after agarose gel electrophoresis.

DNaseI treatment, reverse transcription and (quantitative) PCR
For RT-(q)PCR, 1 μg of total RNA was DNaseI (Ambion) treated according to the manufacturer's instructions. The RNA was subsequently reverse transcribed in a 20 μl reaction using the Taqman reverse transcription kit (Applied Biosystems). Complementary DNA (cDNA) was diluted 5-10 times before PCR amplification. Endpoint PCR was performed using Thermoperfect DNA Polymerase. Quantitative PCR (qPCR) analysis was performed using the GoTaq qPCR SYBR mastermix (Promega) on a LightCycler 480 instrument (Roche). The relative changes in gene expression were calculated using the ΔΔCt method [42] using lysosomal aspartic protease (LAP) as an internal normalization control. Sequences of the PCR primers are indicated in S1 Table. β-elimination Sodium periodate (NaIO 4 ) oxidation and β-elimination of total RNA was performed as described previously [43]. Total RNA (10 μg in 47.5 μl nuclease-free water) was mixed with 12.5 μl 200 mM NaIO 4 and 40 μl 5x borate buffer. As a control, RNA was treated with water instead of NaIO 4 . The reaction was incubated at room temperature for 30 min and 10 μl glycerol was added to the reaction. The reaction was incubated for another 10 min before 10 μl of 500 mM sodium hydroxide (NaOH) was added to induce β-elimination. The reaction was incubated at 45°C for 90 min. After these treatments, total RNA was purified by ethanol precipitation in the presence of 300 mM NaCl and 5 μg of glycogen. Electrophoretic mobility of Aedes aegypti miR-2940-3p and DENV2 piRNAs was then analyzed by small RNA northern blotting as detailed below.

Small RNA northern blotting
Small RNA northern blot was performed as described in [44]. Briefly, total RNA was size-separated on 0.5x TBE, 7 M Urea, 15% Polyacrylamide gels, transferred to Hybond NX nylon membranes (Amersham), and cross-linked using 1-ethyl-3-(3-dimethylaminopropyl) carbodiimide (EDC; Sigma). Individual small RNAs were detected with DNA oligonucleotides that were 5' end-labelled with [ 32 P] γ-adenosine-triphosphate (Perking Elmer) using T4 Polynucleotide kinase (Roche). Hybridization to the oligo-probes was performed overnight at 42°C in Ultrahyb Oligo hybridization buffer (Ambion). Membranes were then washed three times at 42°C in 0.1% SDS with decreasing concentrations of SSC (2x, 1x, 0.1x). Membranes were exposed to Xray films (Carestream) or Phosphorimager screens (BioRad). Sequences of DNA oligonucleotide probes are indicated in S1 Table. Quantification of northern blot panels was performed using ImageJ software. Bands were defined using the rectangular selection tool and the pixel density (area under the curve) was measured and normalized to uninfected dsLuc samples.

Preparation of small RNA libraries
Small RNA libraries were prepared as described previously [37,45]. Briefly, three 25 cm 2 flasks of Aag2 cells were infected in parallel with DENV2 NGC. Three additional flasks were left uninfected. Total RNA was then isolated from these six flasks as specified above and 30 μg of total RNA was size separated on a 15% Polyacrylamide, 7 M urea, 0.5x TBE gel. Subsequently, the small RNAs in the size range from 18 nt to 33 nt were excised from gel using radioactivelylabelled RNA oligos, loaded in the adjacent lanes of the gel, as rulers. The gel was crushed and the small RNAs were eluted in 300 mM sodium acetate overnight at 4°C under constant rotation. The RNA was recovered from the elution buffer using ethanol precipitation and eluted in 10 μl of nuclease-free water. 5 μl of the sample was directly used as input for small RNA deepsequencing library preparation using the TruSeq small RNA library preparation kit (Illumina) following the manufacturer's recommendations. The RNA was ligated to 3' and 5' adapters, reverse transcribed, and PCR amplified. The small RNA libraries were then size-purified from 1x TBE, 6% polyacrylamide gel using overnight elution in 300 mM sodium acetate followed by ethanol precipitation. The individual small RNA libraries were pooled and sequenced on a single sequencing lane on a HighSeq2500 by Baseclear (Leiden, The Netherlands).

Viral small RNA profiling
FASTQ sequence reads were generated using the Casava pipeline (v.1.8.3) and initial quality analysis was performed using the Illumina Chastity filter and an in-house filtering protocol by Baseclear. Subsequent quality assessment was based on the FASTQC quality control tool (v.1.10.0). The individual small RNA sequencing libraries were separated based on the TruSeq indices (no. 1 to 6) that were introduced during PCR amplification. The individual libraries were subsequently analyzed using the Galaxy bioinformatics tool shed [46,47]. For the analysis of viral small RNAs, reads were mapped to the DENV2 NGC genome (GenBank accession: KM204118.1) using Bowtie (v.1.1.2) [48]. Size profiles were obtained from all reads that align to this reference sequence with a maximum of one mismatch. The genome distribution of 21 nt siRNAs, 22-24 nt small RNAs, or 25-30 nt piRNAs was obtained by plotting the number of 5' ends of these reads at each position of the genome. For the pileup plots of UTR-derived miRNA-like small RNAs, the 22-24 nt small RNAs were selected from the initial FASTQ files and mapped to the DENV-NGC genome. From the resulting SAM files the reads mapping to the (+) strand of the virus genome were selected and used as input for the 'Generate pileup from BAM dataset' tool (v. 1.1.2). The values at the nucleotide positions of the DENV2 5'UTR (1-96) and 3' UTR (10273-10723) were selected for display.
Re-analysis of the data published by Hess et al. [36] was performed on the dataset with the accession number SRR921363. The SOLiD-formatted dataset was groomed to fit the requirements for manipulation in Galaxy. After adapter clipping, reads were mapped to the DENV2--JAM1409 genome (GenBank accession: M20558) using Bowtie2 [49]. Size and genome profiles were obtained as described above. Unless specified differently, all small RNA read counts were normalized against the size of the corresponding sequencing library and are expressed as '% of the library' (i.e. reads per hundred).

miRNA analysis and prediction
Analysis of miRNA expression levels was performed using the miRDeep2 tool. Raw data were assessed for quality using FASTQC. Subsequently, adapters were removed from the raw reads, and the reads were quality trimmed using cutadapt software (http://dx.doi.org/10.14806/ej.17.1. 200) with parameters -O 6 -m 17 -n 5 -q 20. Within each library, the resulting reads were collapsed to generate a non-redundant set of FASTA sequences, subsequently processed to the format required for miRNA prediction with miRDeep2 software [50]. Collapsed reads longer than 17 nucleotides were aligned to the Aedes genome (assembly AaegL3, downloaded from vectorbase) and the DENV2 NGC genome using the mapper.pl component of miRDeep2 (parameters: -o 20 -l 19 -r 100 -c). The resulting outputs were parsed to remove alignments that were not full length and perfect match (FLPM). miRDeep2 predictions were generated from the FLPM aligned sequences, with miRBase v21 Arthropoda mature miRNA and pre-miRNA sequences as templates [51]. The 'miRNAs_expressed' output from miRDeep2, which comprises tallies for each known miRNA in each sample, was further processed in the R/Bioconductor environment. Briefly, miRNA read counts were normalized to the number of Aedes-specific genome reads within each sample group, using the lowest number of reads aligning as the baseline. Subsequently, the counts were converted to abundances within each sample, converted to log2 equivalent counts, and all samples quantile normalized prior to linear model fitting with the limma package [52]. MiRNA predictions from the miRDeep2 output were manually curated using the following criteria: i) high-confidence miRNA predictions have reads mapping to both a predicted mature and star sequence, ii) the mature sequences have a homogenous 5' end (80% of reads start at same position), and iii) miRNA-miRNA Ã duplex should resemble a Dicer product on a genomically-encoded hairpin, having a two nt (+/-1 nt) overhang at the 3' end. miRNA predictions supported by >1000 reads as well as miRNA predictions with a seed match to known insect miRNAs were also kept. To be retained, these predictions required a homogenous 5' end, but did not require the presence of reads mapping to the expected star strand. miRNA names and accession numbers were assigned by the miRBase repository.
All deep sequencing libraries have been submitted to NCBI Sequence Read Archive under the accession number SRA303329. All source data are available in S1 Dataset.

DENV2-derived small RNAs in infected Aag2 cells
DENV is a positive (+) strand RNA virus belonging to the Flavivirus genus in the Flaviviridae family. Its RNA genome is approximately 10.7 kilobases in size and encodes a single polypeptide that is processed by proteolytic cleavage events into three structural proteins and seven non-structural proteins (Fig 1A). Since various classes of small RNAs have been implicated in modulating DENV infections in its mosquito vectors, we aimed to characterize the repertoire of virus and host-derived small RNAs in Aedes aegypti Aag2 cells. To this end, we prepared three independent small RNA deep-sequencing libraries from uninfected and DENV2 (NGC strain) infected cells, each. The efficiency of the three infections was comparable as assessed by western blot for the DENV2 NS1 protein ( Fig 1B). As expected, DENV2-derived viral small RNAs showed a clear population of 21 nt vsiRNAs mapping to both the viral positive (+) strand and the negative (-) strand in roughly equal numbers ( Fig 1C). Interestingly, besides siRNAs, a second population of viral small RNAs was produced that resembled vpiRNAs. These were 25-30 nt in length and almost exclusively derived from the viral (+) strand ( Fig  1C). In contrast to vsiRNA, which were distributed along the entire length of the viral genome, these putative vpiRNAs were produced only from few specific positions (Fig 1D). In fact, 85% of all the 25-30 nt reads were derived from four individual vpiRNA sequences, present in the NS5 gene at positions 9180 and 9985, 9989 and 9990 of the DENV2 NGC genome. To test whether these small RNA profiles reflect those from adult mosquitoes, we re-analyzed deep sequencing data from DENV2 infected Aedes aegypti mosquitoes published by Hess et al. [36]. We analyzed the 9 days post infection sample, which showed the highest number of viral siR-NAs and piRNAs. Whereas normalized vsiRNA levels were only 2.2-fold lower in these libraries than in our Aag2 data, vpiRNAs were about forty times lower. Yet, the viral small RNA profiles were strikingly similar, with 21 nt reads being scattered throughout the entire viral genome and piRNA-sized reads being predominantly produced from few positions located towards the 3' end of the viral genome (S1 Fig). These data suggest that similar mechanisms might produce viral piRNAs in Aag2 cells and adult mosquitoes.
vpiRNA production from DENV2 RNA To exclude the possibility that the piRNA-like molecules are sequencing artefacts and to characterize this small RNA population in more detail, we performed small RNA northern blotting for the highly-abundant small RNAs starting at DENV2 genome positions 9180 or 9985-9990. Indeed, small RNAs in the expected size range could readily be detected specifically in DENV2 infected Aag2 cells (Fig 2A). In addition, these sequences were also present in DENV2-infected U4.4 and C6/36 cells derived from Aedes albopictus mosquitoes, which we have previously shown to be competent in producing SINV-derived vpiRNAs [32,37] (Fig 2B). The levels of vpiRNAs did not correlate with the expression of viral genomic RNA. Whereas viral RNA levels were roughly eighteen to nineteen fold higher in U4.4 and C6/36 cells than in Aag2 cells, vpiRNAs were most abundant in Aag2 cells ( Fig 2B). These data suggest that the composition of host factors required for their biogenesis is most favourable in Aag2 cells. As expected, mammalian BHK-21 cells, which lack an active piRNA pathway, did not produce vpiRNAs ( Fig 2B). To exclude that viral piRNA production is an artefact of the use of the specific DENV2 strain, we analyzed piRNA accumulation in Aag2 cells infected with either the DENV NGC or DENV 16681 strain. Small RNA northern blotting revealed that infection with either strain resulted in the production of those viral piRNA sequences that we had found by deep-sequencing ( Fig  2C). Next, we aimed to test whether DENV2-derived piRNAs are methylated at their 3' end. To this end, we performed sodium periodate oxidation followed by beta-elimination, which reveals modifications of the 3' terminal nucleotide of RNA molecules [43]. Unmodified small RNAs, such as animal miRNAs, are susceptible to this treatment and will be shortened by one nucleoside resulting in increased electrophoretic mobility. In contrast to miRNAs, piRNAs are protected against this treatment by methylation of the 2'OH on the ribose of the 3' terminal nucleotide. Indeed, beta-elimination resulted in increased electric mobility of miR-2940-3p. Yet, the migration of DENV2 piRNA bands was not affected by beta-elimination, indicating that their 3' terminal nucleotides are modified, most likely methylated (Fig 2D). Since piRNA methylation occurs after loading into PIWI protein complexes, these data suggest that the identified DENV2 piRNAs are mature piRNAs associated with a PIWI protein.
To identify which PIWI proteins are required for the biogenesis of DENV2 piRNAs in Aag2 cells, we individually knocked down expression of all the eight Aedes PIWI proteins and analyzed the production of vpiRNAs by northern blot. We confirmed knockdown efficiency of roughly 90% for the four PIWI proteins that are detectable by RT-qPCR in Aag2 cells (i.e. Piwi4, Piwi5, Piwi6, Ago3; Fig 2E). Expression levels of Piwi1-3 and Piwi7 were too low to allow reliable quantification. DENV2 piRNAs were almost undetectable upon knockdown of Piwi5 and Ago3 and clearly reduced upon knockdown of Piwi6 ( Fig 2F). These data confirm that the 25-30 nt population of DENV2-derived small RNAs are bona fide piRNAs that require host PIWI proteins for their biogenesis. To test whether knockdown of PIWI expression results in enhanced DENV2 replication, we performed RT-qPCR to compare viral RNA levels in the different knockdown conditions. We found that none of the knockdowns resulted in a significant change in viral RNA levels (S2A Fig). Yet, also knockdown of the well-established antiviral factor Ago2 [53] only resulted in a minor, statistically not-significant, increase of viral RNA replication although knockdown efficiency was higher than 90% (S2B Fig). This suggests that, in our hands, knockdown of small silencing pathway components in Aag2 cells is not suited to uncover robust antiviral activity against DENV2.

DENV2 miRNA-like small RNAs are not expressed in Aag2 cells
The significance of viral miRNA production from DENV genomic RNA is heavily debated [24][25][26]54]. To investigate whether viral miRNA-like molecules are produced in DENV2 infected Aag2 cells, we filtered 22-24 nt small RNA reads that map to the DENV2 genome with a maximum of one mismatch. In general, the number of 22 to 24 nt reads was rather low (~5% of all DENV2 mapping reads) when compared to 21 nt siRNAs (~28%) and 25-30nt piR-NAs (~40%). Furthermore, there were only four outstanding peaks that gave rise to a somewhat higher number of small RNAs (Fig 3A). All of them coincided with the position of a vpiRNA peaks (Fig 1D), suggesting that these small RNAs were by-products of vpiRNA production. Parallel analysis of the virus-derived reads using miRDeep2 did not identify convincing miRNA-like candidates: some reads mapped to two predicted hairpin sequences (genome positions: 9542 (+) strand; 4888 (-) strand), however the mapping patterns showed heterogeneity of the 5' start sites and did not suggest Dicer processing (S3A Fig). Recently, eight miRNA-like small RNAs were computationally predicted based on hairpin structures in the DENV2 genome, but they were not experimentally validated [54]. We specifically looked for small RNA reads in our sequencing data mapping in the proximity of these predicted viral miRNAs, allowing a margin of 3nt around the start site. For each of the predicted miRNAs, we identified only very few (<20) reads in the combined set of DENV2-infected small RNA libraries (total of >3.7Á10 7 reads of which >3.6Á10 5 are DENV specific). In another publication, several 'miRNA-like' RNAs (termed vsRNA-1 to 6) were proposed to be produced from the 5' and 3' UTRs of the DENV2 RNA, based on the analysis of small RNA sequencing data [24]. The DENV UTRs are indeed prone to form RNA structures and hairpins, which were suggested to be processed by the miRNA machinery into specific small RNA species [24]. We specifically looked for the proposed vsRNA sequences in our dataset but could only identify a small RNA population that resembled vsRNA-2 located at the terminal hairpin Small RNAs in DENV2-Infected Mosquito Cells of the DENV2 genome ( Fig 3B). However, small RNAs mapping in that region showed a broad size distribution with the majority ranging in size from 26 to 28 nt, arguing against vsRNA-2 small RNAs being bona fide Dicer products (S3B and S3C Fig). These data suggest that the proposed DENV2 vsRNAs are not an abundant class of small RNAs.

Host miRNA levels are only mildly affected by DENV2 infection
Host miRNAs function as key regulators of gene expression and changes in miRNA expression have been reported during virus infections in various animal hosts, including mosquitoes. To assess host miRNA expression in response to DENV2 infection, we made use of the miRDeep2 toolkit to quantify miRNAs in the uninfected and DENV2-infected Aag2 small RNA libraries. DENV2 infection caused only minimal changes in miRNA levels (Fig 4). The expression of three and seven miRNAs was changed more than 2-fold up or down, respectively in response to DENV infection. Yet, the majority of differentially regulated miRNAs, including all up-regulated miRNAs, were poorly expressed (mean expression levels below twenty reads) making it hard to discriminate these expression changes from experimental noise due to low read counts. Collectively, these data suggest that miRNA expression in Aag2 cells is not heavily affected by DENV infection.

Novel Aedes aegypti miRNAs
The most recent version of miRBase (version 21) contains 101 Aedes aegypti miRNAs, which is considerably less than for other insect species including Drosophila melanogaster (fruitfly; 256 miRNAs), Apis mellifera (honey bee; 254 miRNAs) or Bombyx mori (silkworm; 487 miRNAs). We therefore suspected that the repertoire of published mosquito miRNAs is not yet complete and we set out to identify novel Aedes aegypti miRNAs in our small RNA sequencing data using miRDeep2. We obtained a list of 399 miRNA predictions (S2 Table), 73 of which were known miRNAs annotated in miRBase [51,55,56]. We also confirmed 16 miRNAs that were The remaining 310 predictions of miRNA hairpins were manually inspected for novel miRNAs using a similar approach as described in [59,60], based on three criteria. First, only hairpin predictions that were supported by at least 1000 mature miRNA reads or those with at least one predicted miRNA Ã strand were retained. If a miRNA had an identical seed to a known insect miRNA, it was also retained irrespectively of read count or the presence of a miRNA Ã . In total 140 predictions met these criteria. Second, mature miRNAs were inspected for a homogeneous 5' end of the supporting small RNA reads, defined as having at least 80% of the miRNA reads starting at the same nucleotide. 68 miRNA hairpin predictions met this requirement. Third, the predicted miRNA duplex was checked to resemble a (near) perfect Dicer product, defined as mapping to the stem of a stem-loop structure with one, two or three nucleotide overhangs at the 3' end. Using this approach, we identified 31 unique mature miR-NAs sequences mapping to 39 predicted miRNA hairpins. Nine of the 31 mature miRNA predictions did not have reads mapping to the star strand, but were supported by a homogeneous 5' end in combination with a seed-match to an insect miRNA or having >1000 reads (Table 1, Fig 5). Further inspection identified one of these predictions to be derived from a tRNA which has therefore been removed from the list of predicted miRNAs.

miRNA offset RNAs in Aag2 cells
During the analysis of miRNA predictions, we noted the expression of specific small RNAs adjacent to the mature miRNA and miRNA Ã star strands. These miRNA offset RNAs (moRs) have been detected in small RNA deep-sequencing data from invertebrates, simple chordates, vertebrates and even viruses [61][62][63][64][65][66][67][68]. In total, we identified moRs for 27% (24/89) of previously reported miRNA hairpins. In many of these cases (nine out of 24), the number of moRs per hairpin is below one thousandth of the number of mature miRNA reads. Others were more abundant, with moR reads accumulating up to 3.6% of the number of mature miRNAs. In a single instance, miR-11894b-1, the number of moRs reached 10.8% of the number of mature miRNA reads. In agreement with previous findings in Drosophila [63] 5' moRs were more abundant than 3' moRs ( Fig 6A). moRs have been proposed to be the by-products of Drosha cleavage and in line with this suggestion, we found the 3' end of 5' moRs and the 5' end of 3' moRs to be fixed, reflecting potential Drosha cleavage sites (Fig 6B and 6C). In contrast, the ends of moRs facing the termini of the miRNA stemloop were less well defined, suggesting that they are processed by exonuclease activity (Fig 6B). Why certain miRNA hairpins are prone to accumulation of moRs remains unclear.

Discussion
Small RNA pathways critically influence the outcome of virus infections in many host organisms, including plants, fungi, invertebrates, and vertebrates [3][4][5]15,69]. In plants and invertebrates, siRNA-mediated antiviral immunity is key to the defence against a broad range of virus infections. In Aedes mosquitoes, viral siRNAs were detected from several virus families, including Togaviridae, Flaviviridae, Bunyaviridae, and Reoviridae [5]. In line with previous reports, our analysis identified viral siRNAs derived from the entire genomic RNA of DENV2 in Aedes aegypti cells [35,36]. These vsiRNAs are produced in roughly equal amounts from the (+) strand and the (-) strand of the virus, indicating that the dsRNA replication intermediates serve as substrate for Dcr2. Upon knockdown of either Dcr2 or Ago2 in whole Aedes mosquitoes, DENV2 titres and transmission are enhanced, underlining the pivotal role of RNA interference in limiting DENV2 replication [53]. miRNA predictions above the first thick line have a seed sequence that is present in an insect species, miRNAs between the two thick lines have a seed sequence that is present in a metazoan species, and miRNAs below the second thick line have seed sequences that are not present in any metazoan miRNA. * The mean read count is normalized to the size of the corresponding small RNA sequencing library and presented as reads per million. † miRNA prediction is supported by a seed match to known insect miRNAs or >1000 reads (equivalent to appr. 80 rpm), but not by the presence of a star strand. doi:10.1371/journal.pntd.0004452.t001 Small RNAs in DENV2-Infected Mosquito Cells Novel Aedes aegypti pre-miRNAs. Hairpin structures of novel pre-miRNAs, as predicted using the RNAfold algorithm. Red letters indicate the position of the predicted mature miRNA sequence. Blue letters indicate the miRNA* strand. No miR* strand was found for ten predictions, representing nine eight miRNAs; these predictions are supported either by a seed-sequence known in insects or by high read counts. The mature miRNA sequences of aae-miR11908, aae-miR-11909 and aae-miR-11926 map to multiple different hairpins in the Aedes genome. For aae-miR-11914 and aae-miR-11920, the entire hairpin is encoded at multiple locations in the Aedes genome, as specified in Table 1. doi:10.1371/journal.pntd.0004452.g005 Besides siRNAs, our DENV2 infected small RNA libraries contained a substantial number of virus-mapping reads in the size range of piRNAs. Their expression was confirmed using small RNA northern blotting, validating that these small RNA reads were no sequencing artefacts. Only very few DENV2 genomic locations near the 3' end of the DENV genome give rise to these vpiRNAs but the origin of this spiky pattern remains obscure. We hypothesized that perhaps endogenous, transposon-derived piRNAs would loosely bind the DENV2 genome at these positions triggering the production of secondary vpiRNAs. A similar mechanism has been suggested to initiate piRNA production from specific mRNAs in Drosophila [28]. However, various mapping strategies allowing small RNA alignment with up to six mismatches did not uncover endogenous piRNAs that could trigger vpiRNA production at the observed positions. Previous analyses of DENV2-derived small RNAs identified vpiRNAs in Aedes aegypti mosquitoes or Aag2 cells [35,36]. Also in the viral small RNA population reported by Scott and colleagues [35], a major small RNA spike is located near the 3' end of the DENV genome. Since the entire population of viral RNA is analyzed in this study, it is however hard to assess which type of small RNA contributes to the spike. Interestingly, the spiky genome distribution of vpiRNAs is also recapitulated in adult Aedes mosquitoes. Although in this study the exact location of vpiRNA spikes differs from the positions we found in Aag2 cells, these data suggest that similar mechanisms may be responsible for piRNA biogenesis in Aag2 cells and adult mosquitoes. Yet, genetic evidence for the PIWI protein dependency of vpiRNA production in vivo is lacking. Of note, piRNAs are far less abundant in the small RNA libraries reported by Scott et al. [35] and Hess et al. [36] compared to our data from Aag2 cells. These differences might be due to different experimental conditions, including the chosen MOI, the time point of sampling, or to differences in small RNA library preparation and sequencing methodology. In addition, the specific viral strains may critically influence the accumulation of vpiRNAs. We have tested two laboratory-adapted DENV2 strains which both give rise to vpiRNAs, and it would be interesting to test if pathogenic strains from DENV2 endemic areas would show similar phenotypes.
Using knockdown of PIWI proteins, we identified Piwi5, Ago3 and, to a lesser extent, Piwi6 as responsible for the production of vpiRNAs in Aag2 cells. Therefore, DENV2 vpiRNA biogenesis in Aag2 cells relies on a similar set of PIWI proteins as SINV vpiRNAs, which also depend on Piwi5 and Ago3 [37]. We have recently proposed that Aedes PIWI proteins are specialized in producing piRNAs from various sources. Whereas in mosquito cells piRNA biogenesis from transposons directly and indirectly depends on Piwi4-6 and Ago3, piRNA biogenesis from SINV predominantly requires Piwi5 and Ago3 only [37]. The additional involvement of Piwi6 for piRNA production from DENV2 suggests that Aedes PIWI proteins are even further specialized towards RNA substrates from different viruses. This may be caused by virus-specific sequence elements or structures that are preferentially recognized by certain PIWI proteins. Alternatively, but not mutually exclusive, differences in replication strategies or replication sites might favour recognition of viral RNA by distinct sets of PIWI proteins.
The almost complete loss of DENV2 piRNAs upon knockdown of Piwi5 and Ago3 indicates that both proteins are equally important for vpiRNA biogenesis in Aag2 cells, similar to piRNA biogenesis during SINV infection [37]. SINV piRNAs are produced by a two-step amplification mechanism that resembles ping-pong amplification of transposon piRNAs in Drosophila [27,29]. During this process, a piRNA-loaded PIWI protein (Piwi5 in Aedes or Aubergine in Drosophila) slices a complementary target RNA and transfers the 3' slicer products as the new piRNA precursor to a second PIWI protein (Ago3 in Aedes and Drosophila). From this precursor, an Ago3-bound secondary piRNA is produced that in turn is able to slice a target RNA, giving rise to a new piRNA precursor. This precursor will be matured to generate the same primary piRNA sequence that initiated the amplification. Therefore, this model predicts the presence of piRNAs derived from both strands. Although also during SINV infection (-) strand derived piRNAs are only a minor fraction, they can be identified as the primary piRNAs by a nucleotide bias that is characteristic for Piwi5/Aub bound piRNAs (uridine at position one). During DENV2 infection, 25-30 nt reads from the (-) strand are extremely scarce and they do not have the nucleotide bias that would classify them as primary piRNAs. Therefore, exactly how the production of the secondary, (+) strand vpiRNAs is triggered or whether a different, amplification-independent mechanism is responsible for their production, remains unclear.
The expression of miRNAs from DENV2 genomic RNA is still debated. Based on small RNA sequencing data, Hussain and Asgari have described a set of six viral small RNAs that have miRNA-like properties [24]. Inhibition of one of them, vsRNA-5, by complementary RNA molecules strongly enhances DENV virus replication [24]. However, expression of these viral small RNAs is generally not high and the relevance of such a lowly abundant small RNA during the exponential growth of a virus was therefore questioned [25]. In our dataset, we find a high number of specific vsRNA reads for only vsRNA-2. For all the other predicted vsRNAs we find no or very low numbers of reads. vsRNA-2 is located on a hairpin at the very end of the DENV2 genome. This strongly resembles KUN-miRNA1, a viral small RNA expressed in mosquito cells infected with West Nile virus, a related flavivirus [70]. KUN-miR1 is 21nt in size and its expression is Dcr1-dependent. In contrast, our data demonstrate that vsRNA-2 has a broad size distribution of primarily 26-28 nt, arguing against it being a canonical Dicer-dependent miRNA. Altogether, these data support the notion that miRNA-like small RNAs from DENV2 are extremely lowly abundant, with a questionable role in the regulation of viral replication.
Modulation of host miRNAs after virus infection may be a mechanism that coordinates gene expression during the course of the immune response. Alternatively, it may be a consequence of a viral strategy to manipulate host gene expression. Comparing uninfected with DENV2-infected Aag2 cells showed that the expression of almost all miRNAs was unchanged upon infection. Only a handful of miRNAs were up or down-regulated after exposure to DENV2. The fold changes ranged from approximately 4 fold up to 4 fold down. In whole Aedes aegypti mosquitoes, a total of 31 miRNAs were recently shown to be differentially regulated following infection at three different time points [17]. The set of differentially expressed miRNAs is inconsistent between the different analyzed time points (2, 4 and 9 days post infection, dpi), but the number of differentially expressed miRNAs was higher at nine dpi than at two or four dpi. These data suggest that changes in miRNA expression may be more prominent in a long-term infection and we thus cannot exclude the possibility that prolonged infection of Aag2 cells may result in more pronounced changes in miRNA expression. Alternatively, the observed miRNA changes in adult mosquitoes might not directly happen in infected cells per se, but could reflect an indirect effect of homeostatic or metabolic responses during the infection. It should be remembered that miRNA expression can be highly cell-type specific; if miRNA levels are responsive to DENV2-infection only in selected cell types within the entire mosquito, we may miss those in our Aag2 cell-based assays.
In summary, here we provide an in-depth analysis of small RNAs in DENV2 infected Aag2 cells in comparison to uninfected cells. Aag2 cells provide a powerful model system for studying biochemical details of small RNA biogenesis pathways, as they are fully competent in producing siRNAs, miRNAs and piRNAs originating from both virus and host. Our analyses add both DENV-derived piRNAs and novel Aedes aegypti miRNAs to the small RNA repertoire in this medically important virus-host interaction.  [57,58], and have also been recovered by the miRDeep2 analyses of our libraries. These miRNAs were not yet published in the most recent version of miRBase. The folding of the hairpin was predicted using RNAfold. Red and blue letters indicate the position of the predicted, mature miRNA and miRNA Ã sequences, respectively. The shown miRNA names have been assigned by miRBase; the previous names given by the authors of the indicated studies are presented between brackets. The mature miRNA sequences of aae-miR-11893-1 maps to similar hairpins at two locations in the Aedes genome. Also aae-miR-11894a-1 and aae-miR-11894a-2 are identical sequences that map to two very similar hairpins. Aae-miR-11895-1 and aae-miR-11895-2 are two identical hairpins that map to two locations in the Aedes genome. Similarly, the mature aae-miR-11899 sequence maps to two identical hairpins in the genome. (TIF) S1 Table. Name and sequence of oligonucleotides used in this study. (PDF) S2 Table. Spreadsheet of raw miRDeep2 output. Sheet 1 provides an overview of the selection procedure to reach the final miRNA predictions. Sheets 2 to 6 show the unfiltered miRNA predictions and the filtered miRNA predictions after the individual rounds of selection. (XLSX) S1 Dataset. Data related to Figs 1C, 1D, 2E, 2F, 3, 4, S1, S2, S3B and S3C. (XLSX) and Prevention, Fort Collins, CO, USA) for providing DENV2 16681. We would also like to thank Sam Griffith-Jones and the team from the miRBase repository for assistance with assigning miRNA names.