Local-scale virome depiction in Medellín, Colombia, supports significant differences between Aedes aegypti and Aedes albopictus.

Aedes spp. comprise the primary group of mosquitoes that transmit arboviruses such as dengue, Zika, and chikungunya viruses to humans, and thus these insects pose a significant burden on public health worldwide. Advancements in next-generation sequencing and metagenomics have expanded our knowledge on the richness of RNA viruses harbored by arthropods such as Ae. aegypti and Ae. albopictus. Increasing evidence suggests that vector competence can be modified by the microbiome (comprising both bacteriome and virome) of mosquitoes present in endemic zones. Using an RNA-seq-based metataxonomic approach, this study determined the virome structure, Wolbachia presence and mitochondrial diversity of field-caught Ae. aegypti and Ae. albopictus mosquitoes in Medellín, Colombia, a municipality with a high incidence of mosquito-transmitted arboviruses. The two species are sympatric, but their core viromes differed considerably in richness, diversity, and abundance; although the community of viral species identified was large and complex, the viromes were dominated by few virus species. BLAST searches of assembled contigs suggested that at least 17 virus species (16 of which are insect-specific viruses [ISVs]) infect the Ae. aegypti population. Dengue virus 3 was detected in one sample and it was the only pathogenic virus detected. In Ae. albopictus, up to 11 ISVs and one plant virus were detected. Therefore, the virome composition appears to be species-specific. The bacterial endosymbiont Wolbachia was identified in all Ae. albopictus samples and in some Ae. aegypti samples collected after 2017. The presence of Wolbachia sp. in Ae. aegypti was not related to significant changes in the richness, diversity, or abundance of this mosquito's virome, although it was related to an increase in the abundance of Aedes aegypti To virus 2 (Metaviridae). The mitochondrial diversity of these mosquitoes suggested that the Ae. aegypti population underwent a change that started in the second half of 2017, which coincides with the release of Wolbachia-infected mosquitoes in Medellín, indicating that the population of wMel-infected mosquitoes released has introduced new alleles into the wild Ae. aegypti population of Medellín. However, additional studies are required on the dispersal speed and intergenerational stability of wMel in Medellín and nearby areas as well as on the introgression of genetic variants in the native mosquito population.

Introduction

Members of the mosquito genus Aedes (Diptera; Culicidae) are the primary vectors of arboviruses such as dengue virus (DENV), Zika virus (ZIKV), and chikungunya virus (CHIKV) to humans. Two species widely distributed globally are Ae. (Stegomyia) aegypti (Linnaeus, 1762)—the primary vector for various arboviruses that infect millions of people in tropical and subtropical countries every year [1]—and Ae. (Stegomyia) albopictus (Skuse, 1885)—a native to Asia but has expanded its distribution in the last 40 years, invading not only countries in the tropical and subtropical zones but also those in the temperate zones of North America and Europe, where it has been a primary vector of some arboviruses [2, 3].

Arboviruses affect several countries in the Americas, for example, Colombia; it is highly affected by DENV, ZIKV, and CHIKV [4-6]. The high incidence of arboviruses in Colombia is because of the wide distribution of Ae. aegypti throughout the country, where environmental and social conditions promote viral transmission [7, 8]. In 1998, Ae. albopictus was detected in the southern part of the country, and since then, this species has expanded its distribution to numerous cities such as Medellín, where both Ae. aegypti and Ae. albopictus are infected with DENV and ZIKV [9-11]. Since 2000, the number of dengue cases in Medellín has increased continuously, with periodic outbreaks, and nearly 20,000 dengue cases have been reported in each of the last major outbreaks occurring in 2010 and 2016 [12, 13]. Owing to the significant impact these mosquitoes have had on the health of this city’s residents, multiple vector surveillance and control strategies have been implemented. However, the success of these strategies has been questionable as dengue continues to occur in Medellín. The situation in Colombia reflects those in many other arbovirus-endemic countries, which require new approaches to control these vector-borne diseases.

Currently, there is growing scientific interest regarding the mechanism through which the mosquito microbiome, including insect-specific viruses (ISVs), are involved in arbovirus transmission. Among the microorganisms that can modify the vector competence of mosquitoes, the intracellular endosymbiotic bacterium Wolbachia has been demonstrated in vitro to reduce the replication of multiple arboviruses such as DENV, ZIKV, and CHIKV in Ae. aegypti [14-16]; moreover, it can also alter the native host microbiome in adult mosquitoes [17].

In addition to bacteria, the microbiome of mosquitoes also includes viruses. Recent analyses of virus sequences in metagenomics data have changed our understanding of viral diversity, abundance, evolution, and roles in host biology [18, 19]. In mosquitoes, pathogenic viruses represent only a fraction of the total set of viruses (virome) [20-25]. Furthermore, the identities of ISVs vary between populations and species; some ISVs are acquired from the environment, whereas others may circulate through vertical transmission, forming the virome core of mosquitoes populations [22, 26]. However, metagenomic studies for virome characterization may be conflated by the identification of endogenous viral elements (EVEs) that are not real viruses [27], entailing a challenge in these studies for the correct identification of viruses. Mosquitoes, like many organisms, possess EVEs that are remnants of viral integrations into their genomes [28]. And particularly, a large number of non-retroviral EVEs was detected in Aedes mosquito genomes [29, 30] (ref), increasing the interest in mosquito EVEs due to the hypothesis that they may serve as the source of immunological memory against exogenous viruses in insects [31].

