Skip to main page content
U.S. flag

An official website of the United States government

Dot gov

The .gov means it’s official.
Federal government websites often end in .gov or .mil. Before sharing sensitive information, make sure you’re on a federal government site.

Https

The site is secure.
The https:// ensures that you are connecting to the official website and that any information you provide is encrypted and transmitted securely.

Access keys NCBI Homepage MyNCBI Homepage Main Content Main Navigation
. 2023 Jul 3;9(2):vead041.
doi: 10.1093/ve/vead041. eCollection 2023.

The virome of the invasive Asian bush mosquito Aedes japonicus in Europe

Sandra R Abbo  1 , João P P de Almeida  2 , Roenick P Olmo  3 , Carlijn Balvers  1   4 , Jet S Griep  1   4 , Charlotte Linthout  4 , Constantianus J M Koenraadt  4 , Bruno M Silva  2 , Jelke J Fros  1 , Eric R G R Aguiar  2   5 , Eric Marois  3 , Gorben P Pijlman  1 , João T Marques  2   3
Affiliations

The virome of the invasive Asian bush mosquito Aedes japonicus in Europe

Sandra R Abbo et al. Virus Evol. .

Abstract

The Asian bush mosquito Aedes japonicus is rapidly invading North America and Europe. Due to its potential to transmit multiple pathogenic arthropod-borne (arbo)viruses including Zika virus, West Nile virus, and chikungunya virus, it is important to understand the biology of this vector mosquito in more detail. In addition to arboviruses, mosquitoes can also carry insect-specific viruses that are receiving increasing attention due to their potential effects on host physiology and arbovirus transmission. In this study, we characterized the collection of viruses, referred to as the virome, circulating in Ae. japonicus populations in the Netherlands and France. Applying a small RNA-based metagenomic approach to Ae. japonicus, we uncovered a distinct group of viruses present in samples from both the Netherlands and France. These included one known virus, Ae. japonicus narnavirus 1 (AejapNV1), and three new virus species that we named Ae. japonicus totivirus 1 (AejapTV1), Ae. japonicus anphevirus 1 (AejapAV1) and Ae. japonicus bunyavirus 1 (AejapBV1). We also discovered sequences that were presumably derived from two additional novel viruses: Ae. japonicus bunyavirus 2 (AejapBV2) and Ae. japonicus rhabdovirus 1 (AejapRV1). All six viruses induced strong RNA interference responses, including the production of twenty-one nucleotide-sized small interfering RNAs, a signature of active replication in the host. Notably, AejapBV1 and AejapBV2 belong to different viral families; however, no RNA-dependent RNA polymerase sequence has been found for AejapBV2. Intriguingly, our small RNA-based approach identified an ∼1-kb long ambigrammatic RNA that is associated with AejapNV1 as a secondary segment but showed no similarity to any sequence in public databases. We confirmed the presence of AejapNV1 primary and secondary segments, AejapTV1, AejapAV1, and AejapBV1 by reverse transcriptase polymerase chain reaction (PCR) in wild-caught Ae. japonicus mosquitoes. AejapNV1 and AejapTV1 were found at high prevalence (87-100 per cent) in adult females, adult males, and larvae. Using a small RNA-based, sequence-independent metagenomic strategy, we uncovered a conserved and prevalent virome among Ae. japonicus mosquito populations. The high prevalence of AejapNV1 and AejapTV1 across all tested mosquito life stages suggests that these viruses are intimately associated with Ae. japonicus.

Keywords: Aedes japonicus; RNA interference; anphevirus; bunyavirus; metagenomics; mosquito; narnavirus; rhabdovirus; totivirus; virome.

PubMed Disclaimer

Conflict of interest statement

None declared.

Figures

