The PAZ domain of Aedes aegypti Dicer 2 is critical for accurate and high-fidelity size determination of virus-derived small interfering RNAs
- Melinda Reuter1,2,3,4,12,
- Rhys H. Parry5,12,
- Melanie McFarlane1,
- Rommel J. Gestuveo6,
- Rozeena Arif1,
- Alexander A. Khromykh5,7,
- Benjamin Brennan1,
- Margus Varjak8,
- Alfredo Castello1,
- Lars Redecke9,10,
- Esther Schnettler2,3,4 and
- Alain Kohl1,11
- 1MRC-University of Glasgow Centre for Virus Research, Glasgow G61 1QH, United Kingdom
- 2Bernhard-Nocht-Institute for Tropical Medicine, 20359 Hamburg, Germany
- 3German Centre for Infection Research (DZIF), Partner Site Hamburg-Luebeck-Borstel-Riems, 20359 Hamburg, Germany
- 4Faculty of Mathematics, Informatics and Natural Sciences, University Hamburg, 20148 Hamburg, Germany
- 5School of Chemistry and Molecular Biosciences, The University of Queensland, St. Lucia 4072, Australia
- 6Division of Biological Sciences, University of the Philippines Visayas, Miagao, 5023 Iloilo, Philippines
- 7Australian Infectious Diseases Research Centre, Global Virus Network Centre of Excellence, Brisbane 4072, Queensland, Australia
- 8Institute of Technology, University of Tartu, Tartu 50411, Estonia
- 9University of Lübeck, Institute of Biochemistry, 23562 Lübeck, Germany
- 10Deutsches Elektronen Synchrotron (DESY), Photon Science, 22607 Hamburg, Germany
- 11Centre for Neglected Tropical Diseases, Departments of Tropical Disease Biology and Vector Biology, Liverpool School of Tropical Medicine, Pembroke Place, Liverpool L3 5QA, United Kingdom
- Corresponding authors: schnettler{at}bnitm.de, alain.kohl{at}lstmed.ac.uk
-
↵12 These authors contributed equally to this work.
-
Handling editor: Britt Glaunsinger
Abstract
The exogenous siRNA (exo-siRNA) pathway is a critical RNA interference response involved in controlling arbovirus replication in mosquito cells. It is initiated by the detection of viral long double-stranded RNA (dsRNA) by the RNase III enzyme Dicer 2 (Dcr2), which is processed into predominantly 21 nt virus-derived small interfering RNAs (vsiRNAs) that are taken up by the Argonaute 2 (Ago2) protein to target viral single-stranded RNAs. The detailed understanding of Dicer structure, function and domains owes much to studies outside the context of viral infection and studies in model organisms, and as such how Dcr2 domains contribute to detecting viral dsRNA to mount antiviral responses in infected mosquito cells remains less well understood. Here, we used a Dcr2 reconstitution system in Aedes aegypti derived Dcr2 knockout (KO) cells to assess the contribution of the PIWI–Argonaute–Zwille (PAZ) domain to induction of the exo-siRNA pathway following infection with Semliki Forest virus (SFV; Togaviridae, Alphavirus). Amino acids critical for PAZ activity were identified, and loss of PAZ function affected the production of 21 nt vsiRNAs—with enrichment of 22 nt SFV-derived small RNAs observed—and silencing activity. This study establishes PAZ domain's functional contribution to Dcr2 processing of viral dsRNA to 21 nt vsiRNAs.
Keywords
INTRODUCTION
The emergence and reemergence of arboviruses pose a continuing threat to human and animal health, with mosquito-transmitted pathogens such as dengue virus (DENV; Flaviviridae, Flavivirus), Zika virus (ZIKV; Flaviviridae, Flavivirus), and chikungunya virus (CHIKV; Togaviridae, Alphavirus) as prominent examples of the threats posed (Weaver and Reisen 2010; Weaver et al. 2018; Halstead 2019; Huang et al. 2019; Pierson and Diamond 2020; Zaid et al. 2021). The mosquito Aedes aegypti is an important vector for the transmission of many human-infecting arboviruses, including DENV, ZIKV, and CHIKV. It is found across many areas of the tropics and subtropics where it is well adapted to the human environment, preferentially taking bloodmeals from human hosts, which can result in arbovirus transmission. Prevention of transmission frequently relies on chemical control measures that target vectors, though novel approaches such as Wolbachia endosymbionts and gene drive-based approaches are promising (Carvalho and Moreira 2017; Achee et al. 2019; Brady and Hay 2020; Caragata et al. 2021; Wang et al. 2021, 2022; Facchinelli et al. 2023).
Arboviruses actively replicate in their vectors, which results in the induction of immune responses. Understanding these is relevant for a better comprehension of vector processes that impact arbovirus transmission and thus targets for intervention (Alonso-Palomares et al. 2018; Lambrechts and Saleh 2019). Among host responses that control arbovirus replication in vectors, the RNA interference response—specifically the exogenous siRNA (exo-siRNA) pathway—plays a critical role. Much of our current understanding of this pathway comes from studies in Drosophila melanogaster, including in vitro studies, with critical components and mechanisms also conserved in mosquitoes thus proving the importance of the exo-siRNA pathway across insects. It is initiated following targeting of viral replication induced long double-stranded RNA (dsRNA) by the RNase III enzyme Dicer 2 (Dcr2), which cleaves the dsRNA into virus-derived small interfering RNAs (vsiRNAs). In mosquitoes and D. melanogaster, each vsiRNA strand is predominantly 21 nt in length, with a 19 nt overlap and 2 nt overhangs at the 3′end. In the next step, 21 nt vsiRNAs are taken up by Argonaute 2 (Ago2). Ago2 is part of the RNA-induced silencing complex (RISC) and retains one strand of the siRNA to target and degrade complementary viral RNAs (Bronkhorst and van Rij 2014; Blair and Olson 2015; Olson and Blair 2015; Leggewie and Schnettler 2018; Samuel et al. 2018; Swevers et al. 2018; Tikhe and Dimopoulos 2021; Prince et al. 2023). Indeed, the absence of Dcr2 in mosquito cell lines resulted in the loss of 21 nt vsiRNAs and impaired antiviral responses (Brackney et al. 2010; Scott et al. 2010; Varjak et al. 2017; Gestuveo et al. 2022), whereas Ae. aegypti mosquitoes without functional Dcr2 develop disease phenotypes, increased replication, and dissemination following arbovirus infection (Merkling et al. 2023; Samuel et al. 2023). These data show the importance of Dcr2 as an initiator of the exo-siRNA pathway, but other functional aspects of this effector protein are less understood especially in the context of viral infection and targeting of viral derived dsRNA. The domain structure of Ae. aegypti Dcr2 (Aaeg Dcr2) is similar to that of D. melanogaster Dcr2 (Dm Dcr2), with helicase–DUF–PAZ–RNase 3A–RNase 3B–dsRBD from N to C terminus of the protein (Fig. 1), though other functional domains are present in the protein, too (Paturi and Deshmukh 2021). Structural studies have revealed that Dicer proteins generally have an L shape in which these different domains are arranged (Ciechanowska et al. 2021; Paturi and Deshmukh 2021; Zapletal et al. 2023). Analyses of Dm Dcr2 have confirmed this shape (Sinha et al. 2018; Su et al. 2022; Yamaguchi et al. 2022), which was also predicted for Aaeg Dcr2 using AlphaFold (Gestuveo et al. 2022). The bottom (or base) module of the L-shaped Dcr2 structure contains the helicase domain; the core comprises RNase 3A/B and dsRBD domains, and the cap module the PAZ (PIWI–Argonaute–Zwille) domain. Biochemical and structural studies have investigated the dsRNA-binding properties of Dm Dcr2 domains, with the PAZ domain preferably binding dsRNA termini with 3′ overhangs (Sinha et al. 2018). The structural arrangements of the 3′ nt overhang binding pocket in the PAZ domain and 5′ monophosphate in the adjacent Platform domain have been investigated in detail. Importantly, interactions of dsRNA with PAZ and Platform via 3′ 2 nt overhang and 5′ monophosphate of the dsRNA allow positioning and spacing of the dsRNA for precise dicing activity and siRNA length (Su et al. 2022; Yamaguchi et al. 2022). These structural data consolidate a previous study on the importance of binding of the 5′ monophosphate in dsRNA for length fidelity in 21 nt siRNA production from long dsRNA, though not the efficiency of small RNA production per se, and mutagenesis analysis of the Dm Dcr2 PAZ domain also demonstrated the critical role of 21 nt siRNAs in silencing an inverted repeat transgene in the fly model (Kandasamy and Fukunaga 2016). Interestingly, the presence of 5′ monophosphate and two specific arginine residues in Dm Dcr2 PAZ was found to be critical for cleavage of short dsRNA (Cenik et al. 2011; Fukunaga et al. 2014), though not longer ones where (as stated above) accurate length determination of siRNAs is the critical contribution of 5′ monophosphate binding to the dicing process.
Domain features of Dcr2. (A) Schematic representation of D. melanogaster and Ae. aegypti Dcr2, including functional domains: helicase domain, domain of unknown function (DUF), PAZ domain, RNase IIIA and IIIB domains, and dsRNA-binding domain (dsRBD); aa, amino acid. (B) Multiple sequence alignment of insect Dcr2 referring to Aaeg Dcr2 AAW48725, including the selected mutations of highly conserved amino acids for potential PAZ domain loss of function. Accession numbers for sequences used are indicated in Supplemental Table S1. Alignment was produced using Benchling (https://www.benchling.com/).
A limitation of studies on Dm Dcr2 PAZ domain is their reliance on in vitro studies with dsRNA substrates or endogenous small RNAs, while the role of this domain in antiviral responses during infection when targeting viral derived dsRNA in the cell is, to our knowledge, still unclear. The importance of the Dcr2 helicase domain in antiviral responses has been previously shown (Marques et al. 2013; Donelick et al. 2020; Gestuveo et al. 2022), but there remains a critical gap because the nature of the viral dsRNA substrate is unknown—does Dcr2 detect long, short, or different types of dsRNA in infected cells, and how does recognition occur? Here, we aimed to expand our previous studies on the role of Dcr2 domains in mosquito cell antiviral responses by investigating the role of the Aaeg Dcr2 PAZ domain in mosquito antiviral vsiRNA generation and its requirement for efficient antiviral responses. For this, we used the previously described Aaeg Dcr2 activity reconstitution system in Ae. aegypti derived Dcr2 knockout (KO) cells, with which we demonstrated the importance of the helicase and RNase domains and that 21 nt vsiRNAs are generally produced by a processive production mechanism from viral dsRNA (Gestuveo et al. 2022). To understand the contribution of the PAZ domain to the antiviral activity of Aaeg Dcr2, we identified conserved amino acids in the Aaeg Dcr2 PAZ domain that were mutated to assess their impact on dsRNA processing activity. Mutations abrogated antiviral activity against a positive-sense RNA virus, Semliki Forest virus (SFV), a member of the Togaviridae family, genus Alphavirus. PAZ function inactivating mutations resulted in a loss of 21 nt vsiRNA generation, and enrichment of 22 nt SFV-derived small RNAs. These findings indicate that PAZ domain functionality is critical for the determination of the usually predominant vsiRNA length of 21 nt, but mutations did not affect processive cleavage of viral dsRNA as such.
RESULTS
Identification of amino acids in Aaeg Dcr2 PAZ domain required for antiviral activity
Analysis of the Aaeg Dcr2 sequence (InterProScan, Geneious) indicated that the PAZ domain was located from positions N824 to G941 (Fig. 1A). This was in line with NCBI domain analysis, which also located the PAZ domain between amino acids 824–941. Next, Dcr2 PAZ domain sequences of Ae. aegypti (including recently available sequences, see Gestuveo et al. [2022]), Ae. albopictus, D. melanogaster as well as Anopheles and Culex spec. were aligned to identify highly conserved amino acids in the PAZ domain and thus with potential functional activities. Several positions, as shown in Figure 1B, were fully conserved either between at least Ae. aegypti, Ae. albopictus, and D. melanogaster Dcr2 PAZ domains (amino acids D829, N872, V880, V883, Y909, and V923) or across mosquitoes (amino acid L930) and chosen for mutagenesis in different combinations, which we assumed would give the best chance of success to disrupt PAZ domain functionality. For completeness, corresponding amino acid positions in Dm Dcr2 are indicated in Figure 2A. Using NCBI domain analysis, the PAZ domain of Dm Dcr2 (1722 amino acids in length) was shown as located from amino acids 843–1004.
Mutation and expression of Aaeg Dcr2. (A) Table showing produced Aaeg Dcr2 PAZ domain mutants M1–M4, the respective introduced loss-of-function mutations, as well as corresponding amino acid residues in Dm Dcr2; aa, amino acid. (B) Assessment of (myc-tagged) Aaeg Dcr2 expression by western blot analysis. AF319 cells were transfected with either pPUb plasmids expressing WT Dcr2 or Dcr2 PAZ domain mutants M1–M4, using pPUb-myc-eGFP as control. Anti-myc and anti-α-tubulin (control) antibodies were used. Representative of three independent repeats (other repeats in Supplemental Fig. S2).
AlphaFold was used to determine the location/structure of these conserved amino acids within the structure of the Aaeg Dcr2 PAZ domain, which appeared mainly located (with the exception of L930) around its core (Supplemental Fig. S1) when compared to Dm Dcr2 (Yamaguchi et al. 2022), suggesting structural and/or functional importance.
Combinatorial mutagenesis was carried out as described for mutants M1–M4, with conserved amino acids replaced in the combinations shown as indicated with the aim of achieving disruption of PAZ domain function with a high degree of certainty (Fig. 2A). PAZ domain mutant M1–M4 and WT Dcr2 expression in AF319 Dcr2 KO cells transfected with the corresponding expression constructs was assessed by western blot, showing that all proteins were expressed (Fig. 2B; Supplemental Fig. S2).
Next, we assessed the functionality of the Aaeg Dcr2 mutants for their silencing ability and antiviral activity. WT Dcr2 displayed antiviral activity against SFV-FFLuc in this transfection-based assay (Fig. 3), as previously demonstrated (Gestuveo et al. 2022). M1–M4 Dcr2 mutants lost SFV-FFLuc antiviral activity, compared to WT Dcr2. Similarly, the functional silencing abilities of the mutant Dcr2 were investigated in a reporter plasmid-based assay, where exogenous dsRNA is provided to induce the exo-siRNA pathway against a reporter gene (Gestuveo et al. 2022). Here, M1–M4 Dcr2 mutants only retained residual (though significant compared to control) ability to silence the reporter gene. This indicated that the mutation combinations in M1–M4 all affected Dcr2 PAZ domain functionality, and thus Dcr2's ability as an effector in the exo-siRNA pathway.
Functional analysis of Aaeg Dcr2 PAZ domain mutants. (A) Antiviral activity of Aaeg Dcr2. AF319 cells were transfected with expression plasmids encoding myc-tagged WT Dcr2 or Dcr2 PAZ domain mutants M1–M4 or eGFP (negative control). These cells were infected with SFV-FFLuc (MOI = 0.1) 24 h posttransfection (hpt). FFLuc levels were measured at 48 h postinfection (hpi) and are shown as mean ± SEM relative light units compared to eGFP control set to 1 from three independent repeats, performed in technical triplicate with (*) P < 0.05 according to Student's t-test. (B) Silencing activity of mutant Aaeg Dcr2. AF319 cells were cotransfected with (i) myc-tagged pPUb plasmids expressing WT Dcr2 or Dcr2 PAZ domain mutants M1–M4 (with eGFP as control), (ii) FFLuc and RLuc (internal control) reporter plasmids, and (iii) dsRNA targeting FFLuc (FFLuc dsRNA, blue) or eGFP (eGFP dsRNA; nonsilencing control, gray). At 24 hpt, FFLuc and RLuc levels were measured and are shown as mean ± SEM relative light units (FFLuc/RLuc) normalized to dseGFP from three independent repeats, performed in technical triplicate, with (*) P < 0.05, (**) P < 0.01, (***) P < 0.001 versus controls according to two-way ANOVA.
Aaeg Dcr2 PAZ domain loss-of-function mutations: effects on 21 nt vsiRNA production
The impact of Aaeg Dcr2 PAZ domain mutations on vsiRNA production was investigated using high-throughput small RNA sequencing from AF319 cells transiently expressing Dcr2 PAZ domain mutants or WT Dcr2 and infected with SFV. Raw populations of small RNA read lengths prototypical of miRNAs (22 nt), siRNAs (21 nt), and piRNAs (24–30 nt) (Supplemental Fig. S3) were observed in all cells. Small RNA profiles from Dcr2 PAZ domain mutants M1–M4 and the eGFP control demonstrated a bias toward a read length of 22 nt rather than 21 nt observed in the WT Dcr2 expressing cells. Read mapping profiles of SFV-derived small RNAs indicated that WT Dcr2 efficiently mediated the production of SFV-derived 21 vsiRNAs (Fig. 4A), which are (as expected) (see Gestuveo et al. 2022) the predominant length in the 18–22 nt size range and aligned along the length of the SFV genome and antigenome. Almost no such vsiRNAs were detected in the eGFP control. Importantly, M1–M4 Dcr2 mutants showed reduced production of 21 nt vsiRNAs (though those remaining align along the length of the SFV genome and antigenome), and a more even distribution of read numbers across the 20–22 nt small RNA size range, compared to WT Dcr2, was observed. This indicates that Dcr2 PAZ domain mutants M1–M4 have lost the ability to precisely measure siRNA length and consequently do not predominantly produce 21 nt vsiRNAs, but have retained the ability to process viral dsRNA into small RNAs. Crucially, we examined absolute magnitudes of SFV-derived small RNAs in the sizes of 20 and 22 nt and observed a clear enrichment for 22 nt, ranging from 2.18-fold (Dcr 2 PAZ domain mutant M3) to 3.02-fold (Dcr2 PAZ domain mutant M4) (Fig. 4B). For the 20 nt size, we found modest or no significant change (0.91- to 1.45-fold) consistent with previous results (Kandasamy and Fukunaga 2016). Critically, AlphaFold 3-based modeling of Aaeg Dcr2 PAZ domain interactions with dsRNA suggested that mutations did not abrogate the ability to bind dsRNA (Supplemental Fig. S4). This consolidates the idea that conserved amino acids are required for keeping dsRNA in a precise position for accurate cleavage to occur.
Aaeg Dcr2 PAZ domain mutations affected 21 nt vsiRNA magnitude and length. (A) Small RNA sequencing of AF319 cells, transiently expressing Dcr2 PAZ domain mutants M1–M4 or WT Dcr2, or eGFP (as negative control), and infected with SFV (MOI = 1) at 48 hpi. Histogram of small RNA reads 18–30 nt in length, mapped to the SFV (SFV4 GenBank ID: KP699763) genome (positive numbers) and antigenome (negative numbers), with colors indicating first base nucleotide prevalence per size, shown as mean percentage of mapped reads (y-axis, percentage reads) from two independent experiments ±SD. Mapping of SFV-derived 21 nt vsiRNAs, mapped along the SFV genome (magenta) or antigenome (cyan) (y-axis, vsiRNA reads per million) ±SD. (B) Fold change of magnitude of SFV-derived small RNAs of lengths 20–22 nt of Dcr2 PAZ domain mutants M1–M4 compared to WT Dcr2.
Virus-derived small RNAs retain a 2 nt overhang cleavage pattern independently of Aaeg Dcr2 PAZ domain functionality
The predominant 21 nt vsiRNAs produced by Dcr2 cleavage that are commonly observed following arbovirus infection of mosquito cells are generated as duplexes. Assuming that mechanisms observed for D. melanogaster also apply in mosquitoes, they consist of two 21 nt strands that overlap by 19 nt with a 2 nt overhang at each 3′end. To investigate if and how dsRNA cleavage patterns were affected by the PAZ domain for generating virus-derived small RNAs, we analyzed strand overlapping small RNA patterns (Fig. 5A, heat maps that show mean overlap probability z-scores of SFV-derived small RNAs 18–30 nt; Fig. 5B, showing number of overlapping pairs per million mapped reads of 19–23 nt SFV-derived small RNAs). For WT Dcr2, but also Dcr2 PAZ domain mutants M1–M4, all virus-derived small RNAs in the range of 18–23 nt all showed a clear pattern, with strand complementarity significantly enhanced for overlaps that result in 2 nt overhangs. For example, 21 nt vsiRNAs overlapped with complementary sequences also by 19 nt (thus confirming observations in D. melanogaster); virus-derived small RNAs that are 20 nt in length overlapped with sequences that are 18 nt in length, etc. This showed that PAZ domain mutations affected Aaeg Dcr2's ability to produce the usually predominant 21 vsiRNA over other lengths but did not affect the cleavage mechanism itself. We did, however, note a higher frequency of precise 2 nt overhangs for virus-derived small RNAs in AF319 expressing M1–M4 Dcr2 mutants compared to WT Dcr2. Indeed, further analysis noted a marked z-score drop, for example, for 19 nt overlap (in 21 nt vsiRNA) and +2 nt overhangs for WT Dcr2 compared to Dcr2 PAZ domain mutants M1–M4 (Supplemental Fig. S5, left panel). Interestingly, when we reanalyzed 21 nt vsiRNA data previously published from SFV infection of the Ago2 KO cell line AF525, compared to parental AF5 cells (Scherer et al. 2021), we also observed higher frequencies of virus-derived small RNA duplexes with +2 nt overhangs (Supplemental Fig. S5, right panel). This leaves the possibility that failure to efficiently transfer onto Ago2, or fully process small RNA (and failure to degrade passenger strands) contributes to the increased presence of intact virus-derived small RNAs duplexes that are detected in this overlap analysis.
Aaeg Dcr2 PAZ domain mutations result in dysregulated vsiRNA duplex overlap lengths. (A) Heat maps showing mean overlap probability z-scores of 18–30 nt SFV-derived small RNAs from AF319 transiently expressing Dcr2 PAZ domain mutant M1–M4 or WT Dcr2, or eGFP (control) from two independent repeats, as initially characterized in Figure 4. Differing lengths of vsiRNAs, analyzed individually, are shown horizontally, and nucleotide overlaps are listed vertically. The red arrow labeled Dcr2 indicates the expected 2 nt overlap from dsRNA cleavage with cells boxed in black. The black arrow labeled pp (ping-pong) shows the expected 10 nt overlap from potential ping-pong amplification in the piRNA pathway. (B) Number of overlapping pairs per million mapped reads of 19–23 nt virus-derived small RNAs from AF319 transiently expressing Dcr2 PAZ domain mutants M1–M4 or WT Dcr2, or eGFP (control) infected with SFV. Data are shown as the mean of two independent repeats with the range of values.
Production of virus-derived PIWI-interacting small RNAs (vpiRNAs), which are generally longer (24–30 nt) with a 10 nt overlap of the sense and antisense piRNAs as a result of the so-called ping-pong amplification mechanism in a Dcr2 independent manner (Varjak et al. 2018), was detected across all conditions.
Loss of Aaeg Dcr2 PAZ domain functionality changes the position and magnitude of vsiRNA cleavage
We have previously shown that Aaeg Dcr2 produces a diverse but largely consistent pool of 21 nt vsiRNAs following SFV infection;
a characteristic lost in helicase domain mutants (Gestuveo et al. 2022). To assess the impact of PAZ domain mutations on virus-derived small RNA population, diversity and cleavage bias, we examined
these through two different analyses. First, we defined a position-based “harmony” metric that incorporates the strand bias
and magnitude of the vsiRNAs, and second, using a vsiRNA population approach where each unique siRNA was treated as an individual
gene, we examined differential vsiRNA expression between treatment groups. While visualization of virus-derived small RNA
coverage, as demonstrated in Figure 4, provides a reasonable indicator of the origin and magnitude of vsiRNAs, a more nuanced analysis of the consistency of Dcr2
cleavage over the SFV-derived dsRNA requires consideration of the position, magnitude, and strandedness of 21 nt vsiRNAs.
To develop a unified metric that encapsulates these attributes at each genomic position, we calculated the terminal 5′ vsiRNA
read depth at each position of the SFV genome and antigenome, and normalized the total depth to one. Subsequently, we computed
the normalized differential coverage, di, as the difference between the normalized genome (positive-sense) and antigenome (negative-sense) coverage proportions at
each position. The genome coverage proportion at position i, Gi, is defined as the positive-sense coverage at i divided by the total positive-sense coverage for all N positions. The antigenome coverage proportion, Ai, is similarly calculated as the negative-sense coverage at i normalized by the total negative-sense coverage. The formula for the normalized differential coverage at each position i is given by:
Through this measure, we obtained a ratio reflecting the balance of genome to antigenome virus-derived small RNA coverage, which provides insight into the strand-specific cleavage patterns of WT Dcr2 across the positions of the SFV genome. Using this metric, we generated a series of scatter plots (Fig. 6A) to compare the di values for each position of the SFV genome between WT Dcr2, negative control eGFP and the Dcr2 PAZ domain mutants M1–M4. Importantly, if the di value was 0 for a position for the comparison between treatments it was excluded from the analysis, keeping only nonzero di values for downstream analysis. We used linear regression analysis to assess the relationship between di values for all positions on the SFV genome from the WT Dcr2 and each mutant, with the coefficient of determination (R2) reflecting the extent to which both treatments shared related vsiRNA profile “harmony.” The regression analysis revealed varying degrees of correlation, with the lowest correlation between the eGFP control treatments compared to the WT Dcr2 (R2 = 0.115). This suggests very little positional cleavage harmony between the WT Dcr2 and eGFP control, which is to be expected. Dcr2 PAZ domain mutants M1–M4 showed a more consistent pattern of Dcr2 cleavage to WT with R2 values between 0.39 and 0.41. Finally, we compared di values for all 21 nt vsiRNA reads for all treatments and between mutants and summarized the R2 values as a heat map (Fig. 6B). The results indicated that the PAZ domain mutants M1–M4 were more closely related to each other than WT Dcr2 or eGFP control with R2 between 0.74 and 0.84, indicating that the introduced mutations resulted in more harmonious cleavage patterns across the SFV genome but different from wild type.
Mutations in the PAZ domain of Aaeg Dcr2 change the position and magnitude of cleavage over the SFV genome and virus-derived small RNA populations compared to WT Dcr2. Samples from Figure 4 were analyzed as follows. (A) Scatter plots showing normalized differential coverage (di) for each position in the SFV genome. Each panel compares WT Dcr2 with no Dcr2, an eGFP expressing Dcr2 negative treatment, or Dcr2 PAZ domain mutants M1–M4. The linear regression line is depicted within each plot, with the R2 value indicating the fit of the model. The number of data points (n) with di > 0 over the SFV genome used for the calculation is shown on each graph. (B) Heatmap displaying pairwise R2 values from the linear regressions, comparing the similarity of differential coverage profile between treatments, as indicated. (C) Multidimensional scaling (MDS) plots representing the relationships between different size classes of vsiRNA reads from WT Dcr2, eGFP control, and Dcr2 PAZ domain mutants M1–M4. Each plot corresponds to a specific size class of SFV-derived small RNA: 20, 21, 22, and 23 nt, as well as a combined plot of 20–23 nt. The axes represent the leading two log fold change (logFC) dimensions and the percentage of variation indicated.
In addition to analyzing positional cleavage patterns of 21 nt vsiRNAs (Fig. 6A,B), we examined population-level patterns and diversity of SFV-derived small RNA populations for individual lengths of 20, 21, 22, and 23 nt, as well as combined lengths of 20–23 nt by conducting a differential transcript abundance analysis using edgeR (Fig. 6C). Through the generation of a nonredundant list of unique small RNAs, we generated count tables and normalized for library size using the weighted trimmed mean of M-values (TMM) method. Through MDS plots we were able to visualize the relationships between the treatments. MDS plots are a form of dimensionality reduction; visually, the distances between points represent the dissimilarities or differences between the treatments. Specifically, these distances correspond to the leading log fold changes between the treatments, which reflect the relative expression levels of vsiRNAs across conditions. The MDS plots revealed varying degrees of similarity between the WT Dcr2 and Dcr2 PAZ domain mutant treatments across different virus-derived small RNA size classes. Notably, the 21 nt vsiRNAs displayed distinct clustering patterns, indicating that specific mutations in the Dcr2 PAZ domain contribute to the altered distribution and abundance of these vsiRNA species.
DISCUSSION
While the importance of the helicase domain has been demonstrated for both Dm and Aaeg Dcr2 antiviral responses (Marques et al. 2013; Donelick et al. 2020; Gestuveo et al. 2022), the relevance of the PAZ domain had not yet been investigated for Aaeg Dcr2 and processing of viral dsRNA in the context of arboviral replication in mosquito cells. Here, we identified conserved amino acids in Aaeg Dcr2 that inactivate PAZ domain function and antiviral activity. Critically, we observed that while inactivation did affect the production of 21 nt vsiRNA, it did not affect Aaeg Dcr2's virus-derived small RNA production, in general, but broadened the length of the produced SFV-derived small RNAs and led to an enrichment of specifically 22 nt SFV-derived small RNAs, similar to previous mutational analysis of the Dm Dcr2 PAZ domain (Kandasamy and Fukunaga 2016). Thus, key observations hold true for viral dsRNAs and Aaeg Dcr2. The specific functions that are affected by the mutations we set in Aaeg Dcr2 PAZ will require further investigation. However, structural modeling (Supplemental Fig. S1) suggests positioning mostly in the core of the PAZ domain and thus a role in the binding mechanism of the dsRNA termini. This is further confirmed by modeling Aaeg PAZ domain interactions with dsRNA (Supplemental Fig. S4), suggesting that the location of dsRNA is not affected, but most likely mutations destabilize binding and thus allow flexibility in the interaction, which in turn leads to variable small RNA lengths following cleavage, and in particular enrichment of 22 nt SFV-derived small RNAs. We cannot, of course, fully exclude effects on dsRNA binding or impaired processing that may not be captured by our analysis. Importantly, the critical PAZ domain amino acids we described here in Aaeg Dcr2 are different from those previously identified in Dm Dcr2 (Kandasamy and Fukunaga 2016). Based on our and previous data that demonstrate the importance of the Dcr2 helicase domain in the exo-siRNA pathway, we propose that the helicase binds usually blunt-ended, long viral dsRNA, and cleavage of these termini—which is likely to generate small RNAs of variable length—generates the first dsRNA terminus with 2 nt overhang and 5′ monophosphate that will be used by the PAZ domain to anchor the dsRNA strand and allow Dcr2 to produce mostly, though not exclusively, 21 nt vsiRNAs, and phase variation from occasional unprecise binding that allows the production of a diverse pool of vsiRNAs from dsRNA, may, therefore, start from early cleavage events onward. Importantly, the size of 21 nt is critical and the most efficient siRNA length for silencing to target RNAs in the fly system (Elbashir et al. 2001; Kandasamy and Fukunaga 2016). It is intriguing to note that U is enriched in the first position in libraries (Supplemental Fig. S3) compared to vsiRNAs where a bias toward C or G can be seen. Indeed, this may be due to the presence of host derived sequences, e.g., miRNAs that frequently have a U in the 5′ first base position (Mi et al. 2008; Ghildiyal et al. 2010; Warf et al. 2011). Indeed, in D. melanogaster, it has been shown that siRNAs derived from siRNA pathway (e.g., following processing long dsRNA) give rise to siRNAs with frequently C in position 1, which has been suggested to facilitate sorting into Ago2 (Ghildiyal et al. 2010). A similar mechanism may be at play here in the production of vsiRNAs. Our data also suggest that maturation of vsiRNAs—by cleaving the passenger strand—is less efficient with Dcr2 PAZ domain mutants M1–M4. While the exact processes would need further investigations, the marked increases in duplexes, (with +2 nt overhangs), similar to what is observed in Ago2 KO cells infected SFV where an accumulation of duplexes is seen (Scherer et al. 2021), might be due, for example, to a failure of vsiRNAs to efficiently transfer into Ago2 via the dysfunctional Dcr2 PAZ domain, or a processing issue. This requires further investigation, and there is a degree of “noisiness” in such assays that needs to be taken into account. Indeed, in our previous work (Gestuveo et al. 2022) for WT Dcr2, 19 nt overlaps with 2 nt overhangs, for example, in 21 nt vsiRNAs were more easily observed than in this study. This may be due to technical issues, for example, small RNA sequencing, transfection levels of Dcr2 expression plasmids, etc., and to some extent production of stable Dcr2 expressing cell lines may remediate small differences for future work.
In summary, this study demonstrates the critical role of the Aaeg Dcr2 PAZ domain in the precise production of 21 nt vsiRNAs. Its critical contribution is in accurate vsiRNA length determination, which in turn is critical for the Ae. aegypti exo-siRNA pathway to mount an antiviral response against SFV. Moreover, future work will determine how arboviruses of different families are targeted by Dcr2 and give us further insights into the nature of the viral dsRNA substrate.
MATERIALS AND METHODS
Cells
The Dcr2 KO cell line Aag2-AF319 (subsequently abbreviated to AF319) used here is derived from Ae. aegypti Aag2-AF5 cells, a clone derived from the Aag2 cell line (provided by Kevin Maringer, The Pirbright Institute, UK) (Varjak et al. 2017; Fredericks et al. 2019). As previously described, AF319 cells were cultured in Leibovitz's L-15 medium with GlutaMax (Gibco) with 10% tryptose phosphate broth (TPB; Gibco), 10% fetal bovine serum (FBS; Gibco), and 100 U/mL penicillin–100 µg/mL streptomycin (Pen-Strep; Gibco) at 28°C (Varjak et al. 2017; Fredericks et al. 2019). Cell lines are available as Aag2-AF319 (ECACC 19022602) and Aag2-AF5 (ECACC 19022601) through Public Health England. Baby hamster kidney (BHK-21) cells are a commonly used cell line available at the MRC-University of Glasgow Centre for Virus Research; these cells 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% CO2.
Viruses
Virus production from icDNA and plaque assay titrations were performed as previously described (Ulper et al. 2008; Varjak et al. 2017) to produce SFV4 (abbreviated to SFV; GenBank ID: KP699763) and SFV4(3H)-FFLuc (abbreviated to SFV-FFLuc; reporter virus expressing FFLuc).
Plasmid generation
Mutations in the PAZ domain of Aaeg Dcr2 (GenBank ID: AAW48725) were created by In-Fusion cloning of synthetically generated PAZ domain mutant fragments into the previously generated pPUb-myc-Dcr2 (Varjak et al. 2017) following the manufacturer's guidelines. Briefly, pPUb-myc-Dcr2 was linearized by PCR, removing the WT PAZ domain. Mutant PAZ domain fragments were PCR amplified using primers with extensions that were homologous to the linearized vector. In-Fusion cloning was then performed to create pPUb-myc-Dcr2 constructs containing the various mutants. Mutations were based on multiple sequence alignment, and identification of conserved amino acids in Ae. aegypti, D. melanogaster, and Ae. albopictus.
Virus infections
For investigations of antiviral activity against SFV, AF319 cells were seeded at 2 × 105 cells/well in 24 well plates and transfected with 1 µg of pPUb plasmid expressing WT or mutant Dcr2 or eGFP (as control) after 24 h using DharmaFECT 2 (Horizon Discovery) following the recommended protocol. At 24 hpt, the cells were infected with SFV-FFLuc (MOI = 0.1) and lysed at 48 hpi using 1× passive lysis buffer (PLB; Promega). FFLuc levels were measured using a Luciferase Assay System (Promega) according to the manufacturer's protocol in the GloMax Multi Detection System with dual injectors (Promega). For the production of samples for small RNA sequencing, 2.5 × 105 cells/well AF319 cells were seeded in 8 wells of a 24 well plate, transfected as above after 24 h and infected with SFV (MOI = 1) at 24 hpt for 48 h. For RNA extraction and follow-up small RNA sequencing, cells from 7 wells were combined and lysed using TRIzol (Thermo Fisher Scientific) according to the manufacturer's guidelines, with glycogen as a carrier. The remaining well was lysed in LDS protein lysis buffer and used for confirmatory western blot analysis.
dsRNA production
FFLuc (dsFFLuc) and eGFP (dseGFP) were amplified using T7 RNA polymerase promoter-flanked primers (Varjak et al. 2017; McFarlane et al. 2020) and dsRNAs produced with MEGAscript RNAi kits (Thermo Fisher Scientific) by in vitro transcription as per manufacturer's guidelines. Following in vitro transcription, products were treated with DNase I and RNase A. dsRNA was purified by column purification.
RNAi reporter assay
RNAi reporter assays were performed as previously described (Varjak et al. 2017) to assess the silencing activities of mutant Dcr2 compared to WT Dcr2. For this, 2 × 105 AF319 cells/well were seeded in 24 well plates. After 24 h, cells were cotransfected with 1 µg of pPUb-based expression plasmid, 50 ng of pGL3-PUb (Anderson et al. 2010), and 20 ng of pPUb-RLuc (Alexander et al. 2023), together with either 20 ng of dsRNA targeting FFLuc (dsFFLuc) or 20 ng of dsRNA control (dseGFP), and lysed at 24 hpt in 1× PLB (Promega). To determine FFLuc/RLuc levels, dual-luciferase assays were performed using a dual-luciferase reporter assay system (Promega), according to the manufacturer's protocol in a GloMax Multi Detection System with dual injectors (Promega).
Western blot analysis
Cell lysis was performed in a 1× Bolt sample reducing agent and 1× Bolt LDS sample buffer (Invitrogen). The lysate was separated using a 4%–12% Bis-Tris Plus gel (Invitrogen). For transfer onto a nitrocellulose membrane, a Trans-Blot SD Semi-Dry Transfer Cell (Bio-Rad) was used. Membranes were blocked for a minimum of 1 h using 5% (w/v) nonfat dry milk powder in PBS-Tween (PBS with 0.05% Tween 20, PBS-T) and then washed three times for 10 min using PBS-T. As primary antibodies, mouse anti-myc tag antibody (1:2000; Abcam) and mouse anti-α-tubulin antibody (1:2000; Sigma-Aldrich) in 5% (w/v) nonfat dry milk powder in PBS-T were used. After overnight incubation at 4°C, membranes were again washed three times in PBS-T. This was followed by incubation of 1 h with goat antimouse DyeLight 800 (1:5000; Invitrogen) secondary antibody conjugated with a near-fluorescent dye in 5% (w/v) nonfat dry milk powder in PBS-T. The membranes were again washed three times with PBS-T and a final wash with distilled water and analyzed using an Odyssey CLx and Image Studio Lite v.5.2.5 (LI-COR Biosciences).
Small RNA sequencing
To investigate small RNA production, at least 1 µg isolated total RNA per condition was sent for analysis. Small RNA sequencing was performed at BGI-Tech Solutions (DNBSeq, UMI small library, SE50), UMI being unique molecular identifiers, with a minimum of 10 million clean reads, and with final libraries resulting between 11 million and 18 million reads (Supplemental Table S2).
Small RNA analysis
Using fastp (v0.23.2) (Chen et al. 2018), adapters and low-quality reads from base-called fastq files were trimmed, retaining reads between 18 and 32 nt. The distribution of read sizes and first position bias was calculated using the script “1_fastq_histogram.sh” from the GitHub repository (rhparry/viral_sRNA_tools). Clean, trimmed reads were then mapped to the SFV genome (GenBank ID: KP699763) with Bowtie 2 (v2.4.5) (Langmead and Salzberg 2012), using the sensitive mapping FLAG (—sensitive) as specified in the script “2_mapping_vrnas.sh.” Histograms of mapped read lengths and first base pair bias were generated using samtools (v1.16.1) and “3_bam_sRNA_histogram.sh.” Output BAM files were filtered to include only the 21 nt reads. Coverage statistics across the SFV genome were calculated using “4_viral_srna_coverage.sh,” which uses the bedtools genome coverage tool (v.2.27.1) (Quinlan and Hall 2010). For calculating the harmony metric, coverage at the 5′ end (-5 FLAG) of each read was normalized over the whole genome. Normalized coverage files were processed with awk “{f[NR]=$2;r[NR]=$3;sumf+=$2;sumr+=$3} END {for (i=1;i<=NR;i++) {printf “%d %.6e\n”, i, (f[i]/sumf)-(r[i]/sumr)}}” to calculate di values for each position i of the SFV genome. Linear regression analyses for downstream analysis scatter plots of per position di values of the SFV genome were analyzed in GraphPad Prism (v10.0.2), and R2 values for all pairwise treatments were computed and presented as a heat map. Overlapping virus-derived small RNA pairs and overlap probabilities (z-score) were calculated iteratively from BAM files using the “signature.py” small RNA signatures Python script; both single read lengths from 18 to 30 nt were calculated as per the following conditions (‐‐minquery 18 ‐‐maxquery 18 ‐‐mintarget 18 ‐‐maxtarget 18 ‐‐minscope 1 ‐‐maxscope 18), collective 18–30 nt read lengths were computed as follows: (‐‐minquery 18 ‐‐maxquery 30 ‐‐mintarget 18 ‐‐maxtarget 30 ‐‐minscope 1 ‐‐maxscope 30) (Antoniewski 2014). To examine the diversity and populations of SFV-derived small RNA species, unique 20–23 nt SFV mapped reads were treated as discrete virus-derived small RNAs. Count tables for all extracted SFV-derived small RNAs from extracted fasta files from individual libraries were converted, processed, and enumerated in bash (v.3.2). The dimensions of count tables are as indicated (n) for 20–23 nt virus-derived small RNA species (n = 75,710), 20 nt (n = 13,430), 21 nt (n = 27,592), 22 nt (n = 20,536), and 23 nt (n = 14,155). Differential gene expression analysis was performed using package edgeR (v.3.30.3) (Robinson et al. 2009) in R Studio (v.1.3.1073), with library size normalization using the weighted TMM method in calcNormFactors() and filtered for a minimum count of 10 with filterByExpr(). Relationships between samples were explored using an MDS plot of normalized count tables with plotMDS(top=500) of all SFV-derived small RNA species.
Protein structure predictions
The WT sequence of Aaeg Dcr2 (GenBank ID: AAW48725) was used as a query for structure prediction. AlphaFold 2 algorithm (Jumper et al. 2021) was run in monomer mode using a local install of version 2.3.2 without any restrictions. The resulting models were examined and, in accordance with overall model quality predictions, summarized in mean pLDDT value with the highest quality model selected for further analysis. Structural analyses were performed in PyMOL Molecular Graphics System (version 4.5.0, Schrödinger, LLC). Domain positions were assigned and labeled by color based on GenBank sequence annotations, and residues mutated were highlighted.
Ae. aegypti Dcr2-dsRNA interactions were assessed using AlphaFold 3. The WT sequence of Ae. aegypti Dcr2 (GenBank ID: AAW48725) and derived mutant sequences (M1–M4) were modeled as interacting with dsRNA. The 53 bp dsRNA sequence was obtained from PDB ID: 7W0E for D. melanogaster Dcr2 in the active dicing state (Su et al. 2022). Sequences for WT or M1–M4 Ae. aegypti Dcr2 and each mutant were fed into AlphaFold 3 webserver (Abramson et al. 2024) along with the dsRNA sequence to predict protein structures and their interaction with dsRNA. Molecular graphics and analysis of structures were performed using UCSF ChimeraX (Pettersen et al. 2021).
Data analyses
Virus-derived small RNA metrics, coverage, and overlapping pairs analysis were visualized using GraphPad Prism (v10.0.2).
DATA DEPOSITION
Data underlying figures are available under http://dx.doi.org/10.5525/gla.researchdata.1479. Small RNA sequencing data generated for this study have been deposited in the NCBI Sequence Read Archive (SRA), available under accession number PRJNA1014654. Scripts used for small RNA analysis are available from the viral_sRNA_tools GitHub repository at: https://github.com/rhparry/viral_sRNA_tools.
SUPPLEMENTAL MATERIAL
Supplemental material is available for this article.
ACKNOWLEDGMENTS
For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) license to any Author Accepted Manuscript version arising from this submission. This study was supported by the UK Medical Research Council (MC_UU_12014/8 and MC_UU_00034/4) (A.K.) and (MR/R021562/1 and MC_UU_00034/2) (A.C.), European Research Council (ERC) Consolidator Grant “vRNP-capture” N# 101001634 (A.C.), German Centre for Infection Research (DZIF) (TTU 01.708) (E.S.), DFG (project 497659464) (M.R.), Wellcome Trust/Royal Society Sir Henry Dale Fellowship (210462/Z/18/Z) (B.B.), and Overseas Scholarship Scheme for PhD by Higher Education Commission of Pakistan (R.A.). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Author contributions: Conceptualization: M.R., R.H.P., M.M., R.J.G., A. C., M.V., L.R., E.S., and A.K.; Methodology: M.M., M.R., R.H.P., R.A., R.J.G., and L.R.; Formal analysis: M.R., R.H.P., M.M., R.A., and L.R.; Investigation: M.R., R.H.P., M.M., R.A., and L.R.; Validation: M.R., R.H.P., M.M., R.A., and L.R.; Resources: A.A.K., B.B., A.C., L.R., E.S., and A.K.; Writing—original draft preparation: M.R., R.H.P., M.M., R.J.G., R.A., A.C., L.R., E.S., and A.K.; Writing—review and editing: M.R., R.H.P., M.M., R.J.G., R.A., A.A.K., B.B., M.V, A.C., L.R., E.S., and A.K.; Visualization: M.R., R.H.P, M.M., R.A., and L.R.; Data curation: M.R., R.H.P., M.M., E.S., and A.K.; Supervision: M.M., A.A.K., B.B., A.C., E.S., and A.K.; Project administration: M.M., E.S, and A.K.; Funding acquisition: R.A., B.B., A.C., E.S., and A.K.
Footnotes
-
Article is online at http://www.rnajournal.org/cgi/doi/10.1261/rna.080149.124.
-
Freely available online through the RNA Open Access option.
- Received June 20, 2024.
- Accepted January 20, 2025.
This article, published in RNA, is available under a Creative Commons License (Attribution 4.0 International), as described at http://creativecommons.org/licenses/by/4.0/.
