The effect of many ISVs on the biology of mosquitoes remains unknown, some viruses may either suppress or enhance the replication of medically important arboviruses such as DENV, ZIKV, CHIKV, and West Nile viruses, suggesting that they play a role in modulating vector competence [32-34]. The insect-specific flavivirus Nhumirim virus is an example of an arbovirus suppressor as it can restrict the growth of ZIKV in mosquito cells in vitro and Ae. aegypti mosquitoes, resulting in significantly lower ZIKV infection rates in both orally infected and intrathoracically inoculated mosquitoes [35]. Although ISVs do not replicate in vertebrates, some are phylogenetically related to known pathogenic arboviruses from the families Flaviviridae, Bunyaviridae, Rhabdoviridae, Reoviridae, and Togaviridae. This finding has led to increasing efforts to discover and characterize more ISVs and explore ISVs as models for studying virus restriction or as potential biocontrol agents [33, 34]. A recent virome description of the global populations of Ae. aegypti and Ae. albopictus reveals significant differences in the composition and diversity of ISVs found in these mosquito species. Furthermore, the abundance of some ISVs such as Phasi Charoen-like virus (PCLV) and Humaita-Tubiacanga virus (HTV) may affect the arboviruses infection capacity as well as their transmission dynamics [36].

In Medellín, Colombia, the primary vector of arboviruses is Ae. aegypti, but Ae. albopictus may potentially be a secondary vector as it can be naturally infected with DENV2 [10, 37] and ZIKV [11, 38]. Both Medellín and its adjacent city, Bello, sustain arboviruses at consistent rates. Thus, these cities were selected for pioneering efforts to test an alternative control strategy based on the release of Wolbachia-transfected Ae. aegypti mosquitoes, which started in 2017 [39]. These releases have not yet provided any conclusive epidemiological evidence of arbovirus control in Medellín; however, in Yogyakarta, Indonesia, the introgression of wMel into Ae. aegypti population was related with a decrease in the incidence of symptomatic dengue [40]. Furthermore, other characteristics describing the native mosquitoes’ biology such as virome composition and mosquito population diversity before and after Wolbachia-transfected Ae. aegypti release remain unexplored. Our study used a metataxonomic approach employing RNA-seq to describe the virome composition, temporal stability, and wide mitochondrial diversity of wild sympatric Ae. aegypti and Ae. albopictus captured between 2015 and 2019 in Medellín, Colombia. To our knowledge, the present study is the first to provide the estimates of the diversity and abundance of viromes in wild-caught Ae. aegypti and Ae. albopictus in Colombia. It is also the first study to report on virome diversity in a city in the Americas where Wolbachia-transfected Ae. aegypti release has been explored. This study thus provides evidence showing the local differences in the richness, diversity, and abundance of ISVs between these two mosquito species and discusses the possible impacts of the Wolbachia-transfected mosquitoes on the core composition of the virome and diversity of the native mosquito population on a local scale.

Materials and methods

Mosquito sampling

The municipality of Medellín, Colombia, is located at 75° 34’05” W and 6°13’55” N and covers an area of approximately 376.2 km2. Medellin is the second most populated city in Colombia, with a population of 2.5 million. Its location in the Andes’ central mountains at an altitude between 1500 and 1800 mts gives Medellin a humid subtropical climate with important rainfalls. The precipitation is characterized by a bi-modal distribution, while the temperatures range between 17 °C and 28 °C.

Indoor resting adult mosquitoes were captured randomly from households using mouth aspirators and entomological nets between 2016 and 2019. The captured mosquitoes were transferred alive to the Medical Entomology Laboratory in the University of Antioquia (Medellín), where they were killed and identified using morphological keys [41]. The blood feeding status of female mosquitoes was determined through an external examination of the abdomen, females that had presence of blood were excluded from the study. Subsequently, the mosquitoes were stored at −80°C and sorted according to the sampling periods and species. In total, 430 mosquitoes were thus divided into 14 pools.

RNA isolation and sequencing

Total RNA was extracted from the mosquitoes using the RNeasy Mini Kit (Qiagen) according to the manufacturer’s instructions. The quality of the RNA extracts was evaluated using an Agilent 2100 Bioanalyzer (Agilent Technologies) operated by Macrogen (Seoul, Korea; www.macrogen.org). After confirming the RNA integrity number (> 7), sequencing libraries were constructed using a TruSeq total RNA library preparation kit (Illumina) after first removing the host rRNA with a Ribo-Zero-Gold (Human–Mouse–Rat) kit (Illumina). Each library was sequenced as 100-bp paired-ends on the Novaseq 6000 S4 platform (Illumina).

