# Profile of Small RNAs, vDNA Forms and Viral Integrations in Late Chikungunya Virus Infection of Aedes albopictus Mosquitoes

Michele Marconcini, Elisa Pischedda, Vincent Houé, Umberto Palatini, Nabor

Lozada-Chávez, Davide Sogliani, Anna-Bella Failloux, Mariangela Bonizzoni

► To cite this version:

Michele Marconcini, Elisa Pischedda, Vincent Houé, Umberto Palatini, Nabor Lozada-Chávez, et al. Profile of Small RNAs, vDNA Forms and Viral Integrations in Late Chikungunya Virus Infection of Aedes albopictus Mosquitoes. Viruses, 2021, 13 (4), pp.553. 10.3390/v13040553. pasteur-03323240

HAL Id: pasteur-03323240

https://pasteur.hal.science/pasteur-03323240

Submitted on 20 Aug 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Profile of Small RNAs, vDNA Forms and Viral Integrations in Late Chikungunya Virus Infection of *Aedes albopictus* Mosquitoes

Michele Marconcini<fnref n="1" />, Elisa Pischedda<fnref n="1" />, Vincent Houé<fnref n="2" />, Umberto Palatini<fnref n="1" />, Nabor Lozada-Chávez<fnref n="1" />, Davide Sogliani<fnref n="1" />, Anna-Bella Failloux<fnref n="2" /> and Mariangela Bonizzoni <fnref n="1" />,*

<figure><img src="image_2.png" /></figure>

check for
updates

**Citation:** Marconcini, M.; Pischedda, E.; Houé, V.; Palatini, U.; Lozada-Chávez, N.; Sogliani, D.; Failloux, A.-B.; Bonizzoni, M. Profile of Small RNAs, vDNA Forms and Viral Integrations in Late Chikungunya Virus Infection of *Aedes albopictus* Mosquitoes. *Viruses* 2021, 13, 553. https://doi.org/10.3390/v13040553

Academic Editor: Roy A. Hall

Received: 3 February 2021

Accepted: 23 March 2021

Published: 25 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_3.png" /></figure>

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

**Abstract:** The Asian tiger mosquito *Aedes albopictus* is contributing to the (re)-emergence of Chikungunya virus (CHIKV). To gain insights into the molecular underpinning of viral persistence, which renders a mosquito a life-long vector, we coupled small RNA and whole genome sequencing approaches on carcasses and ovaries of mosquitoes sampled 14 days post CHIKV infection and investigated the profile of small RNAs and the presence of vDNA fragments. Since *Aedes* genomes harbor nonretroviral Endogenous Viral Elements (nrEVEs) which confers tolerance to cognate viral infections in ovaries, we also tested whether nrEVEs are formed after CHIKV infection. We show that while small interfering (si)RNAs are evenly distributed along the full viral genome, PIWI-interacting (pi)RNAs mostly arise from a ~1000 bp window, from which a unique vDNA fragment is identified. CHIKV infection does not result in the formation of new nrEVEs, but piRNAs derived from existing nrEVEs correlate with differential expression of an endogenous transcript. These results demonstrate that all three RNAi pathways contribute to the homeostasis during the late stage of CHIKV infection, but in different ways, ranging from directly targeting the viral sequence to regulating the expression of mosquito transcripts and expand the role of nrEVEs beyond immunity against cognate viruses.

**Keywords:** chikungunya virus; *Aedes albopictus*; RNAi; viral integration; vDNA

## 1. Introduction

Chikungunya virus (CHIKV) was first isolated in Tanzania in 1952 [1] causing sporadic outbreaks. In 2005, it re-emerged in costal Kenya and spread to the islands of the Indian Ocean and India [2]. After 2005, CHIKV moved globally, causing outbreaks not only in tropical regions of the world such as South-East Asia, Central Africa, and Central America, but also in temperate Europe and North America [3–5]. The increasing public heath significance of CHIKV is tightly linked with the global invasion of Asian tiger mosquito *Aedes albopictus* which established in temperate areas of the world [6]. In response to this trend, chikungunya was included in the WHO list of the most relevant neglected tropical diseases in 2017 [7].

Adult female mosquitoes acquire viral pathogens through an infectious blood meal. Ingested viral particles must replicate in the midgut and further disseminate throughout the mosquito general cavity to reach salivary glands. Only if/when infection is established in the salivary glands, mosquitoes can transmit viruses to a new host during a subsequent blood meal, thus continuing the cycle [8,9]. In *Ae. albopictus*, CHIKV infection progresses rapidly and viral particles are detected in the saliva as early as 2 days post infection (dpi) [10]. During the early stages of infection, viral titer increases and mosquitoes activate

<footnote n="1">¹ Department of Biology and Biotechnology, University of Pavia, via Ferrata, 27100 Pavia, Italy; michele.marconcini01@universitadipavia.it (M.M.); elisa.pischedda01@universitadipavia.it (E.P.); umberto.palatini01@universitadipavia.it (U.P.); alejandro.chavez@unipv.it (N.L.-C.); davide.sogliani@gmail.com (D.S.)</footnote>

<footnote n="2">² Arbovirus and Insect Vectors Unit, Department of Virology, Institut Pasteur, 25-28 Rue du Dr Roux, 75015 Paris, France; vincent.houe@pasteur.fr (V.H.); anna-bella.failloux@pasteur.fr (A.-B.F.)</footnote>

<footnote>* Correspondence: m.bonizzoni@unipv.it; Tel.: +39-0382-986324</footnote>
a complex immunity response [11–13]. Rarely, viruses are blocked in the midgut (i.e., midgut infection barrier) or cannot disseminate from the midgut (midgut escape barrier) or the salivary glands (salivary glands barrier) [14,15]. These infection-refractory phenotypes are both viral species and mosquito strain/species specific [16–19]. Most frequently, a balance between viral replication and the mosquito immune system is reached, resulting in the establishment of a non-pathogenic persistent infection which ensures mosquitoes a life-long viral transmission capacity [20,21]. Thus, in the arboviral transmission cycle, there are at least two critical moments: the establishment of viral infection in the mosquito midgut and the shift from acute to persistent infection which renders a mosquito a life-long vector. The understanding of the molecular mechanisms underlying the ability of mosquitoes to acquire, maintain, and transmit arboviruses is expected to aid in development of novel genetic-based control strategies. With the aim of identifying strategies and/or effectors to control viral infection in mosquitoes, most studies focused on biological (i.e., expressional and/or metabolic) changes occurring in mosquitoes in the early time points post infection [22,23]. However, it was recently proposed that unravelling the molecular underpinning of viral persistence could facilitate the design of population replacement strategies of mosquito control which avoid emergence of viral resistance [24].

Both during the initial acute phase of infection and in persistent infection, the cornerstone of mosquito antiviral defense is RNA interference (RNAi) [25–27]. RNAi relies on small single-stranded RNA molecules of 20–30 nucleotides (nt) that recognize and act on target sequences based on sequence-complementarity and through interaction with Argonaute proteins [27,28]. Three classes of small RNAs are known, namely small interfering RNAs (siRNAs), microRNAs (miRNAs), and PIWI-interacting RNAs (piRNAs). These small RNAs differ in length, interacting proteins, mechanisms of biogenesis, targets, and functions ranging from direct target cleavage by Argonaute-2 guided by siRNAs to destabilization of target mRNAs or translation blockage for miRNAs or silencing of transposable elements at the transcriptional or translational levels for piRNAs [28]. Additionally, miRNA-target interaction may also lead to positive regulation of the target by increasing mRNA stability or action on promoter sequences [29]. Despite these differences, increasing experimental evidence reveals crosstalk among RNAi pathways, firstly between the siRNA and piRNA pathways, after arboviral infections of mosquitoes [30,31]. One of the links between the two pathways are viral DNA (vDNAs) fragments, which appear hours after arboviral infections and serve to amplify siRNA-mediated silencing [30,32,33]. In *Drosophila*, vDNA forms are produced by defective viral genomes through the helicase domain of Dicer-2, which also mediates siRNA biogenesis [33]. In *Ae. albopictus*, vDNA forms corresponding to the region of the CHIKV genome encompassing the nsP2 and nsP4 coding sequences have been detected up to 9 days post infection [30]. Moreover, in *Aedes* spp. mosquitoes, DNA fragments from insect-specific virus (ISVs), phylogenetically linked to arboviruses of the *Flaviviridae* and *Rhabdoviridae* families, are integrated into the genome [34–36]. These nonretroviral Endogenous Viral Elements (nrEVEs) are enriched in piRNA clusters, genomic regions from which piRNAs are generated, and produce piRNAs [34–36]. Selected nrEVEs of *Aedes aegypti* limit infection by cognate viruses through EVE-derived piRNAs [37,38], but the genome of neither *Ae. aegypti* nor *Ae. albopictus* has nrEVEs derived from CHIKV [35,36]. Additionally, the landscape of nrEVEs is variable not only across mosquitoes of different geographic populations, but also within laboratory-reared strains [39,40].

The analyses of the profile of small RNA molecules after arboviral infection not only showed differences among mosquito species and viral families [41], but also revealed that small RNA molecules may be produced from infecting viruses, including viral piRNAs (vpiRNAs) and, possibly, miRNAs deriving from pri-miRNAs stem loops which are processed directly in the cytoplasm [29,42–45]. Albeit still debated, arbovirus-encoded miRNA-like molecules appear to either regulate expression of mosquito genes (i.e., the transcription factor GATA4 by the West Nile virus-derived KUN-miR-1) or directly virus
genes, leading to regulation of virus replication [46–48]. The role of vpiRNAs following arboviral infection has not been elucidated yet [45].

Besides changes in the abundance of small RNA molecules following arboviral infections, transcriptional changes occur following infections with CHIKV and dengue virus (DENV) [49].

Here, we combined the analysis of high throughput sequencing data (whole genome sequencing, WGS, and small RNA-seq) and molecular biology techniques to test for vDNA presence and new nrEVE formation and analyze the production of small RNAs, including those from nrEVEs existing in the *Ae. albopictus* genome, 14 days post CHIKV infection, when infection has become persistent [50]. We show that CHIKV infection does not result in new nrEVE formation, but the differential accumulation of piRNAs derived from existing nrEVEs correlates with differential expression of an endogenous transcript. The profile of siRNAs was even along the CHIKV genome. On the contrary, piRNAs mostly arise from a ~1000 bp window, which corresponds to the region encoding the C and E2 proteins. A unique vDNA fragment matching this region was also identified.

