Identification of conserved and functional amino acids in Aaeg Dcr2

Dcr2 proteins are closely related to share functional domains across invertebrate animal species. To determine functionally relevant amino acids in the helicase and RNase III domains, we compared Aaeg Dcr2 to Dmel Dcr2. Multiple sequence alignment of orthologous Dcr2 sequences led to the identification of residues in the helicase domain with K39, Y232, and G488 in Aaeg Dcr2, whose equivalents are critical in Dmel Dcr2 function [35,37,38,40]. Moreover, these residues are part of known helicase motifs [47] that are conserved across divergent species. In addition, D1198, E1341, D1444, and E1548 in Aaeg Dcr2 RNase III A and B domains are well-conserved residues across species that are likely equivalent and have been shown to be functionally important for RNAi activity in Dmel Dcr1 and Dcr2 [37,38,48] (Fig 1A and 1B).

Taking these conserved amino acids, mutations were introduced to Aaeg Dcr2 by mutagenesis cloning of a previously developed expression plasmid, pPUb-myc-Dcr2, encoding a myc-tagged wild-type (WT) Dcr2 (GenBank ID: AAW48725) under the control of a polyubiquitin (PUb) promoter [46] to assess loss of function (Fig 1B, top panel). Mutations individually introduced into the helicase domain were K39N in Motif I, Y232G in Motif IV, and G488R in Motif VI. For the RNase III domains, an expression plasmid with mutations D1198A and E1341A in the RNase IIIA, and D1444A, E1548A in the RNase IIIB domain were collectively introduced and called mutant RNase III (mtR3), with mutations to alanine expected to inactivate RNase III activity. Sorting Intolerant From Tolerant (SIFT) algorithm [49] prediction scores suggested that the individual mutations would affect Dcr2 function (S1 Table). The expression of Dcr2 was verified by immunoblotting lysates obtained from AF319 Dcr2 knockout cells transfected with mutant or WT Dcr2 expression plasmids, with pPUb-myc-eGFP (peGFP) as control (Fig 1B, bottom panel).

To determine the effect of the mutations on Dcr2 silencing activity, a plasmid-based RNAi reporter assay was performed [46]. Briefly, AF319 cells were co-transfected with firefly luciferase (FFLuc) reporter (and Renilla luciferase, RLuc, expressing control plasmid) and WT or mutant Dcr2 plasmids (with peGFP as control) together with dsRNA targeting FFLuc (dsFLuc) or a non-silencing control (dsLacZ). At 24 h post-transfection (hpt), it was observed that all mutants except Y232G mediated no silencing activity against FFLuc and were comparable to dsLacZ-transfected controls (Fig 1C). Specifically, Y232G retained residual silencing that was comparable to WT Dcr2. Its equivalent mutant in Dmel Dcr2 (F225G) also induced vsiRNA production in vivo [40]. The findings suggested that the residues previously identified from studies on Dmel Dcr2 were also important in Aaeg Dcr2 activity, and both helicase and RNase III domain activity appeared to be critical in the overall silencing function and exogenous dsRNA cleavage of Aaeg Dcr2. To verify the functionality of the exo-siRNA pathway downstream of mutant Dcr2, mutations K39N and mtR3 were assessed following transfection of siRNA targeting FFluc in a plasmid-based reporter assay (S1 Fig). Silencing activity was readily detected compared to the control, suggesting that processes post dsRNA cleavage and siRNA loading into Ago2 remained active.

The ability of the loss-of-function mutants to limit arbovirus replication was also examined. The replication of SFV-FFLuc (MOI = 0.1) reporter virus in AF319 cells transiently expressing mutant or WT Dcr2 was assessed (Fig 1D). At 48 h post-infection (hpi), expression of WT Dcr2 resulted in significant reduction of SFV replication compared to eGFP-transfected control cells; no mutant Dcr2 was able to reduce SFV replication as determined by FFLuc levels relative to the eGFP-transfected controls. Including Y232G, which showed some activity in the RNAi reporter assay as mentioned above.

Naturally occurring Dcr2 polymorphisms in field-derived Ae. aegypti mosquitoes

To determine the degree of conservation of the helicase and RNase III domains, Dcr2 cDNA from adult female Ae. aegypti originating from Cameroon (Ben), French Guiana (Cay), Guadeloupe (Guad), Cambodia (KC and Rat), and Gabon (Lope) were sequenced. Thirteen haplotype sequences were identified (S2 Table) and examined further for single nucleotide polymorphisms (SNPs) by comparing them to our WT Dcr2 reference sequence (GenBank ID: AY713296) that was previously cloned to obtain pPUb-myc-Dcr2 [46]. SNP analysis showed that the majority of the polymorphisms were silent transition changes, but it was observed that all field-derived mosquito Dcr2 had a C variant nucleotide at position 1684 (S3 Table). In addition, SNPs were identified spanning the entire length of Dcr2 with a number resulting in amino acid substitutions (Fig 2A).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Naturally occurring Dcr2 variants among field-derived mosquito populations.

