Rapid Intraspecific Evolution of miRNA and siRNA Genes in the Mosquito Aedes aegypti

RNA silencing, or RNA interference (RNAi) in metazoans mediates development, reduces viral infection and limits transposon mobility. RNA silencing involves 21–30 nucleotide RNAs classified into microRNA (miRNA), exogenous and endogenous small interfering RNAs (siRNA), and Piwi-interacting RNA (piRNA). Knock-out, silencing and mutagenesis of genes in the exogenous siRNA (exo-siRNA) regulatory network demonstrate the importance of this RNAi pathway in antiviral immunity in Drosophila and mosquitoes. In Drosophila, genes encoding components for processing exo-siRNAs are among the fastest evolving 3% of all genes, suggesting that infection with pathogenic RNA viruses may drive diversifying selection in their host. In contrast, paralogous miRNA pathway genes do not evolve more rapidly than the genome average. Silencing of exo-siRNA pathway genes in mosquitoes orally infected with arboviruses leads to increased viral replication, but little is known about the comparative patterns of molecular evolution among the exo-siRNA and miRNA pathways genes in mosquitoes. We generated nearly complete sequences of all exons of major miRNA and siRNA pathway genes dicer-1 and dicer-2, argonaute-1 and argonaute-2, and r3d1 and r2d2 in 104 Aedes aegypti mosquitoes collected from six distinct geographic populations and analyzed their genetic diversity. The ratio of replacement to silent amino acid substitutions was 1.4 fold higher in dicer-2 than in dicer-1, 27.4 fold higher in argonaute-2 than in argonaute-1 and similar in r2d2 and r3d1. Positive selection was supported in 32% of non-synonymous sites in dicer-1, in 47% of sites in dicer-2, in 30% of sites in argonaute-1, in all sites in argonaute-2, in 22% of sites in r3d1 and in 55% of sites in r2d2. Unlike Drosophila, in Ae. aegypti, both exo-siRNA and miRNA pathway genes appear to be undergoing rapid, positive, diversifying selection. Furthermore, refractoriness of mosquitoes to infection with dengue virus was significantly positively correlated for nucleotide diversity indices in dicer-2.