Statistically significant differentially accumulated siRNAs were also identified and correlated with differential abundance by qRT-PCR of selected, bioinformatically-deduced, target transcripts. Our results clearly show that RNAi pathways not only target the viral sequence, but host miRNAs and piRNAs are also modulated and contribute to cellular homeostasis during persistent CHIKV infection. The identified correlation between the differential accumulation of nrEVE-piRNAs and the decreased expression of the nrEVE-piRNA target expands the possible role of nrEVEs beyond immunity against cognate viruses.

## 2. Materials and Methods

### 2.1. *Mosquito Infections*

*Aedes albopictus* mosquitoes of the Foshan strain were used in this study. Mosquitoes were reared under constant conditions, at 28 °C and 70–80% relative humidity with a 12/12 h light/dark cycle. One-week old females were infected with CHIKV 06.21 as previously described [51]. CHIKV 06.21 was isolated in La Reunion island and was kindly provided by the French National Reference Center for Arboviruses. Thirty CHIKV-infected mosquitoes were collected 14 days post infection (dpi), along with un-infected controls. Ovaries were dissected from carcasses and the virus was titrated by fluorescence focus assay on C6/36 *Ae. albopictus* cells as previously described [51]. At 14 dpi, mean viral titer was 1515 ± 628 FFU/mL. Dissected ovaries and carcasses were pooled into groups of 15 giving a total of 8 samples, 2 pools for each condition (Figure S1A). Each sample was used for DNA and RNA extraction using the Nucleospin tissue or the Nucleospin miRNA kit (Macherey-Nagel), respectively, following manufacturer's instructions. A second infection experiment was performed as before, resulting in mean viral titer of 1589 ± 50 FFU/mL. Samples from this infection experiment were collected as before and were used for RNA extraction using the standard Trizol protocol.

### 2.2. *Whole Genome and Small RNA Sequencing*

After quality control, DNA of the first infection experiment was sent to the Beijing Genomics Institute (BGI) for Illumina library preparation. Libraries were run on the Illumina Novaseq 6000 by to obtain 150 bp paired end reads. Small RNA libraries were prepared by BGI from the same samples as for WGS and sequenced on the BGI-Seq50 to obtain 40 million reads per sample.

### 2.3. *Detection of vDNA Fragments*

The presence of vDNA fragments was investigated in CHIKV infected ovaries and carcasses by PCR using primers designed to cover the whole viral genome and producing overlapping amplicons of ~300–400 bp each (Table S1). PCR-amplified bands were cloned
using the TA cloning kit (Invitrogen) following manufacturer's instruction and sequenced
by Macrogen Europe (Amsterdam, The Netherlands).

## 2.4. Bioinformatic Analyses of WGS Data

For all genomic analyses regarding the *Ae. albopictus* genome, we used the latest genome assembly (Aalbo_primary.1, Refseq assembly GCF_006496715.1) [34]. WGS data provided by BGI were aligned to the genome of *Ae. albopictus* using BWA-MEM algorithm [52]. Because the landscape of nrEVEs in the genome of *Ae. albopictus* is not fixed, but mosquitoes may have different nrEVE patterns [40], we first assessed which nrEVEs were present in each sample using Genome Analysis ToolKit (GATK) [53] as previously described [40,53] and then assessed an unbiased small RNA profile of nrEVEs after CHIKV infection. Due to sequence similarity among nrEVEs, a nrEVE was considered absent in a sample when the average coverage for that nrEVE was less than 1 read with zero mapping quality. To identify possible viral integrations formed after CHIKV infection, we used the VyPer and ViR pipelines [54,55]. The similarity of nrEVEs sequences to CHIKV 06.21 was calculated by grouping nrEVEs depending on their viral family and then aligning them to the genome of CHIKV 06.21 (NC_004162.2) with MAFFT [56]. We retrieved the identity matrix as implemented by the EMBL-EBI search and sequence analysis tools [57] and plotted the results with custom scripts.

## 2.5. Bioinformatic Analyses of sRNA Data

For the following analyses, sRNA data were aligned using sRNAmapper [58] to either the genome of CHIKV 06.21 (NC_004162.2) or *Ae. albopictus* transcripts (Aalbo_primary.1, Refseq assembly GCF_006496715.1) using the "-best" option. Sequences between 20-22 nt were identified as siRNAs/miRNAs and sequences between 25-30 nt were identified as piRNAs. This classification of sRNAs was based on available literature [58-62], and through bioinformatics analysis as follows. First, we plotted the length distribution of all mature miRNAs reported for *Ae. aegypti*, downloaded from miRBase database v.22 [63], and corroborate that the major range for miRNAs occurs between 20 to 22 nt (Figure S2). Second, we discarded ranges of nts where miRNA and piRNA definitions overlap, as the tail distribution of miRNAs (size range of 23-24 nt) shows (Figure S2). Finally, taking into account that smaller peaks at 20 and 22 nt are often observed as Argonaute-2 products [64], and that thousands of 20 and 22 nt unique reads perfectly mapped to the viral genome (Figure 1), we decided that the broader range of 20-22 nt would provide less biased results than limiting the analyses to 21 nt reads (i.e., exo-siRNA consensus size).

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Size distribution of small RNAs reads mapping to the *Aedes albopictus* genome from ovary and carcass samples. Reads from chikungunya (CHIKV)-infected and uninfected samples are in blue and red, respectively.</figcaption></figure>

### 2.5.1. Small RNA Profile of the CHIKV Genome

sRNA data were mapped against the genome of CHIKV 06.21 (NC_004162.2). To visualize the coverage along the virus genome, alignment (bam) files of biological replicates were merged using bamCompare as implemented in DeepTools2 [59], normalizing the data to reads per kilobase of transcripts (RPKM) values based on total library size, setting bin size to 1 nt to improve resolution, and computing the mean of the number of reads
between the files. Alignment results were exported as bedgraph files and visualized in Integrative Genomics Viewer (IGV) [60]. sRNAs were filtered by length using BBMap reformat.sh [65]. PingPongPro and the small RNA signature tool [61] were implemented in Galaxy [66] to check for ping-pong signature of reads identified as piRNAs; 1U and 10A biases were plotted with weblogo 2.8.2 [64]. siRNAs that had mapped to CHIKV (hereafter called viral-sRNA) were re-mapped to *Ae. albopictus* transcripts and visualized in IGV [60], after merging alignment (bam) files with samtools merge [67]. Candidate sequences of viral miRNA-like precursors were identified using ViralMir [68] and miRNAfold [69] with default parameters by analyzing the secondary structure of the 150 bp viral RNA sequence including candidate viral miRNA-like sequences.

### 2.5.2. Small RNA Profile of Mosquito Transcripts, piRNA Clusters and nrEVEs

sRNA data were aligned to *Ae. albopictus* transcripts. Only sRNAs mapping in reverse orientation with respect to the predicted transcripts were retained. sRNAs were filtered by length using BBMap reformat.sh [65] as indicated above. Counts for each transcript were performed using the subread featureCounts.sh [70]. Transcripts displaying differential abundance of small RNAs between CHIKV infected and control samples were identified and checked against the list of *Ae. albopictus* immunity genes [34] to identify possible regulation of immunity pathways by sRNAs. sRNA data were aligned to the *Ae. albopictus* genome following the same criteria as above and their abundance in each piRNA-cluster was estimated using the subread featureCounts.sh [70]. Comparison between infected vs. non-infected samples in the abundance of sRNA mapping to transcripts (hereafter called T-sRNA) was established based on read counts in edgeR using the weighted trimmed mean of M-values with zero pairing (TMMwzp) normalization method and the likelihood ratio test to assess statistical significance [71]. As implemented in edgeR, genomic regions, or transcripts that had a minimum total read count of 15 across all samples were discarded [72]. Statistical significance was based on a 0.05 FDR threshold. Three normalization strategies were implemented to minimize bias in analysis of differential abundance. Normalizations were based on: (1) the total number of sRNA reads that mapped on the reference genome, (2) the total number of reads analyzed in edgeR, following size-filtering and counts, (3) the fraction of miRNA sized molecules that mapped onto a list of previously identified miRNAs [34], as performed in [73–75]. Results from the three normalization methods were compared and visualized in Venn diagrams using Venny 2.1.0 [76] showing that the normalization based on the total number of sRNA reads has fewer outliers than the other methods and fewer outliers (Figure S3). Transcripts displaying significant differential abundance of small RNAs between infected or uninfected were analyzed using Blast2GO [77]. Deduplicated target transcripts were blasted using default parameters against all currently available Culicidae genome sequences and the top 5 hits were retrieved. GO analyses were performed merging the results from the Interpro classification and the results from the blast database with the subsequent mapping and annotation, as implemented in Blast2GO. GO enrichment was performed based on the annotation of the full sets of transcripts of *Ae. albopictus* with FDR cut-off of 0.05. To find the most represented pathways, target transcripts were analyzed in KOBAS 3.0 [78], selecting *Ae. aegypti* as reference species and searching the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway (K option). Results were considered significant at *p*-value < 0.05.