(A) Dcr2 haplotypes (in grey) were identified from early generation colonies of Ae. aegypti initially collected from Cameroon (Ben), French Guiana (Cay), Guadeloupe (Guad), Cambodia (KC and Rat), and Gabon (Lope) and aligned to a reference WT Dcr2 sequence (GenBank ID: AY713296; thick black line) to determine SNPs. Black bars show disagreements to WT Dcr2 with silent and non-synonymous SNPs highlighted by cyan and magenta bars, respectively. (B) Phylogenetic tree of putative amino acid sequences of haplotypes highlighting the amino acid net charge at positions 246/253 (+: positive, -: negative, =: neutral) and an H or N at position 1514. Numbers on branches indicate % consensus support with the substitution rate measured by the scale bar. (C) Alignment of haplotypes highlighting amino acid variants N246D and Q253K in the helicase, and H1514N in the RNase IIIB domain. (D) Mutagenesis of WT Dcr2 to generate amino acid mutations at position 246/253 (aa 246/253) and a theoretical "D/Q" variant (N246D; negative charge) not observed among the surveyed mosquito populations. To check for Dcr2 expression, AF319 cells transfected with variant or WT Dcr2 (peGFP as control) were lysed and probed for myc-Dcr2 by immunoblot using anti-myc tag and γ-tubulin (loading control) antibodies. Representative blot of n = 2 independent repeats is shown. (E) RNAi reporter assay was performed to determine the effect of the variants on silencing activity. AF319 cells co-transfected with WT or variant Dcr2 plasmids (peGFP as control) and FFLuc and RLuc reporter plasmids with dsRNA targeting FFLuc (dsFLuc) or dsLacZ (non-targeting control). At 24 hpt, luciferase levels were measured and shown as mean±SEM relative light units (FFLuc/RLuc) versus dsLacZ set to 1 compared to eGFP control from n = 3 independent repeats with * = p<0.05 according to two-way ANOVA. (F) Antiviral activity of Dcr2 variants against SFV-FFLuc (MOI = 0.1) was determined in AF319 cells transiently expressing variant or WT Dcr2 (peGFP as control). Cells were infected 24 hpt. At 48 hpi, FFLuc readings were measured and presented as mean±SEM relative light units compared to eGFP control set to 1 from n = 3 independent repeats with * = p<0.05 versus control according to two-way Student’s t-test.

https://doi.org/10.1371/journal.ppat.1010202.g002

A phylogenetic tree of the putative amino acid sequences of the Dcr2 haplotypes showed a higher degree of genetic diversity between Aaeg Dcr2 of African origin compared to those of Asian or American origin (Fig 2B), although pairwise identity between haplotype genetic sequences was fairly high (>99%; S2 Fig). Three major clusters were observed, with all representative haplotypes from Cameroon and Gabon mosquito populations belonging to one clade. Interestingly, Dcr2 amino acid mutations at positions 246 and 253 appear to occur in pairs, with either an N246/Q253 or D246/K253 (Fig 2C). These combinations result in a neutral amino acid net charge except for Ben_Hap11 (N246/K253), resulting in a positive net charge. Another distinct polymorphism was in the RNase IIIB domain at position 1514 where either an H or N was observed. Other polymorphisms were mainly in regions between domains. Importantly, amino acids identified above as critical for helicase and RNase III activities were conserved across Ae. aegypti populations.

Mutations based on the observed net charge conservation at amino acid 246 and 253 were introduced to pPUb-myc-Dcr2 with the addition of a theoretical negative charge mutant (N246D) not observed among the mosquito populations surveyed (Fig 2D, top panel). Site-directed mutagenesis cloning was performed, and the expression of the myc-tagged variants in AF319 cells was verified by immunoblotting (Fig 2D, bottom panel). The predicted effect of the natural haplotype variants on protein function was also assessed by SIFT, with scores indicating tolerated substitutions (S1 Table). To test the silencing efficiency of these natural haplotype variants, the same RNAi reporter assay mentioned above was performed. The assay revealed that all except H1514N could silence FFLuc expression and were comparable to the WT (Fig 2E). The variants were also investigated for their antiviral capacity against SFV. Transient expression of the natural haplotype variants in AF319 significantly reduced SFV-FFLuc (MOI = 0.1) levels at 48 hpi versus eGFP-expressing control cells, with the antiviral activity of the variants not significantly different from WT Dcr2 (Fig 2F). Interestingly, the non-silencing H1514N variant could limit SFV replication, which could be related to how Aaeg Dcr2 processes transfected dsRNA versus viral dsRNA. Although both are considered exogenous dsRNA, it has been shown that Dmel Dcr2 performs termini selective cleavage [37,38]. Overall, these findings indicate that the functionality of the helicase and RNase III domains is important, and thus conserved.

The putative structure of Ae. aegypti Dcr2

The structure of Aaeg Dcr2 was predicted, in order to locate the conserved amino acids important for Dcr2 activity (Fig 3 and S1 Movie). Results of AlphaFold inferencing produced a model with a mean predicted local distance difference test (pLDDT) value of 81.4 (S3 Fig). This indicates a confident prediction, according to guidelines set out on EMBL’s AlphaFold Protein Structure Database [50], available at https://alphafold.ebi.ac.uk. By examining the structure and assigning domain predictions (based on GenBank ID: AAW48725 annotations and InterProScan in Geneious Prime), it can be seen that the putative structure contains two main globular modules: a “foot” module comprising of the predicted helicase domain (yellow) and a “body” domain comprising of the two forecasted RNase IIIA (blue) and IIIB (purple) domains, and the PAZ (violet) domain protruding at the far end. The two domains are bridged by a predicted domain of unknown function (DUF; green). The interface between the two modules forms a cleft, partially occluded by the projected dsRNA binding domain (dsRBD; pink). A previous study [38] yielded cryo-EM maps of moderate resolution of Dmel Dcr2. Fitting the apo-Dcr2 map (EMDB entry 7291) from this publication to our model (S3C Fig; ~50% of atoms inside contour at contour level 0.0066) illustrated that the overall architecture of Dmel Dcr2 and our model is shared (S3D Fig). The location of several conserved amino acids, for example, those in the helicase domain, on the Dcr2 surface suggests that they could be involved in interactions with other proteins, though this remains to be investigated. For completeness, the locations of haplotype variants were also located in the predicted Aaeg Dcr2 structure (Fig 3).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Locations of conserved amino acids critical for Aaeg Dcr2 function can be predicted and are distributed between surface and buried residues.