Bioinformatic analysis

The paired-end reads were processed using the fastp software v.0.20.0 to remove adaptor sequences and perform quality-based filtering [42]. The resulting high-quality reads (Phred quality score > 20) were mapped to their respective reference genomes (Ae. aegypti assembly AaegL5.2 [43] and Ae. albopictus assembly AaloF1.2 [44]) using the BWA-MEM option in VectorBase (https://vectorbase.org/vectorbase/app/) [45]. The unmapped reads were subsequently analyzed as described below.

Viral sequences were evaluated by first creating a custom database of ribosomal RNA sequences using SILVA v.132 LSU, SSU, and 5S rRNA (RF00001), 5.8S rRNA (RF00002), tRNA (RF00005), and ribonuclease P (RF00010, RF00011, and RF00373). Our custom database was then used via the software SortMeRNA v.2.1 [46] to identify, using an e-value cutoff of 10−5, ribosomal sequences in the unmapped reads [47]. Reads with 60% of the read length identical to ribosomal RNA sequences by >60% were excluded from further analysis. The remaining sequences were assembled using SPAdes v.3.12.0 [48] (S1 Table), and contigs > 450 nt were classified using DIAMOND v.2.0.8 [49] and the nr viral database (updated March 2021). The analysis employed the sensitive mode for taxonomic annotation with an e-value cutoff of 10− and a bit score cutoff of 100 [47, 50]. The output file of DIAMOND was parsed using KronaTools v.2.7.1 [51], which found the least common ancestor of the best 25 DIAMOND hits. Only contigs with identity scores of >60% were stored. Contigs identified as viral operational taxonomic units (vOTUs) were confirmed by performing an online BLASTn and BLASTx (https://blast.ncbi.nlm.nih.gov/Blast.cgi) searches, thus eliminating possible false positives. All OTUs were assigned to viral species based on their BLASTx identity. Although we used an RNA-seq strategy and perform a manual review of the viral contigs, allowing to identify and exclude some EVEs, we cannot completely exclude the possibility that some contigs are derived from EVEs. This must be assessed in future studies that evaluate both the identification of EVEs in these mosquito genome populations and their expression.

Virome diversity analysis

The number of viral reads per library was calculated by mapping the libraries against the viral contigs identified, then they were normalized by simply dividing viral reads number by the mosquito mapped reads number, and then multiplying the result by the size of the smaller sample. The alpha diversity index (i.e., diversity within each sample) was evaluated based on the observed virome richness, Shannon index, Simpson index, and evenness of each library at the virus species level using the Rhea script set [52]. Finally, the indices of the host mosquito species were compared using Kruskal–Wallis rank sum and Mann–Whitney tests employing the R software [53].

Beta diversity (i.e., viral diversity between samples) was estimated by first constructing a Bray–Curtis dissimilarity matrix [54], which considers both the shared taxonomic composition and virus abundance in the viromes. The results were plotted using nonmetric multidimensional scaling (NMDS) ordination, and the significance of the resulting clusters was tested using permutational multivariate analysis of variance (PERMANOVA; Adonis tests) employing the R software.

Wide mitochondrial diversity between Ae. aegypti and Ae. albopictus

Paired-end reads were mapped using BWA-MEM v.0.7.17 [45, 55] against the mitochondrial genomes of Ae. aegypti and Ae. albopictus (NC_035159.1 and NC_006817.1, respectively). Only reads with a mapQ score of >20 were retained using SAMtools v.1.9 [56, 57]. Optical polymerase chain reaction (PCR) duplicates were removed, and read grouping was performed in Picard v.2.9.0 (http://broadinstitute.github.io/picard/). Variant calls in each pool were estimated using the “—pooled-discrete” option in the haplotype-based variant detection method implemented in freebayes v.1.2.0 [58]. Multiallelic positions at total depths of >10 x were included and indels were excluded from variant calling. To analyze sequences, we retrieved a single sequence from each sample using FastaAlternateReferenceMaker in GATK v.4.1.2.0 [59]. Then, the overall nucleotide polymorphism per site (Theta-W) was estimated using DnaSP v.6 [60], and a Kimura -2p based distance tree was constructed using MEGA X [61].

DENV PCR assay

Because we detected DENV-3 in one of the RNA-seq samples, we performed a semi-nested reverse transcriptase (RT)-PCR to confirm this finding. The Luna Universal One-Step RT-qPCR kit was used according to manufacturer’s instructions with the primers D3_Fwd1 (5’-GACCCAGAAGGCGGTTATTT-3′) and D3-Rev1 (5’-GCCTCGAACATCTTCCCAATA-3′). These primers amplify a 1,260-bp region of the envelope gene [62]. The second PCR used the Thermo Scientific Taq DNA polymerase according to manufacturer instructions and employed the primers D3-Rev1 (5’-GCCTCGAACATCTTCCCAATA-3′) and D3-MA (5’-ACAAGCCCACGTTGGATATAG-3′). These primers amplify a 1,057-bp fragment of the envelope gene. Amplicons were purified and sequenced (Macrogen). The forward and reverse sequences were joined using BioEdit v7.2 [63], generating the consensus sequences.

Phylogenetic analysis

The obtained consensus sequences were first aligned with reference sequences before phylogenetic analysis. The alignments were performed with MAFFT v.7.475 using the L-INS-i algorithm [64] and 1,000 cycles of iterative refinement. The phylogenetic tree was reconstructed from 1,000 ultrafast bootstrap maximum likelihood-based tree replicates using IQ-TREE v1.6.12 [65]. The best-fitting model was selected using ModelFinder [66]. The phylogenetic tree was drawn using FigTree v.1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/).

Results

In this study, a total of 10 and 4 pooled samples of Ae. aegypti and Ae. albopictus, respectively, were analyzed using mosquitoes captured between 2015 and 2019. Each pooled sample contained 15–36 mosquitoes (Table 1), and these were subjected to RNA extraction followed by RNA sequencing. After processing the raw data, a total of 1,190,309,014 (range 93,199,032–131,645,868) 100-bp paired-end reads were generated from the 10 ribosomal RNA-depleted sequence libraries from Ae. aegypti; a total of 493,223,542 (range 105,740,360–131,645,868) 100-bp paired-end reads were generated from the four Ae. albopictus libraries. Approximately 65.3% of the reads corresponded to the mosquito sequences; therefore, a total of 405,798,910 and 191,058,189 reads were used to characterize the viromes of Ae. aegypti and Ae. albopictus, respectively, collected from Medellín, Colombia.

Table 1

Summary of the sequencing results showing data on pooled samples and the numbers of mosquitoes, mosquito species, read sequences, and mapped mosquito reads.

Sample name	Mosquito species	No. of mosquitoes	Sampling time	No. of clean reads	Mosquito reads mapped
Ae. aeg 2015A	Ae. aegypti	15	January–May, 2015	93,199,032	88,936,067
Ae. aeg 2015B	Ae. aegypti	25	September–November, 2015	113,389,586	104,933,849
Ae. aeg 2016A	Ae. aegypti	35	February–June, 2016	129,886,418	88,083,776
Ae. aeg 2016B	Ae. aegypti	33	July–November, 2016	122,388,802	74,256,344
Ae. aeg 2017A	Ae. aegypti	31	March–June, 2017	105,740,360	56,373,655
Ae. aeg 2017B	Ae. aegypti	34	July–December, 2017	128,342,842	75,574,534
Ae. aeg 2018A	Ae. aegypti	35	March–June, 2018	119,858,346	72,780,110
Ae. aeg 2018B	Ae. aegypti	33	July–November, 2018	131,645,868	85,013,726
Ae. aeg 2019A	Ae. aegypti	31	January–June, 2019	123,578,244	74,904,614
Ae. aeg 2019B	Ae. aegypti	28	July–November, 2019	122,279,516	76,372,131
Ae. alb 2017	Ae. albopictus	25	August–December, 2017	124,650,502	81,419,768
Ae. alb 2018	Ae. albopictus	36	March–August, 2018	114,442,616	63,806,581
Ae. alb 2019A	Ae. albopictus	34	January–June, 2019	130,821,410	71,501,186
Ae. alb 2019B	Ae. albopictus	35	July–November, 2019	123,309,014	85,437,818

Virome characterization

After annotating and analyzing the contigs with DIAMOND BLASTx annotation and KronaTools, respectively, the contigs were verified; 125 of them were identified as viral OTUs (86 from Ae. aegypti and 39 from Ae. albopictus) comprising 21 virus species (Fig 1). These virus species were related to ISVs, arboviruses, and plant viruses belonging to nine different viral families: Flaviviridae, Totiviridae, Reoviridae, Phenuiviridae, Iflaviridae, Bromoviridae, Metaviridae, Xinmoviridae, and Orthomyxoviridae; however, six viruses remained unclassified.

Fig 1

Boxplot of contigs length of each virus species identified.

Sequences from each library were mapped to the identified viral contigs to determine the abundance of each identified virus species. The rarefaction curves shown in Fig 2A are based on the number of viral sequences; these show that the curves describing samples are approaching the saturation plateau, suggesting that most viruses were detected in this experiment. In addition, Ae. albopictus curves showed a tendency to plateau faster, suggesting that the virome of this species is less diverse than that of Ae. aegypti. Furthermore, significantly lower numbers of viral sequences were detected in Ae. albopictus (Fig 2B).

Fig 2

A. Rarefaction curves representing the viruses in Ae. aegypti and Ae. albopictus. B. Proportion of viral sequences relative to the number of mosquito host sequences. Mann–Whitney Test: p < 0.05 (*) and p < 0.01 (**).

A total of 75,219,419 viral sequences (6.3% of libraries) were identified in Ae. aegypti, which were assigned to 17 viral species. The virome of this species was dominated by four viruses, which accounted for approximately 90% of the viral sequences. Of the total viral sequences, 52% corresponded to PCLV (Phenuiviridae); PCLV was present in all libraries and was the most abundant virus in almost all samples, except for Ae. aeg 2015B, where Aedes anphevirus (AaAV; Xinmoviridae) was the most abundant virus detected. The second and third most abundant viruses were Kwale mosquito virus (21.5%) and Guadeloupe mosquito virus (GMV; 12.2%), respectively; both viruses are unclassified but are probably mosquito-specific viruses. Finally, cell fusing agent virus (CFAV; Flaviviridae) accounted for 3.7% of the total viral sequences.

Interestingly, DENV3 (Flaviviridae) was detected in one of the samples from 2016. We confirmed the presence of this virus through RT-PCR targeting a fragment of the envelope gene. Subsequent phylogenetic analysis identified this sequence as DENV3 genotype III (S1 Fig). No other pathogenic arboviruses were identified in this study. The Ae. aegypti virome composition was fairly conserved among the samples throughout the 5 years of sampling, suggesting that the virome of Ae. aegypti in Medellín is stable over time.

The virome of Ae. albopictus had less viral sequences and species than that of Ae. aegypti (Fig 2). In Ae. Albopictus, a total of 12 viral species were identified, which were pooled into five families (84,339 reads, 0.02% of the total). Aedes flavivirus (AeFV; Flaviviridae) was the most abundant virus, accounting for 44.4% of the viral sequences in all samples (Fig 3). Australian Anopheles totivirus (Totiviridae) accounted for 16% of viral sequences, and Katimi river virus (Flaviviridae) accounted for 15.3%. PCLV was predominant in the virome of Ae. aegypti, whereas it accounted for only 1% of the total viral sequence in Ae. albopictus. Of all viruses detected in Ae. albopictus, only eight were present in multiple samples (Fig 3). Tobacco streak virus (Bromoviridae)—a plant virus—and Kwale mosquito virus were detected in sample Ae. alb 2017, whereas GMV (unclassified viruses) and HTV were detected in sample Ae. alb 2019A.

Fig 3

Heatmap showing the normalized abundance of reads on a log2 scale.

The viral species names are from the taxonomic annotation based on DIAMOND and KronaTools. The average BLASTx identity of each virus species to a reference sequence is shown in the blue-green boxes on the left.

Virome comparison between the Aedes spp.

The two Aedes spp. differed in the composition and abundance of virus species in their respective viromes. Fig 3 clearly shows the hierarchical clustering of virus species in each mosquito species based on the Euclidean distance matrix. In total, 21 virus species were identified, of which only eight were shared between the mosquito species studied (Australian Anopheles totivirus, CFAV, PCLV, GMV, HTV, Kwale mosquito virus, Kaiowa virus, and Aedes aegypti To virus 2). Similarly, viruses in both mosquito species had significantly different levels of alpha diversity, as evident from the indexes of diversity and richness (Fig 4A). Although the Ae. aegypti virome is richer in virus species, the Ae. albopictus virome is more diverse, based on Shannon and Simpson indices. This higher diversity may be because of the higher evenness of the Ae. albopictus viral community.

Fig 4

Diversity analysis of the Ae. aegypti and Ae. albopictus viromes.

A. Alpha diversity of viruses in Ae. aegypti and Ae. albopictus at the species level. Mann–Whitney Test: p < 0.05 (*), p < 0.01 (**), and p < 0.001 (***). B. Nonmetric multidimensional scaling (NMDS) of viruses associated with the two mosquito species. STRESS = 0.023, PERMANOVA test on mosquito species: p = 0.004 and R2 = 0.59. C. Comparison of normalized and relative abundance levels of Aedes aegypti To virus 2 in the samples of Ae. aegypti with and without Wolbachia sp.

Beta diversity is based on Bray–Curtis dissimilarities that were calculated from the normalized abundance of viral species and then analyzed via unconstrained ordination using the PERMANOVA test on mosquito species (p = 0.004 and R2 = 0.59) and NMDS. Fig 4B shows a clear separation of viral communities according to mosquito species, suggesting that the viromes of these two mosquito species differ significantly in their compositions.

Effects of Wolbachia on the Ae. aegypti and Ae. albopictus viromes

The presence of the intracellular symbiotic bacterium, Wolbachia, was detected in all Ae. albopictus samples, which might have been because this species is naturally infected by the bacterium. In contrast, Ae. aegypti is not a natural host for Wolbachia; however, Ae. aegypti infected with Wolbachia have been released in Medellín since 2017. Wolbachia sequences were detected in five Ae. aegypti samples that were collected from mid- 2017 to 2019. The consensus contigs of the gene encoding the surface protein (wsp) were used for phylogenetic analysis (Fig 5), which shows that Ae. albopictus was infected by Wolbachia spp. from clades A and B. In Ae. aegypti, Wolbachia sp. belonging to clade A, which is related to Drosophila melanogaster isolates, was detected. This result was expected because the released mosquitoes were transfected with the wMel strain, which originated from D. melanogaster.

Fig 5

Phylogenetic reconstruction of the Wolbachia consensus sequences associated with Ae. albopictus (blue) and Ae. aegypti (red).

Whether the Wolbachia infection in Ae. aegypti could influence the richness or diversity of the Ae. aegypti virome was investigated. For this purpose, Ae. aegypti samples with and without Wolbachia sp. were compared. No significant differences in richness or diversity were found between the samples; however, Aedes aegypti To virus 2, anErrantivirus, was significantly more abundant in samples with Wolbachia (Fig 4C).

Mitochondrial diversity of Aedes mosquitoes

On average, 459,117 (± 222,695) reads were mapped in Ae. aegypti, whereas 443,887 (± 195,269) reads were mapped in Ae. albopictus (Table 2). A total of 285 and 150 single-nucleotide polymorphisms (SNPs) were detected in Ae. aegypti and Ae. albopictus, respectively (Table 2).

Table 2

Results of mitochondrial diversity analysis showing the numbers of mitochondrial mapped reads and SNPs.

Sample name	No. of mitochondrial mapped reads	No. of SNPs
Ae. aeg 2015A	859,461	163
Ae. aeg 2015B	207,324	143
Ae. aeg 2016A	352,912	155
Ae. aeg 2016B	328,784	150
Ae. aeg 2017A	109,532	122
Ae. aeg 2017B	435,763	197
Ae. aeg 2018A	565,002	175
Ae. aeg 2018B	682,894	177
Ae. aeg 2019A	495,284	187
Ae. aeg 2019B	554,218	195
Average	459,117	285
Ae. alb 2017	308,580	112
Ae. alb 2018	591,312	128
Ae. alb 2019A	630,538	142
Ae. alb 2019B	245,118	106
Average	443,887	150

In general, the number of SNPs per sample was higher in Ae. aegypti than in Ae. albopictus, except for sample Ae. aeg 2015A. However, variable sites were observed throughout the entire mitochondrial genome in both species. SNP density peaked at around position 7,000 of the mitochondrial genome, a region involved with coding for the ND5 gene in both species (Fig 6). In addition, a second high-density SNP region located between positions 12,400 and 14,600 of the mitochondrial genome was observed in Ae. aegypti samples. This location harbors coding regions of the ND1, tRNA-Leu, tRNA-Val, and 16s, and 12s rRNA genes (Fig 6).

Fig 6

Mitochondrial SNP density per site in Ae. aegypti and Ae. albopictus samples collected between 2015 and 2019 in Medellín, Colombia.

The overall nucleotide polymorphism (Theta-W) was 0.004 for Ae. aegypti, whereas it was 0.002 for Ae. albopictus. Thus, the neighbor-joining tree showed two well-supported clades harboring both species. A secondary clade including samples collected after June 2017 (i.e., sample Ae. aeg 2017A) was also revealed in Ae. aegypti (Fig 7).

Fig 7

Distance-based tree based on mitochondrial sequences in Ae. aegypti and Ae. albopictus collected between 2015 and 2019 in Medellín, Colombia.

Discussion

Recent studies on mosquito microbiomes have begun to elucidate the complex interactions between microorganisms and their hosts as well as the effects the former has on the latter. Most studies have focused on the bacteriome, thus the role and transmission of the virome, which comprises a highly diverse community of viruses, are less understood. Metagenomic studies rely on database analysis, which is a problem in virome studies because sequence information for many viral families and genera is still limited. Therefore, this study used stringent restrictions such as OTUs identity and cover greater than 60% to reduce the possibility of reporting false positives. Although false positive errors are minimized through this approach, it can also result in the underestimation of the true diversity of viruses in these mosquito populations.

This study, to the best of our knowledge, is the first to characterize the viromes of Ae. aegypti and Ae. albopictus populations in Colombia, where mosquito-borne diseases such as dengue, Zika, and chikungunya are major public health concerns [5]. The findings of the present study add significant information to the existing knowledge regarding the mosquito viromes in the Americas and the world. The study results reveal that the compositions of the viromes of sympatric, field-caught Ae. aegypti and Ae. albopictus differ significantly, suggesting that virome composition is species-specific and reflects differences in the host evolutionary history, host immunological response, virus–mosquito interactions, and perhaps vector competence. The two mosquito species also differ largely in the richness of virus species associated with them and the abundances of these viruses. This result was consistent across all comparisons (Fig 3).

A larger proportion of viral sequences were associated with Ae. aegypti; these sequences represented a greater richness of viruses than those associated with Ae. albopictus. These results support the hypothesis that Ae. aegypti has a higher capacity to carry viruses and, therefore, to disperse these viruses more widely [36]. This also suggests that this mosquito species is more susceptible to infection by arboviruses such as DENV, ZIKV, and CHIKV [36]. Globally, the Ae. aegypti virome is highly variable in terms of the richness of viral species [22, 36]. In the present study, the Ae. aegypti virome comprised least 17 virus species predominated by PCLV, which is consistent with the results of previous studies showing that this virus superinfects Ae. aegypti, for example, PCLV has been reported to be a predominant ISV in Ae. aegypti populations in Australia [24], Guadeloupe [67], South China [68], Thailand [69], Brazil [36], and Grenada [70]. Interestingly, PCLV is less abundant in African populations such as those in Senegal, and it is absent in Gabon, where the dominant subspecies is related to Ae. aegypti formosus, which is considered less susceptible to arbovirus infection [71].

GMV and Kwale mosquito virus were also highly abundant in Ae. aegypti samples. GMV has been identified recently in Ae. aegypti from Guadeloupe, where it shows high abundance in the mosquito virome [67, 72]. Kwale mosquito virus was detected in a metagenomic analysis of Ae. aegypti mosquitoes from Kenya; it is phylogenetically closely related to Hubei mosquito virus 2 and GMV related viruses [22]. GMV and Kwale mosquito virus are currently considered unclassified viruses, and their effects on the biology of Aedes remain unknown. Another abundant member of the Ae. aegypti virome is CFAV, a flavivirus that was the first insect-specific flavivirus found in Aedes [73]. Since then, it has been detected in multiple regions [22, 24, 74–76]. Furthermore, CFAV infection can significantly enhance DENV replication in Ae. aegypti cells [77].

In this study, Flaviviridae was the most diverse family, and three ISVs as well as one arbovirus (CFAV, Menghai flavivirus, Xishuangbanna AeFV, and DENV3) were detected in most samples. Menghai flavivirus has been isolated from Ae. albopictus mosquitoes in China, and it is phylogenetically related to Xishuangbanna AeFV, which belongs to the Aedes-associated flaviviruses cluster [78]. Although ISVs are abundant and highly prevalent, arbovirus infections in wild mosquitoes are rare. Nonetheless, DENV3 was detected in sample Ae. aeg. 2016A. Results of the phylogenetic analysis revealed that this virus belongs to genotype III, the same genotype that has been circulating in Colombia and other South American countries [61]. A dengue epidemic occurred in Medellín in 2016 [13]; however, the Ae. aegypti virome from that year did not differ from those from other years.

ISVs from the Totiviridae family (Aedes aegypti toti-like virus and Australian Anopheles totivirus) and one from the Iflaviridae family (iflavirus) were detected in sample Ae. aeg 2016B. The latter is related to slow bee paralysis virus (72.2% BLASTx identity), although it may also be a new noncharacterized iflavirus. Similarly, several viruses from Aedes mosquitoes have been tentatively identified as Iflavirus-like viruses [72, 79, 80]. These results indicate that this virus family needs to be expanded and revised. Furthermore, AaAV (Xinmoviridae) was detected in all samples. This Ae. aegypti-associated virus has been detected in field-caught and colony mosquitoes as well as in cells lines. AaAV can interact with Wolbachia and, in the process, influence DENV replication in cell lines [81].

Virus species richness was lower in Ae. albopictus samples, where 12 virus species were detected. The core virome of Ae. albopictus mainly comprised eight viruses that were found in multiple samples. The most diverse set of viruses belonged to the Flaviviridae family and included three ISVs, namely, AeFV, CFAV, and Kamiti river virus. AeFV is a flavivirus associated with both Ae. albopictus and Ae. flavopictus [82]. A temporal study of AeFV in Ae. albopictus in China suggests that it is a seasonal virus [83]. Because Colombia is a tropical country, AeFV may be perennially present in Ae. albopictus.

Australian Anopheles totivirus was the second most abundant species in Ae. albopictus, and it was present in all samples in addition to Kaiwá virus, Aedes aegypti To virus 2, Lampyris noctiluca errantivirus 1 (Metaviridae), and PCLV. In contrast to the findings obtained from Ae. aegypti, PCLV was rare in Ae. albopictus samples. Notably, both PCLV and CFAV have been reported to inhibit ZIKV, DENV, and La Crosse virus (Peribunyaviridae) in the Ae. albopictus cell line Aa23 [84].

More unique viruses per sample were detected in Ae. albopictus than in Ae. aegypti; furthermore, Tobacco streak virus (Bromoviridae), HTV, Kwale mosquito virus, and GMV were each detected in only one sample, suggesting that these viruses have been acquired from the environment. The higher level of variability in the virome composition of Ae. albopictus may be because of this mosquito’s tendencies to breed in natural habitats and feed on different types of vertebrates [85, 86], thus exposing itself to more viruses in the environment.

Mosquitoes are not known to be vectors of plant viruses; however, they can acquire these viruses from nectar during feeding; thus, some plant viruses have been detected in mosquito viromes [24, 25]. For example, in China, a plant-specific Tymoviridae-like virus that was isolated from Culex spp. was able to replicate in mosquito cells, suggesting that it had the potential to be transmitted by mosquitoes [87].

Despite the lower level of virus richness in Ae. albopictus than that in Ae. aegypti, the Alpha diversity of the Ae. albopictus virome was higher. This is because the abundance levels of viruses in Ae. albopictus were more evenly distributed. Multiple virus sequences identified in this study are consistent with other recently described sequences that currently lack a formal taxonomic classification. Therefore, the mosquito virome needs to be better characterized and considered in relationship with vector competence. However, studies relating the mosquito viroma profile and vector competence must take into consideration that the virome composition of mosquitoes is easily modified when they are used for establishing laboratory colonies from field-caught mosquitoes [22].

An important difference between Ae. albopictus and Ae. aegypti is the high prevalence of the intracellular bacterium Wolbachia in the former [88, 89] and its absence in the latter [90]. The effect of Wolbachia on mosquitoes viromes is still unclear, however, it has been hypothesized that Wolbachia colonization can reduce the abundance and richness of its associated viruses, which may explain the differences reported between the Ae. aegypti and Cx. quinquefasciatus viromes [67], and also the differences observed into this work between the Ae. albopictus and Ae. aegypti viromes [17]. The present study investigated whether differences between the microbiomes of Ae. aegypti and Ae. albopictus influence the compositions of their respective viromes. Samples collected between 2015 and the first half of 2017 that were negative for Wolbachia were compared with those that were positive for this bacterium. The results suggest that the presence of wMel Wolbachia did not influence the diversity or richness of the Ae. aegypti virome. However, the number of Wolbachia-infected mosquitoes in each of the five pools of mosquitoes was not quantified. In the wild populations of D. melanogaster, Wolbachia infection does not influence the diversity of native viruses in the insect [91]. Nevertheless, the results of the present study suggest that the presence of Wolbachia may enhance the level of infection of mosquitoes with Aedes aegypti To virus 2. A similar result was observed in other populations of field-caught Ae. aegypti mosquitoes, where Wolbachia enhanced ISF infection rates and loads, demonstrating that Wolbachia does not act as an antiviral against all flaviviruses [24, 92]. Therefore, we consider that Wolbachia presence in Ae. aegypti can alter the infection levels of low frequency viruses, but does not affect the mosquito core virome. However, this needs to be better addressed in future studies.

Although Wolbachia did not influence the Ae. aegypti virome diversity, a higher number of mitochondrial SNPs clustered apart from those of previous periods were observed in Ae. aegypti samples collected after the second half of 2017 (when the first Wolbachia-transfected mosquitoes were released; Table 2). Therefore, the release of these Wolbachia-transfected mosquitoes appeared to affect the native Ae. aegypti population (Figs 6 and 7). The transfected mosquitoes were the product of a cross between the native population and Wolbachia-bearing mosquitoes from Australia [93]. Therefore, a temporal interpopulation substructure might have formed after transfected mosquitoes were released in Bello and Medellín. In particular, some genetic variants of the Australian strain may have remained after the crosses, and these were introduced into the Medellín mosquito population, increasing its diversity. Evidence for this can be obtained by analyzing biparental genetic markers and determining the levels of introgression between the native and introduced populations. Future Wolbachia-infected mosquitoes release can achieve faster rates of invasion at a lower cost by determining whether these variants become fixed in the native population and by assessing their effects on maintaining Wolbachia infections and its local spread capacities [94, 95].

The analysis of mitochondrial diversity during the 3 years of Ae. albopictus sampling revealed lower numbers of mitochondrial variable sites and SNPs than those observed in Ae. aegypti; moreover, Ae. albopictus diversity was similar between the sampling periods. Ae. albopictus is an invasive species detected in Medellín in 2011 [96], and its relatively low diversity suggests that its population size is smaller than that of Ae. aegypti and reflects a more recent founder event. To the best of our knowledge, however, no population genetic studies of this species have been conducted in either Medellín or Colombia. Therefore, it is not possible to establish how the diversity of this species has changed locally.

Conclusions

In conclusion, the levels of abundance and diversity between the viromes of Ae. aegypti and Ae. albopictus differ strikingly, indicating that virome profiles are species-specific and can be determined based on their evolutionary history. Most of the viruses associated with both species were detected throughout different sampling periods, suggesting that the viromes are temporally stable. Wolbachia spp. from clades A and B were prevalent in Ae. albopictus, whereas one strain of Wolbachia was present in post-release samples of Ae. aegypti. The presence of Wolbachia in Ae. aegypti was not related to any changes in the composition of its virome, except for an increase in the abundance of one ISV, Aedes aegypti To virus 2, which could be related to the presence of Wolbachia. Furthermore, mitochondrial genetic diversity was lower in Ae. albopictus than in Ae. aegypti, which suggests that the population size of the former is smaller and reflects a more recent founder event. The genetic diversity of Ae. aegypti increased after the second half of 2017, likely related to the release of mosquitoes transfected with Wolbachia. This study thus provides a baseline for future virome studies on Ae. aegypti and Ae. albopictus in Colombia and other countries in the Americas.

Phylogenetic reconstruction of dengue virus sequences from Aedes aegypti mosquitoes in Medellín (Colombia) using a 1,057-bp region of the envelope gene.

The contig was assembled from RNA-seq data, and the fragment was amplified and sequenced using the Sanger method. The sequences generated in this study are indicated in red. The scale bar indicates the number of substitutions per site.

(TIF)

Click here for additional data file.

Summary of assembly metrics.

(DOCX)

Click here for additional data file.