Small RNA coverage of nrEVEs was tested after having verified nrEVE presence in the sample under study as described above. Despite the high stringency implemented in sRNAmapper (18 nt perfect seed match, 1 mismatch allowed after the seed and 2 mismatches allowed on the 3' end tail of the read), single sRNA molecules could be mapped to several nrEVEs with equal confidence because of the high sequence similarity among nrEVEs [34]. To identify sRNAs mapping to nrEVEs and accumulate differentially across samples, we used two strategies to reduce any bias due to nrEVE sequence similarity: (1) we focused on single small RNAs (level 1) or (2) we took in consideration the overall small RNAs mapping to a nrEVE (level 2) (Figure S1B). For the first strategy, reads from
each condition that mapped to nrEVEs were collapsed based on their sequence identity with fastx-toolkit [79] and each sequence was used to count the amount of each piRNA and siRNA sequence present in every fastq file, using custom scripts. Individual sRNA names were given following fastx-toolkit results in the format "rank—number" where rank starts at 1 for sRNA with the highest read count and number is the actual number of reads. For the second strategy, piRNAs and siRNAs abundance from each nrEVE was estimated using subread featureCounts.sh [70]. Differential abundance of small RNAs mapping to nrEVEs (hereafter defined as EVE-sRNA) between infected vs. uninfected samples was assessed in edgeR as described in the previous paragraph. EVE-sRNAs that were found statistically differentially abundant between infected and not infected samples were used to blast (BLASTn, default settings) [80] the collection of *Ae. albopictus* transcripts [34]. When found, the top 10 blast hits of differentially-abundant EVE-sRNAs were retrieved and thermodynamic analyses on the RNA-RNA pairing were performed using the software IntaRNA locally with heuristic mode under the model B [81,82], setting the seed between positions 2–7 [83] and the seed energy threshold < −3.2 kcal/mol [81]. Candidate target transcripts were analyzed using Blast2GO [77] to find their Gene Ontology terms as previously described. Functional annotation was also performed in Argot 2.5 [84] searching against Culicidae with default parameters and a cut-off ≥200.

## 2.6. *Differential Expression of Ae. albopictus Annotated miRNAs*

sRNA data were mapped to the consensus sequences of 229 loci corresponding to 121 validated *Ae. albopictus* pre-miRNAs [34]. Counting and differential expression analysis were performed in edgeR as described above. Results of identical miRNAs from different loci were averaged. For each differentially expressed miRNA, both the consensus star sequence and the consensus mature sequence were retrieved, and the possible targets investigated using blast and IntaRNA as described for EVE-sRNAs. Reads that mapped on differentially expressed miRNAs were also aligned to the CHIKV genome to investigate their potential as antiviral miRNAs.

## 2.7. *Quantitative PCR*

Selected target transcripts from all previous analyses (i.e., from the list of differentially abundant T- sRNAs and transcripts targeted by differentially abundant EVE-sRNAs or viral-sRNAs) were tested for differential expression by quantitative PCR (qPCR) on samples from the second infection experiment. RNA from these samples was extracted using standard Trizol protocol and reverse transcribed with the GoScript reverse transcriptase kit (Promega), following manufacturer's instructions. qRT-PCR reactions were performed using the QuantiNova SYBR Green PCR Kit (Qiagen) following the manufacturer's instructions on an Eppendorf Mastercycler RealPlex4 (Table S1). Estimates of relative quantification were performed with the delta-delta-Ct method implemented in the software qBase+ (Biogazelle), using RPL34 gene as housekeeping [85]. Statistical significance was assessed by ANOVA and unpaired t-tests as implemented in qBASE+. Selected piRNAs were also investigated using qPCR (Table S1). Small RNA molecules from the same samples used for transcript quantification were reverse transcribed using the miScript II RT Kit (Qiagen) following manufacturer's instruction with the HiSpec buffer to ensure only mature sRNAs would be selected. qPCR reactions were set up using the miScript SYBR Green PCR Kit (Qiagen), following manufacturer's instruction. piRNA expression was normalized to the expression of U6B small nuclear RNA (RNU6B).

# 3. Results and Discussion

*Aedes albopictus* mosquitoes of the reference Foshan strain were infected with CHIKV. At 14 dpi, when persistence has established, saliva and two pools of each 15 female carcasses and ovaries were collected. Pools of carcasses and ovaries were processed to extract both DNA and RNA for WGS and small RNA-seq, respectively. WGS data were used to derive the pattern of nrEVEs of the sampled mosquitoes and verify whether CHIKV infection
results in new nrEVEs formation. Small RNA-seq was used to quantify the changes in small RNAs against CHIKV, the *Ae. albopictus* transcriptome and nrEVEs (Figure S1A).

WGS libraires resulted in a total of 1,324,733,838 and 1,713,328,214 clean reads in CHIKV-infected carcasses and ovaries, respectively, and 1,467,049,196 and 774,603,540 clean reads in uninfected carcasses and ovaries, respectively, after grouping biological replicates. Sequencing of small-RNA libraries resulted in a total of 1,551,205,753 and 1,010,736,055 reads, in the range of 18–30 nt, which could be mapped to the reference *Ae. albopictus* genome from infected ovaries and carcasses, respectively, after grouping biological replicates. Similarly, 1,568,888,349 and 1,651,306,629 reads were mapped for un-infected control ovaries and carcasses, respectively. Reads size distribution of ovaries samples shows a peak at 21 nt in both infected and uninfected samples, as well as marked piRNAs production peaking at 28 nt (Figure 1). In uninfected and CHIKV- infected carcasses samples, only a peak at 21 nt was seen, indicating a lower production of piRNAs in carcasses (Figure 1).

### 3.1. Small RNA Profile of CHIKV, vDNA Fragments, and New nrEVEs

The CHIKV genome consists of a single positive strand RNA of ~11–12 Kb. Upon entry into a host cell and uncoating, the CHIKV genome is used both as mRNA for translation of nonstructural proteins and as the template for the synthesis of a complementary negative strand RNA molecule, which is the replicative form of the virus. A genomic length positive strand molecule is synthesized from the negative strand RNA, along with a shorter RNA molecule from an internal promoter. This subgenomic positive strand RNA molecule encodes for nonstructural proteins [86]. Small RNA sequencing data from CHIKV-infected and control samples were mapped against the CHIKV genome. Mapped sRNAs were filtered by size to assess siRNA- (20–22 nt) and piRNA- (25–30 nt) responses. A total of 19,008 and 316,527 siRNA size reads mapped to the CHIKV genome from infected ovaries and carcasses, respectively (Figure 2A). Average strand-specific RPKM normalized counts show that small RNAs map predominantly on the positive strand of CHIKV, with the only exception of siRNAs in CHIKV-infected ovaries. This result indicates that the RNAi activity is mostly focused on the negative strand (Table 1).

<table><caption>**Table 1.** RPKM values of strand-specific piRNA and siRNA reads which mapped to CHIKV.</caption><thead><tr><th rowspan="2">Sample</th><th colspan="2">piRNAs</th><th colspan="2">siRNA</th></tr><tr><th>+RPKM</th><th>-RPKM</th><th>+RPKM</th><th>-RPKM</th></tr></thead><tbody><tr><td>Ovaries</td><td>14,025 ± 1007</td><td>3036 ± 61</td><td>28,453 ± 1308</td><td>33,094 ± 766</td></tr><tr><td>Carcasses</td><td>11,610 ± 2993</td><td>1852 ± 613</td><td>38,731 ± 2230</td><td>26,837 ± 127</td></tr></tbody></table>

Consistently, piRNA size reads were found to map against CHIKV (Figure 2B), with a predominance of piRNA size reads on the positive strand of the virus (Figure 2C). We found a total of 5510 piRNA size reads in ovaries and 60,824 in carcasses of CHIKV-infected samples with a clear pick between positions ~7600–8700 of the viral genome, corresponding to the region encoding the C and E2 proteins. A pick of piRNAs at the 5' terminus of the CHIKV genome, around nucleotide position 8000, was previously identified in both *Ae. aegypti* and *Ae. albopictus* mosquitoes 4 dpi with CHIKV [87]. Accordingly, we saw a single vDNA form between positions 7964–8742 of the viral genome in infected ovaries (Figure 2B). A smaller fragment in the same region was found in carcasses. In a laboratory colony of *Ae. albopictus* derived from Vietnam, vDNA forms from nonstructural proteins (i.e., nsP2 and nsP4) were identified as early as 2 dpi and up to 9 dpi with CHIKV [30]. These vDNA forms were shown to influence the profile of siRNAs during the early phase of infection. Taken together, these results indicate that the subgenomic mRNA is a general hot spot for piRNAs in *Aedes* spp. mosquitoes both during acute and persistent CHIKV infection [42]. On the contrary, vDNA forms can be generated from different parts of the CHIKV genome, possibly depending on the combination of viral and mosquito strains. Additionally, while during the early phases of arboviral infection, vDNAs modulate siRNAs
in mosquitoes [30,33], during persistent infection, we could identify a unique vDNA form corresponding to a pick of piRNAs. Overlap analysis of piRNA reads revealed that the ping-pong amplification loop is active at this stage of CHIKV infection both in ovaries and carcasses (Figure 2D).

Using WGS data from the same samples, the presence of novel nrEVEs was tested, implementing VyPER [54] followed by ViR [55], resulting in no signals. This result indicates that no integration of vDNA fragments occurred following CHIKV infection in either soma (carcass) or germline (ovary) cells. The absence of nrEVEs following CHIKV infection may be due to the fact that integrations of viral sequences is a continuous but rare event in *Aedes* genomes and may depend on the viral load and the frequency of infection [36]. Additionally, no nrEVE with similarity to CHIKV has been characterized in arthropod genomes so far [34,35,88], suggesting that vDNA forms of CHIKV may remain episomal.

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Small RNA response against the infecting CHIKV. (A) siRNA profile along the CHIKV genome. Values on the y-axis represent quantification of the signal as Reads Per Kilobase Million (RPKM). (B) piRNA profile along the CHIKV genome and PCR amplification of the vDNA fragment. Wells 1 and 2 show the PCR results from the DNA of the two biological replicates of CHIKV-infected ovaries; wells 3–4 show the PCR results from the DNA of the corresponding uninfected controls. The red box indicates the position in which the vDNA fragment was identified. (C) Read size distribution of piRNA reads in forward (black), reverse (grey), and displaying 1U bias (dots). (D) Nucleotide bias at each position of piRNAs against CHIKV mapping to the sense (+) strand (left panel) and antisense (−) strand (right panel) in ovaries and carcasses, and corresponding PingPong overlap probabilities.</figcaption></figure>

## 3.2. *Viral SmallRNAs on Ae. albopictus Transcripts*

After having identified small RNAs that mapped to CHIKV, we checked whether any of these viral small RNA molecules also mapped to *Ae. albopictus* transcripts. A total of 68 siRNA reads from the CHIKV-infected samples, deriving primarily from two hotspots on the nsP3 coding region (Figure 3A), mapped to *Ae. albopictus* transcripts with an enrichment for afadin-like protein (LOC115259774), glucose transporter type1-like proteins (LOC109428702—LOC109424279), and an uncharacterized protein (LOC115263499) (Table S2). These siRNA reads and the corresponding CHIKV genomic sequences were further analyzed to check for hairpins using viralMir [68] and miRNAfold [69]. These small RNAs may originate from hairpin structures within the CHIKV genomic sequence
(Figure S4). The expression levels of the above-mentioned three transcripts were compared by qPCR between CHIKV-infected samples and un-infected controls (Figure 3B). Results showed a significant decrease in expression for transcripts encoding the glucose transporter and the uncharacterized protein (p-values < 0.05), but not that of afadin-like protein, in infected carcasses. Overall, these results suggest that CHIKV encodes for small RNAs which may target mosquito transcripts. However, whether these small RNAs are canonical miRNAs or are part of the heterogenous population of siRNAs produced following infection has to be further investigated. Modulation of pathways linked to energy production, including carbohydrate metabolism, has already been observed in mosquitoes following arboviral infection [89–91]. It was proposed that induction of glucose metabolism following viral infection may compensate the lack of energy, which derive from the breakdown of mitochondrial membranes resulting in a lack of ATP [90,92].

A

<figure><img src="image_6.png" /></figure>

B

<figure><img src="image_7.png" /><figcaption>**Figure 3.** (A) Graphical representation of the hotspots on the viral genome for small RNA molecules that map to both CHIKV and *Ae. albopictus* transcripts. (B) Bar plot of the qPCR results on three candidate target transcripts of viral small RNA molecules: afadin (XM_029860564.1), glucose- transporter (XM_029879742.1), and uncharacterized transcript (XM_029866737.1). ANOVA significant results at *p* < 0.05 are indicated by an asterisk.</figcaption></figure>

### 3.3. *The Profile of SmallRNAs on Ae. albopictus Transcripts*

Modulation of mosquito small RNAs has been shown after arboviral infection not only as a direct response to the virus, but also as a way to maintain cellular homeostasis, by regulating gene expression [28,47,48]. On these bases, we mapped small RNA data to *Ae. albopictus* transcripts and kept only reads that mapped on antisense orientation to assess differential abundance between infected and uninfected samples. A total of 376,119 siRNA reads and 5,684,597 piRNA reads mapped to mosquito transcripts in CHIKV-infected ovaries; 1,084,531 siRNA reads and 690,008 piRNA reads in CHIKV-infected carcasses. Regarding the corresponding un-infected controls, 302,702 siRNA reads and 5,780,782 piRNA reads mapped to transcripts in ovaries, and 244,103 siRNA reads and 243,947 piRNA reads in carcasses. A total of 30 and 534 transcripts displayed significant differential abundance (SDA) of siRNAs in CHIKV-infected ovaries and carcasses, respectively, representing six over-represented GO terms, including “DNA helicase activity”, “DNA recombination”, and “DNA repair” (Table S4). Interestingly, transcriptional regulation of genes associated
with DNA repair and chromatin function has been recently linked to ovaries and carcasses, respectively [93] (Figure 4A,B, Table S3).

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Differential abundance of small RNA molecules mapping to *Ae. albopictus* transcripts represented by average log counts-per-million (CPM) on the x-axis and log fold-change on the y-axis. Red dots indicate significantly over-expressed sRNAs and blue dots significantly under-expressed sRNAs upon infection compared to un-infected controls. Statistical significance of the data was assessed by likelihood ratio test. (A) siRNA abundance of transcripts in CHIKV-infected ovaries; (B) siRNA abundance of transcripts in CHIKV-infected carcasses; (C) piRNA abundance of transcripts in CHIKV-infected ovaries; (D) piRNA abundance of transcripts in CHIKV-infected carcasses. Triangles indicate transcripts tested by qPCR. (E) Bar plot of the qPCR results on transcripts to which SDA small RNAs are mapped to. Significant results at *p* < 0.05 are indicated by an asterisk. Significance was established by ANOVA. (F) Normalized coverage of the *Ae. albopictus* genomic region NW_021838576.1:117401889-117402177. The two shades of grey represent the biological replicates. Query bar indicates alignment similarity to XM_029878700.1 as computed by BLASTn online.</figcaption></figure>

GO enrichment analyses reveled an over representation representation of 6 terms, including "DNA helicase activity", "DNA recombination", and "DNA repair" (Table S4). Interestingly, transcriptional regulation of genes associated with DNA repair and chromatin function has been recently linked to transgenerational immunopriming (TGIP) in *D. melanogaster* after infection with Sindbis virus, an Alphavirus-like CHIKV, which was clearly demonstrated only in females [93]. Thus, it is tempting to speculate that during persistent infection of female mosquitoes, there is transcriptional regulation of functions that not only helps with maintenance of infection, but also prepares to TGIP. Because piRNAs have been shown to contribute in regulating mosquito gene expression [94], we also checked for transcripts with significantly different piRNA profile between infected and uninfected samples. A total of 299 and 81 transcripts displayed SDA of piRNAs in CHIKV-infected ovaries and carcasses, respectively (Figure 4C,D, Table S5). In carcasses, the majority of these transcripts showed lower abundance of piRNAs during infection. KOBAS gene-list enrichment analysis of these transcripts revealed an enrichment for several pathways, including oxidative phosphorylation, purine metabolism, and DNA damage-related pathways (Figure S5, Table S6).

Two transcripts displayed differential abundance for piRNAs in CHIKV-infected ovaries and carcasses as well as for siRNAs in CHIKV-infected ovaries; four transcripts displayed differential abundance for piRNAs and siRNAs in CHIKV-infected ovaries; 23 transcripts displayed differential abundance for piRNAs and siRNAs in CHIKV-infected carcasses (Figure S6). The expression levels of 3 transcripts (XM_029853235.1, uncharacterized; XM_029877572.1 encoding for a neurofilament heavy polypeptide-like protein and XM_029878700.1 for a hydrolase), which had consistent differential abundance of small RNAs across conditions and the relatively highest number of reads, were compared by qPCR between infected and un-infected samples. qPCR results show all three transcripts
are significantly under-expressed in carcasses during infection compared to un-infected controls (Figure 4E).

The hydrolase XM_029878700.1 was the only tested transcript with SDA of only piRNAs. These piRNAs were extracted and re-mapped to the *Ae. albopictus* genome after masking the genomic region encoding for XM_029878700.1, leading to the identification of genomic region NW_021838576.1:117401889-117402177. This region, which maps outside annotated piRNA clusters [34], shows high sequence homology to XM_029878700.1, but without annotated ORFs (Figure 4F). Annotation as "hydrolase" is a generic term as hydrolysis can be performed on different chemical bonds, including sugar glucosyl bonds and lipids. It is interesting to notice that sugar and lipid homeostasis is a common feature of *Aedes* spp. mosquitoes after infection not only with CHIKV, but also Flaviviruses [22,90–92,95–97], suggesting that the candidate transcript identified here should be further functionally investigated.

The database of *Ae. albopictus* miRNAs was recently published [34], so we took advantage of this resource to check for specific changes in mosquito miRNA expression levels during persistent CHIKV infection. A total of 14 miRNAs were detected as differentially modulated (Table 2), with the majority being upregulated (8) in ovaries but depleted in carcasses. All the DE miRNAs (except for miR-new6) have a homolog in *Ae. aegypti*. Differential expression levels ranged between ~0.20-fold for aal-miR-new6 and ~20-fold for aal-miR-210.

<table><caption>**Table 2.** Differentially expressed miRNAs in ovaries and carcasses. Total counts are calculated across sequencing libraries of infected replicates. Expression values are expressed as log2 transformed fold-changes (infected/control).</caption><thead><tr><th>miRNA</th><th>Total Count</th><th>Log2FC</th><th>FDR</th></tr></thead><tbody><tr><td><strong>Ovaries</strong></td><td></td><td></td><td></td></tr><tr><td>aal-miR-210</td><td>791</td><td>4.396</td><td>0.003</td></tr><tr><td>aal-miR-124</td><td>2592</td><td>4.209</td><td>0.001</td></tr><tr><td>aal-miR-1000</td><td>1147</td><td>3.177</td><td>0.014</td></tr><tr><td>aal-miR-219</td><td>115</td><td>2.980</td><td>0.014</td></tr><tr><td>aal-miR-932</td><td>2858</td><td>2.890</td><td>0.001</td></tr><tr><td>aal-miR-981b</td><td>260</td><td>2.864</td><td>0.026</td></tr><tr><td>aal-miR-193</td><td>93</td><td>2.725</td><td>0.028</td></tr><tr><td>aal-miR-285</td><td>2257</td><td>2.179</td><td>0.012</td></tr><tr><td>aal-miR-2941</td><td>94,365</td><td>-1.026</td><td>0.026</td></tr><tr><td>aal-miR-7</td><td>6849</td><td>-1.054</td><td>0.012</td></tr><tr><td>aal-miR-316</td><td>2558</td><td>-1.200</td><td>0.012</td></tr><tr><td>aal-miR-9b</td><td>4903</td><td>-1.247</td><td>0.006</td></tr><tr><td><strong>Carcasses</strong></td><td></td><td></td><td></td></tr><tr><td>aal-miR-new6</td><td>235</td><td>-2.399</td><td>0.016</td></tr><tr><td>aal-miR-1891</td><td>15,269</td><td>-1.929</td><td>0.016</td></tr></tbody></table>

We compared these results to a previous report of differentially expressed miRNAs in *Ae. albopictus* cells [98] and found only one miRNA present in both analyses (aal-miR-1000). Interestingly, miR-1000 was reported to be under-expressed in cells 24 h post infection, whereas we found it over-expressed in adult mosquitoes 14 dpi. We searched the *Ae. albopictus* transcriptome for potential target sites of the DE miRNAs, as well as the corresponding star miRNAs, performing thermodynamic analyses on the RNA-RNA pairing of the top 10 blast hits for each miRNA. We found a total of 229 possible targets (140 from over-expressed miRNAs, 89 from under-expressed miRNAs) (Table S7). KOBAS analysis revealed 12 significantly enriched pathways (Table 3), including the Toll and Imd pathways, suggesting the under regulation of these immunity pathways during persistent infection. "Lysine degradation" was also identified among the enriched pathways when looking at differentially abundant piRNAs (Figure S5).
<table><caption>**Table 3.** Biological pathways affected by differentially accumulated miRNAs. Biological pathways based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) database of transcripts identified as targets of differentially abundant miRNAs. The list of target transcripts was obtained by KEGG Orthology Based Annotation System (KOBAS) gene-list enrichment analysis.</caption><thead><tr><th>Term</th><th>Total Count</th><th>KEGG ID</th><th>p-Value</th></tr></thead><tbody><tr><td><strong>Over-Expressed miRNAs Targets</strong></td><td></td><td></td><td></td></tr><tr><td>Lysine degradation</td><td>9</td><td>aag00310</td><td>0.000</td></tr><tr><td>mTOR signaling pathway</td><td>9</td><td>aag04150</td><td>0.000</td></tr><tr><td>Toll and Imd signaling pathway</td><td>3</td><td>aag04624</td><td>0.003</td></tr><tr><td>Phosphatidylinositol signaling system</td><td>3</td><td>aag04070</td><td>0.006</td></tr><tr><td>Apoptosis</td><td>3</td><td>aag04214</td><td>0.007</td></tr><tr><td>Terpenoid backbone biosynthesis</td><td>2</td><td>aag00900</td><td>0.014</td></tr><tr><td>AGE-RAGE signaling pathway</td><td>2</td><td>aag04933</td><td>0.019</td></tr><tr><td>MAPK signaling pathway</td><td>3</td><td>aag04013</td><td>0.022</td></tr><tr><td>Endocytosis</td><td>3</td><td>aag04144</td><td>0.045</td></tr><tr><td><strong>Under-Expressed miRNAs Targets</strong></td><td></td><td></td><td></td></tr><tr><td>SNARE interactions in vesicular transport</td><td>3</td><td>aag04130</td><td>0.000</td></tr><tr><td>Proteasome</td><td>3</td><td>aag03050</td><td>0.001</td></tr><tr><td>Phagosome</td><td>3</td><td>aag04145</td><td>0.003</td></tr></tbody></table>