Center, AlphaFold prediction (mean pLDDT score = 81.4) for the structure of WT Aaeg Dcr2 (GenBank ID: AAW48725), colored by predicted domains from annotations and Interproscan. Residues representing conserved amino acids involved in Dcr2 activity and natural haplotype variants are highlighted in red and purple, respectively. Inset panels show enlarged representations of conserved amino acid residues targeted in mutagenesis experiments, with dashed lines showing contacts in red, hydrogen bonds in cyan, and distances in yellow. Clockwise from upper right: loss-of-function mutations in the RNase III domains line a predicted cleft in the core of the protein; G488 and K39 are forecast to be buried near the surface of the helicase domain; Y232 appears to form hydrogen bonds and extensive contacts with a neighbouring α-helix (residues 360–374) and β-strand (residues 496–502); natural haplotype variant residues N246 and Q253 (distance = 10.71Å) are likely positioned on the exposed surface of the helicase domain; H1514 is seen at the surface exposed tip of an α-helix.

https://doi.org/10.1371/journal.ppat.1010202.g003

Role of Dcr2 helicase and RNase III domains in viral dsRNA targeting and vsiRNA production

Next, we assessed whether the loss of silencing ability of the helicase or RNase III mutants were linked to changes in vsiRNA production. In mosquito cells, the Dcr2-dependent exo-siRNA response to arbovirus infection results in the production of predominantly 21 nt vsiRNAs that map to the viral genome and anti-genome [26–28]. Similar to our previous data [46], expression of WT Dcr2 in AF319 cells followed by SFV infection resulted in the production of 21 nt vsiRNAs that mapped to the SFV genome and anti-genome (Fig 4A). An obvious difference between samples from cells with WT and mutant Dcr2 was the depletion of 21 nt vsiRNAs that could be mapped to the SFV genome and anti-genome (Fig 4B). Mutations K39N and G488R in the helicase domain and mtR3 in the RNase IIIA and B domains resulted in a substantial decrease of SFV-derived vsiRNAs. In contrast, mutant Y232G showed better production of vsiRNAs compared to the other mutants but fewer in number than the WT. This may explain the silencing ability observed in the RNAi reporter assay but is insufficient for the silencing of SFV replication in our experimental set-up.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Domain loss-of-function mutations in Aaeg Dcr2 affect vsiRNA production.

(A) Distribution of 18–35 nt vsRNAs mapped to the SFV (GenBank ID: KP699763) genome (sense; magenta) or anti-genome (anti-sense; cyan) shown as mean % mapped reads in bars with replicates shown as circles from n = 2 independent repeats from AF319 cells transiently expressing mutant or WT Dcr2 and infected with SFV (MOI = 1) at 48 hpi. (B) SFV-derived 21 nt vsiRNAs from (A), mapped to SFV genome (magenta) or anti-genome (cyan) with the mean number of mapped reads shown from n = 2 independent repeats. The SFV genome organization (top panel) is shown for reference.

https://doi.org/10.1371/journal.ppat.1010202.g004

Generally, vsRNAs in the canonical vpiRNAs size range (24–30 nt) mapping to the SFV genome were abundant, as we previously observed [46], though their contribution to controlling virus replication remains unclear. These exo-siRNA pathway-independent vpiRNAs are linked to the piRNA pathway. Our vsRNA data also detected vpiRNAs that were abundant regardless of Dcr2 sequence and functionality (S4A Fig). Importantly, vpiRNAs display a nucleotide bias of an A in position 10 on the genome strand and U in position 1 of the anti-genome strand (A₁₀/U₁); with a 10 nt overlap between genome and anti-genome strands, which is due to the ping-pong amplification mechanism of the piRNA pathway [34]. Indeed, this characteristic signature was detected here in 27 nt vpiRNAs (representative for vpiRNAs; S4B Fig). Summaries of numbers and percentages of all small RNAs mapped to SFV are shown in the S4 Table.

Characterisation of 21 nt vsiRNAs with piRNA signature: vsi-piRNAs