Much of what we know about RNAi in insects has been elucidated in Drosophila melanogaster, where the biogenesis and regulatory functions of each of the small RNA classes have been separated into distinct pathways [7]. The exo-siRNA pathway has a central role in Drosophila antiviral immunity [8], [9] and is initiated by Dicer-2 (Dcr2). Dcr2 is an RNase III family protein that recognizes cytoplasmic long dsRNA and cleaves it into ,21 bp siRNAs [10,11]. The siRNAs, in association with Dcr2 and the dsRNA-binding protein R2D2, are loaded into a multiprotein RNA-induced silencing complex (RISC), which contains Argonaute-2 (Ago2) [12,13,14]. In the effector stage of the pathway, the RISC unwinds and degrades one of the siRNA strands and retains the other strand as a guide for recognition and sequence-specific cleavage of viral mRNA, mediated by the ''slicer'' endonuclease activity of Ago2 [15,16,17,18].
MicroRNAs (miRNAs) are 22-23 nt RNA guides that regulate cellular functions such as differentiation and development and metabolic homeostasis. Although only invertebrates have siRNAs, both vertebrates and invertebrates have miRNAs, which are transcribed from the cellular genome as primary miRNAs by RNA polymerase II and are processed sequentially by two distinct endonucleases in the RNase III family, nuclear Drosha and cytoplasmic Dicer 1 (Dcr1), the only ortholog of the dcr gene family in mammals. Dcr1 processes pre-miRNA to imperfectly basepaired duplex miRNA with ,23 nt strands and acts with the dsRNA-binding protein R3D1 to load the miRNA guide strand into an Argonaute-1 (Ago1)-containing RISC [13,19]. Typically miRNAs recognize targets in the 39 non-coding region of cellular mRNAs by imperfect complementarity and inhibit their translation [20].
There is currently a great deal of interest in identifying genes that condition the ability of arthropods to transmit RNA viruses that are pathogenic to humans and domestic animals. Of particular interest is the mosquito Aedes aegypti, which is an important vector of a number of pathogenic arthropod borne viruses (arboviruses), including the dengue viruses (DENV1-4), yellow fever virus and chikungunya virus, and is also a tractable genetic system with which to identify candidate genes [21]. Aedes aegypti is distributed in all subtropical and tropical regions of the world. Most importantly, Ae. aegypti populations demonstrate a great deal of variation in their susceptibility to arboviral infection [22].
Several lines of evidence suggest the importance of the exo-siRNA pathway in antiviral immunity in Drosophila and mosquitoes. Drosophila with mutations in or depletion of known exo-siRNA pathway components are hypersensitive to RNA virus infections and develop a dramatically increased viral load [9,23,24]. Increases in arboviral replication occur after knockdown of one or more genes in the exo-siRNA pathway [25,26]. siRNAs derived from the infecting virus genome (viRNAs) have been discovered and characterized in infected insects [27,28,29,30]. Many insect pathogenic viruses encode suppressors of RNAi that counteract insect immunity [31].
Noting that interaction between RNA viruses that encode suppressors of RNAi and their insect hosts may lead to a coevolutionary ''arms race'' and directional selection on RNAi genes, Obbard et al. [32] undertook a comparative study of the rates of amino acid evolution in exo-siRNA and miRNA pathway components in three species of Drosophila. They showed that among Drosophila species, the ratio of replacement to silent amino acid substitutions (w = K A /K S ) among the exo-siRNA genes dcr2, r2d2, and ago2 is much greater than w among their miRNApathway counterparts dcr1, r3d1, and ago1. In fact it was shown that Dcr2, R2D2, and Ago2 are among the fastest evolving 3% of all Drosophila proteins [32]. Recent selective sweeps in ago2 have reduced genetic variation across a region of more than 50 kb in the genomes of Drosophila melanogaster, D. simulans, and D. yakuba, and it was estimated that selection has fixed adaptive substitutions in this gene every 30-100 thousand years [33]. The rapid evolution of exo-siRNA pathway genes compared with miRNA pathway genes supported a hypothesis for directional selection on host antiviral RNAi genes driven by evolution of virus-encoded suppressors of RNAi (VSRs) that enable viruses to evade RNA silencing [32,34].
More than 50 VSRs encoded by plant and insect pathogenic viruses have been described. The Flock House virus (Nodaviridae) protein B2 is one of the best characterized animal VSRs [35]. B2 binds both siRNA duplexes and long dsRNA, thereby preventing dsRNA binding by proteins in the exo-siRNA pathway. No arbovirus VSRs have been identified [31,36,37]; however, mosquito-borne alphaviruses were engineered to express B2 and then used to infect mosquitoes orally or by injection [28] [38]. These mosquitoes had reduced pools of viRNAs, increased infectious virus titers and, most importantly, greatly decreased survival rates relative to mosquitoes infected with non-recombinant viruses.
The ability of arboviruses to cause persistent, non-cytopathic infections in both mosquito cells and mosquitoes despite the RNAi response has led to speculation about arboviral mechanisms of immune suppression or evasion in insect cells. Entomopathogenic VSRs are generally virulence factors that increase the likelihood of virus transmission by killing the insect hosts; however, pathogenic rather than persistent infections of mosquitoes by arboviruses would be detrimental to transmission and maintenance in nature. Thus, arboviral mechanisms of evasion are unlikely to involve VSRs that increase virulence, but could involve strategies such as rapid evolution of genome 'decoy' regions [39] or RNAi-escape mutations [29].
A recent review concluded that mosquito RNAi is the major innate immune pathway controlling arbovirus infection and transmission in mosquitoes in a similar way to Drosophila antiviral immunity [40]. However, there has been no characterization and comparison of the molecular evolution of miRNA and exo-siRNA genes in a mosquito. Herein we describe the intraspecific patterns of molecular evolution of ago1, ago2, dcr1, dcr2, r3d1, and r2d2 within and among collections of Ae. aegypti from throughout its geographic range. We tested whether the interspecific evolutionary patterns of small RNA pathway genes among Drosophila species [32] are also apparent intraspecifically in Ae. aegypti. We compared nearly complete sequences of all exons in each of the six genes from 104 individual Ae. aegypti collected in six geographically distinct sites throughout the mosquito's range. We determined if amino acids encoded by the six major genes in the two pathways appear to be under positive selection by performing a fixed-site phylogenetic analysis by maximum likelihood (PAML) [41] on each of the six genes. Identified variable sites were further characterized as to the likelihood that they would alter protein structure or function using the program SIFT (Sorting Intolerant From Tolerant) [42]. Finally, we sought to determine if sequence diversity in the six genes is correlated with susceptibility to infection with DENV2 by calculating Pearson's correlation coefficients between vector competence data gathered in previous studies for four of the six mosquito collections [43,44,45] and various measures of genetic diversity.