We also tested for potential targets site of the DE miRNAs on the CHIKV genome, but no significant alignment was found, indicating that these miRNAs cannot directly target CHIKV. Overall, these results show a different modulation of miRNAs in ovaries and carcasses. Depletion of miRNAs in carcasses is a trend consistent with what has been observed following Flavivirus infection (i.e., Zika and DENV) [99,100].

### 3.4. Differential Abundance of nrEVEs-sRNAs Following Infection and Their Targets on the Ae. albopictus Transcriptome

The genome of *Ae. albopictus* harbors hundreds of nrEVEs, which are variably distributed among the genomes of Foshan mosquitoes [40]. None of these nrEVEs derive from Alphaviruses and the maximum sequence similarity between annotated nrEVE and CHIKV is 70% (Figure 5A). Nevertheless, nrEVEs produce piRNAs [35] and responder and trailer piRNAs are produced from viral RNAs guided by endogenous piRNAs, indicating that piRNA biogenesis can spread outside a putative nrEVE targeted region [31,101]. On this basis, we first verified which nrEVEs our mosquitoes had, then tested whether the piRNA profile of these nrEVEs was different in infected vs. un-infected-fed samples and finally checked for putative targets of the differentially accumulated nrEVEs-derived piRNAs.

Differences in the piRNA profile of nrEVEs were observed only in CHIKV-infected ovaries. A total of 24 nrEVEs mapping in different clusters showed a piRNA profile that was significantly different between infected and un-infected samples (Figure S7, Table S8). piRNAs from these nrEVEs peaked between 27 and 29 nt (Figure 5C). When looking at single piRNAs (Figure S1B, level 1), a total of 34 and 4 piRNA molecules were significantly differentially expressed in infected vs un-infected ovaries and carcasses (Figure 5B), respectively, and peaked at 27 nt (Figure 5C). Of these 38 piRNAs, 30 originated from nrEVEs embedded in piRNA clusters. We focused on CHIKV-infected ovaries. Intersecting the genomic origins of the DE piRNAs with the SDA nrEVEs, we identified 24 DE piRNAs from 8 SDA nrEVEs (Table S9). These 24 piRNAs do not blast on the CHIKV genome indicating that they do not target the viral mRNA, but they blast against *Ae. albopictus* transcripts. The free energy of the seed pairing between these piRNAs and their candidate target transcripts was evaluated with IntaRNA, showing that these piRNAs could interfere with the detected transcripts. Two piRNAs (1049-820 and 3135-266), both originating from nrEVE Rhabdo22 (Figure 5D, red bar), which is embedded in piRNA cluster 461, have a common target, transcript XM_029880513.1, an uncharacterized protein (Figure 5E). In accordance with bioinformatic data, qPCR analysis confirmed a higher expression of 3135-266 in infected vs. un-infected ovaries and a lower expression of the candidate target
transcript (Figure 5F). Taken together, these results suggest that nrEVEs-derived piRNAs
can be differentially accumulated during CHIKV infection and regulate the expression
profile of endogenous transcripts.