Mutations in the helicase and RNase III domains clearly depleted vsiRNA production, with only a few 21 nt reads still detectable, predominantly in the region matching the start of the subgenomic RNA. To investigate the origins of these 21 nt vsiRNAs, their sequences were characterized based on overlap probabilities using a previously developed algorithm that maps vsRNA duplexes by generating overlapping pairs between sequences [51]. This analysis revealed a strong bias towards a 10 nt overlap for 21 nt vsiRNAs produced in cells transfected with mutant Dcr2, except in the Y232G mutant and WT Dcr2 (Fig 5A). This 10 nt bias is reminiscent of the Dcr2-independent ping-pong amplification observed in piRNAs [34,51], which is apparent in the longer vsRNAs (24–30 nt). This suggested a ping-pong biogenesis of these 21 nt vsiRNAs detected in cells with a dysfunctional Dcr2. Sequence analysis of these 21 nt vsiRNAs also showed the classic A₁₀/U₁ bias for all Aaeg Dcr2 mutants except Y232G (Figs 5B and S5) that was also observed in the 27 nt vpiRNAs (S4B Fig) mapping to the viral genome and anti-genome. These 21 nt vsiRNAs with piRNA signature are not true vsiRNAs and will be referred to as “vsi-piRNAs”, separating the 21 nt vsRNAs into two groups. Such 21 nt vsi-piRNAs may be produced in presence of WT Dcr2 or mutant Y232G (which functions similarly to WT in terms of dsRNA cleavage but less efficiently). Still, the signature signal is likely outnumbered by “true”, canonical 21 nt vsiRNAs.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Characterising 21 nt vsiRNAs with vpiRNA signature: vsi-piRNAs.

(A) Heat maps showing mean overlap probabilities of z-scores of 18–30 nt SFV-derived vsRNAs from AF319 transiently expressing mutant or WT Dcr2 from n = 2 independent repeats. vsRNA lengths are shown horizontally, and the number of nucleotide overlaps listed vertically. Red arrow labelled Dcr2 indicates expected 2 nt overlap from dsRNA cleavage with cells boxed in black. Black arrow labeled pp (ping-pong) shows expected 10 nt overlap from potential ping-pong amplification. (B) Representative sequence logos of n = 2 independent repeats (other replicate in S5 Fig) of 21 nt vsRNAs mapping to the genome (left panel) or anti-genome (right panel) of SFV from AF319 transiently expressing mutant or WT Dcr2. Overall height of the stack indicates conservation presented as bits, with the height of each symbol reflecting the relative frequency of the corresponding nucleotide at a given position. Error bars indicate an approximate Bayesian 95% credible interval.

https://doi.org/10.1371/journal.ppat.1010202.g005

Previous work in Aag2 cells had led to the identification of Ae. aegypti orthologs of the endonuclease Zucchini and exonuclease Nibbler, both involved in piRNA 3’ end formation in D. melanogaster. Zucchini-mediated piRNA 3’ end formation depends on sequence motif recognition and preferably a U downstream of the cleavage site. This results in piRNA sharp ends and a so-called +1U bias. Importantly, analysis of vpiRNAs from Aag2 cells infected with the alphavirus Sindbis (SINV) had suggested a Zucchini-like biogenesis mechanism in vpiRNA 3’ termini formation, as sharp ends and +1U bias were observed [52]. This led us to further investigate the origins of 21 nt vsi-piRNAs by subjecting vpiRNAs (25–30 nt) and 21 nt vsi-piRNAs sequences to a similar analysis. Strikingly, following SFV infection of AF319 cells, we observed +1U bias in 25–30 nt vpiRNAs but not across 21 nt vsRNAs (S6 Fig), suggesting common vpiRNA production mechanisms across alphaviruses in Ae. aegypti-derived cells, regardless of the presence of mutant or WT Dcr2.

Analysis of true 21 nt vsiRNAs revealed predominant patterns of sequence overlaps

Virus-derived dsRNA cleavage by Dmel Dcr2 has been shown to result in predominant vsiRNA sequences with a distinct 19 nt overlap between small RNA duplexes with 2 nt overhangs, presumably vsiRNA duplexes [51,53,54]. Other previous work has since led to a model for Dmel Dcr2 where the helicase domain recognizes blunt end viral dsRNA, and processive dsRNA cleavage leading to diverse vsiRNAs [55]. The distribution of helicase domain-dependent, true 21 nt vsiRNAs produced by WT Aaeg Dcr2 in AF319 cells along the SFV genome and anti-genome also suggests a processive cleavage mechanism covering the entire length of the viral RNA (Fig 4B). By using the overlapping pairs analysis mentioned earlier, the extent of overlap of 21 nt vsiRNA sequences was determined to investigate vsiRNA production patterns and phasing. A pattern was observed in the number of overlapping pairs wherein sequences of true 21 nt vsiRNAs produced by WT Dcr2 predominantly overlapped by either short stretches of 5–7 nt or longer stretches of 15–18 nt as confirmed by the z-scores of the overlaps and the corresponding probability z-scores (Figs 6 and S7A).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Viral dsRNA cleavage pattern of Aaeg Dcr2 during SFV infection.

Number of overlapping pairs (top panel) and probability z-scores (bottom panel) of 21 nt virus-derived small RNAs from AF319 transiently expressing mutant or WT Dcr2 infected with SFV. Mean of n = 2 independent repeats is presented as bars or red line (grey region indicates the range).

https://doi.org/10.1371/journal.ppat.1010202.g006