Mosquito strains and DNA extraction
Six geographically distinct Ae. aegypti populations were analyzed in this study ( Figure 1). DNA was analyzed from 18 Ae. aegypti individuals collected from Poza Rica, 18 from Lerdo de Tejada, and 18 from Chetumal in Mexico. Details on these collection sites and determination of their vector competence for DENV2 were published previously [43,44]. DNA was analyzed from 20 mosquitoes from PK-10 (near Kedouguo) and ten from Mindin, Senegal [45]. All mosquitoes from PK-10 and five from Mindin were the formosus subspecies of Ae. aegypti as determined by the absence of silver scales on the first abdominal tergite [46]. Vector competence in Ae. aegypti formosus for flaviviruses tends to be lower than in subspecies Ae. aegypti aegypti [47,48,49]. The 20 Ae. aegypti from Pai Lom, Thailand were from a laboratory colony provided by Dr. L. C. Harrington at Cornell University. DNA was extracted from all 104 mosquitoes using the Qiagen DNeasy Blood and Tissue Kit (Valencia, CA).
Single strand conformational polymorphism (SSCP) analysis was performed on all PCR products [50]. Gels were scored based on the SSCP banding patterns observed for each of the 56 amplicons for each individual mosquito in the study. Unique haplotypes were identified and recorded from each SSCP gel. Each unique haplotype for each of the 56 amplicons was sequenced on both strands on the ABI 31306L Genetic Analyzer at the Proteomics and Metabolomics facility at Colorado State University. Sequences were obtained from at least two mosquitoes to test the sensitivity of the SSCP technique when a haplotype pattern occurred more than once.