<figure><img src="image_9.png" /></figure>

**Figure 5.** (A) Box plot representing nonretroviral Endogenous Viral Element (nrEVE) sequence identity with respect to the CHIKV genome. Chu stands for *Chuviridae* (dark blue); Chulike stands for *Chuviridae*-like (red); Flavi stands for *Flaviviridae* (light green); Meso stands for *Mesoniviridae* (gray); Phasma stands for *Phasmaviridae* (light blue); Phenui stands for *Phenuiviridae* (light orange); Rhabdo stands for *Rhabdoviridae* (purple); Toti stands for *Totiviridae* (dark green); Un stands for unclassified (dark red); Virgo stands for *Virgoviridae* (pink); Xinmo stands for *Xinmoviridae* (dark orange). (B) Differential abundance of piRNAs from nrEVEs in CHIKV-infected samples represented by average log counts-per-million (CPM) on the x-axis and log fold-change on the y-axis. Red dots indicate significantly over-expressed piRNAs and blue dots significantly under-expressed piRNAs upon infection compared to un-infected controls. Significance was assessed by likelihood ratio test. Triangles indicate piRNAs tested by qPCR. (C) Size distribution of piRNAs in CHIKV-infected ovaries. Left: size distribution of all piRNAs derived from nrEVEs. Right: size distribution of differentially abundant piRNAs derived from nrEVEs. Values are normalized according to the TMMwzp normalization method implemented in edgeR. (D) piRNA expression levels during infection in cluster visualized in IGV; nrEVEs in this region are highlighted (Rhabdo22 and Rhabdo23) as well as the peak representing the piRNAs used for qPCR analysis. (E) Graphical sequence panel (NCBI) representation of transcript XM_029880513.1; the piRNA target sites are represented on the bottom along with interaction energy as computed by IntaRNA. (F) Bar plot of the qPCR results on selected piRNAs and their putative target transcript. Significant results at *p* < 0.05 are indicated by an asterisk.
## 4. Conclusions

Here, we describe the profile of small RNAs and the presence of vDNA forms and nrEVEs during persistent infection of *Ae. albopictus* with CHIKV. Because RNAi is the main antiviral pathway of *Aedes* spp. mosquitoes, a number of studies have described the profile of small RNAs following infection with arboviruses, including CHIKV [87,98–100,102–105]. In line with previous reports, we found a strong exogenous siRNA response throughout the whole CHIKV genome in both ovaries and carcasses [30,88]. On the contrary piRNAs aligned primarily to a region encompassing the C and E2 coding sequences from which a unique vDNA form was also identified, in both analyzed tissues. This result differs from previous studies in which various vDNAs from multiple regions of the genome of the flaviviruses Zika, dengue, and West Nile virus were identified 24 h post infection and up to 21 dpi in mosquito cells and adults, when infection has become persistent [30,32,106]. Whether this discrepancy is dependent on the specific viral species and mosquito species/strains studied or it indicates that only the vDNA molecule producing the most piRNAs is retained once infection has become persistent requires further investigations.

CHIKV infection did not result in the formation of novel viral integrations in the genome of *Ae. albopictus*. However, the piRNA profile of nrEVEs in the genomes of the tested mosquitoes changed and two nrEVE-derived piRNAs were further validated by qRT-PCR to be more abundant in CHIKV-infected than uninfected ovaries. Accordingly, the transcript to which these piRNAs map showed decreased in expression in CHIKV-infected ovaries. Thus, nrEVEs may have other functions than antagonistic action against cognate viruses, which was demonstrated in ovaries [38]. This result may explain previous findings describing that the presence/absence of certain nrEVEs correlates with differences in dissemination efficiency for CHIKV and dengue viruses, even if the genome of neither *Ae. aegypti* not *Ae. albopictus* has nrEVEs derived from these arboviruses [34,35,39]. To conclude, here we provided an in-depth analysis of small RNAs in CHIKV-infected *Ae. albopictus* mosquitoes at 14 dpi, when CHIKV has established persistent infection and identify target transcripts of differentially abundant small RNAs that merits further investigations. These results clearly show that the activity of the RNAi pathway modulates sRNAs and piRNAs, which contribute to cellular homeostasis, affecting primarily functions related to lipid and sugar metabolisms and DNA repair/binding. Finally, the identified correlation between the differential accumulation nrEVE-piRNAs and the decreased expression of the piRNA-target expands the possible role of nrEVEs beyond immunity against cognate viruses.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/1999-4911/5/13/4/553/s1, Figure S1: (A) Scheme of the infection experiments. (B) Strategy used to study the small RNA profile of nrEVEs, Figure S2: Length distribution of the full collection of *Aedes aegypti* miRNAs downloaded from miRbase, Figure S3: Venn diagram representing the intersection of the sRNA differential abundance across each transcript based on the three different library sizes for normalization, Figure S4: Hairpin predictions computed with miRNAfold, Figure S5: Bar plot of significantly enriched pathways based on differential sRNA abundance in CHIKV-infected ovaries and carcasses, Figure S6: Venn diagram of the transcripts displaying differential abundance across conditions, Figure S7: Venn diagram representing the intersection of the sRNA differential abundance across each nrEVE based on the three different library sizes for normalization, Table S1: List of primers used in this study, Table S2: List of candidate targets of potential viral-sRNAs, Table S3: List of transcripts displaying significant differential abundance for siRNAs, Table S4: GO enrichment analysis of siRNA SDA transcripts, Table S5: List of transcripts displaying significant differential abundance for piRNAs, Table S6: KOBAS gene-list enrichment analysis of SDA transcripts, Table S7: List of candidate targets of DE miRNAs, Table S8: Differential abundance analysis of nrEVEs, Table S9: DE piRNAs originating from SDA nrEVEs.

**Author Contributions:** Conceptualization, M.B. and M.M.; methodology, M.M., E.P., U.P., V.H., D.S. and N.L.-C.; resources, M.B. and A.-B.F.; data curation, M.M.; writing—original draft preparation, M.M.; writing—review and editing, M.B.; supervision, M.B. and A.-B.F.; project administration, M.M.
and M.B.; funding acquisition, M.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Italian Ministry of Education, University and Research FARE MIUR project R1623HZAH5 to M.B.; by a European Research Council Consolidator Grant (ERC-CoG) under the European Union's Horizon 2020 Programme (Grant Number ERC-CoG 682394) to M.B.; and by the Italian Ministry of Education, University and Research (MIUR): Dipartimenti Ecellenza Program (2018–2022) to the Dept. of Biology and Biotechnology “L. Spallanzani”, University of Pavia.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** WGS data were deposited in the Sequence Read Archive (SRA) under Bioproject: PRJNA659388. Small RNA data were deposited in SRA under Bioproject: PRJNA607026.