This suggests that viral dsRNA cleavage by Aaeg Dcr2 results in diverse vsiRNA sequences with a wide range of overlaps. However, 19 nt overlaps were also detected, which may denote the original vsiRNA duplexes. Whether this finding represents critical differences to observations with Dmel Dcr2 [51,53,54] requires further investigation. In addition, the overlap pattern generated by mutant Y232G showed a striking resemblance to WT Dcr2, corroborating further that this mutant can still produce true 21 nt vsiRNAs. In contrast, for the other mutant Dcr2 we identified more conspicuous patterns with 10 nt overlaps as mentioned earlier, which had low z-score probabilities for 21 nt vsRNAs with 5–7 and 15–18 nt overlap, unlike what was observed for WT and Y232G Dcr2. This further confirms that these mutant Dcr2 had lost their functionality, and presumably, only a few vsiRNAs are produced and masked by vsi-piRNAs. To verify these findings, we re-analyzed previously published small RNA data of SFV infected AF319 cells, AF319 cells transiently expressing WT Dcr2, and the parental AF5 cells that naturally express Dcr2 [46]. The re-analysis gave similar results for true 21 nt vsiRNAs (S7B Fig) following a comparable pattern of overlapping pairs (S7C Fig), suggesting a degree of consistency by which functional Dcr2 processes SFV-derived dsRNA in Ae. aegypti cells.

Functional Aaeg Dcr2 generated diverse but consistent true 21 nt vsiRNAs during SFV infection

To explore the possibility that Aaeg Dcr2 loss-of-function mutations affected population-level diversity of SFV vsRNA we conducted a differential transcript abundance analysis using edgeR (Fig 7). To determine this, we treated all the SFV-mapped small RNA reads (16–46 nt) as unique transcripts and generated count tables with normalized reads for all vsRNAs or 21 nt vsRNA subset (S7D Fig). Counts were then normalized for library size using the weighted trimmed mean of M-values (TMM) method and visualized in a multi-dimensional scaling (MDS) plot of the samples. MDS plots can be viewed as a type of unsupervised clustering where distances correspond to average log-fold changes between each pair of samples. Interpretations of the MDS plot indicated that the distances (and therefore similarity) between each sample mutant or WT Dcr2 correlate. This indicated a similar population of vsRNAs (16–46 nt) within treatments, except for SFV derived vsRNAs from cells transiently expressing WT and Y232G Dcr2. Indeed, WT and Y232G showed the most variation in vsRNA populations, suggesting greater diversity of vsRNAs compared to dysfunctional Dcr2 mutants K39N, G488R, and mtR3 (Fig 7A).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Functional Aaeg Dcr2 generates diverse but highly consistent true 21 nt vsiRNAs.

(A) Multi-dimensional scaling (MDS) plot of vsRNAs (16–46 nt) and (B) 21 nt vsRNAs from AF319 cells transiently expressing functional Dcr2 (WT or Y232G; mainly true 21 nt vsiRNAs) or mutant Dcr2 (G488R, K39N, mtR3; mainly vsi-piRNAs) per replicate with axes presented as log fold change (LogFC). (C) Heatmap showing the relationship of 21 nt vsRNAs generated from AF319 cells transiently expressing WT or mutant Dcr2 (as indicated above) during SFV infection. Hierarchical clustering of treatment replicates (columns) and 21 nt vsRNAs (rows) performed by first calculating a distance matrix using Pearson correlation distance and clustered using Ward’s minimum variance method. Each library is shown, and the scale corresponds to log-transformed TMM normalized counts per million.

https://doi.org/10.1371/journal.ppat.1010202.g007

When 21 nt vsRNAs were isolated independently, the pattern was comparable for mutant Dcr2. However, true 21 nt vsiRNAs from WT and mutant Y232G were highly similar (Fig 7B) and grouped together. The 4 samples from the WT and Y232G clustered closely across both leading log dimensions. In contrast, the remaining dysfunctional mutants had greater inter-treatment diversity and grouped separately within the MDS space. Importantly it could also be observed that the greatest dissimilarity between replicates of any sample was 21 nt vsRNAs generated by mtR3, as shown by the plot versus K39N and G488R. This inter-treatment grouping is further illustrated in the clustered heatmap (Fig 7C). One striking difference between treatments is that select pools of 21 nt vsRNAs generated by WT and Y232G were predominantly absent in the other mutants. Moreover, there are 21 nt vsRNAs that were generated across treatments but differed greatly in the number of reads. It can also be observed that the two pools of 21 nt vsRNAs were markedly different between functional and dysfunctional Aaeg Dcr2. Overall, this shows that in the presence of functional Aaeg Dcr2 (WT and Y232G), diverse pools of vsRNA populations are produced during SFV infection, though pools of true 21 nt vsiRNAs were highly consistent.

Cells

The Dcr2 knockout cell line AF319 used in this study was previously developed; it is derived from Ae. aegypti AF5 cells (a clone derived from the Aag2 cell line; provided by Kevin Maringer, The Pirbright Institute, UK) [46]. AF319 cells were grown in Leibovitz’s L-15 medium plus GlutaMax (Gibco) with 10% tryptose phosphate broth (TPB; Gibco), 10% fetal bovine serum (Gibco), and penicillin-streptomycin (pen-strep; 100 U/mL-100 μg/mL; Gibco) at 28°C [46,70]_. These cell lines are available as Aag2-AF319 (ECACC 19022602) and Aag2-AF5 (ECACC 19022601) through Public Health England. Baby hamster kidney cells (BHK-21; a commonly used cell line available at the MRC-University of Glasgow Centre for Virus Research, UK) were grown in Glasgow Minimum Essential Medium (GMEM; Gibco) supplemented with 10% TPB, 10% newborn calf serum (Gibco), and pen-strep at 37°C with 5% CO₂.

Viruses

SFV4 (abbreviated to SFV; GenBank ID: KP699763) and SFV4(3H)-FFLuc (abbreviated to SFV-FFLuc; reporter virus expressing firefly luciferase, FFLuc) have been previously described and were produced from icDNA in BHK-21 cells with virus plaque assays done as previously described [46,71].