Sequence analysis
Forward and reverse strand sequences were assembled for each unique haplotype with SeqMan 2 (DNAStar, Madison, WI). Amino acid sequences for each locus were aligned in MAFFT ver. 6.606 [51] using the iterative refinement method ($1,000 iterations) with consistency and weighted sum-of pairs scores (G-INS-i). The corresponding nucleotide sequence alignments were derived from amino acid alignments in MacClade ver. 4.03 [52]. No indels were inferred for any of the six sequence sets. DNAsp 5.10 (http://www.ub.es/dnasp) was used to determine the number of haplotypes (h), numbers of segregating sites that were synonymous (S s ) or caused amino acid replacements (S A ), and the average number of nucleotide differences per site between all sequence pairs (pi) [53]. DNAsp also estimated K, the average number of nucleotide differences between all sequence pairs [54] (equation A3) and among replacement sites (K A ) and synonymous sites (K S ) and their ratio (w = K A /K S ). The effective recombination rate between adjacent sites (R), the PHI test for recombination and Fu and Li's F* test were also calculated. The w ratio was used to infer the action of natural selection from comparative sequence data [41]. If replacements/substitutions are deleterious and therefore subject to purifying selection, K A ,K S and w,1. Alternatively, when replacement substitutions are favored by natural selection K A .K S and w.1.

Phylogenetic analyses
Maximum likelihood (ML) analyses of nucleotide characters from each of the molecular data matrices were performed using RAxML ver. 7.0.3 [55]. Because RAxML only implements general time reversible (GTR) Q-matrices for nucleotide characters, more restrictive variants of the GTR matrix were not used when selected by the Akaike Information Criterion. Optimal likelihood trees were searched for using 1,000 independent searches starting from randomized parsimony trees with the GTR-GAMMA model and four discrete rate categories. Likelihood bootstrap (BS) analyses [56] were conducted with 2,000 replicates with ten searches per replicate using the ''-f i'' option, which ''refine[s] the final BS tree under GAMMA and a more exhaustive algorithm'' [57].

Tests for positive selection
CodeML [41] in the PAML package was used to perform a maximum likelihood analysis of protein-coding DNA sequences using codon substitution models [58]. CodeML estimated K A and K S and performed likelihood ratio tests (LRTs) of positive selection along lineages based on w to identify amino acid sites potentially under positive selection. We loaded the ML tree topology derived from RAxML and enforced that as a constraint. PAML was set to explicitly account for the nucleotide content at each codon position when calculating K A and K S . PAML optimized the branch lengths based on the particular model that it employed for each analysis. Five models (M0, M1a, M2a, M7 and M8) were compared. These five were reported to be the most effective models for detecting positive selection based on both simulations and empirical data [41]. Model M0 assumes and calculates one v for all codons. Model M1a is a neutral model that assumes two site classes: w 0 ,1 (estimated empirically from the data) and w 1 = 1. Model M2a is a selection model that is compared with M1a by a LRT. It adds a third site class to M1a, with w 2 .1 estimated empirically. Model M7 (beta) is a flexible null model, in which v for a codon    [59]. A star tree is commonly used as a null model in phylogenetic comparative methods. All branches emerge from a single common ancestor (no topology) rather than emerging internally from one another and all branches are of equal length (no differential stasis). All branches in a star tree evolve independently of one another. Recombination among nuclear genes can homogenize haplotypes and thus create a star phylogeny. If a phylogeny derived from a dataset is resolved as a star phylogeny or if the data fit a star phylogeny as well as they fit any alternative phylogeny then the branches are said to evolve independently. Independence is desirable when testing for evidence of selection because specific hierarchies in which branches evolve in a dependent manner may lead to false detection of sites under selection because the hierarchy may be incorrect for sites under selection [59]. Results obtained from the analysis of star trees were highly similar to those obtained from the estimated gene tree (results not shown), suggesting that positive selection results were not seriously affected by possible recombination events.

Phenotype correlation analysis
Phenotypic data consisted of DENV2 midgut infection barriers (MIB = proportion of orally exposed female mosquitoes that fail to develop a midgut infection) and midgut escape barriers (MEB = proportion of females with a midgut infection that fail to develop a disseminated infection). Phenotypic data were not available for the Thailand or Mindin collections. Pearson correlation coefficients and Fisher exact tests were performed to compare the average number of nucleotide differences per site between all sequence pairs (pi), the average number of nucleotide differences between all sequence pairs (K), among replacement sites (K A ), among synonymous sites (K S ) and w between each phenotype (i.e. MIB, MEB) using R v2.11.1 (http://cran.r -project.org/).

Results
Intraspecific patterns of siRNA and miRNA gene variation in Ae. aegypti Table 2 lists, for each gene and each mosquito collection, sample sizes (N), numbers of segregating sites encoding synonymous (S s ) and replacement (S a ) substitutions, numbers of haplotypes (h), average numbers of nucleotide differences per site between all sequence pairs (pi), average numbers of nucleotide differences between all sequence pairs (k), average numbers of nucleotide differences among replacement sites (K A ), average numbers of nucleotide differences among synonymous sites (K S ), K A /K S = w, effective recombination rates between adjacent sites (R), and the PHI tests for recombination alongside Fu and Li's F* test. Graphs of w for all six genes appear in Figure 3. In Ae. aegypti, w was 1.4-fold higher in dcr2 than in dcr1, as compared to 5.4-fold higher in an interspecific study in Drosophila [32]. The Ae. aegypti v ratio in ago2 was 27.4 fold higher than in ago1 (the K A for ago1 in all Drosophila spp. was zero [32]) . The nearly six-fold higher w ratios found in r2d2 as compared to r3d1 in Drosophila spp. [32] were not found among Ae. aegypti collections. Instead w was approximately the same in both genes. The probability values from Fisher's exact test of S a and S s appear above the bar graphs for the three gene pairs in Figure 3; only the ago1 vs. ago2 comparison in mosquitoes was significantly different. In contrast, among Drosophila spp. all three comparisons were significant.
Obbard et al. [32] found no replacement substitutions in ago1 among Drosophila spp. In contrast, we found 10 nucleotide substitutions that caused amino acid replacements in Ae. aegypti ago1. In ago2, v was 2.5 fold higher among Drosophila spp. than among Ae. aegypti populations. In dcr1, w was approximately the same among Drosophila spp. as among Ae. aegypti populations, while w in dcr2 was 3.5 fold higher among Drosophila spp. than among Ae. aegypti populations. In r3d1, w was 2.2 fold higher among Drosophila spp. than among Ae. aegypti populations while in r2d2, v was 13.2 fold higher among Drosophila spp. The same trends in w occurred in all six Ae. aegypti collections ( Table 2).
While not significant, comparison of w between Dcr1 and Dcr2 among Ae. aegypti populations reveals the same trend as among Primer locations are numbered with respect to supercontigs in the Ae. aegypti genome project in Vector Base ( Figure 2). An * indicates that all or part of that primer is located in an intron. doi:10.1371/journal.pone.0044198.t001 Drosophila spp, with a higher v in the siRNA pathway gene. In contrast, replacement substitutions in ago1 were detected in four of the six Ae. aegypti collections while no replacement substitutions were found among four species of Drosophila. The same large disparity in v between Drosophila Ago1 and Ago2 was evident among Ae. aegypti populations, but the same was not true of the Ae. aegypti dsRNA-binding proteins R3D1 and R2D2, amongst which w was approximately identical.  Table 3). The average numbers of nucleotide differences per site between all sequence pairs (pi) were compared for each gene. The proportion of segregating sites in regions of known function versus regions where no function has been assigned to date were not significantly different for any of the six genes as determined by Fisher's Exact Test. Average values of pi were compared between regions of known versus unassigned function using a Student's ttest and a significant difference was only seen for R3D1 in which pi was greater in regions of known function (Table 3).

Tests for positive selection
Maximum likelihood analysis of protein-coding DNA sequences was used to estimate the average number of nucleotide differences among replacement sites (K A ) and synonymous sites (K S ) and  Sample size (N), numbers of segregating sites that were synonymous (S s ) or led to amino acid replacements (S a ), the number of haplotypes (h), the overall pi(the average number of nucleotide differences per site between two sequences) are listed for each gene and collection. Also listed are the average number of nucleotide differences, k and K among synonymous sites (K s ) and among replacement sites (K a ) and their ratio (K a /K s ). The effective recombination rate between adjacent sites (R) and the PHI test for recombination are listed alongside Fu and Li's F* test and associated significance tests (*P,0.05, **P,0.01). doi:10.1371/journal.pone.0044198.t002 perform LRTs for positive selection. The ML tree topologies derived from RAxML were enforced as a constraint in CodeML from the PAML package [41]. Five models (M0, M1a, M2a, M7 and M8) were compared and results for each of the five models and each of the six genes are shown in Table S1. For each gene the model with the greatest likelihood score is highlighted in grey. In each case the M2a model had a significantly better fit than the M1a model and the M8 model had a significantly better fit than the M7 model, implying that models of positive selection had a better fit than neutral models. All of the amino acids identified using the naïve empirical Bayes (NEB) method were also identified using the Bayes empirical Bayes (BEB) method, while a few additional sites were identified with BEB. The alternative amino acids found to be under positive selection using BEB are listed in Table S2 alongside w from the M8 model and its standard error. All replacement to synonymous substitution ratios (w) are mapped across dicer genes in Figure 4, across argonaute genes in Figure 5 and across genes encoding dsRNA-binding proteins in Figure 6. Even though clusters are visually evident, the distances between positively selected sites (PSSs) were formally subjected to a chi-square goodness-of-fit test to determine if they were distributed in a random (Poisson) fashion. Only ago2, dcr1, and dcr2 had enough sites to perform this test and the PSSs all three of these genes were clustered rather than randomly distributed (analyses available on request).
An additional test was conducted on PSSs using the program Sorting Intolerant from Tolerant (SIFT) [42] to determine whether the amino acid changes in Table S2 are predicted to affect protein function. SIFT scores replacement substitutions on a scale from 0-1, where scores at or below 0.05 are likely to change protein function, whereas higher scores are not. Thirty PSSs were detected in Dcr1, eight (27%) of which were likely to change protein function; in Dcr2, 16 of 38 (42%) PSSs were likely to change protein function (Figure 4). In Dcr1, a single cluster consisting of positively-selected sites 384, 385, 388 and 395 was detected in a region without assigned function between the DExD/H-like helicase and helicase superfamily C-terminal domains. Replacements at PSSs 497, 502 and 592 were in the helicase superfamily C-terminal domain and two of these were predicted to change function. Replacements at PSSs 795, 817 and 978 were in the dsRNA-binding domain but were not predicted to change function. In Dcr2, one cluster of PSSs occurred in the ribonuclease III carboxyl-terminal domain and included amino acid sites 1446, 1450, and 1454, which are among the residues responsible for dimerization. Missense mutations in an RNase IIIencoding domain of Drosophila dcr2 resulted in a profound loss of dsRNA processing activity and destabilization of the protein [24]. A second Dcr2 PSS cluster consisted of 8 PSSs, six with significant SIFT scores, in a region of unassigned function between the ribonuclease III C terminal domain and the dsRNA binding domain. A cluster also occurred in a region of unassigned function amino-terminal to the PAZ domain.
Replacement substitutions in ago1 were detected in four of the six Ae. aegypti collections even though no replacement substitutions were found in ago1 among four species of Drosophila. It is unclear why diversifying selection would occur within and among collections of Ae. aegypti while only purifying selection is evident among Drosophila spp. Three PSSs were identified in Ago1 but none of these was predicted to change protein function. In strong contrast, 16 of 38 PSSs were likely to change protein function in Ago2 ( Figure 5). Four clusters of PSSs were found in Ago2; one of the clusters occurred in the amino-terminal region of the PAZ domain at residues 372, 375, 378, 383, 385, and 390, and the 383 and 390 replacements had significant SIFT scores. However, the oligonucleotide binding domain of PAZ occurs in the carboxylterminus in residues 418-472 according to GenBank annotation ACR56327. The amino-terminus of ago2 encodes a series of polyglutamines, however the function(s) of these residues are unknown. One intriguing PSS occurred in the 59 RNA guide strand anchoring site in the amino-end of the Piwi domain. This corresponds to amino acid residue 700 and involves a threoninemethionine replacement. Of two PSSs identified in R3D1, one was . The ratio (w = K a /K s ) for miRNA (ago1, dcr1, and r3d1) and siRNA (ago2, dcr2, and r2d2) pathway genes among Ae. aegypti collections. The average number of nucleotide differences among replacement sites (K a ) relative to the average number of nucleotide differences among synonymous sites (K a ). doi:10.1371/journal.pone.0044198.g003 likely to change protein function. In R2D2, four of five PSSs were predicted to alter function.

Phenotype correlation analysis
Phenotype correlation analysis was performed to identify potential correlations of mosquitoes with DENV2 midgut infection barrier (MIB) and escape barrier (MEB) phenotypes with the number of nucleotide differences between all sequence pairs (k), nucleotide differences per site between all sequence pairs (pi) and and/or the number of nucleotide differences among replacement sites (K A ) and among synonymous sites (K S ) and/or the K A /K S ratio (w). Phenotype data were available for Poza Rica (MIB: 21%, MEB: 18%), Lerdo de Tejada (MIB: 25%, MEB: 36%), and Chetumal, Mexico (MIB: 9%, MEB: 8%) [43], [44] and PK10, Senegal (MIB: 8%, MEB: 76%) [45] collections. Pearson correlation coefficients are displayed in Table 4. No significant MIB correlations were found. No significant MEB correlations with these variables were observed for the argonaute, r3d1 or r2d2 genes. However, w in dcr1 was significantly positively correlated with MEB and pi, k, and K S were significantly correlated with MEB for dcr2. This suggests the possibility that diversity in Dicers may increase the frequency of mosquitoes with MEBs.

Discussion
Intraspecific patterns of variation between miRNA and siRNA pathway genes in Ae. aegypti In contrast to Drosophila spp., we found that in Ae. aegypti mosquitoes both exo-siRNA and miRNA pathway genes appear to be undergoing rapid, positive, diversifying selection. However, in similarity to Drosophila, most positively selected sites occurred in protein regions without defined functions (Table S2 and Figs. [3][4][5]. Our observations are consistent with a hypothesis that diversifying selection acts on both dcr1 and dcr2 to maintain the intraspecific w ratios among dcr1 genes of Ae. aegypti populations at the same level as among Drosophila species. Mandibulate arthropods are unique in having two Dicers [34]. With the exception of Cnidarians and Porifera, the remainder of animals so far examined, including vertebrates, possess a single Dicer that produces miRNAs and, in the case of invertebrates such as C. elegans, also siRNAs. It is possible that Dcr1 retains some role in antiviral RNAi in mosquitoes but not in Drosophila. Mosquitoes (Culicidae) are members of the primitive fly suborder Nematocera while Drosophilidae are one of the most recently evolved of fly families [60]. In addition, replication of some mammalian viruses has been shown to be indirectly inhibited by host miRNAs and some viruses exploit host miRNAs during replication [61]; however, potential roles of mosquito miRNAs in arbovirus replication have not been explored. It is possible that components that have been assigned functions in distinct RNA silencing pathways, including the miRNA pathway, interact with or serve as redundant or backup antiviral mechanisms for the exo-siRNA pathway in insects. Evidence of interactions between components of RNA silencing pathways in Drosophila was provided by the demonstration that R3D1-Dcr2 heterodimers, rather than the canonical R2D2-Dcr2 complexes, are involved in Ago2-RISC-loading of endo-siRNAs [5,6]. The endo-siRNA pathway, which has been shown to function in suppressing transposon activity in somatic cells of Drosophila, can also be triggered by transient transfection of exogenous dsRNA, suggesting a potential role in antiviral defense [62]. Potential interactions between siRNA and miRNA pathways in mosquitoes, particularly in antiviral defense or control of transposon activity, remain to be examined.
Evidence for a role of the piRNA pathway in insect antiviral defense also has emerged recently. In our examination of antiviral RNAi in Anopheles gambiae mosquitoes we found that dsRNAmediated silencing of the ago3 gene concurrent with O'nyongnyong virus infection resulted in increased virus titers, hinting at a possible role for Ago3 in antiviral immunity [25]. Furthermore, RNA virus-specific piRNAs, in addition to viral siRNAs, were recently described in Drosophila ovary cells [63] and other studies showed that piwi-family mutants of Drosophila were more susceptible to Drosophila virus X infection than wild-type flies [8]. We also found that in cultured C6/36 (Aedes albopictus) cells a singlenucleotide deletion in dcr2 causes a defect in the exo-siRNA pathway-mediated antiviral defense that was apparently compensated by the piRNA pathway [30,64]. The redundant role of the piRNA pathway in antiviral defense in mosquito somatic cells and  its particular relevance in dcr2-null cell culture lines has recently been confirmed [65].

Sources of diversifying selection
Obbard et al. [34] suggested that a ''molecular arms race'' with pathogenic virus suppressors of RNAi drives siRNA pathway gene diversity in Drosophila spp. There are various reasons to believe that infections with arboviruses are not the drivers of diversifying selection that we have documented in the miRNA and siRNA pathways in Ae aegypti. A hallmark of arboviruses is that they have little, if any, effect on survivorship or fecundity in their insect hosts [66]. Two studies further suggest that selection has, in fact, prevented the evolution of potent VSRs in arboviruses [28,38]. Even if infection by arboviruses imposed a minimal fitness effect, a number of field studies have demonstrated that very few mosquitoes (10 23 to 10 24 ) collected in areas endemic for certain arboviruses are actually infected at any given time [67,68,69] thus providing only rare opportunities for selection. An interesting alternative possibility for selection may involve the recently discovered high prevalence of persistent insect-only flaviviruses in natural populations [70,71,72]. These viruses are maintained through vertical transmission from one generation to the next without obvious pathogenesis and without requiring horizontal transmission through infected vertebrates.
We suggest that infections with entomopathogenic viruses or transposon invasion and movement are more likely causes of the diversifying selection detected in this study. Unfortunately, few mosquito-pathogenic RNA viruses have been identified or wellstudied. One such virus, Nodamura virus (Nodaviridae), was isolated from Culex tritaeniorhynchus in Japan [73] and can experimentally infect and produce pathogenesis in Ae. aegypti [74]. It has a bipartite, positive-strand RNA genome that replicates through a dsRNA intermediate. Unique features of viruses in this family are that their replication complexes are sequestered in membraneenclosed spherules within the mitochondrial outer membrane during replication and they encode B2-type VSRs [75]. Small RNAs that appeared to be derived from the RNA genome of a previously-undescribed nodavirus were recently discovered by deep sequencing analysis of a small RNA library from Ae. aegypti [63], suggesting that other potentially pathogenic mosquito viruses remain to be found.

Phenotype correlation analysis
Our a priori hypothesis was that Ae. aegypti collections that were more refractory for DENV2 disseminated infection would also have higher rates of evolution in genes encoding components of the siRNA pathway compared to DENV2 susceptible populations. The trends shown in correlation of vector competence with certain measures of genetic diversity in RNAi pathway genes in Table 4 need to be tested in more Ae. aegypti populations and possibly in other Aedes species. Furthermore, we need to perform quantitative trait loci mapping and association studies to test for a correlation between miRNA and siRNA pathway alleles and arbovirus susceptibility. If a correlation is detected it could suggest that RNA silencing evolved in mosquitoes as a means to combat entomopathogenic virus infection or genome invasion by transposons but that this evolution may have indirectly provided a regulatory mechanism for replication and transmission of arboviruses. The proportions of segregating sites in regions of known function versus regions where functions have not been assigned were compared by Fisher's Exact Test (FET) and the resulting probability is listed. Average values of pi(average number of nucleotide differences per site between sequence pairs) were compared between regions of known versus unassigned function using a student's t-test and a significant difference was only seen for R3D1, in which pi was greater in regions of known function. doi:10.1371/journal.pone.0044198.t003 Table 4. Phenotype correlation analysis. To identify potential correlations of Ae aegypti susceptibility to dengue virus infection with measures of nucleotide and amino acid diversity in RNAi genes, Pearson correlation analyses were performed between DENV-2 MIB (the proportion of orally virus-exposed female mosquitoes that fail to develop a midgut infection) and MEB (the proportion of females with a midgut infection that fail to develop a disseminated infection) with p (number of nucleotide differences per site between all sequence pairs), k (number of nucleotide differences between all sequence pairs), K A (number of nucleotide differences among replacement sites), K S (number of nucleotide differences among synonymous sites) and K A /K S . Pearson correlation coefficients are shown. *P#0.05, **P#0.01. doi:10.1371/journal.pone.0044198.t004