**Acknowledgments:** We acknowledge the technical assistance of Francesca Scolari and Davide Sacco for mosquito maintenance.

**Conflicts of Interest:** The authors declare no conflict of interest.

## References

1. Mason, P.J.; Haddow, A.J. An epidemic of virus disease in Southern Province, Tanganyika Territory, in 1952–1953; An additional note on Chikungunya virus isolations and serum antibodies. *Trans. R. Soc. Trop. Med. Hyg.* **1957**, *51*, 238–240. [CrossRef]

2. Njenga, K.M.; Nderitu, L.; Ledermann, J.P.; Ndirangu, A.; Logue, C.H.; Kelly, C.H.L.; Sang, R.; Sergon, K.; Breiman, R.; Powers, A.M. Tracking epidemic Chikungunya virus into the Indian Ocean from East Africa. *J. Gen. Virol.* **2008**, *89*, 2754–2760. [CrossRef]

3. Zeller, H.; Van Bortel, W.; Sudre, B. Chikungunya: Its history in Africa and Asia and its spread to new regions in 2013–2014. *J. Infect. Dis.* **2016**, *214*, S436–S440. [CrossRef] [PubMed]

4. Leparc-Goffart, I.; Nougairede, A.; Cassadou, S.; Prat, C.; de Lamballerie, X. Chikungunya in the Americas. *Lancet* **2014**, *383*, 514. [CrossRef]

5. Amraoui, F.; Failloux, A.-B. Chikungunya: An unexpected emergence in Europe. *Curr. Opin. Virol.* **2016**, *21*, 146–150. [CrossRef] [PubMed]

6. Weaver, S.C.; Forrester, N.L. Chikungunya: Evolutionary history and recent epidemic spread. *Antivir. Res.* **2015**, *120*, 32–39. [CrossRef] [PubMed]

7. World Health Organization. Neglected Diseases. Available online: https://www.who.int/neglected_diseases/diseases/en/ (accessed on 7 September 2020).

8. Cox, J.; Mota, J.; Sukupolvi-Petty, S.; Diamond, M.S.; Rico-Hesse, R. Mosquito bite delivery of dengue virus enhances immunogenicity and pathogenesis in humanized mice. *J. Virol.* **2012**. [CrossRef]

9. Liu, J.; Swevers, L.; Kolliopoulou, A.; Smagghe, G. Arboviruses and the challenge to establish systemic and persistent infections in competent mosquito vectors: The interaction with the RNAi mechanism. *Front. Physiol.* **2019**, *10*, 890. [CrossRef]

10. Dubrulle, M.; Mousson, L.; Moutailler, S.; Vazeille, M.; Failloux, A.-B. Chikungunya virus and *Aedes* mosquitoes: Saliva is infectious as soon as two days after oral infection. *PLoS ONE* **2009**, *4*, e5895. [CrossRef]

11. Merkling, S.H.; van Rij, R.P. Beyond RNAi: Antiviral defense strategies in Drosophila and mosquito. *J. Insect Physiol.* **2013**, *59*, 159–170. [CrossRef]

12. Cheng, G.; Liu, Y.; Wang, P.; Xiao, X. Mosquito defense strategies against viral infection. *Trends Parasitol.* **2016**, *32*, 177–186. [CrossRef] [PubMed]

13. Lee, W.S.; Webster, J.A.; Madzokere, E.T.; Stephenson, E.B.; Herrero, L.J. Mosquito antiviral defense mechanisms: A delicate balance between innate immunity and persistent viral infection. *Parasites Vectors* **2019**, *12*, 1–12. [CrossRef] [PubMed]

14. Bosio, C.F.; Fulton, R.E.; Salasek, M.L.; Beaty, B.J.; Black, W.C. 4th Quantitative trait loci that control vector competence for dengue-2 virus in the mosquito *Aedes aegypti*. *Genetics* **2000**, *156*, 687–698.

15. Bennett, K.E.; Flick, D.; Fleming, K.H.; Jochim, R.; Beaty, B.J.; Black, W.C., 4th. Quantitative trait loci that control dengue-2 virus dissemination in the mosquito *Aedes aegypti*. *Genetics* **2005**, *170*, 185–194. [CrossRef] [PubMed]

16. Franz, A.W.E.; Kantor, A.M.; Passarelli, A.L.; Clem, R.J. Tissue barriers to arbovirus infection in mosquitoes. *Viruses* **2015**, *7*, 3741–3767. [CrossRef]

17. Palmer, W.H.; Varghese, F.S.; van Rij, R.P. Natural variation in resistance to virus infection in Dipteran insects. *Viruses* **2018**, *10*, 118. [CrossRef]

18. Garcia-Luna, S.M.; Weger-Lucarelli, J.; Rückert, C.; Murrieta, R.A.; Young, M.C.; Byas, A.D.; Fauver, J.R.; Perera, R.; Flores-Suarez, A.E.; Ponce-Garcia, G.; et al. Variation in competence for ZIKV transmission by *Aedes aegypti* and *Aedes albopictus* in Mexico. *PLoS Negl. Trop. Dis.* **2018**, *12*, e0006599. [CrossRef]

19. Dickson, L.B.; Sanchez-Vargas, I.; Sylla, M.; Fleming, K.; Black, W.C., 4th. Vector competence in West African *Aedes aegypti* is flavivirus species and genotype dependent. *PLoS Negl. Trop. Dis.* **2014**, *8*, e3153. [CrossRef] [PubMed]
20. Lambrechts, L.; Scott, T.W. Mode of transmission and the evolution of arbovirus virulence in mosquito vectors. *Proc. R. Soc. B Biol. Sci.* **2009**. [CrossRef]

21. Salas-Benito, J.S.; De Nova-Ocampo, M. Viral interference and persistence in mosquito-borne flaviviruses. *J. Immunol. Res.* **2015**, *2015*, 873404. [CrossRef]

22. Chotivan, N.; Andre, B.G.; Sanchez-Vargas, I.; Islam, M.N.; Grabowski, J.M.; Hopf-Jannasch, A.; Gough, E.; Nakayasu, E.; Blair, C.D.; Belisle, J.T.; et al. Dynamic remodeling of lipids coincides with dengue virus replication in the midgut of *Aedes aegypti* mosquitoes. *PLoS Pathog.* **2018**, *14*, e1006853. [CrossRef]

23. Da Cruz-Nunes, E.; Canuto, G.A.B. Metabolomics applied in the study of emerging arboviruses caused by *Aedes aegypti* mosquitoes: A review. *Electrophoresis* **2020**, *41*, 2102–2113. [CrossRef] [PubMed]

24. Lambrechts, L.; Saleh, M.-C. Manipulating mosquito tolerance for arbovirus control. *Cell Host Microbe* **2019**, *26*, 309–313. [CrossRef] [PubMed]

25. Olson, K.E.; Blair, C.D. Arbovirus-mosquito interactions: RNAi pathway. *Curr. Opin. Virol.* **2015**, *15*, 119–126. [CrossRef]

26. Samuel, G.H.; Adelman, Z.N.; Myles, K.M. Antiviral immunity and virus-mediated antagonism in disease vector mosquitoes. *Trends Microbiol.* **2018**, *26*, 447–461. [CrossRef] [PubMed]

27. Mongelli, V.; Saleh, M.-C. Bugs are not to be silenced: Small RNA pathways and antiviral responses in insects. *Annu. Rev. Virol.* **2016**, *3*, 573–589. [CrossRef]

28. Ozata, D.M.; Gainetdinov, I.; Zoch, A.; O'Carroll, D.; Zamore, P.D. PIWI-interacting RNAs: Small RNAs with big functions. *Nat. Rev. Genet.* **2019**, *20*, 89–108. [CrossRef] [PubMed]

29. Asgari, S. Role of microRNAs in arbovirus/vector interactions. *Viruses* **2014**, *6*, 3514–3534. [CrossRef]

30. Goic, B.; Stapleford, K.A.; Frangeul, L.; Doucet, A.J.; Gausson, V.; Blanc, H.; Schemmel-Jofre, N.; Cristofari, G.; Lambrechts, L.; Vignuzzi, M.; et al. Virus-derived DNA drives mosquito vector tolerance to arboviral infection. *Nat. Commun.* **2016**, *7*, 1–10. [CrossRef]

31. Joosten, J.; Van Rij, R.P.; Miesen, P. Slicing of viral RNA guided by endogenous piRNAs triggers the production of responder and trailer piRNAs in *Aedes* mosquitoes. *bioRxiv* **2020**. [CrossRef]

32. Nag, D.K.; Brecher, M.; Kramer, L.D. DNA forms of arboviral RNA genomes are generated following infection in mosquito cell cultures. *Virology* **2016**, *498*, 164–171. [CrossRef]

33. Poirier, E.Z.; Goic, B.; Tomé-Poderti, L.; Frangeul, L.; Boussier, J.; Gausson, V.; Blanc, H.; Vallet, T.; Loyd, H.; Levi, L.I.; et al. Dicer-2-dependent generation of viral DNA from defective genomes of RNA viruses modulates antiviral immunity in insects. *Cell Host Microbe* **2018**, *23*, 353–365.e8. [CrossRef]

34. Palatini, U.; Masri, R.A.; Cosme, L.V.; Koren, S.; Thibaud-Nissen, F.; Biedler, J.K.; Krsticevic, F.; Johnston, J.S.; Halbach, R.; Crawford, J.E.; et al. Improved reference genome of the arboviral vector *Aedes albopictus*. *Genome Biol.* **2020**, *21*, 215. [CrossRef]

35. Palatini, U.; Miesen, P.; Carballar-Lejarazu, R.; Ometto, L.; Rizzo, E.; Tu, Z.; van Rij, R.P.; Bonizzoni, M. Comparative genomics shows that viral integrations are abundant and express piRNAs in the arboviral vectors *Aedes aegypti* and *Aedes albopictus*. *BMC Genom.* **2017**, *18*, 1–15. [CrossRef] [PubMed]

36. Crava, C.M.; Varghese, F.S.; Pischedda, E.; Halbach, R.; Palatini, U.; Marconcini, M.; Gasmi, L.; Redmond, S.; Afrane, Y.; Ayala, D.; et al. Population genomics in the arboviral vector *Aedes aegypti* reveals the genomic architecture and evolution of endogenous viral elements. *Mol. Ecol.* **2021**. [CrossRef]