Plasmids and mutagenesis

Site-directed mutagenesis of a previously developed pPUb-myc-Dcr2 [46] expression plasmid of a myc-tagged WT Aaeg Dcr2 (GenBank ID: AAW48725; coding sequence GenBank ID: AY713296) under the control of PUb promoter [72] was performed. In-Fusion cloning (Takara Bio) was performed following the manufacturer’s suggestions using mutagenesis primers designed with 15 bp extensions homologous to target vector ends (S5 Table). The mutations were based on previous studies on Dmel Dcr2 that disrupted its activity. The helicase domain mutations (K39N, Y232G, and G448R; Fig 1B) were part of known conserved amino acid motifs I, IV, and VI [47] across helicase sequences that have been identified in Dmel Dcr2 [35,37,38,40]. The RNase IIIA and IIB mutations (D1198A, E1341A, D1444A, and E1548A; Fig 1B) were based on conserved residues observed during domain IIIA and IIIB dimerization that have been investigated in Dmel Dcr2 and Dcr1 [37,48].

To identify the conserved residues for mutagenesis in Aaeg Dcr2, multiple sequence alignment of reference sequences of Dmel Dcr2 (GenBank ID: NP_523778), H. sapiens Dcr1 (GenBank ID: NP_085124), and Ae. albopictus Dcr2 (GenBank ID: NP_001339903) together with Aaeg Dcr2 WT was performed (Fig 1B, top panel) using Clustal Omega in Geneious Prime v.2021.1.1. For the helicase mutants, PCR with mutagenesis primers was performed on pPUb-myc-Dcr2 and allowed to recombine after In-Fusion reaction. For the RNase III mutants, the mutations were collectively introduced using a single construct (mtR3) from a synthetic dsDNA fragment (Eurofins). The pPUb-myc-Dcr2 plasmid was first linearized by inverse PCR, and the PCR amplified mtR3 construct mixed in an In-Fusion reaction (Takara Bio) to allow recombination.

For mutations based on Dcr2 haplotypes (S2 Table) identified from Ae. aegypti mosquitoes collected from various geographical origins, three mutations in the helicase domain were generated, N246D (negative net charge), Q253K (positive net charge), and N246D/Q253K (neutral net charge). One mutation in the RNase IIIB domain H1514N was also generated. These were based on identified SNPs (S3 Table) between the haplotype sequences and a reference WT Dcr2 (GenBank ID: AY713296). In making the constructs, pPUb-myc-Dcr2 was first subjected to double restriction enzyme digest using NruI and Esp3I (New England Biolabs) for the helicase mutations or SmaI and FseI (New England Biolabs) for the RNase IIIB mutant. The double restriction enzyme digest linearized the plasmid and excised the fragments of Dcr2 intended for PCR mutagenesis. The purified PCR products were then used in different combinations with the linearized backbone during In-Fusion cloning. All plasmid constructs were transformed in Stellar Competent Cells (Takara Bio) and purified using PureLink HiPure Plasmid Maxiprep kit (Thermo Fisher Scientific) following manufacturer’s protocol. All PCRs were carried out using CloneAmp HiFi PCR Premix (Takara Bio) as per manufacturer’s instructions. The mutations in all the constructs were verified through Sanger sequencing. The firefly (FFLuc; pIZ-Fluc) and Renilla luciferase (RLuc; pAcIE1-Rluc) expressing plasmids used for the RNAi reporter assay have been previously described [73]. A pPUb-myc-eGFP plasmid encoding an enhanced Green Fluorescent Protein (eGFP) was used as control [46].

dsRNA production

Primers flanked with T7 RNA polymerase promoters previously described [46,69] were used to amplify FFLuc (dsFLuc) or LacZ (dsLacZ) for in vitro transcription using MEGAscript RNAi kit (Thermo Fisher Scientific) as per manufacturer’s guidelines. PCR products were treated with Dnase I and Rnase A followed by column purification of dsRNA.

RNAi reporter assay

To determine the efficiency of Aaeg Dcr2 activity, an RNAi reporter assay was performed as previously described [46]. Briefly, AF319 cells at 2.5 x 10⁵ cells/well were seeded in 24-well plates and allowed to settle for 24 h. The cells were then co-transfected with 1 μg WT or mutant Dcr2 expression plasmids (with pPUb-myc-eGFP as control), 150 ng pIZ-Fluc, and 50 ng pAcIE1-Rluc plasmids together with 1 ng dsRNA targeting FFLuc (dsFLuc) or dsRNA control (dsLacZ) and 2 μL DharmaFECT 2 (Horizon Discovery). To verify activity of the RNAi pathway at the post-dsRNA cleavage step, the cells were co-transfected with the expression plasmids as above and 1 ng siRNA targeting FFLuc (siFLuc) or siRNA control (siHyg). Cells were lysed at 24 hpt with 1X Passive Lysis Buffer (PLB; Promega), and FFLuc/RLuc levels measured using a Dual-Luciferase Reporter Assay System (Promega) observing the manufacturer’s protocols in a luminometer (GloMax Multi-Detection System with Dual Injectors; Promega).

Virus infection