Figure 1.
Figure 1.
Analysis of the Ae. japonicus virome using a small RNA-based metagenomic approach. (A) Map of Europe indicating mosquito collection sites: Strasbourg, France and Lelystad, the Netherlands . Pools with the number of sampled mosquitoes from Strasbourg are indicated inside the circles FR_01, two mosquitoes and FR_02, four mosquitoes; and from Lelystad inside the circles NL_01, four mosquitoes, and NL_02, six mosquitoes. Captured mosquitoes were morphologically identified by species. Samples were used to prepare small RNA libraries for high-throughput sequencing. Sequencing results were analysed using our metagenomic pipeline. Assembled contigs were classified into non-viral, viral, and unknown sequences based on sequence similarity against reference databases. (B) Individual results from our sequence similarity analysis for each of the four small RNA libraries in this study. The total number of contigs larger than or equal to 200 nt (n) and the proportion of non-viral, viral and unknown contigs are shown.
Figure 2.
Figure 2.
Co-occurrence of viral and unknown contigs. Hierarchical clustering of viral and unknown contigs assembled from small RNAs derived from Ae. japonicus. Clustering was based on the Euclidean distance of RPKM values of small RNA counts with size from 20 to 22 nt in each library applying average method. Contig clusters were defined using the dendrogram. Contigs inferred to be from the same virus were coloured equally. Heatmap on the left represents the small RNA abundance for each curated contig in Ae. japonicus libraries, but this was not considered for the hierarchical clustering. We plotted the Log2 of the RPKM values of small RNA counts with sizes from 20 to 22 nt (maximum value: 10; minimum value: 0). Heatmap on the right represents Z-score values for small RNAs from each size from 15 to 35 nt divided by strand polarity in the library in which the contig was originally assembled (maximum value: 7; minimum value: −1).
Figure 3.
Figure 3.
Phylogeny of viruses identified in Ae. japonicus mosquitoes. Phylogenetic trees were generated using the multiple sequence alignments of RdRp amino acid sequences. The trees were inferred by using the maximum likelihood method. The tree with the highest log likelihood is shown for each virus. Number of conserved sites and the substitution models used for each tree: (A) AejapNV1, 1446 sites, LG + G + F; (B) AejapTV1, 1269 sites, LG + G; (C) AejapBV1, 616 sites, LG + G + F; and (D) AejapAV1, 1188 sites, LG + G + I + F. Node bootstraps were calculated with 1,000 replicates and are shown close to each clade, and values < 60 were omitted. Trees were midpoint-rooted, and RdRp sequences from distinct viral families were included in the alignments as outgroups. The trees are drawn to scale, and branch lengths represent expected numbers of substitutions per amino acid site. Accession numbers for the nucleotide sequences from which the corresponding protein sequences were derived or the direct protein sequences are shown with the virus names. Viruses identified in this study are in bold.
Figure 4.
Figure 4.
Phylogeny of glycoprotein and nucleocapsid of AejapBV1 and AejapBV2. Phylogenetic trees were generated using the glycoprotein and nucleocapsid amino acid sequences. The trees were inferred by using the maximum likelihood method. The trees with the highest log likelihood are shown. Number of conserved sites and the substitution models used for each tree: glycoprotein, 796 sites, WAG + G + F, and nucleocapsid, 573 sites, LG + G. Node bootstraps were calculated with 1,000 replicates and are shown close to each clade, and values < 60 were omitted. Trees were midpoint-rooted, and sequences from distinct viral families were included in the alignments. The trees are drawn to scale, and branch lengths represent expected numbers of substitutions per amino acid site. Accession numbers for the nucleotide sequences from which the corresponding protein sequences were derived or the direct protein sequences are shown with the virus names. Viruses identified in this study are in bold.
Figure 5.
Figure 5.
CpG dinucleotide usage and GC content of viral and unknown sequences discovered in Ae. japonicus. Data points indicate the ratio of GC content (X-axis) and CpG dinucleotide frequency (Y-axis) of individual viral or unknown sequences. The GC content of 0.5 (vertical dotted line) is the expected frequency if the genome would consist of 50 per cent GC and 50 per cent AT. The observed/expected (O/E) CpG ratio of 1.0 (horizontal dotted line) is the expected frequency of CpG occurrence when all mononucleotides in a given RNA sequence would be randomly distributed. Besides the viral and unknown sequences obtained from Ae. japonicus, the genome sequences of previously discovered narnaviruses in other mosquito species were also included in the analysis. GenBank accession numbers of these mosquito-associated narnaviruses are MW226855.1, MW226856.1, NC_035120.1, KF298284.2, MF176344.1, KF298275.2, MW520409.1, MK285331.1, MK285333.1, and MK285336.1.
Figure 6.
Figure 6.
Small RNA profiles and genome organization of AejapNV1. (A) Left: size distribution and 5ʹ base preference of small RNAs derived from AejapNV1 S1 and S2. Middle and right: coverage of 21 and 24-29 nt sized small RNAs across S1 and S2. Viral reads mapping to the positive strand are shown in blue, whereas viral reads mapping to the negative strand are shown in brown. (B) Genome organisation of AejapNV1. The ambigrammatic coding strategy of S1 and S2 is shown. RdRp ORF and fORF on the positive strand are shown in blue, whereas reverse ORFs (rORFs) on the negative strand are shown in brown. Untranslated regions are indicated by black lines. Predicted SL structures at the 3ʹ terminus of the positive-sense RNA strand are also shown for both segments. The locations of start and stop codons are indicated by arrows and coloured blue and red, respectively. (C) Multiple sequence alignment of the 5ʹ and 3ʹ termini of the positive strand for indicated narnaviruses. Asterisks (*) indicate conserved, complementary runs of G or C nucleotides at the 5ʹ or 3ʹ end, respectively. Dots represent the remainder of the viral genome. Start codons of RdRp ORF are in green, stop codons of RdRp ORF/fORF are in red, and start codons of rORF are in blue. The start codons of the fORFs of AejapNV1 S2 and CxNV1 S2, as well as the start codon of the RdRp ORF of Scer20SNV, are located more than 10 nt downstream of the 5ʹ end and therefore not shown in the alignment. The nucleotides involved in 3ʹ SL formation are indicated.
Figure 7.
Figure 7.
Small RNA profiles and genome organisation of AejapTV1. (A) Left: size distribution and 5ʹ base preference of small RNAs derived from AejapTV1. Middle and right: coverage of 21 and 24-29 nt sized small RNAs across the genome of AejapTV1. Viral reads mapping to the positive strand are shown in blue, whereas viral reads mapping to the negative strand are shown in brown. (B) Genome organization of AejapTV1. Untranslated regions are indicated by black lines. Capsid and RdRp ORFs on the positive strand are shown in blue. These ORFs are encoded in different frames, and a putative −1 ribosomal frameshift area was observed in between the two ORFs. This area consisted of a slippery heptamer, a spacer region, and a predicted three-stemmed pseudoknot with a free energy of −33.97 kcal/mol.
Figure 8.
Figure 8.
Small RNA profiles of AejapBV1 and AejapBV2. Size distribution and 5ʹ base preference of small RNAs derived from (A) AejapBV1 segments L, M and S and (B) AejapBV2 segments M and S are shown on the left, whereas coverage of 21 and 24-29 nt sized small RNAs across the same segments of the respective viruses is shown in the middle and on the right. For the coverage profiles, viral reads mapping to the positive strand are indicated in blue, whereas viral reads mapping to the negative strand are indicated in brown. Inset figures in the middle figures show a zoomed in 21 nt small RNA coverage. Inset figures in the right figures show features of 24-29 nt reads aligned to its respective viral references. Left inset figures show a histogram with normalized frequency of overlap sizes of reads aligned to positive and negative strands. Right inset figures show sequence logo representations of reads aligned to positive and negative strand separately.
Figure 9.
Figure 9.
Small RNA profiles of AejapAV1 and AejapRV1. Size distribution and 5ʹ base preference of small RNAs derived from (A) AejapAV1 and (B) AejapRV1 are shown on the left. In the middle and on the right, the coverage of 21 and 24-29 nt sized small RNAs across the assembled contigs of the same respective viruses is shown. For the coverage profiles, viral reads mapping to the positive strand are indicated in blue, whereas viral reads mapping to the negative strand are indicated in brown. The inset figure in the middle figure of (A) shows a zoomed in 21 nt small RNA coverage. Inset figures in the right figures show features of 24-29 nt reads aligned to its respective viral references. Left inset figures show a histogram with normalized frequency of overlap sizes of reads aligned to positive and negative strands. Right inset figures show sequence logo representations of reads aligned to positive and negative strand separately.
Figure 10.
Figure 10.
Detection of AejapNV1 S1 and S2, AejapTV1, AejapAV1, and AejapBV1 in field-collected Ae. japonicus from the Netherlands by RT-PCR. (A) Individual adult female of Ae. japonicus (Aejap individual females, 1 and 2) or Cx. pipiens (Cpp individual females, 1 and 2), a pool of four Ae. japonicus adult females (indicated by ‘Aejap pool’), and a water sample (no RNA; negative control, indicated by ‘-’) were tested for the presence of the RdRp segment of the four viruses we identified. All samples were also tested for mosquito RPS7 to check the quality of the RNA. Expected amplicon sizes were 1060 bp (AejapNV1 S1), 774 bp (AejapNV1 S2), 998 bp (AejapTV1), 1043 bp (AejapAV1), 1095 bp (AejapBV1), and 175 bp (RPS7). Aejap individual Female 1 tested positive for AejapNV1 S1 and S2 and also for AejapTV1. Aejap individual Female 2 tested negative for AejapAV1, but positive for AejapBV1. The pool of Ae. japonicus females was positive for all tested viruses. Culex pipiens females tested negative for all viruses. The lanes indicated with ‘M’ contain the DNA marker. (B) Prevalence of AejapNV1 S1 and S2 in Ae. japonicus adult females, adult males, larvae, and pools of twenty-five eggs. The number of samples tested is indicated by ‘n’. Individual samples were tested twice for each virus by RT-PCR. Samples with different results in the first and second tests were considered ‘inconclusive’. ‘Not determined’ refers to samples teste for S1 that we did not have enough RNA for additional testing for the presence of S2. (C) Prevalence of AejapTV1 in Ae. japonicus adult females, adult males, and larvae. Individual samples were screened by RT-PCR. The number of samples tested is indicated by ‘n’.
See this image and copyright information in PMC

Similar articles

See all similar articles

Cited by

See all "Cited by" articles

References

    1. Abbo S. R. et al. (2020) ‘The Invasive Asian Bush Mosquito Aedes japonicus Found in the Netherlands Can Experimentally Transmit Zika Virus and Usutu Virus’, PLoS Neglected Tropical Diseases, 14: e0008217. - PMC - PubMed
    1. Agboli E. et al. (2019) ‘Mosquito-Specific Viruses—Transmission and Interaction’, Viruses, 11: 873. - PMC - PubMed
    1. Aguiar E. R. et al. (2015) ‘Sequence-independent Characterization of Viruses Based on the Pattern of Viral Small RNAs Produced by the Host’, Nucleic Acids Research, 43: 6191–206. - PMC - PubMed
    1. Aguiar E. R., Olmo R. P., and Marques J. T. (2016) ‘Virus-derived Small RNAs: Molecular Footprints of Host–Pathogen Interactions’, WIREs RNA, 7: 824–37. - PMC - PubMed
    1. Auewarakul P. (2005) ‘Composition Bias and Genome Polarity of RNA Viruses’, Virus Research, 109: 33–7. - PMC - PubMed