37. Tassetto, M.; Kunitomi, M.; Whitfield, Z.J.; Dolan, P.T.; Sánchez-Vargas, I.; Garcia-Knight, M.; Ribiero, I.; Chen, T.; Olson, K.E.; Andino, R. Control of RNA viruses in mosquito cells through the acquisition of vDNA and endogenous viral elements. *eLife* **2019**, *8*, e41244. [CrossRef] [PubMed]

38. Suzuki, Y.; Baidaliuk, A.; Miesen, P.; Van Rij, R.P.; Lambrechts, L.; Saleh, M.-C. Non-retroviral endogenous viral element limits cognate virus replication in *Aedes aegypti* ovaries CFAV-eve efficient viral replication CFAV. *Curr. Biol.* **2020**, *30*, 3495–3506.e6. [CrossRef]

39. Houé, V.; Gabiane, G.; Dauga, C.; Suez, M.; Madec, Y.; Mousson, L.; Marconcini, M.; Yen, P.-S.; de Lamballerie, X.; Bonizzoni, M.; et al. Evolution and biological significance of flaviviral elements in the genome of the arboviral vector *Aedes albopictus*. *Emerg. Microbes Infect.* **2019**, *8*, 1265–1279. [CrossRef] [PubMed]

40. Pischedda, E.; Scolari, F.; Valerio, F.; Carballar-Lejarazu, R.; Catapano, P.L.; Waterhouse, R.M.; Bonizzoni, M. Insights into an unexplored component of the mosquito repeatome: Distribution and variability of viral sequences integrated into the genome of the arboviral vector *Aedes albopictus*. *Front. Genet.* **2019**, *10*, 93. [CrossRef] [PubMed]

41. Göertz, G.P.; Miesen, P.; Overheul, G.J.; Van Rij, R.P.; Van Oers, M.M.; Pijlman, G.P. Mosquito small RNA responses to West Nile and insect-specific virus infections in *Aedes* and *Culex* mosquito cells. *Viruses* **2019**, *11*, 271. [CrossRef]

42. Miesen, P.; Joosten, J.; van Rij, R.P. PIWIs go viral: Arbovirus-derived piRNAs in vector mosquitoes. *PLoS Pathog.* **2016**, *12*, e1006017. [CrossRef] [PubMed]

43. Cardin, S.-E.; Borchert, G.M. Viral microRNAs, host microRNAs regulating viruses, and bacterial microRNA-like RNAs. *Methods Mol. Biol.* **2017**, *1617*, 39–56. [CrossRef] [PubMed]

44. Frangeul, L.; Blanc, H.; Saleh, M. Differential Small RNA Responses against co-infecting insect-specific viruses in *Aedes albopictus* mosquitoes. *Viruses* **2020**, *12*, 468. [CrossRef] [PubMed]

45. Blair, C.D. Deducing the role of virus genome-derived PIWI-associated RNAs in the mosquito-arbovirus arms race. *Front. Genet.* **2019**, *10*, 1114. [CrossRef] [PubMed]
46. Skalsky, R.L.; Olson, K.E.; Blair, C.D.; Garcia-Blanco, M.A.; Cullen, B.R. A "microRNA-like" small RNA expressed by Dengue virus? *Proc. Natl. Acad. Sci. USA* **2014**, *111*, E2359. [CrossRef] [PubMed]

47. Hussain, M.; Asgari, S. MicroRNA-like viral small RNA from Dengue virus 2 autoregulates its replication in mosquito cells. *Proc. Natl. Acad. Sci. USA* **2014**, *111*, 2746–2751. [CrossRef]

48. Hussain, M.; Torres, S.; Schnettler, E.; Funk, A.; Grundhoff, A.; Pijlman, G.P.; Khromykh, A.A.; Asgari, S. West Nile virus encodes a microRNA-like small RNA in the 3' untranslated region which up-regulates GATA4 mRNA and facilitates virus replication in mosquito cells. *Nucleic Acids Res.* **2012**, *40*, 2210–2223. [CrossRef]

49. Sigle, L.T.; McGraw, E.A. Expanding the canon: Non-classical mosquito genes at the interface of arboviral infection. *Insect Biochem. Mol. Biol.* **2019**, *109*, 72–80. [CrossRef] [PubMed]

50. Mousson, L.; Martin, E.; Zouache, K.; Madec, Y.; Mavingui, P.; Failloux, A.B. Wolbachia modulates Chikungunya replication in *Aedes albopictus*. *Mol. Ecol.* **2010**, *19*, 1953–1964. [CrossRef]

51. Marconcini, M.; Hernandez, L.; Iovino, G.; Houé, V.; Valerio, F.; Palatini, U.; Pischedda, E.; Crawford, J.E.; White, B.J.; Lin, T.; et al. Polymorphism analyses and protein modelling inform on functional specialization of Piwi clade genes in the arboviral vector *Aedes albopictus*. *PLoS Negl. Trop. Dis.* **2019**, *13*, e0007919. [CrossRef]

52. Li, H. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. *Broad Inst. Harv. MIT* **2013**, *00*, 1–3.

53. McKenna, A.; Hanna, M.; Banks, E.; Sivachenko, A.; Cibulskis, K.; Kerytsky, A.; Garimella, K.; Altshuler, D.; Gabriel, S.; Daly, M.; et al. The Genome Analysis Toolkit: A MapReduce framework for analyzing next-generation DNA sequencing data. *Genome Res.* **2010**, *20*, 1297–1303. [CrossRef] [PubMed]

54. Forster, M.; Szymczak, S.; Ellinghaus, D.; Hemmrich, G.; Rühlemann, M.; Kraemer, L.; Mucha, S.; Wienbrandt, L.; Stanulla, M.; Franke, A.; et al. Vy-PER: Eliminating false positive detection of virus integration events in next generation sequencing data. *Sci. Rep.* **2015**, *5*, 11534. [CrossRef] [PubMed]

55. Pischedda, E.; Crava, C.; Carlassara, M.; Zucca, S.; Gasmi, L.; Bonizzoni, M. ViR: A tool to solve intrasample variability in the prediction of viral integration sites using whole genome sequencing data. *BMC Bioinform.* **2021**, *7*, 1–15. [CrossRef]

56. Katoh, K.; Standley, D.M. MAFFT multiple sequence alignment software version 7: Improvements in performance and usability. *Mol. Biol. Evol.* **2013**, *30*, 772–780. [CrossRef]

57. Madeira, F.; Park, Y.M.; Lee, J.; Buso, N.; Gur, T.; Madhusoodanan, N.; Basutkar, P.; Tivey, A.R.N.; Potter, S.C.; Finn, R.D.; et al. The EMBL-EBI search and sequence analysis tools APIs in 2019. *Nucleic Acids Res.* **2019**, *47*, W636–W641. [CrossRef]

58. Rosenkranz, D.; Han, C.T.; Roovers, E.F.; Zischler, H.; Ketting, R.F. Piwi proteins and piRNAs in mammalian oocytes and early embryos: From sample to sequence. *Genom. Data* **2015**, *5*, 309–313. [CrossRef]

59. Ramírez, F.; Ryan, D.P.; Grüning, B.; Bhardwaj, V.; Kilpert, F.; Richter, A.S.; Heyne, S.; Dündar, F.; Manke, T. deepTools2: A next generation web server for deep-sequencing data analysis. *Nucleic Acids Res.* **2016**, *44*, W160–W165. [CrossRef]

60. Robinson, J.T.; Thorvaldsdóttir, H.; Winckler, W.; Guttman, M.; Lander, E.S.; Getz, G.; Mesirov, J.P. Integrative genomics viewer. *Nat. Biotechnol.* **2011**, *29*, 24–26. [CrossRef]

61. Uhrig, S.; Klein, H. PingPongPro: A tool for the detection of piRNA-mediated transposon-silencing in small RNA-Seq data. *Bioinformatics* **2018**, *35*, 335–336. [CrossRef]

62. Antoniewski, C. Computing siRNA and piRNA overlap signatures. *Methods Mol. Biol.* **2014**, *6*, 135–146.

63. Kozomara, A.; Griffiths-Jones, S. miRBase: Integrating microRNA annotation and deep-sequencing data. *Nucleic Acids Res.* **2010**, *39*, D152–D157. [CrossRef]

64. Crooks, G.E.; Hon, G.; Chandonia, J.-M.; Brenner, S.E. WebLogo: A sequence logo generator. *Genome Res.* **2004**, *14*, 1188–1190. [CrossRef]

65. BBMAP. Available online: https://sourceforge.net/projects/bmap/ (accessed on 15 September 2020).

66. Afgan, E.; Baker, D.; Batut, B.; van den Beek, M.; Bouvier, D.; Čech, M.; Chilton, J.; Clements, D.; Coraor, N.; Grüning, B.A.; et al. The Galaxy platform for accessible, reproducible and collaborative biomedical analyses: 2018 update. *Nucleic Acids Res.* **2018**, *46*, W537–W544. [CrossRef]

67. Li, H.; Handsaker, B.; Wysoker, A.; Fennell, T.; Ruan, J.; Homer, N.; Marth, G.; Abecasis, G.; Durbin, R.; 1000 Genome Project Data Processing Subgroup. The sequence alignment/map format and SAMtools. *Bioinformatics* **2009**, *25*, 2078–2079. [CrossRef]

68. Huang, K.-Y.; Lee, T.-Y.; Teng, Y.-C.; Chang, T.-H. ViralmiR: A support-vector-machine-based method for predicting viral microRNA precursors. *BMC Bioinform.* **2015**, *16* (Suppl. S1), S1–S9. [CrossRef]

69. Tav, C.; Tempel, S.; Poligny, L.; Tahi, F. miRNAFold: A web server for fast miRNA precursor prediction in genomes. *Nucleic Acids Res.* **2016**, *44*, W181–W184. [CrossRef]

70. Liao, Y.; Smyth, G.K.; Shi, W. featureCounts: An efficient general purpose program for assigning sequence reads to genomic features. *Bioinformatics* **2013**, *30*, 923–930. [CrossRef]

71. Robinson, M.D.; McCarthy, D.J.; Smyth, G.K. edgeR: A Bioconductor package for differential expression analysis of digital gene expression data. *Bioinformatics* **2010**, *26*, 139–140. [CrossRef] [PubMed]

72. Chen, Y.; Lun, A.T.L.; Smyth, G.K. From reads to genes to pathways: Differential expression analysis of RNA-Seq experiments using Rsubread and the edgeR quasi-likelihood pipeline. *F1000Research* **2016**, *5*, 1438. [CrossRef] [PubMed]

