Virome and nrEVEome diversity of Aedes albopictus mosquitoes from La Reunion Island and China
Virology Journal volume 19, Article number: 190 (2022)
Aedes albopictus is a public health threat for its worldwide spread and ability to transmit arboviruses. Understanding mechanisms of mosquito immunity can provide new tools to control arbovirus spread. The genomes of Aedes mosquitoes contain hundreds of nonretroviral endogenous viral elements (nrEVEs), which are enriched in piRNA clusters and produce piRNAs, with the potential to target cognate viruses. Recently, one nrEVE was shown to limit cognate viral infection through nrEVE-derived piRNAs. These findings suggest that nrEVEs constitute an archive of past viral infection and that the landscape of viral integrations may be variable across populations depending on their viral exposure.
We used bioinformatics and molecular approaches to identify known and novel (i.e. absent in the reference genome) viral integrations in the genome of wild collected Aedes albopictus mosquitoes and characterize their virome.
We showed that the landscape of viral integrations is dynamic with seven novel viral integrations being characterized, but does not correlate with the virome, which includes both viral species known and unknown to infect mosquitoes. However, the small RNA coverage profile of nrEVEs and the viral genomic contigs we identified confirmed an interaction among these elements and the piRNA and siRNA pathways in mosquitoes.
Mosquitoes nrEVEs have been recently described as a new form of heritable, sequence-specific mechanism of antiviral immunity. Our results contribute to understanding the dynamic distribution of nrEVEs in the genomes of wild Ae. albopictus and their interaction with mosquito viruses.
Following its rapid global spread, the Asian tiger mosquito Aedes albopictus has emerged as a serious public health threat for its ability to transmit over twenty human pathogenic viruses, including dengue, chikungunya and Zika viruses, and the nematode Dirofilaria immitis [1, 2]. Besides arboviruses, mosquitoes can be infected with viruses that are not able to replicate in vertebrate cells and are thus called insect-specific viruses (ISVs) [3, 4]. ISVs are phylogenetically related to arboviruses, for instance, the Flaviviridae family contains two ISVs groups called classical Insect-Specific Flaviviruses (cISFVs) and dual-host ISFVs, respectively. Dual-host ISFVs are phylogenetically nested within mosquito-borne viruses, whereas cISFs are a monophyletic group thought to be an ancestral lineage of Flaviviruses . ISVs are also abundant in the Rhabdoviridae, Bunyaviridae and Reoviridae families and are present, but have not been so frequently identified yet, in the Togaviridae family . Phylogenetic analyses within the Bunyaviridae family support the hypothesis that arboviruses of this viral family evolved from ISVs by overcoming the barriers preventing infection of vertebrates. ISVs have also been proposed as new biological control agents against arboviruses based on the observation that some ISVs, such as Eilat Virus in Ae. albopictus, alter mosquito vector competence by upregulating the antiviral immune responses or through superinfection exclusion . Despite the recognized significance of ISVs, knowledge on the virome diversity of wild-collected Ae. albopictus mosquitoes is still limited [8,9,10,11]. Aedes albopictus populations are usually classified into native, old and invasive . Native populations are endemic in the species home range in Southeast Asia; old populations can be found in the islands of the Indian Ocean, which were colonized in the ‘800 century, and invasive are populations resulting from overlapping and chaotic incursions which took place in both tropical and temperate regions of the world in the past 50–60 years .
Both ISVs and arboviral infections elicit RNA interference (RNAi) pathways in mosquitoes. RNAi depends on the production of small (s)RNAs, which are classified on the basis of their size, biogenesis and interacting proteins into microRNAs of ~ 20–22 bp; small interfering (si)RNAs of ~ 20–22 bp, and PIWI interacting (pi)RNAs of ~ 25–30 bp . miRNAs derive from hairpin-like structure forming in long single-stranded RNA molecule, siRNAs are processed from long double stranded RNAs, while piRNAs are produced both in the nucleus (primary piRNAs) from ad hoc genomic regions called piRNA clusters and in the cytoplasm (secondary piRNAs) from self-amplification of primary piRNAs after target recognition and processing . While during viral infection of mosquitoes, siRNAs are consistently found and homogeneously cover infecting viral genomes, production of virus-derived miRNAs is controversial and the profile of piRNAs appears to be tissue-specific and both viral and host-specific [15,16,17]. On this basis, the pattern of small RNAs, in primis siRNAs, has been used to reconstruct the sequences of viruses infecting mosquitoes and other animals [16, 18, 19].
piRNA clusters of Aedes mosquitoes, but not those of Drosophila melanogaster nor Anophelinae, are enriched in non-retroviral RNA virus endogenous elements (nrEVEs) (20,21,22,23,24). The genome of Ae. albopictus hosts a total of 456 nrEVEs . The most represented viral families from which nrEVEs originate are Rhabdoviridae (Rhabdo-nrEVEs), Flaviviridae (Flavi-nrEVEs) and the newly discovered Chuviridae [20, 23, 25], besides unclassified viruses identified in arthropods but whose host range is still unknown. Aedes albopictus nrEVEs produce piRNAs which are antisense with respect to cognate viruses suggesting antiviral activity . Reduction of cognate viral infection due to nrEVE-derived piRNAs has been shown for selected nrEVEs of both Ae. aegypti and Aag2 cells [6, 26, 27]. Despite an antiviral activity of nrEVEs have not been proven in Ae. albopictus, results in Ae. aegypti suggest that nrEVEs are involved in shaping the relationships among mosquitoes, ISVs and, possibly, arboviruses. However, information on the co-occurrence between nrEVEs and ISVs in geographic mosquito populations is lacking. To address this knowledge gap, we sampled Ae. albopictus mosquitoes from the native home range in Southern China and in La Reunion Island  and combined whole genome (WG) and small RNA sequencing approaches to characterize the concomitant patterns of viral integrations and their the small RNA profile while also detecting the mosquito virome.
Materials and methods
Wild samples collection and sequencing
Adult Ae. albopictus individuals were collected in public areas in La Reunion Island (France) in the Guangzhou prefecture (China), which did not require ethical approval (Fig. 1). In La Reunion Island, sampling sites included sites from the Eastern (Bras-Panon [RE1, RE2] and Saint-Rose [RE3, RE4]) and the western (Tampon [RW1, RW2] and SaintPierre [RW3, RW4]) regions, which experience low and high incidence of arboviral disease transmission, respectively . In the Guangzhou prefecture, we sampled mosquitoes from the Baiyun district (GZU1 sample), which is located at the northern outskirts of prefecture and is a forested area known for its natural attractions and from the campus of the Southern Medical University (GZU2 sample), an urban and densely populated area of the Guangzhou city, which is subjected to intense vector control interventions .
Mosquitoes collected in La Reunion Island were immediately processed in pools of 10 females and homogenized in DNA/RNA-shield (Zymo Research, Irvine CA, USA) with a Kimble electric pestle (DWK Life Sciences, Mainz, Germany). Samples were stored at 4 °C before being transferred to the University of Pavia where each pool was divided into two volumes that were used for either DNA or RNA extraction. DNA was extracted from mosquito pools using the DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany), according to manufacturer’s protocol and stored at − 20 °C. RNA was extracted using a standard protocol based on Trizol (Ambion/Invitrogen, Waltham MA, USA), treated with DNase I (Sigma-Aldrich, St. Louis MO, USA) and stored at − 80 °C. Nucleic acids from three pools were further pooled to generate a sample. For each locality we processed two samples, which were sent to Biodiversa (Rovereto, Italy) for library preparation and sequencing resulting in the analyses of a total of 60 mosquitoes per site (Additional file 1: Table S1). WG and Small RNA libraries were prepared using the Nextera DNA Library Preparation Kit for paired end reads (Illumina, San Diego CA, USA) and the NEBNext Small RNA Library Prep Set (New England Biolabs, Ipswich MA, USA), respectively. DNA and small-RNA Libraries were sequenced on an Illumina HiSeq 2500 instrument.
Samples from China were also processed in pools of 10 mosquitoes each. DNA was extracted using the Universal Genomic DNA extraction kit (Takara, Kusatsu, Japan) according to the supplier protocol and stored at − 20 °C. RNA was extracted using a custom protocol based on Trizol (Ambion/Invitrogen, Waltham MA, USA), treated with DNase I (Takara, Kusatsu, Japan) and stored at − 80 °C. DNA and RNA. As for samples from La Reunion, 3 pools of the nucleic acids of 10 mosquitoes each were further pooled and sent to the Beijing Genomics Institute (BGI), China, for library preparation and sequencing (Additional file 1: Table S1). DNA WGS libraries were prepared with the Nextera DNA Library Preparation Kit (Illumina, San Diego CA, USA) and were sequenced on an Illumina HiSeq X instrument. Small RNA libraries were prepared with the TruSeq Small RNA Library Prep Kit (Illumina, San Diego CA, USA) and sequenced on an Illumina HiSeq 4500 sequencer.
WGS data were analyzed to identify the landscape of nrEVEs of wild-collected mosquitoes using as reference the nrEVEs that we previously annotated in the Ae. albopictus reference genome (AalbF2) . Reads were checked and cleaned with FastQC v. 0.11.9 (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and trimmomatic v. 0.40 , respectively, and mapped on the AalbF2 genome using the SVD pipeline as previously described. Given that we sequenced pools of 30 individuals, we estimated the frequency of nrEVEs in each pool based on read coverage. Briefly, single copy orthologs were identified in the AalbF2 genome assembly using BUSCO v5.3.1 , configured with the diptera_odb10 database. Normalized depth values of a selected set of 359 orthologs, which satisfied a minimum BUSCO score of 1000 (Additional file 2: Table S2), and of all the reference nrEVEs were extracted using FeatureCounts in the SubRead package . Average normalized depth for single copy orthologs was calculated for each sample and used as a reference to estimate nrEVE allelic frequency in pooled samples (Additional file 3: Fig. S1). nrEVEs were considered “present” in a sample when the estimated allelic frequency was higher than 0.2. The Principal Component Analysis (PCA) classifying our samples based on their pattern of nrEVEs was computed using the Exponential family PCA method in the LogisticPCA R  package and plotted with ggplot2.
Characterization of novel viral integration sites in wild collected mosquitoes
The Vy-PER and ViR pipelines were run as previously described  to identify in the DNA of wild-collected mosquitoes nrEVEs which are absent in the AalbF2 assembly. Bioinformatic predictions of each of the newly identified viral integrations were molecularly confirmed by PCR using the DreamTaq Green PCR Master Mix 2X (Thermo Fisher) and nrEVE specific primers (Additional file 4: Table S3). PCR products were cloned into E. Coli OneShot TOP10 chemically competent cells (Invitrogen) using the TOPO TA Cloning Kit for Sequencing (Thermo Fisher, USA). Plasmids containing the PCR product were extracted with the QIAprep Spin Miniprep Kit (Qiagen, Germany). Plasmid products were Sanger sequenced by Macrogen Europe (Amsterdam, The Netherlands) using the M13F and M13R plasmid primers (Additional file 5).
The presence of the newly identified nrEVEs was also tested by PCR on single mosquitoes of laboratory strains derived from mosquitoes sampled in Canton (China), Chiang Mai (Thailand), La Reunion Island, Crema (Italy) and Tapachula (Mexico). DNA was extracted from single mosquitoes using the Wizard genomic DNA purification kit (Promega) and PCRs were performed as described above.
Raw small-RNA (sRNA) sequencing data were checked for quality with FastQC v. 0.11.9 (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/); the presence of adapters was confirmed with DNApi v. 1.1 . Adapters and bases with baseQ inferior to 20 were trimmed from the reads using Cutadapt v. 2.9 , also removing clean reads smaller than 10 bp and longer than 50 bp. Low complexity reads were excluded with the Duster.pl script from the NGS-toolbox package (https://sourceforge.net/p/ngs-toolbox/wiki/Home/). Briefly, reads were mapped to the Ae. albopictus genome (AalbF2 assembly) after masking nrEVEs using BEDtools maskfasta v. 2.28 . Small-RNA reads were mapped on the masked genome with Bowtie v. 1.1.2  optimizing the parameters for small-RNA reads (-n 1 -l 18 -best). The reads that did not map to the Ae. albopictus genome were extracted with SAMtools view v. 1.4  and converted to fastq with SAMtools fastq. Contigs were reconstructed from unmapped sRNA reads using the Oases pipeline  with Velvet v.1.2.10 (https://www.ebi.ac.uk/~zerbino/velvet/) testing k-mer lengths from 13 to 31 bp. Sequences assembled by Oases and longer than 100 bp were filtered to remove nrEVEs-derived assemblies using BLASTn v. 2.6.0  by comparing all the assembled sequences against nrEVEs annotated in AalbF2 . Contigs with a percentage of identity higher than 90% and covering at least 50% of an existing nrEVE were considered nrEVEs assemblies rather than viral contigs. Redundant assemblies were clustered using CD-HIT v. 4.8.0  with the following options: -c 0.9 -n 8 -d 60 -g 1. The clustered assemblies were screened against the entire NCBI nr peptides database (https://www.ncbi.nlm.nih.gov/refseq/about/nonredundantproteins/, downloaded in march 2020) using DIAMOND v. 0.9.31  BLASTx algorithm with an e-value cutoff of 1e−5. DIAMOND hits were loaded on MEGAN6 community edition v. 6.19.9 . MEGAN6 was run with the naive LCA algorithm. Assemblies assigned to viruses by MEGAN6 were manually validated against the nr database using BLASTx .
Small RNA profile of viral contigs and novel nrEVEs
Selected assembled viral contigs, including viral contigs of at least 1500 bp and Flavivirus-like viral contigs shorter than 500 bp, were used for the realignment of sRNA reads from the samples in which the contigs were detected. Clean sRNA reads were mapped using Bowtie v. 1.1.2  as described above. Sequences between 20–22 bp and 25–30 bp were filtered as siRNAs/miRNAs and piRNAs, respectively, using BBMap reformat (sourceforge.net/projects/bbmap/). Small RNAs that mapped to viral contigs (hereafter called viral sRNA) were then re-mapped to the AalbF2 genome assembly, then split by orientation with samtools view and counted with BBMap pileup.sh (sourceforge.net/projects/bbmap/). The reads were categorized as either mapping exclusively on viral contigs or mapping on both viral contigs and nrEVEs. Reads size distributions were plotted on Python 3 using the Plotly v. 5.7.0 graphing library .
To understand the siRNA coverage of the novel viral integrations identified in wild-collected mosquitoes, the sRNA reads from each mosquito pool in which the viral integration had been identified were re-mapped to the viral integration using the methodology described above. Coverage profile of mapped piRNA and siRNA was plotted with GraphPad Prism 8 (www.graphpad.com).
Adult Ae. albopictus mosquitoes were collected from four localities from La Reunion Island and two sites in the Guangzhou prefecture (China) (Fig. 1). From all sites, WGS and small RNA data were concomitantly produced to characterize nrEVE patterns, test for the presence of new viral integrations, estimate the sRNA profile of nrEVEs and assemble viral contigs.
Patterns of nrEVEs across samples
Out of the 456 nrEVEs annotated in the Ae. albopictus genome, 222 nrEVEs appear to be absent in all samples, while 234 nrEVEs were estimated to be present based on read coverage with respect to mosquito single-copy ortholog genes. The majority of nrEVEs, which were not detected in any samples, included Chuviridae-like nrEVEs and nrEVEs from unclassified viruses followed by Rhabdovirus-like nrEVE. To determine if nrEVE patterns differ across samples based on their geographic origin, we used PCA. While results showed a clear separation between samples from the Guangzhou prefecture in China and La Reunion Island (Fig. 2a), within La Reunion, samples from the East and West coast did not group separately, consistent with high gene flow among populations of the island .
We further used the Vy-PER and ViR pipelines to align WGS data to the Ae. albopictus genome (AalbF2 assembly), identify unmapped reads containing viral sequences and characterize novel integration sites. A total of 7 new integrations were identified: 6 were detected in samples from La Reunion Island and one in mosquitoes from the Guangzhou prefecture (Table 1, Fig. 2b, Additional file 1: Table S1). All newly identified nrEVEs derived from ISVs of the Flaviviridae family, including Kamiti River Virus (KRV), Cell Fusing Agent Virus (CFAV) and Aedes Flavivirus, with a percentage of nucleotide identity ranging from 67 to 78%. Additionally, these nrEVEs, which range in size between 814 and 4705 bp, encompass one or more viral proteins without breaks or rearrangement with respect to the original viral genome. Each of the newly identified nrEVE was molecularly validated by PCR using specific primers followed by Sanger sequencing (Additional file 4: Table S3). Additional wild-caught mosquitoes from La Reunion Island and Canton (China), besides samples from Chiang Mai (Thailand), Crema (Italy) and Tapachula (Mexico) were used to assess the frequency of the newly identified nrEVE (Fig. 2c). Apart from KRV-1 and KRV-3 which were detected only in RE2 and GZU1, respectively, all other newly identified nrEVEs were found in mosquitoes from both native and invasive populations, at frequency between 4 and 90% (Fig. 2c).
Since it is known that viral integrations are enriched in piRNA clusters and may produce piRNAs, we checked if, and where, piRNAs and siRNAs mapped on the newly detected nrEVE. Small RNAs in the size of piRNA, but not siRNAs were seen mapping across the novel viral integrations; coverage was not homogenous across the whole viral integration, but the overall profile was similar across samples (Fig. 3).
The virome of Ae. albopictus mosquitoes
A total of 122 viral contigs corresponding to 25 viral species were assembled from small RNAs of mosquitoes collected in the Guangzhou prefecture. A total of 126 viral contigs corresponding to 26 viral species were assembled from small RNAseq data of samples from La Reunion (Fig. 4). Viral contigs ranged in size between 101 and 5847 bp and had a highly variable percentage of nucleotide identity (between 33 and 100%, with a median of 70%) with respect to known circulating viruses (Fig. 4, Additional file 6). Most assembled contigs had similarities to bona fide arthropod viruses and included ISVs of the Flaviviridae family (i.e. CFAV, KRV, La Tina virus, and Aedes Flavivirus [AeFV]) [5, 7, 46], the Mononegavirales order (Culex tritaeniorhynchus rhabdovirus [CTRV]  and Merida virus  and Aedes Anphenvirus  from the Rhabdoviridae and Xinmonoviridae families, respectively); Yongsang bunyavirus 1, previously identified in Aedes vexans mosquitoes from Korea , Phasi-Charoen Like Virus (family Phenuiviridae), which is highly prevalent in both wild-caught Ae. aegypti and Ae. albopictus mosquitoes , along with less known viruses, which have been recently characterized from insect samples, such as Wuhan Mosquito Virus, Wenzhou sobemo-like virus, Guadalupe mosquito virus, Hubei mosquito virus, Chuvirus Mos8Chu0, Croada virus and Kaiowa virus [11, 52, 53]. We did not assemble any contigs from arboviruses, except for a contig of 168 bp with sequence similarity to a member of the Vesicular stomatitis virus serogroup of the Rhabdoviridae family, the Yug Bogdanovac virus (YBV), which was initially isolated from sandflies , and we assembled from small RNAs of the RE1 sample from the Eastern side of La Reunion.
Because viral replication elicits RNAi and results in a small RNAs profile skewed towards siRNAs that target both strands of a viral genome, we checked this size signature on the viral species for which we assembled contigs of at least 1500 bp, including Sarawak virus; Aedes aegypti toti-like virus; Aedes aegypti virga-like virus; Guadeloupe mosquito mononega-like virus,; Chuvirus Mos8Chu0, Point-Douro narna-like virus, Aedes albopictus negev-like virus and Shinobi tetravirus with a contig of 3450 bp (Fig. 5). Importantly, none of the assembled viral contigs matched to Ae. albopictus nrEVEs or newly identified nrEVEs. However, because endogenous piRNAs can trigger production of trailer and responder piRNAs from viral RNAs , we realigned piRNAs from the viral contigs to nrEVEs and tagged them as possible nrEVEs-derived small RNAs. A profile indicative of active viral infection was observed for the Sarawak, Shinobi tetravirus, the Aedes albopictus negev-like and, possibly, the Guadalupe mosquito mononega-like virus (Fig. 5a). On the contrary the smallRNA profile of the Aedes aegypti toti-like virus was skewed towards piRNAs, some of which may also derive from viral integrations despite the absence of sequence complementarity between the viral contig and nrEVEs . Furthermore, the profiles of contigs classified as Chuvirus Mos8Chu0, Aedes aegypti toti-like virus and Aedes aegypti virga-like virus were skewed to small RNAs of exclusively one strand and included smallRNAs which also map to nrEVEs. Because piRNAs against a single orientation of the viral genome with no siRNAs is a classic signature of nrEVEs, whereas siRNAs of positive and negative strand polarity are typical for viral DNA (vDNA) fragments, which have been shown to be formed after mosquito infection by a wide range of RNA viruses, possibly from their defective viral genomes [56,57,58,59], we cannot exclude that the smallRNA pattern observed for the Aedes aegypti toti-like virus contig represent an undetected novel nrEVE, whereas the profiles of the Chuvirus Mos8Chu0 and the Aedes aegypti virga-like virus are difficult to interpret. We also verified the profile of assembled flaviviral contigs shorter than 500 bp, which may represent vDNA fragments of flaviviruses or viral integrations which we were not able to identify from WGS data. The observed profiles, with a prevalence of small RNAs in the size of piRNAs of only one orientation (Fig. 5b), point to nrEVEs, which are probably rare in the pool samples analyzed thus resulting in DNA reads below the threshold for novel nrEVE detection.
As our analyses were being completed, independent groups of investigators used a metagenomic approach to describe the virome of mosquitoes collected between 2017 and 2018 from urban Guangzhou  and that of Aedes mosquitoes using published RNA-seq datasets [9, 11]. A total of 42% of the viral contigs we identified had already been detected in world-wide Ae. albopictus samples (Fig. 3), further validating small RNA seq based virome identification and suggesting the existence of a core virome.
Eukaryotic genomes have integrated viral sequences, most of which derive from retroviruses encoding for retrotranscriptase and integrases. Previous work conducted by us, and others, revealed the presence of nrEVEs with similarity to, mostly, flaviviruses and rhabdoviruses in Ae. albopictus and Ae. aegypti, but not in the genomes of Anopheles species and Culex quinquefasciatus . Additionally, we have shown that the pattern of nrEVEs varies across geographic populations, with new integrations being detectable in the genome of wild-collected mosquitoes [20, 23]. Here we report the identification of 7 novel nrEVEs in wild-collected Ae. albopictus mosquitoes. On average, these nrEVEs have a higher sequence identity to corresponding viral genomes in comparison to reference nrEVEs, suggesting more recent integration events. Newly identified Ae. albopictus nrEVEs have similarities to a range of cISFs, including CFAV, AeFV and KRV, unlike the new nrEVEs identified in Ae. aegypti wild mosquitoes . Additionally, newly identified nrEVEs appear broadly distributed across mosquitoes from different locations, probably because of the quick and chaotic global invasion of Ae. albopictus . Of note CFAV-2 and KRV-1 were exclusive of mosquitoes from Italy and La Reunion Island, among the sites we tested, suggesting that these integrations occurred after Ae. albopictus established on La Reunion island, which is a primary source of mosquito introductions to Italy . Small RNA profiles of newly identified nrEVEs confirm the presence of low abundance piRNAs, and not siRNAs, which appear in hotspots that are conserved across samples similarly to what is observed for reference nrEVEs of both Ae. aegypti and Ae. albopictus [20, 23]. Overall, these results indicate that, independently of the timing of the integration event, piRNAs will derive from nrEVEs, but not throughout the nrEVE sequence. Thus, endeavors aiming at using nrEVEs in strategies for vector control either as landmarks for effector integration or taking advantage of their inherent antiviral activity against cognate viruses might require experimental validation of their piRNA profile prior transgenic applications .
Given the abundance of nrEVEs in the genome of Aedes mosquitoes, we further asked whether the presence of newly identified nrEVEs correlates with on-going viral infections and analyzed the virome of the same mosquitoes in which we tested for novel nrEVEs. While our samples could be grouped similarly based both on their nrEVE landscape and their virome (Fig. 2), indicating that sampling location is the main factor influencing both virome and nrEVEome, we did not assemble viral contigs from the same viruses for which novel nrEVEs were identified in any sample (Figs. 3, 4). These results along with the number of newly identified nrEVEs with respect to the number of wild mosquitoes we analyzed show that integration events are rare. Whether the load of infecting viruses or other ecological and/or biological circumstances such as thermal stress, co-infections with bacteria or fungi and/or the presence of active transposable elements are needed for an integration event to occur requires further investigations.
Overall, our small-RNA-based virome analysis led to the identification of 43 viruses; 18 of these viral species have already been described in Aedes spp. mosquitoes, with Phasi Charoen-like phasivirus, Guadeloupe mosquito monomega-like virus, Aedes anphevirus, toti-like viruses and Chuvirus Mos8Chu0 being the most prevalent [9, 11, 60]. Detailed analysis of the smallRNA profiles of the assembled viral contigs showed a profile with 21 bp peaks in sense and antisense orientations indicative of active replication only for few viruses, such as Sarawak virus, Shinobi tetravirus and the Aedes albopictus negev-like virus. Consequently, we cannot exclude that some of the assembled viral contigs are simply related to the diet or are associated with mosquito microbiota. For instance, we assembled contigs with similarities to Virgaviridae, which are plant viruses, and Totiviridae for which fungi and protozoans are natural hosts, in both mosquitoes from China and La Reunion Island. Sarawak virus was first identified in 2017 from Ae. albopictus males collected in 2013 on the Sarawak island of Borneo and its genome sequence has been produced showing it is an Alphatetraviridae with three open reading frames (ORFs) . Shinobi tetravirus was discovered in 2018 in Ae. albopictus C6/36 cells and its presence was shown to suppress replication of Zika virus. Shinobi tetravirus readily infects Ae. albopictus strains in laboratory conditions , but its prevalence in wild mosquitoes is uncertain. Aedes albopictus negev-like virus was discovered in 2020 in the Ae. albopictus Aa23 cell line, where Wolbachia wAlbB was shown to negatively affect its replication .
Overall, our results highlight the need not only to continue this type of studies to further understand the diversity of viruses infecting mosquitoes, but also to start introducing a functional characterization of these viruses and the effects of their infections on mosquitoes’ fitness and metabolism.
Distinct small RNA profiles were observed for additional viral contigs, including Chuvirus Mos8Chu0, which warrant further investigations. Chuviridae are a poorly characterized group of viruses, which have been recently identified through metagenomic analyses in several insects from different orders including Hemiptera, Diptera and Coleoptera [67, 68]. Insects in which Chuviridae have been identified tend to also harbor Chuviridae-like nrEVEs [25, 67], thus these viruses may represent an ideal model to study the arm-race with insects resulting in integration events.
Availability of data and materials
Raw sequences will be available under NIH BioProject number PRJNA866864, upon acceptance of this manuscript for publication. Biosamples numbers are listed in Additional file 1: Table S1.
Bonizzoni M, Gasperi G, Chen X, James AA. The invasive mosquito species Aedes albopictus: current knowledge and future perspectives. Trends Parasitol. 2013;29:460–8.
Little EAH, Harriott OT, Akaratovic KI, Kiser JP, Abadam CF, Shepard JJ, Molaei G. Host interactions of Aedes albopictus, an invasive vector of arboviruses, in Virginia, USA. PLoS Negl Trop Dis. 2021;15:e0009173.
Huhtamo E, Cook S, Moureau G, Uzcátegui NY, Sironen T, Kuivanen S, Putkuri N, Kurkela S, Harbach RE, Firth AE, Vapalahti O, Gould EA, de Lamballerie X. Novel flaviviruses from mosquitoes: mosquito-specific evolutionary lineages within the phylogenetic group of mosquito-borne flaviviruses. Virology. 2014;464–465:320–9.
Agboli L, Altinli S. Mosquito-specific viruses—transmission and interaction. Viruses. 2019;11:873.
Cook S, Moureau G, Kitchen A, Gould EA, de Lamballerie X, Holmes EC, Harbach RE. Molecular evolution of the insect-specific flaviviruses. J Gen Virol. 2012;93:223–34.
Olson KE, Bonizzoni M. Nonretroviral integrated RNA viruses in arthropod vectors: an occasional event or something more? Curr Opin Insect Sci. 2017;22:45–53.
Carvalho VL, Long MT. Insect-specific viruses: an overview and their relationship to arboviruses of concern to humans and animals. Virology. 2021;557:34–43.
de Almeida JP, Aguiar ER, Armache JN, Olmo RP, Marques JT. The virome of vector mosquitoes. Curr Opin Virol. 2021;49:7–12.
Parry R, James ME, Asgari S. Uncovering the worldwide diversity and evolution of the virome of the mosquitoes Aedes aegypti and Aedes albopictus. Microorganisms. 2021;9:1653.
Kubacki J, Flacio E, Qi W, Guidi V, Tonolla M, Fraefel C. Viral metagenomic analysis of Aedes albopictus mosquitos from Southern Switzerland. Viruses. 2020;12:E929.
Shi C, Zhao L, Atoni E, Zeng W, Hu X, Matthijnssens J, Yuan Z, Xia H. Stability of the virome in lab- and field-collected Aedes albopictus mosquitoes across different developmental stages and possible core viruses in the publicly available virome data of Aedes mosquitoes. Systems. 2020;5:e00640-20.
Manni M, Guglielmino CR, Scolari F, Vega-Rúa A, Failloux A-B, Somboon P, Lisa A, Savini G, Bonizzoni M, Gomulski LM, Malacrida AR, Gasperi G. Genetic evidence for a worldwide chaotic dispersion pattern of the arbovirus vector, Aedes albopictus. PLoS Negl Trop Dis. 2017;11:e0005332.
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.
Miesen P, Joosten J, van Rij RP. PIWIs go viral: arbovirus-derived piRNAs in vector mosquitoes. PLoS Pathog. 2016;12:e1006017.
Goic B, Vodovar N, Mondotte JA, Monot C, Frangeul L, Blanc H, Gausson V, Vera-Otarola J, Cristofari G, Saleh M-C. RNA-mediated interference and reverse transcription control the persistence of RNA viruses in the insect model Drosophila. Nat Immunol. 2013;14:396–403.
Aguiar ERGR, Olmo RP, Marques JT. Virus-derived small RNAs: molecular footprints of host–pathogen interactions. WIREs RNA. 2016;7:824–37.
Parameswaran P, Sklan E, Wilkins C, Burgon T, Samuel MA, Lu R, Ansel KM, Heissmeyer V, Einav S, Jackson W, Doukas T, Paranjape S, Polacek C, dos Santos FB, Jalili R, Babrzadeh F, Gharizadeh B, Grimm D, Kay M, Koike S, Sarnow P, Ronaghi M, Ding S-W, Harris E, Chow M, Diamond MS, Kirkegaard K, Glenn JS, Fire AZ. Six RNA viruses and forty-one hosts: viral small RNAs and modulation of small RNA repertoires in vertebrate and invertebrate systems. PLoS Pathog. 2010;6:e1000764.
Webster CL, Waldron FM, Robertson S, Crowson D, Ferrari G, Quintana JF, Brouqui J-M, Bayne EH, Longdon B, Buck AH, Lazzaro BP, Akorli J, Haddrill PR, Obbard DJ. The discovery, distribution, and evolution of viruses associated with Drosophila melanogaster. PLoS Biol. 2015;13:e1002210.
Kreuze JF, Perez A, Untiveros M, Quispe D, Fuentes S, Barker I, Simon R. Complete viral genome sequence and discovery of novel viruses by deep sequencing of small RNAs: a generic method for diagnosis, discovery and sequencing of viruses. Virology. 2009;388:1–7.
Palatini U, Masri RA, Cosme LV, Koren S, Thibaud-Nissen F, Biedler JK, Krsticevic F, Johnston JS, Halbach R, Crawford JE, Antoshechkin I, Failloux A-B, Pischedda E, Marconcini M, Ghurye J, Rhie A, Sharma A, Karagodin DA, Jenrette J, Gamez S, Miesen P, Masterson P, Caccone A, Sharakhova MV, Tu Z, Papathanos PA, Van Rij RP, Akbari OS, Powell J, Phillippy AM, Bonizzoni M. Improved reference genome of the arboviral vector Aedes albopictus. Genome Biol. 2020;21:215.
Palatini U, Miesen P, Carballar-Lejarazu R, Ometto L, Rizzo E, Tu Z, van Rij RP, 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:512.
Suzuki Y, Frangeul L, Dickson LB, Blanc H, Verdier Y, Vinh J, Lambrechts L, Saleh M-C. Uncovering the repertoire of endogenous flaviviral elements in Aedes mosquito genomes. J Virol. 2017;91:e00571-e617.
Crava CM, Varghese FS, Pischedda E, Halbach R, Palatini U, Marconcini M, Gasmi L, Redmond S, Afrane Y, Ayala D, Paupy C, Carballar-Lejarazu R, Miesen P, van Rij RP, Bonizzoni M. Population genomics in the arboviral vector Aedes aegypti reveals the genomic architecture and evolution of endogenous viral elements. Mol Ecol. 2021;30:1594–611.
Whitfield ZJ, Dolan PT, Kunitomi M, Tassetto M, Seetin MG, Oh S, Heiner C, Paxinos E, Andino R. The diversity, structure, and function of heritable adaptive immunity sequences in the Aedes aegypti genome. Curr Biol. 2017;27:3511-3519.e7.
Dezordi FZ, Vasconcelos CRDS, Rezende AM, Wallau GL. In and outs of chuviridae endogenous viral elements: origin of a potentially new retrovirus and signature of ancient and ongoing arms race in mosquito genomes. Front Genet. 2020;11:542437.
Tassetto M, Kunitomi M, Whitfield ZJ, Dolan PT, Sánchez-Vargas I, Garcia-Knight M, Ribiero I, Chen T, Olson KE, Andino R. Control of RNA viruses in mosquito cells through the acquisition of vDNA and endogenous viral elements. Elife. 2019;8:e41244.
Ophinni Y, Palatini U, Hayashi Y, Parrish NF. piRNA-guided CRISPR-like immunity in eukaryotes. Trends Immunol. 2019;40:998–1010.
Latreille AC, Milesi P, Magalon H, Mavingui P, Atyame CM. High genetic diversity but no geographical structure of Aedes albopictus populations in Réunion Island. Parasites Vectors. 2019;12:597.
Li L, Liu W-H, Zhang Z-B, Liu Y, Chen X-G, Luo L, Ou C-Q. The effectiveness of early start of Grade III response to dengue in Guangzhou, China: a population-based interrupted time-series study. PLoS Negl Trop Dis. 2020;14:e0008541.
Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.
Simão FA, Waterhouse RM, Ioannidis P, Kriventseva EV, Zdobnov EM. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics. 2015;31:3210–2.
Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;30:923–30.
Landgraf AJ, Lee Y. Dimensionality reduction for binary data through the projection of natural parameters. J Multivar Anal. 2020;180:104668.
Palatini U, Pischedda E, Bonizzoni M. Computational methods for the discovery and annotation of viral integrations. Methods Mol Biol. 2022;2509:293–313.
Tsuji J, Weng Z. DNApi: a de novo adapter prediction algorithm for small RNA sequencing data. PLoS ONE. 2016;11:e0164228.
Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. 1. EMBnet J. 2011;17:10–2.
Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26:841–2.
Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009;10:R25.
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–9.
Schulz MH, Zerbino DR, Vingron M, Birney E. Oases: robust de novo RNA-seq assembly across the dynamic range of expression levels. Bioinformatics. 2012;28:1086–92.
Camacho C, Coulouris G, Avagyan V, Ma N, Papadopoulos J, Bealer K, Madden TL. BLAST+: architecture and applications. BMC Bioinform. 2009;10:421.
Li W, Godzik A. Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics. 2006;22:1658–9.
Buchfink B, Xie C, Huson DH. Fast and sensitive protein alignment using DIAMOND. Nat Methods. 2015;12:59–60.
Bağcı C, Beier S, Górska A, Huson DH. Introduction to the analysis of environmental sequences: metagenomics with MEGAN. Methods Mol Biol. 2019;1910:591–604.
Sievert C. Interactive web-based data visualization with R, plotly, and shiny. https://plotly-r.com/. Retrieved 1 Aug 2022.
Hoshino K, Isawa H, Tsuda Y, Sawabe K, Kobayashi M. Isolation and characterization of a new insect flavivirus from Aedes albopictus and Aedes flavopictus mosquitoes in Japan. Virology. 2009;391:119–29.
Kuwata R, Isawa H, Hoshino K, Tsuda Y, Yanase T, Sasaki T, Kobayashi M, Sawabe K. RNA splicing in a new rhabdovirus from Culex mosquitoes. J Virol. 2011;85:6185–96.
Charles J, Firth AE, Loroño-Pino MA, Garcia-Rejon JE, Farfan-Ale JA, Lipkin WI, Blitvich BJ, Briese T. Merida virus, a putative novel rhabdovirus discovered in Culex and Ochlerotatus spp. mosquitoes in the Yucatan Peninsula of Mexico. J Gen Virol. 2016;97:977–87.
Parry R, Asgari S. Aedes anphevirus: an insect-specific virus distributed worldwide in Aedes aegypti mosquitoes that has complex interplays with Wolbachia and dengue virus infection in cells. J Virol. 2018;92:e00224-e318.
Sanborn MA, Klein TA, Kim H-C, Fung CK, Figueroa KL, Yang Y, Asafo-adjei EA, Jarman RG, Hang J. Metagenomic analysis reveals three novel and prevalent mosquito viruses from a single pool of Aedes vexans nipponii collected in the Republic of Korea. Viruses. 2019;11:222.
Zhang X, Huang S, Jin T, Lin P, Huang Y, Wu C, Peng B, Wei L, Chu H, Wang M, Jia Z, Zhang S, Xie J, Cheng J, Wan C, Zhang R. Discovery and high prevalence of Phasi Charoen-like virus in field-captured Aedes aegypti in South China. Virology. 2018;523:35–40.
de Lara Pinto AZ, Santos de Carvalho M, de Melo FL, Ribeiro ALM, Morais Ribeiro B, Dezengrini SR. Novel viruses in salivary glands of mosquitoes from sylvatic Cerrado, Midwestern Brazil. PLoS ONE. 2017;12:e0187429.
Atoni E, Zhao L, Hu C, Ren N, Wang X, Liang M, Mwaliko C, Yuan Z, Xia H. A dataset of distribution and diversity of mosquito-associated viruses and their mosquito vectors in China. Sci Data. 2020;7:342.
Pfeffer M, Dilcher M, Tesh RB, Hufert FT, Weidmann M. Genetic characterization of Yug Bogdanovac virus. Virus Genes. 2013;46:201–2.
Joosten J, Overheul GJ, Van Rij RP, Miesen P. Endogenous piRNA-guided slicing triggers responder and trailer piRNA production from viral RNA in Aedes aegypti mosquitoes. Nucleic Acids Res. 2021;49:8886–99.
Poirier EZ, Goic B, Tomé-Poderti L, Frangeul L, Boussier J, Gausson V, Blanc H, Vallet T, Loyd H, Levi LI, Lanciano S, Baron C, Merkling SH, Lambrechts L, Mirouze M, Carpenter S, Vignuzzi M, Saleh M-C. 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.
Nag DK, Kramer LD. 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–7.
Nag DK, Brecher M, Kramer LD. DNA forms of arboviral RNA genomes are generated following infection in mosquito cell cultures. Virology. 2016;498:164–71.
Goic B, Stapleford KA, Frangeul L, Doucet AJ, Gausson V, Blanc H, Schemmel-Jofre N, Cristofari G, Lambrechts L, Vignuzzi M, Saleh M-C. Virus-derived DNA drives mosquito vector tolerance to arboviral infection. Nat Commun. 2016;7:12410.
He W, Chen Y, Zhang X, Peng M, Xu D, He H, Gao Y, Chen J, Zhang J, Li Z, Chen Q. Virome in adult Aedes albopictus captured during different seasons in Guangzhou City, China. Parasit Vectors. 2021;14:415.
Palatini U, Contreras CA, Gasmi L, Bonizzoni M. Endogenous viral elements in mosquito genomes: current knowledge and outstanding questions. Curr Opin Insect Sci. 2022;49:22–30.
Macias VM, Palatini U, Bonizzoni M, Rasgon JL. Leaning into the bite: the piRNA pathway as an exemplar for the genetic engineering need in mosquitoes. Front Cell Infect Microbiol. 2020;10:614342.
Sadeghi M, Popov V, Guzman H, Phan TG, Vasilakis N, Tesh R, Delwart E. Genomes of viral isolates derived from different mosquitos species. Virus Res. 2017;242:49–57.
Fujita R, Kato F, Kobayashi D, Murota K, Takasaki T, Tajima S, Lim C-K, Saijo M, Isawa H, Sawabe K. Persistent viruses in mosquito cultured cell line suppress multiplication of flaviviruses. Heliyon. 2018;4:e00736.
Frangeul L, Blanc H, Saleh M-C, Suzuki Y. Differential small RNA responses against co-infecting insect-specific viruses in Aedes albopictus mosquitoes. 4. Viruses. 2020;12:468.
Bishop C, Parry R, Asgari S. Effect of Wolbachia wAlbB on a positive-sense RNA negev-like virus: a novel virus persistently infecting Aedes albopictus mosquitoes and cells. J Gen Virol. 2020;101:216–25.
Li C-X, Shi M, Tian J-H, Lin X-D, Kang Y-J, Chen L-J, Qin X-C, Xu J, Holmes EC, Zhang Y-Z. Unprecedented genomic diversity of RNA viruses in arthropods reveals the ancestry of negative-sense RNA viruses. Elife. 2015;4:e05378.
Käfer S, Paraskevopoulou S, Zirkel F, Wieseke N, Donath A, Petersen M, Jones TC, Liu S, Zhou X, Middendorf M, Junglen S, Misof B, Drosten C. Re-assessing the diversity of negative strand RNA viruses in insects. PLoS Pathog. 2019;15: e1008224.
We thank all members of the Bonizzoni’s lab. for fruitful discussions.
This research was funded by a Human Frontier Science Program Research grant (RGP0007/2017) to M.B.; 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 Eccellenza Program (2018–2022) to the Dept. of Biology and Biotechnology “L. Spallanzani”, University of Pavia.
Ethics approval and consent to participate
Adult Ae. albopictus individuals were collected in public areas and immediately homogenized in DNA/RNA-shield (Zymo Research, Irvine CA, USA). Mosquitoes collected in La Reunion Island were transferred to the University of Pavia upon a Material Transfer Agreement between the University of Pavia and CIRAD. Mosquitoes collected in China were processed in loco.
Consent for publication
All authors approved this manuscript and consent to its publication.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
The original version of this article was revised: the authors corrected the author name of Rebeca Carballar-Lejarazu.
. Table S1.
. Table S2
. Presence or absence of AalbF2 annotated nrEVEs in samples from China and La Reunion. (A) Average distribution of reads coverage for 359 BUSCO genes calculated from the WGS mosquito data from La Reunion island and the Guangzhou prefecture in China described in the study. (B) WGS reads coverage for flaviviridaei- and rhabdoviridae-derived annotated nrEVEs in samples from Guangzhou (GZU) and the Eastern and Western sides of La Reunion Island (RE and RW, respectively).
. Table S3
. Dataset Fasta 1.
. Dataset Fasta 2.
About this article
Cite this article
Palatini, U., Alfano, N., Carballar-Lejarazu, R. et al. Virome and nrEVEome diversity of Aedes albopictus mosquitoes from La Reunion Island and China. Virol J 19, 190 (2022). https://doi.org/10.1186/s12985-022-01918-8