As above, 2.5 x 10⁵ AF319 cells/well in a 24-well plate were seeded, and the same amount of WT or mutant Dcr2 and eGFP (control) expression plasmids were transfected 24 h later. To test antiviral function against SFV, cells were infected at 24 hpt with SFV-FFLuc (MOI = 0.1) for 48 h at 28°C. To measure FFLuc levels from SFV-FFLuc infection, the Luciferase Assay System (Promega) was used following the manufacturer’s protocol. For small RNA sequencing samples, AF319 cells as mentioned above were transfected with 1 μg WT or mutant Dcr2 plasmid and followed by infection at 48 hpt with SFV (MOI = 1) for 48 h. Seven corresponding wells per WT or mutant Dcr2-transfected cells were pooled for total RNA extraction using TRIzol (Thermo Fisher Scientific) with glycogen as carrier, according to the manufacturer’s protocol.

Western blot analysis

Cells were lysed with 1X Bolt sample reducing agent and 1X Bolt LDS sample buffer (Invitrogen), and lysate separated on Bolt 4–12% Bis-Tris Plus gels (Invitrogen) and transferred to 0.45 μm nitrocellulose membrane (Thermo Fisher Scientific) using a Trans-Blot SD Semi-Dry Transfer Cell (BioRad). Membranes were then blocked with 5% (w/v) non-fat dry milk in Tween-PBS (PBS with 0.1% Tween 20) for at least 1 h, and proteins of interest probed with mouse anti-myc tag antibody (1:2000; Abcam) and rabbit anti-γ tubulin antibody (1:2000; Sigma-Aldrich) in 2% (w/v) non-fat dry milk in Tween-PBS for 1 h. Three 10 minute Tween-PBS washes were carried out before incubation with goat anti-mouse DyLight 800 (1:5000; Invitrogen) and goat anti-rabbit DyLight 680 (1:5000; Invitrogen) secondary antibody conjugated with a near fluorescent dye in 2% (w/v) non-fat dry milk in Tween-PBS for 1 h. Membranes were again washed three times with Tween-PBS and viewed using Odyssey CLx with images analyzed using Image Studio Lite v.5.2.5 (LI-COR Biosciences).

Sequencing of Dcr2 and analysis

Natural genetic polymorphisms in Aaeg Dcr2 gene were surveyed in early-generation colonies (<10 laboratory generations) of mosquitoes collected from Cameroon (Ben; Bénoué National Park, 2014), Gabon (Lope; Lopé National Park, 2014), Cambodia (Rat; Ratanakiri province, 2015), Cambodia (KC; Kampong Cham province, 2014), Guadeloupe (Guad; Saint François, 2015) and French Guiana (Cay; Cayenne, 2015). Before Dcr2 sequencing, mosquito colonies were reared under controlled insectary conditions (28°C, 12 h:12 h light: dark cycle, and 70% relative humidity). Total RNA was purified from 10–15 adult females per colony using NucleoSpin RNA kit (Macherey-Nagel) and reverse transcribed into cDNA with oligo(dT) using M-MLV reverse transcriptase (Invitrogen) following manufacturer’s instructions. Full-length Dcr2 cDNA was amplified by PCR using Q5 High Fidelity DNA polymerase (New England Biolabs) and custom forward and reverse primers (S5 Table). PCR products were purified with ChargeSwitch PCR Clean-Up Kit (Invitrogen), and the final concentration adjusted to 2 ng/μL following Qubit fluorometric quantification (Invitrogen). Multiplexed libraries were prepared using Nextera XT DNA Library Preparation Kit (Illumina) according to manufacturer’s instructions and sequenced on an Illumina NextSeq 500 (2x75 bp cycles, paired ends). Raw reads were demultiplexed using bcl2fastq v.2.15.0 (Illumina) and trimmed to remove Illumina adaptor sequences using Trimmomatic v.0.33 [74]. Clean reads were aligned to the VectorBase reference Dcr2 sequence, AAEL006794, to extract unphased genotypes for each sample. Phased haplotypes were reconstructed using the PHASE algorithm implemented in DnaSP v.5 [75]. Thirteen Dcr2 haplotype sequences were identified (GenBank ID: MW924847 to MW924859; S2 Table). To determine SNPS (Fig 2A and S3 Table), the 13 haplotype sequences were aligned to a WT Dcr2 sequence (GenBank ID: AY713296, previously cloned into by pPUb-myc-Dcr2) using MUSCLE [76]. Putative amino acid sequences of haplotypes were used to generate unrooted phylogenetic trees constructed using the Neighbor-Joining method following the Jukes-Cantor model with 1000 bootstrap replicates [76–78]. The consensus tree (Fig 2B) is presented with corresponding % consensus support for each node and substitution rates. All sequences were analyzed and visualized in Geneious Prime (v.2020.1). In silico analysis prediction on the effect of the different mutations on Aaeg Dcr2 function (S1 Table) was performed in SIFT [49] available at https://sift.bii.a-star.edu.sg/.

Generating a putative Aaeg Dcr2 structure

The WT sequence of Dcr2 (GenBank ID: GenBank record AAW48725) was issued as a query to DeepMind’s AlphaFold structure inference algorithm v.2.0.1 [79], using template models 1 through 5 and not restricting available templates. The resulting models were examined and in accordance with overall model quality predictions summarised in mean pLDDT value with the highest quality model selected for further analysis. Structural analyses were performed in ChimeraX [80]. Domain positions were assigned and labeled by color based on GenBank sequence annotations and predictions from InterProScan in Geneious Prime v.2021.1.1. Residues of interest were also labeled and contacts, hydrogen bonds, and distances estimated using built-in functions. Plots of error metrics (predicted local distance difference test and predicted aligned error) were produced using AlphaPickle [81].