73. McCormick, K.P.; Willmann, M.R.; Meyers, B.C. Experimental design, preprocessing, normalization and differential expression analysis of small RNA sequencing experiments. *Silence* **2011**, *2*, 2. [CrossRef]
74. Locati, M.D.; Terpstra, I.; de Leeuw, W.C.; Kuzak, M.; Rauwerda, H.; Ensink, W.A.; van Leeuwen, S.; Nehrdich, U.; Spaink, H.P.; Jonker, M.J.; et al. Improving small RNA-seq by using a synthetic spike-in set for size-range quality control together with a set for data normalization. *Nucleic Acids Res.* **2015**, *43*, e89. [CrossRef]

75. Akkouche, A.; Mugat, B.; Barckmann, B.; Varela-Chavez, C.; Li, B.; Raffel, R.; Pélisson, A.; Chambeyron, S. Piwi is required during drosophila embryogenesis to license dual-strand piRNA clusters for transposon repression in adult ovaries. *Mol. Cell* **2017**, *66*, 411–419.e4. [CrossRef]

76. Venny. Available online: https://bioinfogp.cnb.csic.es/tools/venny/ (accessed on 3 August 2020).

77. Götz, S.; García-Gómez, J.M.; Terol, J.; Williams, T.D.; Nagaraj, S.H.; Nueda, M.J.; Robles, M.; Talón, M.; Dopazo, J.; Conesa, A. High-throughput functional annotation and data mining with the Blast2GO suite. *Nucleic Acids Res.* **2008**, *36*, 3420–3435. [CrossRef] [PubMed]

78. Xie, C.; Mao, X.; Huang, J.; Ding, Y.; Wu, J.; Dong, S.; Kong, L.; Gao, G.; Li, C.-Y.; Wei, L. KOBAS 2.0: A web server for annotation and identification of enriched pathways and diseases. *Nucleic Acids Res.* **2011**, *39*, W316–W322. [CrossRef]

79. FASTX. Available online: http://hannonlab.cshl.edu/fastx_toolkit/ (accessed on 3 August 2020).

80. Altschul, S.F.; Gish, W.; Miller, W.; Myers, E.W.; Lipman, D.J. Basic local alignment search tool. *J. Mol. Biol.* **1990**, *215*, 403–410. [CrossRef]

81. Mann, M.; Wright, P.R.; Backofen, R. IntaRNA 2.0: Enhanced and customizable prediction of RNA-RNA interactions. *Nucleic Acids Res.* **2017**, *45*, W435–W439. [CrossRef]

82. Batz, Z.A.; Goff, A.C.; Armbruster, P.A. MicroRNAs are differentially abundant during *Aedes albopictus* diapause maintenance but not diapause induction. *Insect Mol. Biol.* **2017**, *26*, 721–733. [CrossRef] [PubMed]

83. Zhang, D.; Tu, S.; Stubna, M.; Wu, W.-S.; Huang, W.-C.; Weng, Z.; Lee, H.-C. The piRNA targeting rules and the resistance to piRNA silencing in endogenous genes. *Science* **2018**, *359*, 587–592. [CrossRef]

84. Lavezzo, E.; Falda, M.; Fontana, P.; Bianco, L.; Toppo, S. Enhancing protein function prediction with taxonomic constraints—The Argot2.5 web server. *Methods* **2016**, *93*, 15–23. [CrossRef]

85. Reynolds, J.A.; Poelchau, M.F.; Rahman, Z.; Armbruster, P.A.; Denlinger, D.L. Transcript profiling reveals mechanisms for lipid conservation during diapause in the mosquito, *Aedes albopictus*. *J. Insect Physiol.* **2012**, *58*, 966–973. [CrossRef]

86. Pietilä, M.K.; Hellström, K.; Ahola, T. Alphavirus polymerase and RNA replication. *Virus Res.* **2017**, *234*, 44–57. [CrossRef]

87. Morazzani, E.M.; Wiley, M.R.; Murreddu, M.G.; Adelman, Z.N.; Myles, K.M. Production of virus-derived ping-pong-dependent piRNA-like small RNAs in the mosquito soma. *PLoS Pathog.* **2012**, *8*, e1002470. [CrossRef]

88. ter Horst, A.M.; Nigg, J.C.; Dekker, F.M.; Falk, B.W. Endogenous viral elements are widespread in arthropod genomes and commonly give rise to PIWI-interacting RNAs. *J. Virol.* **2019**, *93*, e02124-18. [CrossRef]

89. Lee, R.C.H.; Chu, J.J.H. Proteomics profiling of chikungunya-infected *Aedes albopictus* C6/36 cells reveal important mosquito cell factors in virus replication. *PLoS Negl. Trop. Dis.* **2015**, *9*, e0003544. [CrossRef]

90. Patramool, S.; Surasombatpattana, P.; Luplertlop, N.; Sévéno, M.; Choumet, V.; Thomas, F.; Missé, D. Proteomic analysis of an *Aedes albopictus* cell line infected with Dengue serotypes 1 and 3 viruses. *Parasit. Vectors* **2011**, *4*, 138. [CrossRef]

91. Tchankouo-Nguetcheu, S.; Khun, H.; Pincet, L.; Roux, P.; Bahut, M.; Huerre, M.; Guette, C.; Choumet, V. Differential protein modulation in midguts of *Aedes aegypti* infected with chikungunya and dengue 2 viruses. *PLoS ONE* **2010**, *5*, e13149. [CrossRef] [PubMed]

92. Ritter, J.B.; Wahl, A.S.; Freund, S.; Genzel, Y.; Reichl, U. Metabolic effects of influenza virus infection in cultured animal cells: Intra- and extracellular metabolite profiling. *BMC Syst. Biol.* **2010**, *4*, 61. [CrossRef] [PubMed]

93. Mondotte, J.A.; Gausson, V.; Frangeul, L.; Suzuki, Y.; Vazeille, M.; Mongelli, V.; Blanc, H.; Failloux, A.-B.; Saleh, M.-C. Evidence for long-lasting transgenerational antiviral immunity in insects. *Cell Rep.* **2020**, *33*, 108506. [CrossRef] [PubMed]

94. Halbach, R.; Miesen, P.; Joosten, J.; Taşköprü, E.; Rondeel, I.; Pennings, B.; Vogels, C.B.F.; Merkling, S.H.; Koenraadt, C.J.; Lambrechts, L.; et al. A satellite repeat-derived piRNA controls embryonic development of *Aedes*. *Nature* **2020**, *580*, 274–277. [CrossRef] [PubMed]

95. Perera, R.; Riley, C.; Isaac, G.; Hopf-Jannasch, A.S.; Moore, R.J.; Weitz, K.W.; Pasa-Tolic, L.; Metz, T.O.; Adamec, J.; Kuhn, R.J. Dengue virus infection perturbs lipid homeostasis in infected mosquito cells. *PLoS Pathog.* **2012**, *8*, e1002584. [CrossRef]

96. Vial, T.; Tan, W.-L.; Deharo, E.; Missé, D.; Marti, G.; Pompon, J. Mosquito metabolomics reveal that dengue virus replication requires phospholipid reconfiguration via the remodeling cycle. *Proc. Natl. Acad. Sci. USA* **2020**, *117*, 27627–27636. [CrossRef]

97. Gottipati, K.; Woodson, M.; Choi, K.H. Membrane binding and rearrangement by chikungunya virus capping enzyme nsP1. *Virology* **2020**, *544*, 31–41. [CrossRef] [PubMed]

98. Shrinet, J.; Jain, S.; Jain, J.; Bhatnagar, R.K.; Sunil, S. Next generation sequencing reveals regulation of distinct *Aedes* microRNAs during chikungunya virus development. *PLoS Negl. Trop. Dis.* **2014**, *8*, e2616. [CrossRef] [PubMed]

99. Campbell, C.L.; Harrison, T.; Hess, A.M.; Ebel, G.D. MicroRNA levels are modulated in *Aedes aegypti* after exposure to Dengue-2. *Insect Mol. Biol.* **2014**, *23*, 132–139. [CrossRef]

100. Saldaña, M.A.; Etebari, K.; Hart, C.E.; Widen, S.G.; Wood, T.G.; Thangamani, S.; Asgari, S.; Hughes, G.L. Zika virus alters the microRNA expression profile and elicits an RNAi response in *Aedes aegypti* mosquitoes. *PLoS Negl. Trop. Dis.* **2017**, *11*, e0005760. [CrossRef]
101. Gainetdinov, I.; Colpan, C.; Arif, A.; Cecchini, K.; Zamore, P.D. A single mechanism of biogenesis, initiated and directed by PIWI proteins, explains piRNA production in most animals unified piRNA biogenesis model. *Mol. Cell* **2018**, *71*, 775–790.e5. [CrossRef] [PubMed]

102. Miesen, P.; Ivens, A.; Buck, A.H.; van Rij, R.P. Small RNA profiling in dengue virus 2-infected *Aedes* mosquito cells reveals viral piRNAs and novel host miRNAs. *PLoS Negl. Trop. Dis.* **2016**, *10*, 1–22. [CrossRef]

103. Sánchez-Vargas, I.; Scott, J.C.; Poole-Smith, B.K.; Franz, A.W.E.; Barbosa-Solomieu, V.; Wilusz, J.; Olson, K.E.; Blair, C.D. Dengue virus type 2 infections of *Aedes aegypti* are modulated by the mosquito's RNA interference pathway. *PLoS Pathog.* **2009**, *5*, e1000299. [CrossRef] [PubMed]

104. Scott, J.C.; Brackney, D.E.; Campbell, C.L.; Bondu-Hawkins, V.; Hjelle, B.; Ebel, G.D.; Olson, K.E.; Blair, C.D. Comparison of dengue virus type 2-specific small RNAs from RNA interference-competent and -incompetent mosquito cells. *PLoS Negl. Trop. Dis.* **2010**, *4*, e848. [CrossRef] [PubMed]

105. Hess, A.M.; Prasad, A.N.; Ptitsyn, A.; Ebel, G.D.; Olson, K.E.; Barbacioru, C.; Monighetti, C.; Campbell, C.L. Small RNA profiling of Dengue virus-mosquito interactions implicates the PIWI RNA pathway in anti-viral defense. *BMC Microbiol.* **2011**, *11*, 24–30. [CrossRef] [PubMed]

106. Nag, D.K.; Kramer, L.D. Patchy DNA forms of the Zika virus RNA genome are generated following infection in mosquito cell cultures and in mosquitoes. *J. Gen. Virol.* **2017**, *98*, 2731–2737. [CrossRef] [PubMed]