Assessing the Potential Interactions between Cellular miRNA and Arboviral Genomic RNA in the Yellow Fever Mosquito, Aedes aegypti

Although the role of exogenous small interfering RNA (siRNA) and P-element induced wimpy testis (PIWI)-interacting RNA (piRNA) pathways in mosquito antiviral immunity is increasingly better understood, there is still little knowledge regarding the role of mosquito cellular microRNA (miRNA). Identifying direct interactions between the mosquito miRNAs and the RNA genome of arboviruses and choosing the relevant miRNA candidates to explore resulting antiviral mechanisms are critical. Here, we carried out genomic analyses to identify Aedes aegypti miRNAs that potentially interact with various lineages and genotypes of chikungunya, dengue, and Zika viruses. By using prediction tools with distinct algorithms, several miRNA binding sites were commonly found within different genotypes/and or lineages of each arbovirus. We further analyzed those miRNAs that could target more than one arbovirus, required a low energy threshold to form miRNA-viralRNA (vRNA) complexes, and predicted potential RNA structures using RNAhybrid software. We predicted miRNA candidates that might participate in regulating arboviral replication in Ae. aegypti. Even without any experimental validation, which should be done as a next step, this study can shed further light on the role of miRNA in mosquito innate immunity and targets for future studies.


Importance
The role of the small interfering RNA (siRNA) and P-element induced wimpy testis (PIWI)-interacting RNA (piRNA) pathways in controlling virus infections in mosquitoes has been extensively studied and are considered to be a major part of the antiviral innate immune response.
Although several studies have indicated that cellular microRNAs (miRNA) are involved in mosquito antiviral immunity, the miRNA-viralRNA (vRNA) direct interaction and its effect on virus replication in mosquitoes are still unclear. The cellular miRNAs of the mosquito Aedes aegypti, a vector for many arboviral diseases, may directly interact with three major arboviruses: chikungunya, dengue, and Zika viruses. By using the miRanda and TargetSpy tools (http://bioinfo5.ugr.es/srnatoolbox), several miRNAs were predicted to have potential binding sites that are common to multiple viral genotypes or lineages. Further analysis was carried out on miRNA-vRNA interactions that required a low energy threshold to form a complex. This study shows a broad picture of possible interactions between the miRNA transcriptomes were also changed in response to chikungunya virus (CHIKV), DENV, and Zika virus (ZIKV) infection in Aedes spp [29,31,[47][48][49][50].
Mosquito miRNAs are involved in viral replication by regulating host factors or direct miRNA-vRNA interaction. Some miRNAs reported in Aedes spp. mosquitoes can limit viral replication as a consequence of the down-regulation of mosquito genes. A miRNA, miR-2940 in Ae. aegypti and Aedes albopictus could be triggered by virus infection leading to a down-regulation of metalloprotease, which is essential for virus replication. This antiviral miRNA inhibits the replication of CHIKV, DENV, WNV, and even Palm Creek virus (PCV, an insect-specific flavivirus) in Ae. aegypti or Ae. Albopictus [30,48,51,52]. Notably, even though bioinformatics approaches suggested a target site of miR-2940 in the 3'UTR of WNV, it has been experimentally demonstrated that this miRNA has no significant effect on viral replication [52]. Another example is aae-miR-2b-3p in Ae. aegypti; aae-miR-2b-3p is up-regulated in response to CHIKV infection. Knocking down an ubiquitin-related modifier (URM) [53] impedes the thiolation of tRNA and consequently, limits CHIKV infection [54].
Surprisingly, only a few studies have shown that mosquito miRNAs interfere with virus replication via complementarity with viral RNA genome in mosquitoes. Two binding sites of aal-miR-252 and aal-miR-281 (from Ae. albopictus) were found in the genomic RNA of DENV at the envelope protein and 5 UTR regions, respectively [55,56]. However, aal-miR-252 can limit the viral replication, whereas aal-miR-281 facilitates DENV replication via a direct interaction at the 5 UTR. Besides, expressing synthetic miRNAs complementary to arboviral RNA genome was shown to reduce the transmission of CHIKV, DENV, and ZIKV [57,58], thus indirectly supporting the hypothesis that cellular miRNA can target viral nucleic acids and replication.
There is still a remarkable lack of information about the direct miRNA-vRNA interaction in mosquitoes. However, miRNAs are more easily identified by next-generation sequencing (NGS) techniques and are now available in databases while improved bioinformatics tools have been developed. To predict potential mosquito miRNA binding sites in the viral RNA, we compared published small RNA sequencing data from Ae. aegypti to virus sequence data of three major arboviruses, e.g., CHIKV, DENV (serotypes 1-4), and ZIKV. We predicted and analyzed the potential target sites on each virus genome to reveal practicable miRNA-vRNA interactions by combining thermodynamics and miRNA expression profiles. This approach can underpin future studies on the role of miRNAs in regulating arbovirus replication in mosquito cells. Notably, the biological meaning of the prediction results was only important after experimental validations.

Identification Strategies for miRNA and vRNA Interactions
Key human pathogenic arboviruses (flaviviruses, DENV1-4, ZIKV, and the alphavirus, CHIKV) were chosen for analyzing the relationship between miRNAs and viral genomes (vRNA). The genome sequences for each virus were collected from the virus database Virus Pathogen Resource (ViPR, www.ViPRbrc.org) [59,60], while the miRNA sequences of Ae. aegypti were retrieved from the miRNA database, miRBase, and published results of small RNA sequencing [61]. Predictions of miRNA-vRNA interactions were carried out mainly using miRanda software [62] and in coordination with TargetSpy [63] via the online tool sRNAtoolbox [64] with default settings (http://bioinfo5.ugr.es/ srnatoolbox). The consensus binding sites predicted by both software were extracted by BEDtools (version 2.25.0) [65]. Only the prediction sites shared by the two prediction algorithms were chosen for a further case study and evaluation of the affinity of each miRNA-vRNA complex. In addition, the structures of these complexes were predicted using the tool RNAhybrid [66] via BiBiServ2 (https://bibiserv2.cebitec.uni-bielefeld.de). Total of 261, 1671, 1244, 884, 164, and 157 complete viral genomes of CHIKV, DENV1-4, and ZIKV, respectively, were retrieved from ViPR as data input for miRNA-vRNA prediction; the viruses from each genotype were chosen for further analysis, as they showed the highest number of potential miRNA binding sites.

Flowchart Validation
Experiment data of a luciferase reporter of Ae. aegypti miRNA-mRNA interactions published by Zhang, et al. [67] were applied to validate the workflow for predicting miRNA binding sites we have adopted in this study. AAEL013070, AAEL006834, AAEL000577, and AAEL010015 of Ae. aegypti were shown to be regulated by aae-miR-11-3p (AAEL013070), aae-miR-275-3p (AAEL006834 and AAEL000577), and aae-miR-286b-3p (AAEL010015) in the 3 UTR [67]. Among the four most significant miRNA-mediated reductions validated by Zhang, Aksoy, Girke, Raikhel and Karginov [67], three of them could be correctly identified with our workflow. These transcripts and the miRNA database of Ae. aegypti were used as input data. Using the default setting, the 3 UTR binding sites for aae-miR-11-3p in AAEL013070 (−13.04 kcal/mol) and aae-miR-275-3p for AAEL006834 (−18.7 kcal/mol) and AAEL000577 (−15.89 kcal/mol), were commonly predicted by both miRanda and TargetSpy. Even though no consensus binding site at the 3 UTR for aae-miR-286b-3p could be found in AAEL000577 with both algorithms, a potential binding site for aae-miR-286b-3p at the 3 UTR could be detected at two distinct positions by miRanda and TargetSpy (http://bioinfo5.ugr.es/srnatoolbox).

Several Potential miRNA Binding Sites Were Predicted in Viral RNA Genomes
With our prediction flowchart, a total of 674 miRNA-vRNA interactions were predicted consensually by both algorithms for each virus (Figure 1). Among them, 93 potential binding sites could be found in CHIKV genome, 151, 130, 123, and 98 potential binding sites found in the genomes of DENV1-4, respectively, and 79 potential binding sites in the ZIKV genome.

miRNA Binding Sites in the CHIKV Genome
CHIKV belongs to the family of Togaviridae and the genus Alphavirus, with three genotypes circulating worldwide: East/Central/South African (ECSA), West African (WA), and Asian. One additional CHIKV lineage, Indian Ocean lineage (IOL), emerged in 2004-5 from the ECSA phylogroup and has spread throughout many tropical regions [68]. The IOL lineage predominates in regions where the vector Ae. albopictus is present/dominant, in part due to the selection of an Ae. albopictus-adaptive substitution in the CHIKV E1 envelope glycoprotein (E1-A226V). This substitution confers efficient infection and dissemination in Ae. albopictus for IOL of CHIKV [69,70]. The Asian and IOL/ECSA genotypes were responsible for the most recent outbreaks [71].
When comparing the binding sites previously identified among all four genotypes, the highest number of shared binding sites was between ECSA and IOL genotypes (n = 18), whereas WA and Asian genotypes were relatively independent to each other with no common binding sites (Figure 2a). Among them, six miRNA binding sites were common among Asian, ECSA, and IOL genotypes, which could be targeted by aae-miR-263a-5p, aae-miR-279-3p, aae-miR-305-5p, aae-miR-34-3p, and aae-miR-996-3p (Supplementary Table S1). And among three potential miRNA binding sites that are commonly found among WA, ECSA, and IOL genotypes, two might potentially be targeted by aae-miR-285-5p and aae-miR-989-3p in the E1 and Capsid coding regions, respectively. The other miRNA binding site shared among WA, ECSA, and IOL genotypes may potentially be targeted by aae-miR-iab-4-5p in the 3 UTR. Notably, this miRNA might have an additional binding site on the 5 UTR of ECSA and IOL genotypes. Any activity of this miRNA might be increased by multiple binding sites on the 5 -and potential binding site that could be consensually found in the minus strand of CHIKV RNA genome; for aae-miR-282-5p (Supplementary Figure), the predicted binding sites were shown in Supplementary  Table S2.
TargetSpy. Even though no consensus binding site at the 3′UTR for aae-miR-286b-3p could be found in AAEL000577 with both algorithms, a potential binding site for aae-miR-286b-3p at the 3′UTR could be detected at two distinct positions by miRanda and TargetSpy (http://bioinfo5.ugr.es/srnatoolbox).

Several Potential miRNA Binding Sites were Predicted in Viral RNA Genomes
With our prediction flowchart, a total of 674 miRNA-vRNA interactions were predicted consensually by both algorithms for each virus ( Figure 1). Among them, 93 potential binding sites could be found in CHIKV genome, 151, 130, 123, and 98 potential binding sites found in the genomes of DENV1-4, respectively, and 79 potential binding sites in the ZIKV genome.

miRNA Binding Sites in DENV Genomes
DENV has evolved independently into four serotypes from distinct sylvatic progenitors and then into several genotypes [73]. They only share 60-75% sequence similarity and need to be analyzed separately to find potential miRNA-vRNA interactions [74]. Several conserved miRNA binding sites could be found in each serotype which might participate in viral regulation, as outlined below.

DENV-1
Five genotypes of DENV-1 were selected for analysis: genotype I (AF298808), genotype II (JQ922547), genotype III (DQ285562), genotype IV (EF025110), and genotype V (JX669462). A total of Viruses 2019, 11, 540 6 of 16 28,30,34,33, and 26 potential binding sites were found in genotypes I, II, III, IV, and V, respectively ( Figure 2b). Among them, only two potential miRNA binding sites were common to all five genotypes of DENV-1, which could be targeted respectively by aae-miR-1-3p and aae-miR-282-5p on the capsid and NS3 protein coding regions (Supplementary Table S1). Furthermore, six potential miRNA binding sites were commonly found in four of the five genotypes (Supplementary Table S1). The NS5 region of genotypes I, II, III, and V could potentially be targeted by aae-miR-316-5p, aae-miR-92a-3p, and aae-miR-92b-3p. In addition, the miRNA aae-miR-316-5p might have a further potential binding site within the capsid region (genotypes I, III, IV, and V), which may increase the probability of miRNA-vRNA interactions. The only miRNA that was found with a common binding site in genotypes I, II, IV, and V within NS4B is aae-miR-11-5p (Supplementary Table S1). Two potential miRNA binding sites for aae-miR-263a-3p and aae-miR-998-3p were found conserved on the 3 UTR of genotypes I, III, and V (Supplementary Table S1).
Three other miRNA binding sites were shared between at least four genotypes. Two binding sites were common between Asian I, Asian II, Asian American, and Cosmopolitan genotypes: aae-miR-281-3p within NS2B, and aae-miR-998-3p within the 3 UTR region of DENV-2. Another miRNA, aae-miR-315-5p, was also predicted to target the NS5 region of all DENV-2 genotypes except the Asian I genotype. The only miRNA predicted to target Asian I, Asian American, American, and Cosmopolitan genotypes of DENV-2 is aae-miR-263a-3p, which has a potential binding site on the 3 UTR. In addition, the same miRNA binding site within the 3 UTR of DENV-2 was also found within the same region of DENV-1.

miRNA Targeting the ZIKV Genome
ZIKV was first isolated from Aedes africanus mosquitoes in 1948 [75] although serological evidence has shown a broader geographic distribution of human infections including North/East Africa and South/Southeast Asia [76][77][78][79][80][81][82][83]. Since the first human case reported in Nigeria in 1952 [81], only 13 cases of mild, febrile illness were reported until the outbreak in the State of Yap (Federated States of Micronesia) in 2007 [84][85][86][87], where more than 70% of the population were infected. [88] Later, cases  [89,90] The first ZIKV case reported in America was in 2015 in Bahia (Brazil) [91,92]. ZIKV caused a total of 51,473 suspected cases and more than 4300 cases of microcephaly in Brazil by March 2016 [93,94], and the virus spread to at least 33 countries or areas in the Americas [93,95]. Related to DENV, ZIKV belongs to the Flaviviridae family and the genus Flavivirus. The three distinct genotypes East Africa (EA), West Africa (WA), and Asian were likely to be originated in East Africa [96,97]. ZIKV has a conserved genome with less than 12% divergence at a nucleotide level among all virus strains, and with 99% nucleotide similarity for strains from the Americas [98]. Thus, the interactions between mosquito miRNA and viral RNA could be more relevant than between highly divergent viruses and might provide a new insight for evaluating the antiviral immunity of mosquitoes against newly emergent viruses. Three genotypes of ZIKV were selected for analysis: EA (KF268949), WA (JU955592), and Asian (KU365778). According to our results, 30, 22, and 27 potential binding sites were found in EA, WA, and Asian genotypes respectively (Figure 2f). Brazil by March 2016 [93,94], and the virus spread to at least 33 countries or areas in the Americas [93,95]. Related to DENV, ZIKV belongs to the Flaviviridae family and the genus Flavivirus. The three distinct genotypes East Africa (EA), West Africa (WA), and Asian were likely to be originated in East Africa [96,97]. ZIKV has a conserved genome with less than 12% divergence at a nucleotide level among all virus strains, and with 99% nucleotide similarity for strains from the Americas [98]. Thus, the interactions between mosquito miRNA and viral RNA could be more relevant than between highly divergent viruses and might provide a new insight for evaluating the antiviral immunity of mosquitoes against newly emergent viruses. Three genotypes of ZIKV were selected for analysis: EA (KF268949), WA (JU955592), and Asian (KU365778). According to our results, 30, 22, and 27 potential binding sites were found in EA, WA, and Asian genotypes respectively (Figure 2f).  Five potential miRNA binding sites within the 5′UTR, NS2A, NS5, and 3′UTR were common among all three genotypes (Supplementary Table 1): aae-miR-286a-3p, aae-miR-286b-3p, aae-miR-34-3p, and aae-miR-new8. Among them, aae-miR-286a-3p, aae-miR-286b-3p, and aae-miR-34-3p were predicted to target the 5′ or 3′UTR region of ZIKV. Moreover, aae-miR-286a-3p has an additional potential binding site on the NS5 region.

Selection of miRNAs and Thermodynamic Analysis
Sequence complementary is an important but not the only factor that modulates the miRNA-vRNA interactions. The thermodynamic between miRNA and vRNA is another crucial factor that has an effect on the formation of miRNA-vRNA complex. Several miRNAs were predicted to have a low minimum free energy (MFE), indicating a relatively high affinity to form miRNA-vRNA Five potential miRNA binding sites within the 5 UTR, NS2A, NS5, and 3 UTR were common among all three genotypes (Supplementary Table S1): aae-miR-286a-3p, aae-miR-286b-3p, aae-miR-34-3p, and aae-miR-new8. Among them, aae-miR-286a-3p, aae-miR-286b-3p, and aae-miR-34-3p were predicted to target the 5 or 3 UTR region of ZIKV. Moreover, aae-miR-286a-3p has an additional potential binding site on the NS5 region.

Discussion
As poikilotherms, the body temperature of mosquitoes is influenced by ambient temperature, resulting in a different miRNA thermodynamic stability profile in insects than vertebrates [100]. Using adequate algorithms that have high prediction accuracy in insects is then important for predicting the potential miRNA-vRNA interactions. Here, two algorithms, miRanda and TargetSpy, were chosen for analyzing the mosquito potential miRNA-vRNA interactions. Several potential mosquito miRNA-vRNA interactions could be consensually predicted by both algorithms. miRanda and TargetSpy are two powerful algorithms for predicting the miRNA potential binding sites on target sequences in insects. miRanda is one of the most extensively used miRNA target prediction tools and was applied originally for identifying miRNA binding sites in Drosophila and mosquitoes. The miRanda algorithm works in three phases. Firstly, the complementarity matches between the input miRNAs and RNA sequences were identified based on dynamic programming algorithm alignment. Secondly, a thermodynamic calculation was made to rule out the matches with an MFE value above the threshold. Finally, the remaining results were filtered by checking the sequence conservation with Drosophila melanogaster, Drosophila pseudoobscura, and Anopheles gambiae. In addition, miRanda can also weigh the matches between the 2 nd to the 8 th nucleotide from the 5 arm of miRNA (seed-region), to evaluate the potential for that miRNA-vRNA complex to form [62]. For a more stringent selection, the same database was filtered with the other algorithm, TargetSpy. Unlike miRanda, TargetSpy is an algorithm based on machine learning and automatic feature selection with a broad spectrum of compositional, structural, and base pairing of each miRNA to the targeted sequence [63]. TargetSpy is able to predict species-specific targets, as the miRNA-vRNA interactions are not extensively studied. Besides, even though trained on mouse data, TargetSpy is an algorithm which shows the highest accuracy among several other algorithms on experimentally proven datasets (including insects) [63]. Using a combination of miRanda and TargetSpy, we obtained consensus results which should be more robust for detecting miRNAs that are interacting with viral genomes. Moreover, a more reliable miRNA-vRNA interaction was examined by comparing these results with the structure information predicted by RNAhybrid.
Despite the amount of false positive miRNA-target interactions predicted by several algorithms due to the nature of miRNA target recognition [101], other factors like MFE that affect miRNA-target interaction were applied for results sorting. The lower the MFE, the higher the potential to form the miRNA-target complex [102,103]; in order to narrow down the possible miRNA candidates, a relatively stringent cut-off, -20 kcal/mol, was applied in this study [62]. Even though MFE is an important factor affecting miRNA-vRNA complex formation, it does not certify that the interaction will lead to functional changes.
In this study, we found eight miRNA binding sites that both have low MFE and are also relatively conserved in each genotype/lineage. In theory, these miRNAs have a higher chance to form a miRNA-vRNA complex with more than one genotype/lineage of virus, and thus might subsequently be involved in virus replication. Moreover, according to the predicted binding sites, depending on the binding regions, some of them might have an even higher chance to participate in virus replication; the predicted binding sites for aae-mir-11-5p, aae-miR-286a/b-3p, aae-miR-316-5p, and aae-miR-34-3p are located in important viral regulatory regions such as the 3 UTR, 5 UTR, or the subgenomic RNA untranslated region of the three tested viruses. Aae-mir-286a/b-3p was predicted to target all genotypes of ZIKV at the 5 UTR; although not highly expressed, this miRNA is an ovary specific miRNA and could be deposited into embryos [104]. It might be worth exploring the potential relationship between this miRNA and the vertical transmission of ZIKV. Except for DENV-4, most DENV genotypes were predicted to have a binding site for aae-miR-316-5p at the 5 UTR region; aae-miR-316-5p being a miRNA that could be detectable in the fat body and slightly up-regulated in midguts 24 h after receiving a blood meal [67,105]. However, the effectiveness of aae-miR-316-5p might be reduced due to the relatively low expression level [105]. Conserved aae-miR-34-3p binding sites could be found in the 3 UTR of ZIKV, DENV-1, and DENV-3; this miRNA was reported to be up-regulated in response to CHIKV and Wolbachia infection [53,106]. Finally, the only miRNA that has low MFE and is relatively conserved among lineages of CHIKV is aee-miR-11-5p, which was predicted to target the end of the subgenomic untranslated RNA region, which also covers the start codon of the viral structural polyprotein of CHIKV. Therefore, the expression of CHIKV structural proteins might be influenced by the formation of a miRNA-vRNA complex at the end of subgenomic RNA untranslated regions, which might antagonize the binding region with host translational factors. In addition, aae-miR-11-5p is also one of the most abundantly expressed miRNAs in Ae. Aegypti [104], which could increase the chance of complex formation. Notably, the replication of the three arboviruses relies on the duplication of viral polyproteins. Thus, a single effective miRNA-vRNA interaction in the viral genome is sufficient to regulate viral replication. In mammals, more evidence is available on the role of virus replication regulated by direct miRNA-vRNA interactions [107][108][109][110][111][112][113][114][115], whether facilitation or reduction.
The critical roles of cellular miRNAs in host viral immunity have received more attention elsewhere with mammalian host miRNA-vRNA interactions [107,108]. In humans, the interaction between the liver specific miR-122 and the 5'UTR of Hepatitis C virus (HCV) stabilizes the viral RNA facilitating viral replication [116,117]. On the contrary, the replication of the arbovirus Eastern equine encephalitis virus (EEEV) could be reduced by the human hematopoietic cell specific miR-142-3p that has four binding sites in the 3'UTR [109]. In contrast, there is still a lack of information on the miRNA targetome for Ae. aegypti. However, few studies have demonstrated that the replication of some genotypes of DENV-2 could be affected by direct miRNA-vRNA interaction in Ae. Albopictus [55,56], and CHIKV, DENV-3, and ZIKV transmission could be reduced in synthetic antiviral miRNA expressing genetically modified Ae. Aegypti [57,58], suggesting the cellular miRNA can target viral RNA, and such interactions could potentially affect viral RNA replication.
In this study, we used the consensus results from two adequate algorithms to predict the possible interactions between Ae. aegypti cellular miRNAs and three arboviruses. We further sorted the results according to the sequence conservation among genotypes or the miRNA-vRNA interactions with low MFE and provided several candidates for further investigation. Even though we could merely provide a list of candidates, it remains necessary to evaluate the role of each potential miRNA on virus replication, and the subsequent reduction or facilitation of virus replication has to be demonstrated experimentally. Inhibitor studies should shed light on individual miRNA-vRNA interactions. However, miRNA-vRNA interactions can also be proven by techniques such as Argonaute-crosslinking immunoprecipitation (AGO-CLIP), and mutagenic analysis of virus genomes is possible for many of the arboviruses investigated here.
Author Contributions: P.-S.Y. and A.-B.F. contributed to conceiving, performing the analysis and writing the paper. C.-H.C., V.S., and A.K. helped in interpreting the data. A.K. helped editing the paper. All authors reviewed the paper.