Small RNA sequencing and analysis

RNA library preparations and sequencing reactions were carried out by GENEWIZ, LLC (South Plainfield, NJ, USA). Initial RNA and final DNA samples were quantified using Qubit 2.0 Fluorometer (Life Technologies), and RNA integrity was checked with 4200 TapeStation (Agilent Technologies). The small RNA sequencing library was prepared using Illumina TruSeq Small RNA library Prep Kit as per the manufacturer’s specifications. The final DNA library was also validated by qPCR (KAPA Biosystems). Sequencing libraries were multiplexed and clustered to a flowcell. After clustering, the flowcell was loaded to an Illumina HiSeq 4000 according to the manufacturer’s instructions. The samples were sequenced using a 2x150bp paired-end configuration, with data trimmed to 2x50bp paired-end reads. The HiSeq Control Software conducted image analysis and base calling. The generated raw sequence data (.bcl files) from Illumina HiSeq was converted to fastq files and de-multiplexed using Illumina bcl2fastq 2.17 software with one mismatch allowed for index sequence identification. The Illumina TruSeq small RNA 3’ adapter (5’-TGGAATTCTCGGGTGCCAAGG-3’) was trimmed from fastq basecalled libraries using the fastx toolkit (Galaxy Version 1.0.3; retrieved from http://hannonlab.cshl.edu/fastx_toolkit/). Only reads containing the adapters that were 16 nt in length or greater and without ambiguous nucleotides were kept for downstream analysis. Clean reads were then mapped against the SFV genome (GenBank ID: KP699763) using Bowtie2 (Galaxy Version 2.3.4.3) [82], with the sensitive mapping flag (—sensitive). Aligned reads outputted using the–al flag for downstream analysis. The trimming and mapping statistics summary is available in S4 Table. Only forward reads were used for the analysis, and data were deposited in the NCBI Biosample database (BioProject ID: PRJNA705585).

For the analysis of host RNAi responses to SFV, individual binary alignment files for small RNA sizes for SFV derived vsiRNA (21 nt), and vpiRNA (24–30 nt) were made using BAM filter (Galaxy version 0.5.9) from Bowtie2 output BAM files. For visualization of vsiRNA and vpiRNA mapping profiles, coverage statistics for each position of the SFV genome were extracted using bedtools (v.2.27.1) genome coverage tool [83] and visualized using GraphPad Prism (v.7). To obtain nucleotide position biases in vsi/vpiRNAs, extracted positive-sense and negative-sense mapped reads were converted to fasta and visualized using WebLogo3 [84]. For determination of overlapping vsi/vpiRNAs, overlapping pairs and overlap probabilities were determined from Bowtie2 BAM files using the small RNA signatures tool (v.3.2.1) hosted on the Mississippi Galaxy server (https://mississippi.sorbonne-universite.fr/) [51].

For the analysis of downstream SFV genomic nucleotide bias (+1–4) of SFV small RNA reads, 21 nt, representing true vsiRNAs or vsi-piRNAs, and 25–30 nt, representing vpiRNAs, were extracted from each genome orientation of BAM files and converted to fasta format. The downstream +1–4 nucleotides in the SFV genome of exact matches were iteratively extracted using grep (v2.20). Biases in resultant downstream +1–4 nucleotides were visualised using WebLogo3 [84].

To examine the broad diversity of SFV-derived small RNA species, we treated all unique 16–46 nt mapped reads as discrete vsRNA and produced counts tables for all extracted SFV small RNAs. Mapped fastq reads from individual libraries were converted, processed, and enumerated in bash (v.3.2). The resulting dimensions of the counts table for all 16–46 nt vsRNA species and 21 nt vsRNAs were 398,279 and 23,934, respectively. Count tables were then imported to the differential gene expression package edgeR (v.3.30.3) [85] in the R Studio (v.1.3.1073) environment and normalized for library size using the weighted TMM method in calcNormFactors() and filtered for a minimum count of 5 with filterByExpr(). To explore the relationship between samples, an MDS plot of the edgeR normalized count tables with plotMDS(top = 500) of both the SFV vsRNAs and 21 nt vsRNAs was constructed. A heatmap of the log-transformed TMM normalized 21 nt vsRNAs was produced using pheatmap (v.1.0.12; https://CRAN.R-project.org/package=pheatmap). Hierarchical clustering was performed in pheatmap by calculating a distance matrix using Pearson correlation distance ("correlation"), and clustered using Ward’s minimum variance method (“ward.d2”).

Data analyses

All experiments were performed in n = 3 independent biological repeats with 4 technical repeats except for the small RNA sequencing and immunoblot analyses which were performed with duplicate samples. Luciferase readings for RNAi reporter assay were expressed as FFLuc/RLuc ratio to account for transfection efficiency and normalized to the negative control (dsLacZ) readings set to 1, and relative light units were presented as mean±SEM. SFV-FFLuc replication measured as FFLuc readings and presented as mean±SEM light units relative to the negative control (peGFP) set to 1. Statistical significance was assessed by a two-tailed Students t-test or two-way ANOVA with Tukey’s Multiple Comparison, when appropriate at alpha = 0.05. All data analyses and visualizations were performed in GraphPad Prism (v.7).