Front. Cell. Infect. Microbiol.Frontiers in Cellular and Infection MicrobiologyFront. Cell. Infect. Microbiol.2235-2988Frontiers Media S.A.10.3389/fcimb.2023.1252656Cellular and Infection MicrobiologyOriginal ResearchGenome-wide detection of Wolbachia in natural Aedes aegypti populations using ddRAD-SeqMuharromahAtikah Fitria12ReyesJerica Isabel L.1KagiaNgure1WatanabeKozo1*1Molecular Ecology and Health Laboratory, Center for Marine Environmental Studies (CMES), Ehime University, Matsuyama, Japan2Entomology Laboratory, Department of Tropical Biology, Faculty of Biology, Universitas Gadjah Mada, Yogyakarta, Indonesia
Edited by: Joseph James Gillespie, University of Maryland, United States
Reviewed by: A. Alwin Prem Anand, University of Tübingen, Germany
S. Gowri Sankar, Vector Control Research Centre (ICMR), India
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
Background
Wolbachia, an endosymbiotic bacterium, is globally used to control arboviruses because of its ability to block arboviral replication and manipulate the reproduction of Wolbachia host, Aedes aegypti. Polymerase chain reaction (PCR)-based Wolbachia detection has been recently reported from natural Ae. aegypti populations. However, due to the technical limitations of PCR, such as primer incompatibility, PCR-based assays are not sufficiently reliable or accurate. In this study, we examined double digestion restriction site-associated DNA sequencing (ddRAD-Seq) efficiency and limitations in Wolbachia detection and quantification in field-collected Ae. aegypti natural populations in Metro Manila, the Philippines, compared with PCR-based assays.
Methods
A total of 217 individuals Ae. aegypti were collected from Metropolitan Manila, Philippines. We separated it into 14 populations consisting of 7 female and male populations. We constructed a library for pool ddRAD-Seq per population and also screened for Wolbachia by PCR assays using wsp and 16S rRNA. Wolbachia density per population were measured using RPS17 as the housekeeping gene.
Results
From 146,239,637 sequence reads obtained, 26,299 and 43,778 reads were mapped across the entire Wolbachia genome (with the wAlbA and wAlbB strains, respectively), suggesting that ddRAD-Seq complements PCR assays and supports more reliable Wolbachia detection from a genome-wide perspective. The number of reads mapped to the Wolbachia genome per population positively correlated with the number of Wolbachia-infected individuals per population based on PCR assays and the relative density of Wolbachia in the Ae. aegypti populations based on qPCR, suggesting ddRAD-Seq-based semi-quantification of Wolbachia by ddRAD-Seq. Male Ae. aegypti exhibited more reads mapped to the Wolbachia genome than females, suggesting higher Wolbachia prevalence rates in their case. We detected 150 single nucleotide polymorphism loci across the Wolbachia genome, allowing for more accurate the detection of four strains: wPip, wRi, TRS of Brugia malayi, and wMel.
Conclusions
Taken together, our results demonstrate the feasibility of ddRAD-Seq-based Wolbachia detection from field-collected Ae. aegypti mosquitoes.
WolbachiaAedes aegyptiddRAD-seqPhilippinesgenome-wideJapan Society for the Promotion of Science10.13039/501100001691Sumitomo Electric Industries10.13039/100016690section-in-acceptanceBacteria and HostIntroduction
Dengue, Zika, and Chikungunya represent major public health concerns worldwide (Silva et al., 2020). These arboviral diseases are transmitted by the vector mosquito Aedes aegypti. A novel approach for combating these mosquito-borne diseases using Wolbachia bacteria has been established in various countries (Hoffmann et al., 2011; Nguyen et al., 2015; Schmidt et al., 2017; Nazni et al., 2019; Zheng et al., 2019; Crawford et al., 2020; Beebe et al., 2021; Pinto et al., 2021; Utarini et al., 2021; Ahmad et al., 2021). Exploiting the cytoplasmic incompatibility of Wolbachia on host reproduction could suppress mosquito populations by resulting in an unviable embryo, thereby replacing natural mosquito populations with Wolbachia-infected ones. For instance, dengue incidences in Wolbachia-treated areas in Australia and Indonesia were reduced by 96% (Ryan et al., 2020) and 77% (Utarini et al., 2021), respectively.
Wolbachia inhibits arboviral replication in host mosquitoes. The wMel Wolbachia strain reduces CHIKV (Aliota et al., 2016b), ZIKV (Aliota et al., 2016a), and DENV (Walker et al., 2011) transmission from Ae. aegypti to other hosts, including humans. Another Wolbachia strain, wAlbA, blocks ZIKV (Chouin-Carneiro et al., 2019) and the wAlbB strain might inhibit dengue and Zika virus transmission in Ae. aegypti (Hugo et al., 2022). For an effective Wolbachia-based arbovirus control, important information such as Wolbachia infection prevalence and that related to natural mosquito population strains should be available, as Wolbachia strain co-infection in Ae. aegypti could potentially induce inter-strain competition in the host mosquito. For example, triple-strain infection (wMel, wAlbA, and wAlbB) in Ae. albopictus and Ae. aegypti inhibited cytoplasmic incompatibility expression and showed low maternal transmission fidelity (Ant and Sinkins, 2018; Liang et al., 2020). Prior knowledge of Wolbachia strains in native mosquito populations in the deployment area could help identify the most suitable Wolbachia strain to use.
Natural Wolbachia infection in Ae. aegypti remains controversial. For instance, using PCR assays, Gloria-Soria et al. (2018) did not detect Wolbachia in any Ae. aegypti collected from 27 countries. Other studies also confirmed the lack of Wolbachia detection in Ae. aegypti from Cape Verde islands (da Moura et al., 2023), Singapore (Ding et al., 2020) and California (Torres et al., 2020). In contrast, natural Wolbachia infection could be identified in Ae. aegypti in the USA (Coon et al., 2016; Kulkarni et al., 2019), Malaysia (Teo et al., 2017), Thailand (Thongsripong et al., 2017), India (Balaji et al., 2019), the Philippines (Carvajal et al., 2019), Panama (Bennett et al., 2019), and China (Zhang et al., 2022). A possible explanation for these results is the different susceptibility to Wolbachia between host individuals, potentially influenced by host genotype and environmental conditions (Mouton et al., 2007). False positive detections due to Wolbachia contamination from other mosquito host species during the larval stage could also be suspected. False negative detections could be due to PCR primer incompatibility or low bacterial concentration of Wolbachia in the mosquito. For example, a previous study described host age- and sex-related Wolbachia density variances in mosquito bodies (Tortosa et al., 2010). Further data would be required to validate natural Wolbachia infection in Ae. aegypti.
Polymerase chain reaction (PCR) is frequently used to diagnose Wolbachia infection in insects (de Oliveira et al., 2015). Independent PCR tests for multiple Wolbachia genes would be encouraged to reduce the possibility of false negative results. In addition, the authors of previous studies involving PCR assays did not use locally designed primers specific to the tested local Wolbachia populations, potentially resulting in false negative detection due to primer incompatibility if large genetic variations were present in the target genes among local populations. To design such local primers for multiple genes, prior genomic information on the local Wolbachia population would be highly desirable, but such prior information is usually not available.
High-throughput sequencing technologies, such as double digestion restriction site-associated DNA sequencing (ddRAD-Seq), could serve as a powerful alternative to address these limitations. DdRAD-Seq uses two different restriction enzymes (RE) to segment the whole genome of organisms into short fragments (Peterson et al., 2012), and sequences numerous randomly selected DNA fragments in parallel. When using DNA extracted from Ae. aegypti individuals, the genome of the microorganisms present in the mosquito would also be sequenced. The ddRAD-Seq provides genome-wide information without requiring prior knowledge of the local target populations, potentially reducing PCR primer incompatibility-related false negative detection. So far, few studies have used ddRAD-Seq to detect Wolbachia. Lee et al. (2020) used ddRAD-Seq to observe the coevolution of Wolbachia and its host, Anoplolepis gracilipes. Yang et al. (2022) characterized wAlbA and wAlbB Wolbachia strain infection in Ae. albopictus using ddRAD-Seq. However, to the best of our knowledge, no study has applied ddRAD-Seq to detect Wolbachia in Ae. aegypti.
Sequencing the microorganism endosymbiont DNA extracted from Ae. aegypti might allow the detection of Wolbachia DNA sequences. In this study, we examined the feasibility of using ddRAD-Seq for detecting and quantifying Wolbachia from field-collected female and male Ae. aegypti populations in Metropolitan Manila, the Philippines. We assessed the accuracy of Wolbachia detection using ddRAD-Seq compared to the results of PCR assays (both conventional and quantitative PCR) and explored the advantages and limitations of these methods in Wolbachia detection and quantification. We also estimated the genetic diversity of Wolbachia in field-collected Ae. aegypti samples using ddRAD-Seq-derived reads.
Materials and methodsMosquito sampling
Ae. aegypti mosquitoes were collected from Metropolitan Manila, the Philippines. A total of 217 Ae. aegypti individuals (93 males and 124 females) that had been previously used by Carvajal et al. (2020) and Regilme et al. (2021) were used in this study. We assessed male and female populations from seven different regions in Metropolitan Manila (Figure 1). Adult mosquitoes were collected using a UV light trap (Mosquito Trap, Jocanima Corporation, Las Pinas City, Philippines) from May 2014 to January 2015 (Carvajal et al., 2020) and from September to October 2017 (Regilme et al., 2021). Mosquito identification was conducted using pictorial keys from Rueda (2004) and the molecular method using species-specific microsatellite markers undertaken by previous studies (Carvajal et al., 2020; Regilme et al., 2021). The current study used the same DNA samples as the two previous studies (Carvajal et al., 2020; Regilme et al., 2021) to detect and quantify Wolbachia.
Sampling site locations of Ae. aegypti in Manila City, Metropolitan Manila, Philippines. Big circles in red, yellow, dark green, light green, gray, dark blue, and purple indicate the geographical midpoints of Ae. aegypti populations per location; small circles near each big circle indicate the households in the sampling locations. F and M indicate the total number of female and male individuals per population, respectively (Muharromah et al., 2023).
DNA preparation, ddRAD-sequencing, and data processing
The ddRAD-Sequencing data were obtained from (Muharromah et al., 2023). For the library construction, firstly, the DNA sequences of each Ae. aegypti mosquito was determined using a Quantus Fluorometer (Promega, USA). Individual DNA samples of all 14 populations were pooled in equimolar DNA amounts (Pool-Seq) (Schlötterer et al., 2014) based on the sex (female and male) and location (Central, East, West, North, South, Manila North and Manila South). Prior to the library preparation, we optimized in selecting the restriction enzymes for ddRAD-Seq. The restriction enzymes for ddRAD-Seq were selected using two approaches: in silico and empirical approach. We compared seven restriction enzyme combinations (DraI-NlaIII, MluCI-NlaIII (Rašić et al., 2014), DraI-MluCI, SbfI HF-MspI (Sherpa et al., 2018), EcoRI-NlaIII (Rašić et al., 2014) SbfI HF-HaeIII (Gamboa and Watanabe, 2019), and SspI-NlaIII) to produce desired sequenceable DNA fragments of c.a. 100–500 bp, which following the addition of adapters and sequencing primers will result in an acceptable library size for sequencing (c.a. 200–700 bp). Double digestion with in silico analysis allows the prediction of the number of sequenceable DNA fragments using restriction-site information from the enzymes and reference genome of Ae. aegypti via the DDsilico program (Rašić et al., 2014). Empirical digestion analysis is an experimental method for observing DNA fragment distribution using the actual DNA of Ae. aegypti and restriction enzymes visualized with a High-Sensitivity DNA Assay 2100 Bioanalyzer (Agilent, USA). We selected MluCI and NlaIII (New England Biolabs, Beverly MA, USA) as the optimal combination because it generated the highest number of potential ddRAD loci using the in silico and empirical approaches (Supplementary Figures 1, 2). The ddRAD-Sequencing library preparation was performed using the Restriction Enzymes (REs) NlaIII and MluCI (New England Biolabs, USA) (Rašić et al., 2014) to digest the Ae. aegypti DNA for 3 h at 37°C. Next, the REs were inactivated at 65°C for 20 min and purified using a QiaQuick PCR Purification Kit (Qiagen, Hilden, Germany). The digested DNA was ligated to Illumina P1 and P2 adapters using a T4 Ligation mix containing 0.5 µl of 4 nM/µl P1 Adapter, 0.5 µl of 6 nM/µl P2 Adapter, T4 DNA ligase (Takara Bio, Japan), T4 ligase buffer and H2O at 16°C for 16 h with the total volume 15 µl, after which the ligase was inactivated at 65°C for 20 min. The adapter-ligated DNA was amplified in a 10-µl PCR reaction mix containing 5 µl of Phusion High Fidelity Master Mix (New England Biolabs, USA), 2 µl of P1 primer (5’-AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACG-3’), and 2 µl of P2 primer (5’-CAAGCAGAAGACGGCATACGAGATCGTGATGTGACTGGAGTTCAGACGTGTGC- 3’) with the PCR cycling conditions as follows: 98°C for 30 s; 12 cycles of 98°C for 10 s, 60°C for 30 s, and 72°C for 90 s; and elongation at 72°C for 5 min. The final library was formed by pooling seven PCR replicates and purified using a Qiaquick PCR Purification Kit (Qiagen, Hilden, Germany). The library was checked for quality and quantity using Bioanalyzer (Agilent Technologies, USA) and KAPA Quantification kits (Roche, USA). After that, the library was sequenced using a HiSeq X Ten Illumina sequencer (paired-end, 2 × 150 bp) at the Beijing Genomics Institute, China.
The raw sequence data were verified for quality using FASTQC v0.11.8 (Andrews, 2010). The reads were trimmed and filtered to remove the adapters and the barcodes using Trimmomatics 0.39 (Bolger et al., 2014), retaining 100 bp of read length. The reads were mapped to the Wolbachia reference genome (Sinha et al., 2019; Martinez et al., 2022) using the bwa mem algorithm in BWA (Li and Durbin, 2009), generating a SAM format file per population. The ambiguously mapped reads from the mapping were filtered with a minimum MAPQ score of 20. This MAPQ score indicated the possibility that less than 1 out of 100 mappings were incorrect. The SAM files were converted to BAM files using SAMTOOLS 1.9 (Danecek et al., 2021) to have a memory-efficient file form. The reads mapped to the Wolbachia genome were extracted using the samtools view command in SAMTOOLS 1.9. The extracted reads were sorted toward the reference coordinates using SAMTOOLS 1.9. Calling single nucleotide polymorphisms (SNPs) was conducted using bcftools (Li, 2011; Danecek et al., 2021). First, all of population files were merged using bcftools mpileup command. We converted the bcf file to vcf file form using bcftools filter command. The SNPs were filtered using bcftools for minimum quality 20 and minimum read depth of 10. After that, the nucleotide diversity were calculated using vcftools (Danecek et al., 2011) over 10 kb windows of the genome (nucleotide diversity value was estimated for every 10,000 bases across the genome). To identify the species and strain obtained in each mosquito population, the sorted BAM file per population was converted to FASTA file format, then identified using MMSeqs2 version 13.45111 (Steinegger and Söding, 2017) using the UniProtKB/SwissProt database (Bairoch and Apweiler, 2000) released in 27 April 2022 by comparing the amino acids from the query sequence with those from the database and high sensitivity value (-s 9) to improve accuracy. To visualize the mapped read in the Wolbachia genome, we used Proksee (Stothard et al., 2019). The gene annotation was performed using Prokka 1.14.6 (Seemann, 2014) provided in Proksee.
In identifying other bacteria from the samples, we classified the raw reads using Kaiju (Menzel et al., 2016) with the NCBI non-redundant (NR) database and minimum occurrence percentage ≥0.0001%. We calculated the number of Wolbachia contigs by de novo assembly from the filtered data using metaSPAdes (Nurk et al., 2017) and then we classified it using Kaiju with NCBI non-redundant database.
PCR for <italic>Wolbachia</italic> detection and quantification
The PCR-based Wolbachia detection/non-detection data (targeting the 16S rRNA and wsp genes of each individual in the ddRAD analysis) were obtained from Reyes et al. (2022) for samples from 10 populations in Metropolitan Manila. For this, a 16S rRNA and wsp gene marker has been used (Table 1). In addition, the PCR data results of four populations in Manila City were obtained from Regilme et al. (2022) using the same target gene markers (Table 1). In addition, we also performed quantitative PCR (qPCR) assays targeting the wsp gene in the pooled DNA samples of each population for ddRAD using Real-Time Quantitative PCR (Bio-Rad, USA). In this pool-based qPCR, we used primers designed by Reyes et al. (2022) for detection of the wsp gene from each mosquito individual comprising the 14 populations of Metropolitan Manila. Reyes et al. (2022) designed the primers for 118 wsp sequences extracted from Ae. aegypti samples and sequenced by Carvajal et al. (2019) (GenBank popset 1712729902). Next, Multiple Sequence Comparison by Log-Expectation was used for multiple sequence alignment and Codon Code Aligner version 1.2.4 (available at https://www.codoncode.com/aligner/) to display the outcomes. The consensus sequence of the alignment was then used to create wsp primers for the Ae. aegypti samples using Primer-BLAST. Five primer pairs were produced using Primer-BLAST and they were confirmed using a known positive sample of Cx. quinquefasciatus. Two of five primer pairs (wspAAML 01 and wspAAML 05) were chosen from the group for further optimization as they yielded the proper band size of the target markers in the sample without any nonspecific binding. Then, Reyes et al. (2022) established the optimal annealing temperature and primer concentration for both pairings to select the best wspAAML primer pair for further investigation. After careful consideration, wsp 05 was chosen since its PCR efficiency was within the typical MIQE criterion of ≥90%. This approach enabled the design of primers capable of detecting the variable sequences of wsp genes present in the local populations of Wolbachia in this region.
PCR primers used for Wolbachia detection.
Name
Gene
Oligonucleotide sequence (5’-3’)
Probe
Reference
16SF
16S rRNA
5′-AGTGAAGAAGGCCTTTGGG-3′
5′TET-CTGTGAGTACCGTCATTATCTTCCTCACT-BHQ13′
Fraser et al. (2020)
16SR
16S rRNA
5′-CACGGAGTTAGCCAGGACTTC-3′
wspAAML F
wsp
5′-AGCATCTTTTATGGCTGGTGG-3′
5′FAM-ACGACGTTGGTGGTGCAACATTTGC-TAMRA3′
Reyes et al. (2022)
wspAAML R
wsp
5′- AATGCTGCCACACTGTTTGC-3′
WolbF
16S rRNA
5′-GAAGATAATGACGGTACTCAC-3′
Zhou et al. (1998)
Wspecr
16S rRNA
5′-AGCTTC GAGTGAAACCAATTC-3′
wsp 81F
wsp
5′-TGGTCCAATAAGTGATGAAGAAAC-3′
Simoes et al. (2011)
wsp 691R
wsp
5′-AAAAATTAAACGCTACTCCA-3′
17SF
RPS17
5′-TCCGTGGTATCTCCATCAAGCT-3′
5′HEX-CAGGAGGAGGAACGTGAGCGCAG-BHQ13′
Frentiu et al., 2014
17SR
RPS17
5′-CACTTCCGGCACGTAGTTGTC-3′
The relative density of Wolbachia was calculated using the Ae. aegypti ribosomal S17 (RPS17) gene as a housekeeping/reference gene (Table 1). The relative density of Wolbachia was assessed using the delta CT calculation method (where CT refers to the qPCR threshold cycle) as follows: 2 CT(using the RPS17 reference gene)/ 2 CT(target 16S rRNA or wsp genes) (Fraser et al., 2017). For the 16S rRNA gene amplification, we used 5 µl of 1X iTaq mix (Bio-Rad) with 0.2 µl of 0.2 µM 16S rRNA primers. The 16S rRNA reaction was performed in a volume of 10 µl, using 0.15µl of the probe at a concentration of 0.15 µM and completing the reaction with 3.45 µl of H2O. Concerning the wsp amplification, we also used 5 µl of iTaq mix (Bio-Rad) at 1x concentration with 0.5 μM of wsp forward and reverse primers in a volume of 0.5 μl per primer, 0.3 µM of wsp probe at a volume of 0.3 µl, and 2.7 µl of H2O to obtain a total reaction volume of 10 µl. The RPS17 gene was amplified using 5µl of iTaq mix (Bio-Rad) at 1x concentration, 0.3 µl of 0.3 µM forward and reverse RPS17 primers, 0.2 µl of RPS17 gene probe at 0.2 µM, and 3.2 µl of H2O were added to complete the reaction volume to a total of 10 µl. The PCR cycling conditions for the RPS 17 and 16S rRNA genes were as follows: 95°C for 30 s; 95°C for 5 s; 60°C for 10 s with 40 cycles. The PCR cycling conditions for the wsp gene were as follows: 95°C for 2 min; 95°C for 30 s; 58.8°C for 30 s with 40 cycles.
Data analysis
The total number of PCR-detected Wolbachia-positive individuals for the 16S rRNA and wsp genes per population was standardized by dividing it by the total number of analyzed individuals per population. The total number of reads mapped to the Wolbachia genome per population was also standardized by dividing it by the total number of reads remaining after quality filtering using Trimmomatics per population. The correlation between the percentage of ddRAD-Seq reads mapped to the Wolbachia genome per population and the percentage of Wolbachia-positive individuals detected by PCR per population was examined using Spearman’s correlation test in RStudio version 1.4.1106. Similarly, we tested the correlation between the percentage of reads mapped to the Wolbachia genome per population, measured the relative Wolbachia density per population by qPCR assays, and analyzed the total number of individuals by ddRAD per population using Spearman’s correlation test.
Results
<italic>Wolbachia</italic> detection using pooled ddRAD-sequencing
ddRAD-Seq produced a total of 377,047,648 raw reads with an average of 26,931,975 reads per population. After quality filtering and trimming, we obtained a total of 146,239,637 reads with a minimum length of 100 bp. The ddRAD-Seq data showed varying numbers of reads mapped to the Wolbachia wAlbA and wAlbB genomes among the 14 Ae. aegypti populations in Metropolitan Manila (Table 2). Ten populations (Female and Male Central, East, North, South, and West) were confirmed to display Wolbachia genome sequences with more than 100 reads mapped to the Wolbachia genome. However, four populations from Manila City were detected to exhibit a few reads mapped to Wolbachia (< 100 reads). The highest read number mapped to the Wolbachia genome was found in the female South population, followed by the male West and male South populations. We observed a higher total number of reads mapped to the wAlbB genome compared to that to the wAlbA genome. In the following analysis, we only used the reads mapped to the wAlbB genome.
The number of reads obtained by ddRAD-Seq analysis in 14 populations of Ae. aegypti and the number of the reads that mapped to the wAlbA or wAlbB reference genomes and all bacteria classified using Kaiju.
No
Population
N
Reads
wAlbA
wAlbB
All bacteria classified reads using Kaiju
1
F_Central
24
13,278,566
502
553
247,145 (0.74%)
2
F_East
19
10,486,484
86
88
33,201 (0.11%)
3
F_North
12
8,938,915
1,016
1,643
68,807 (0.33%)
4
F_South
28
8,938,153
2,480
14,459
34,386 (0.14%)
5
F_West
18
10,250,794
192
199
68,524 (0.25%)
6
F_North_Manila
12
7,469,088
17
17
2,224 (0.012%)
7
F_South Manila
16
9,624,898
75
77
5,658 (0.02%)
8
M_Central
20
10,250,794
941
1,925
117,544 (0.59%)
9
M_East
12
8,056,812
374
473
57,652 (0.26%)
10
M_North
7
8,572,293
215
211
68,895 (0.38%)
11
M_South
20
14,003,836
8,588
11,509
72,053 (0.24%)
12
M_West
12
11,627,817
11,747
12,547
41,987 (0.14%)
13
M_North_Manila
15
10,280,531
13
13
3,485 (0.013%)
14
M_South_Manila
9
14,460,656
53
64
12,296 (0.03%)
Total
217
146,239,637
26,299
43,778
833,857
N = total number of analyzed individuals, Reads = number of reads obtained after trimming and filtering, wAlbA and wAlbB = total number of reads mapped to the wAlbA and wAlbB reference genomes, respectively.
Detection accuracy of ddRAD-sequencing compared with the PCR assays
A total of 217 samples used in this study were tested using PCR assays both on individual and pooled data. Table 2 shows the results of Wolbachia detection and relative density using the PCR assays. The PCR assay results indicated that 85 and 49 individuals from the 217 samples could be positively detected with Wolbachia using the wsp and 16S rRNA genes, respectively. The relative density of Wolbachia per individual was in the range of 0.0002–147.03 (for 16S rRNA) and 0.0002–64.44 (for wsp). The qPCR results on the pooled DNA samples showed positive results for seven populations (Female Central, North, and South; Male Central, East, South, and West) with the relative density of Wolbachia in the range of 0.010–1.51 (Table 3).
Individual-based detection of Wolbachia in Ae. aegypti using PCR and relative density Wolbachia using individual-based and pool-based estimation.
Population
N
wsp
16S rRNA
wsp, 16S rRNA
X̄ Relative Density (wsp)- Ind
X̄ Relative Density (16S rRNA)- Ind
Relative Density (wsp)-Pool
F_Central
24
15
8
7
0.158
0.015
0.669
F_East
19
6
2
1
0.026
0.002
0
F_North
12
3
2
2
0.194
0.160
0.146
F_South
28
13
13
11
0.019
2.184
0.028
F_West
18
7
4
3
0.003
0.001
0
F_North_Manila
12
0
0
0
0
0
0
F_South Manila
16
0
0
0
0
0
0
M_Central
20
14
6
5
0.039
3.702
0.170
M_East
12
5
2
2
0.036
0.003
0.010
M_North
7
1
1
1
0.002
0.015
0
M_South
20
11
9
9
1.470
3.046
1.505
M_West
12
8
2
2
8.145
73.776
0.551
M_North_Manila
15
0
0
0
0
0
0
M_South_Manila
9
0
0
0
0
0
0
* Total
217
83
49
43
N = total number of individuals, x̄ = mean value.
The percentage of Wolbachia-positive individuals detected by PCR per population showed a positive correlation with the percentage of ddRAD-Seq reads mapped to the Wolbachia genome per population both for the 16S rRNA (Figure 2A, p = 0.0001) and wsp (Figure 2B, p = 0.0002) gene. Furthermore, the percentage of Wolbachia-positive individuals detected using both the 16S rRNA and wsp genes also positively correlated with the percentage of reads mapped to the Wolbachia genome (Figure 2C, p < 0.0001). The relative Wolbachia density per population measured by qPCR also showed a positive correlation with the percentage of ddRAD-Seq reads mapped to the Wolbachia genome both for the 16S rRNA (Supplementary Figure 1A, p < 0.0001) and wsp (Supplementary Figure 1B, p = 0.0004) gene. For the pooled data, the percentage of mapped reads toward the Wolbachia genome showed a positive correlation with the relative Wolbachia density estimated by qPCR (Figure 2D, p = 0.0005). Finally, the percentage of reads mapped to the Wolbachia genome per population did not correlate with the total number of individuals analyzed with ddRAD per population (p > 0.05) (Supplementary Figure 1C).
Correlation plots between the percentage of mapped reads to the Wolbachia genome and that of Wolbachia-positive mosquitoes with 16S rRNA gene (A), wsp gene (B), both 16S rRNA and wsp genes (C), and the relative Wolbachia density in the pooled data (D). Percentage of reads mapped in the Wolbachia genome = total number of reads mapped in the Wolbachia genome divided by the total number of reads after trimming and filtering, percentage of Wolbachia-positive mosquitoes = total number of Wolbachia-positive individuals divided by the total number of individuals per population.
The genetic diversity of <italic>Wolbachia</italic>
A total of 150 SNPs were detected from the cumulative reads mapped to the Wolbachia genome from the 14 Ae. aegypti populations. We observed low nucleotide diversity from these SNPs (π = 0.00000651). A total of 21 regions in the 10,000 bp of sliding windows in the Wolbachia genome showed an SNP number in the range of 2–18. Only two SNPs were found in the conserved 16S rRNA gene region (Figure 3). A high number of SNPs (SNPs > 10) could be observed in 6 of the 21 regions (Supplementary Table 2). These regions were further annotated and the regions of 98,001–99,000 bp (11 SNPs), 634,001–635,000 bp (11 SNPs), 881,001–882,000 bp (10 SNPs), 1,211,001–1,212,000 bp (11 SNPs), 1,370,001–1,371,000 bp (15 SNPs), and 1,431,001–1,432,000 bp (18 SNPs) were located in the RCSc_1, trxB, gph, hypothetical protein, IS982 family transposase ISWpi16, IS481 family transposase ISWpi2 genes, respectively.
The visualization of sequence reads generated based on the ddRAD-Seq approach (green) and gap regions (white) on the Wolbachia pipientis wAlbB complete genome. The green arrow indicates the region of single nucleotide polymorphisms (SNPs) located across the Wolbachia pipientis wAlbB complete genome. Supplementary Table 2 includes detailed information about the 150 SNPs.
From the reads mapped to the Wolbachia genome, we identified the following four Wolbachia strains using MMSeqs2 (wPip, wRi, TRS of Brugia malayi, and wMel). One of the four strains was from two Wolbachia species (Wolbachia endosymbiont Culex quinquefasciatus and Wolbachia pipientis), while another strain (wRi) could not be identified at the species level as the DNA sequence searched from the database did not indicate the species name (Supplementary Table 1).
Diversity of bacteria identified from the ddRAD-Seq reads and number of <italic>Wolbachia</italic> contigs using kaiju
The bacteria diversity in the Ae. aegypti populations were different from each other (Supplementary Table 3). The percentage of all bacteria reads from the raw reads of all 14 population samples is 0.22% with the Wolbachia reads percentage is 0.00327%. We observed the percentage of bacteria in each Ae. aegypti per population is less than 1%. The higher bacteria diversity is in F Central population with the percentage 0.74% and the lowest bacteria diversity is in F North Manila (0.012%) (Table 2). The Wolbachia contigs from our data showed high number in the population of F Central, F North, F South, M Central, M South and M West (Supplementary Table 4).
Discussion
In this study, we tested the feasibility of using ddRAD-Seq for Wolbachia detection and quantification from field-collected Ae. aegypti populations in Metropolitan Manila, the Philippines. Overall, the number of ddRAD-Seq reads mapped to the Wolbachia genome in each population showed a consistent pattern with the results of PCR- and qPCR-based Wolbachia detections and quantifications. As expected, the ddRAD-Seq reads revealed numerous Wolbachia genes across the entire genome. This result suggests that ddRAD-Seq might complement conventional Wolbachia detection PCR assays that rely only on a few DNA markers, thereby providing stronger and more reliable support for genome-wide Wolbachia detection. The ddRAD-Seq approach enabled us to obtain information on a large number of genes randomly sampled from the Wolbachia genome without using Wolbachia-specific primers. Therefore, theoretically, it could be expected to reduce false negative detections due to primer incompatibility for genetically diverse populations using PCR. The Wolbachia detection based on the PCR assay targets only a limited number of loci (e.g., 16S rRNA or wsp). ddRAD-Seq targets a large number of loci randomly selected from the Wolbachia genome, which increases the possibility of detection at any of the loci. In support of this theory, we found Wolbachia sequences from ddRAD-Seq reads in a population in Manila City, where the PCR assays did not detect Wolbachia, although the number of mapped reads was small (<100 reads) (Table 2).In addition, while PCR-based assays did not detect Wolbachia from the populations of Female West and Male North (Table 3), the ddRAD-Seq detected sequences of Culex quinquefasciatus Pel Wolbachia endosymbiont and unclassified Wolbachia from these populations (Supplementary Table 1), which are not 16S rRNA and wsp genes in the Wolbachia genome.
However, the ddRAD-Seq approach also has limitations in detecting and quantifying Wolbachia. One is the possibility of false positive detection. The ddRAD-Seq randomly generates DNA sequence fragments of the organisms. Therefore, ddRAD-Seq reads from other bacteria evolutionarily close to Wolbachia could be mistakenly identified as Wolbachia, especially if the evolutionary rate of that sequence region is low (i.e., no/small interspecific variation). Although we tried to remove such ambiguously mapped reads after mapping the reference genome and increasing the sensitivity criteria in the identification using MMSeqs2, this possibility still cannot be completely excluded. Furthermore, mechanical errors could still occur related to the ddRAD-Seq data generated by the Illumina platform even after quality filtering, and thus the occurrence of erroneous sequences unintentionally identified as Wolbachia cannot be completely prevented. In addition, it is known that Wolbachia transfer its genes to the host genome (Kondo et al., 2002; Sieber et al., 2017). Klasson et al. (2009) investigated the horizontal gene transfer between Wolbachia and the Ae. aegypti genome and concluded that the gene transfer most likely occurs from Wolbachia to the host genome. It is not possible to determine whether the Wolbachia genome sequence detected in this study originated from Wolbachia-infected mosquitoes or from the host-integrated Wolbachia genome. However, it is worth mentioning that this study identified contigs not only from a limited portion of the genomic regions of Wolbachia but also from numerous genome-wide regions (Supplementary Table 4). This finding supports the former possibility because, if the detected contigs were sequences of Wolbachia integrated into the mosquito genome, we would expect to observe fewer number of contigs from a narrower range of genomic regions that are integrated into the mosquito genome. Furthermore, ddRAD-Seq cannot completely eliminate the possibility of Wolbachia contamination from the environment or from the commensal or parasitic species such as nematodes within the mosquito’s body. This method detects Wolbachia based on the presence of DNA fragment from Wolbachia that are sequenced alongside DNA fragments from the host. It is impossible to determine whether the detected DNA fragments are from an authentic Wolbachia infection, contamination from other host species at larval stage or derived from the parasites within the mosquito’s body.
Another limitation of ddRAD-Seq is the less quantitative nature of the data. Using ddRAD-Seq makes accurate estimation of the relative Wolbachia gene concentrations per individual or population theoretically difficult, while it is possible using qPCR. However, in this study, we observed an interesting phenomenon: the number of reads mapped to the Wolbachia genome positively correlated with the number of Wolbachia-infected individuals (Figures 2A–C) or the relative density of Wolbachia in the Ae. aegypti population (Figure 2D; Supplementary Figures 3A, B). This result suggests the possibility of Wolbachia semi-quantification using ddRAD-Seq. However, ddRAD-Seq is theoretically unlikely to reflect the amount of Wolbachia in the template DNA due to PCR bias that might occur during library preparation. Future studies should continue to explore the possibility of using ddRAD-Seq data for the quantification of Wolbachia or other host organism-infecting bacteria (e.g., mosquitoes). For example, in this study, we found a pattern indicating that 3 of the 14 populations exhibited notably high numbers of reads (> 10,000 reads) mapped to the Wolbachia wAlbB genome (Table 2). Therefore, ddRAD-Seq could be potentially used as a tool to explore which local Ae. aegypti populations of Wolbachia could be potentially infected with a high prevalence rate.
Using ddRAD-Seq reads mapped to Wolbachia enabled us to discover significant differences in the levels of genetic diversity or evolutionary rates among different regions across the Wolbachia genome. This is a new discovery that would have likely remained undetected using simple PCR amplification and sequencing of only a part of the genome. The 150 SNP loci were expected to have a relatively fast evolutionary rate in the genome. Six of the 37 regions in the genome exhibited higher SNP numbers and genetic diversity than the other 31 regions (STable 2), suggesting higher mutation rates of the 6 regions (RCSc_1 (catalytic activity), trxB (catalytic activity), gph (phosphoglycolate phosphatase activity), hypothetical protein (function not determined), IS982 family transposase ISWpi16 (function not determined), and IS481 family transposase ISWpi2 (nucleic acid binding) genes). To capture a wide range of genetically diverse Wolbachia strains, it is recommended to analyze a large number of loci, including those with high evolutionary rates (Held and Leese, 2007).
Wolbachia strains detected in this study were as follows: wPip, wRi, TRS of Brugia malayi, and wMel. The ddRAD-Seq reads mapped in multiple genes across the Wolbachia genome, which might contribute to more accurate Wolbachia strain classification. The high evolutionary gene markers are useful for classifying Wolbachia into strains with phylogenetically high resolutions, such as subpopulations within populations. Such phylogenetically finer classification could contribute to unraveling how the arbovirus-blocking effect of different strains (e.g., DENV) might function (Flores et al., 2020), forecast the potential competition among different strains (Liang et al., 2020), and guide mass release programs. We could identify only two SNPs from the ddRAD-Seq reads in the conserved region of 16S rRNA with 3.18E-05 nucleotide diversity, which is known for its low mutation rate and might not be appropriate for phylogenetic analysis (Held and Leese, 2007; Rodrigues and Silva, 2016; ColwelL and Haig, 2019). Only one read hit the wsp gene region, and no SNP was detected in the wsp gene. The wsp gene appears to be a fast-evolving gene marker and an informative gene for discriminating Wolbachia strains. The limited number of SNPs found in the 16S rRNA gene could be attributed to its highly conserved nature. On the other hand, the absence of SNPs detected in the wsp gene may be due to the small number of ddRAD reads mapped to this region, resulting in an insufficient number of sequences for detecting sequence variations. In order to obtain a sufficient number of reads for mapping to the high number of genes across the Wolbachia genome, an optimization in the ddRAD-Seq library preparation (e.g., restriction enzyme selection) should be performed both for the host and the Wolbachia genomes. An alternative to the utilization of 16S rRNA gene marker for Wolbachia detection, a study of Sankar et al. (2021) detected low prevalence of natural Wolbachia from Anopheles culicifacies and Anopheles stephensi using a nested PCR method of 16S rRNA. Regarding the evolutionary rates of the 16S rRNA and wsp markers, the detection of Wolbachia using PCR assays showed different results between targets genes where wsp gene were detected from higher number of individuals (83 individuals) than 16S rRNA (49 individuals) (Table 3). 16S rRNA gene is known to evolve at a slow rate, while the wsp gene is known to evolve at a fast rate. The use of fast evolving gene as a marker increases the possibility of false negative detection due to primer incompatibility in PCR. The different evolutionary rates of the two markers may have caused the difference in detection rates. Other published studies also reported different percentage of positive results for each marker (e.g., Wong et al., 2020) while defining true positive as positive for both. In general, PCR results require the use of at least two markers for detection to avoid bias and increase reliability.
In Metropolitan Manila, male Ae. aegypti yielded more reads that could be mapped to the Wolbachia genome than females (Table 2). This result indicates higher Wolbachia prevalence rates in males. The different feeding behavior and dispersal capabilities of female and male mosquitoes might affect the occurrence of Wolbachia. Female mosquitoes are hematophagous insects and tend to be anthropophagic (feeding on human blood), which might influence the composition of their gut microbiome. Sarma et al. (2022) reported that human blood-fed female Ae. aegypti exhibited higher microbiome species diversity than Ae. aegypti not fed human blood. The increased species diversity in the microbiome might increase competition among microorganisms, making it more difficult for Wolbachia to persist or multiply within the mosquito body. The fact that male mosquitoes prefer to stay close to mating sites and disperse less than females could be another factor limiting bacterial diversity (Minard et al., 2013). A high Wolbachia prevalence rate in four species of male mosquitoes was also identified in the study of Yang et al. (2021) in China using PCR. This study showed higher infection rates of Wolbachia in males than in females in other mosquito species: Ae. albopictus (male = 98.8%; female = 96.5%), Armigeres subalbatus (male = 98.1%; female = 93.2%), Culex pipiens (male = 95.7%; female = 80.4%), and Culex tritaeniorhynchus (male = 100%; female = 5.6%). Kulkarni et al. (2019) reported that male Ae. aegypti exhibited a higher infection rate than female mosquitoes in Florida (female = 3.6%, male = 5.5%) but a lower infection rate than female mosquitoes in New Mexico (female = 58.8%, male = 54.9%).
Conclusions
In this study, we demonstrated that ddRAD-Seq could be applied efficiently for detecting, quantifying, and assessing the genetic diversity of Wolbachia strains in Ae. aegypti populations in Metropolitan Manila, the Philippines. As expected, the ddRAD-Seq reads revealed various Wolbachia genes across the genome. This result suggests that ddRAD-Seq might complement the conventional PCR assays that detect Wolbachia relying only on a few DNA markers and provide more reliable support for genome-wide Wolbachia detection. Moreover, we demonstrated that the number of ddRAD-Seq reads mapped to the Wolbachia genome in each population tended to be consistent with the conventional PCR- and qPCR-based Wolbachia detection results. These results suggest the significance of further validating the quantitative assessment of Wolbachia infection by ddRAD-Seq in future studies. The prevalence and genetic diversity of the Wolbachia strains infecting the mosquito populations, as revealed by ddRAD-Seq, might provide useful insights into the design of a mass release program of Ae. aegypti artificially infected with Wolbachia for mosquito-borne disease control.
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/, PRJNA954465.
Ethics statement
The manuscript presents research on animals that do not require ethical approval for their study.
Author contributions
AM was involved in conceptualizing, designing and performing the experiments, analyzing the data, and writing and editing the original draft. KW contributed to the conceptualization of the research, reviewed and edited the draft, and supervised the work. JR performed the qPCR experiments for the individual-based, reviewed, and edited the draft. NK was involved in supervising the data analysis. All authors contributed to the article and approved the submitted version.
Funding
The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This work was funded by grants from JSPS Core-to-Core Program B. Asia-Africa science platforms (JPJSCCB20190008) and the Sumitomo Electric Industries Group Corporate Social Responsibility (CSR) Foundation, Endowed Chair program.
Acknowledgments
We are grateful to Dr. Thaddeus Carvajal and Dr. Maria Angenica F Regilme for their contribution in collecting Ae. aegypti samples from Metropolitan Manila, Philippines. We also grateful to Dr. Christopher Robinson for the help in editing the manuscript. Finally, we would like to thank Micanaldo Ernesto Francisco for contributing to the Metropolitan Manila map.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcimb.2023.1252656/full#supplementary-material
ReferencesAhmadN. A.ManciniM. V.AntT. H.MartinezJ.KamarulG. M. R.NazniW. A.. (2021). Wolbachia strain wAlbB maintains high density and dengue inhibition following introduction into a field population of Aedes aEgypti. Philos. Trans. R Soc. Lond B Biol. Sci.376 (1818), 20190809. doi: 10.1098/rstb.2019.0809AliotaM. T.PeinadoS. A.VelezI. D.OsorioJ. E. (2016a). The wMel strain of Wolbachia Reduces Transmission of Zika virus by. Aedes aEgypti. Sci. Rep.6, 28792. doi: 10.1038/srep28792AliotaM. T.WalkerE. C.Uribe YepesA.Dario VelezI.ChristensenB. M.OsorioJ. E. (2016b). The wMel Strain of Wolbachia Reduces Transmission of Chikungunya Virus in Aedes aEgypti. PloS Negl. Trop. Dis.10 (4), e0004677. doi: 10.1371/journal.pntd.0004677AndrewsS. (2010). FastQC: A quality control tool for high throughput sequence data. Available at: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/.AntT. H.SinkinsS. P. (2018). A Wolbachia triple-strain infection generates self-incompatibility in Aedes albopictus and transmission instability in Aedes aEgypti. Parasites Vectors11, 295. doi: 10.1186/s13071-018-2870-0BairochA.ApweilerR. (2000). The SWISS-PROT protein sequence database and its supplement TrEMBL in 2000. Nucleic Acids Res.28 (1), 45–48. doi: 10.1093/nar/28.1.45BalajiS.JayachandranS.PrabagaranS. R. (2019). Evidence for the natural occurrence of Wolbachia in Aedes aEgypti mosquitoes. FEMS Microbiol. Lett.366 (6), fnz055. doi: 10.1093/femsle/fnz055BeebeN. W.PagendamD.TrewinB. J.BoomerA.BradfordM.FordA.. (2021). Releasing incompatible males drives strong suppression across populations of wild and Wolbachia-carrying Aedes aEgypti in Australia. Proc. Natl. Acad. Sci. U S A.118 (41), e2106828118. doi: 10.1073/pnas.2106828118BennettK. L.Gómez-MartínezC.ChinY.SaltonstallK.McMillanW. O.RoviraJ. R.. (2019). Dynamics and diversity of bacteria associated with the disease vectors Aedes aEgypti and Aedes albopictus. Sci. Rep.9 (1), 12160. doi: 10.1038/s41598-019-48414-8BolgerA. M.LohseM.UsadelB. (2014). Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics.30 (15), 2114–2120. doi: 10.1093/bioinformatics/btu170CarvajalT. M.HashimotoK.HarnandikaR. K.AmalinD. M.WatanabeK. (2019). Detection of Wolbachia in field-collected Aedes aEgypti mosquitoes in metropolitan Manila, Philippines. Parasit Vectors.12 (1), 361. doi: 10.1186/s13071-019-3629-yCarvajalT. M.OgishiK.YaegeshiS.HernandezL. F. T.ViacrusisK. M.HoH. T.. (2020). Fine-scale population genetic structure of dengue mosquito vector, Aedes aEgypti, in Metropolitan Manila, Philippines. PloS Negl. Trop. Dis.14 (5), e0008279. doi: 10.1371/journal.pntd.0008279Chouin-CarneiroT.AntT. H.HerdC.LouisF.FaillouxA. B.SinkinsS. P. (2019). Wolbachia strain wAlbA blocks Zika virus transmission in Aedes aEgypti. Med. Vet. Entomol.34 (1), 116–119. doi: 10.1111/mve.12384ColwelLM. A.HaigS. M. (2019). The Population Ecology and Conservation of Charadrius Plovers (United Kingdom: CRC Press).CoonK. L.BrownM. R.StrandM. R. (2016). Mosquitoes host communities of bacteria that are essential for development but vary greatly between local habitats. Mol. Ecol.25 (22), 5806–5826. doi: 10.1111/mec.13877CrawfordJ. E.ClarkeD. W.CriswellV.DesnoyerM.CornelD.DeeganB.. (2020). Efficient production of male Wolbachia-infected Aedes aEgypti mosquitoes enables large-scale suppression of wild populations. Nat. Biotechnol.38 (4), 482–492. doi: 10.1038/s41587-020-0471-xda MouraA. J. F.ValadasV.Da Veiga LealS.Montalvo SabinoE.SousaC. A.PintoJ. (2023). Screening of natural Wolbachia infection in mosquitoes (Diptera: Culicidae) from the Cape Verde islands. Parasit Vectors.16 (1), 142. doi: 10.1186/s13071-023-05745-wDanecekP.AutonA.AbecasisG.AlbersC. A.BanksE.DePristoM. A.. (2011). The variant call format and VCFtools. Bioinformatics.27 (15), 2156–2158. doi: 10.1093/bioinformatics/btr330DanecekP.BonfieldJ. K.LiddleJ.MarshallJ.OhanV.PollardM. O.. (2021). Twelve years of SAMtools and BCFtools. Gigascience10 (2), giab008. doi: 10.1093/gigascience/giab008de OliveiraC. D.GonçalvesD. S.BatonL. A.ShimabukuroP. H.CarvalhoF. D.MoreiraL. A. (2015). Broader prevalence of Wolbachia in insects including potential human disease vectors. Bull. Entomol Res.05 (3), 305–315. doi: 10.1017/S0007485315000085DingH.YeoH.PuniamoorthyN. (2020). Wolbachia infection in wild mosquitoes (Diptera: Culicidae): implications for transmission modes and host-endosymbiont associations in Singapore. Parasit Vectors13 (1), 612. doi: 10.1186/s13071-020-04466-8FloresH. A.Taneja de BruyneJ.O’DonnellT. B.Tuyet NhuV.Thi GiangN.Thi Xuan TrangH.. (2020). Multiple Wolbachia strains provide comparative levels of protection against dengue virus infection in Aedes aEgypti. PloS Pathog.16 (4), e1008433. doi: 10.1371/journal.ppat.1008433FraserJ. E.De BruyneJ. T.Iturbe-OrmaetxeI.StepnellJ.BurnsR. L.FloresH. A.. (2017). Novel Wolbachia-transinfected Aedes aEgypti mosquitoes possess diverse fitness and vector competence phenotypes. PloS Pathog.13 (12), e1006751. doi: 10.1371/journal.ppat.1006751FraserJ. E.O’DonnellT. B.DuyvestynJ. M.O’NeillS. L.SimmonsC. P.FloresH. A. (2020). Novel phenotype of Wolbachia strain wPip in Aedes aEgypti challenges assumptions on mechanisms of Wolbachia-mediated dengue virus inhibition. PloS Pathog.16 (7), e1008410. doi: 10.1371/journal.ppat.1008410FrentiuF. D.ZakirT.WalkerT.PopoviciJ.PykeA. T.van den HurkA.. (2014). Limited dengue virus replication in field-collected Aedes aEgypti mosquitoes infected with Wolbachia. PloS Negl. Trop. Dis.8 (2), e2688. doi: 10.1371/journal.pntd.0002688GamboaM.WatanabeK. (2019). Genome-wide signatures of local adaptation among seven stoneflies species along a nationwide latitudinal gradient in Japan. BMC Genomics20 (1), 84. doi: 10.1186/s12864-019-5453-3Gloria-SoriaA.ChiodoT. G.PowellJ. R. (2018). Lack of evidence for natural wolbachia infections in aedes aEgypti (Diptera: culicidae). J. Med. Entomol.55 (5), 1354–1356. doi: 10.1093/jme/tjy084HeldC.LeeseF. (2007). The utility of fast evolving molecular markers for studying speciation in the Antarctic benthos. Polar Biol.30 (4), 513–521. doi: 10.1007/s00300-006-0210-xHoffmannA. A.MontgomeryB. L.PopoviciJ.Iturbe-OrmaetxeI.JohnsonP. H.MuzziF.. (2011). Successful establishment of Wolbachia in Aedes populations to suppress dengue transmission. Nature.476 (7361), 454–457. doi: 10.1038/nature10356HugoL. E.RašićG.MaynardA. J.AmbroseL.LiddingtonC.ThomasC. J. E.. (2022). Wolbachia wAlbB inhibit dengue and Zika infection in the mosquito Aedes aEgypti with an Australian background. PloS Negl. Trop. Dis.16 (10), e0010786. doi: 10.1371/journal.pntd.0010786KlassonL.KambrisZ.CookP. E.WalkerT.SinkinsS. P. (2009). Horizontal gene transfer between Wolbachia and the mosquito Aedes aegypti. BMC Genomics10, 33. doi: 10.1186/1471-2164-10-33KondoN.NikohN.IjichiN.ShimadaM.FukatsuT. (2002). Genome fragment of Wolbachia endosymbiont transferred to X chromosome of host insect. Proc. Natl. Acad. Sci. U. S. A.99 (22), 14280–14285. doi: 10.1073/pnas.222228199KulkarniA.YuW.JiangJ.SanchezC.KarnaA. K.MartinezK. J. L.. (2019). Wolbachia pipientis occurs in Aedes aEgypti populations in New Mexico and Florida, USA. Ecol. Evol.9 (10), 6148–6156. doi: 10.1002/ece3.5198LeeC. C.LinC. Y.TsengS. P.MatsuuraK.YangC. S. (2020). Ongoing coevolution of wolbachia and a widespread invasive ant. Anoplolepis gracilipes. Microorganisms.8 (10), 1569. doi: 10.3390/microorganisms8101569LiH. (2011). A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data. Bioinformatics.27 (21), 2987–2993. doi: 10.1093/bioinformatics/btr509LiH.DurbinR. (2009). Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics.25 (14), 1754–1760. doi: 10.1093/bioinformatics/btp324LiangX.LiuJ.BianG.XiZ. (2020). Wolbachia inter-strain competition and inhibition of expression of cytoplasmic incompatibility in mosquito. Front. Microbiol.11. doi: 10.3389/fmicb.2020.01638MartinezJ.AntT. H.MurdochyS. M.TongL.da Silva FilipeA.SinkinsS. P. (2022). Genome sequencing and comparative analysis of Wolbachia strain wAlbA reveals Wolbachia-associated plasmids are common. PloS Genet.18 (9), e1010406. doi: 10.1371/journal.pgen.1010406MenzelP.NgK. L.KroghA. (2016). Fast and sensitive taxonomic classification for metagenomics with Kaiju. Nat. Commun.7, 11257. doi: 10.1038/ncomms11257MinardG.MavinguiP.MoroC. V. (2013). Diversity and function of bacterial microbiota in the mosquito holobiont. Parasites Vectors6, 146. doi: 10.1186/1756-3305-6-146MoutonL.HenriH.CharifD.BoulétreauM.VavreF. (2007). Interaction between host genotype and environmental conditions affects bacterial density in Wolbachia symbiosis. Biol. Lett.3 (2), 210–213. doi: 10.1098/rsbl.2006.0590MuharromahA. F.CarvajalT. M.RegilmeM. A.WatanabeK. (2023). Finescale adaptive divergence of Aedes aEgypti in heterogeneous landscapes and among climatic conditions inMetropolitan Manila, Philippines. bioRxiv (Under Review). doi: 10.1101/2023.04.12.536525NazniW. A.HoffmannA. A.NoorAfizahA.CheongY. L.ManciniM. V.GoldingN.. (2019). Establishment of Wolbachia Strain wAlbB in Malaysian Populations of Aedes aEgypti for Dengue Control. Curr. Biol.29 (24), 4241–4248.e5. doi: 10.1016/j.cub.2019.11.007NguyenT. H.NguyenH. L.NguyenT. Y.VuS. N.TranN. D.LeT. N.. (2015). Field evaluation of the establishment potential of wMelPop Wolbachia in Australia and Vietnam for dengue control. Parasit Vectors.8, 563. doi: 10.1186/s13071-015-1174-xNurkS.MeleshkoD.KorobeynikovA.PevznerP. A. (2017). metaSPAdes: a new versatile metagenomic assembler. Genome Res.27 (5), 824–834. doi: 10.1101/gr.213959.116PetersonB. K.WeberJ. N.KayE. H.FisherH. S.HoekstraH. E. (2012). Double digest RADseq: an inexpensive method for De Novo SNP discovery and genotyping in model and non-model species. PloS One7 (5), e37135. doi: 10.1371/journal.pone.0037135PintoS. B.RibackT. I. S.SylvestreG.CostaG.PeixotoJ.DiasF. B. S.. (2021). Effectiveness of Wolbachia-infected mosquito deployments in reducing the incidence of dengue and other Aedes-borne diseases in Niterói, Brazil: A quasi-experimental study. PloS Negl. Trop. Dis.15 (7), e0009556. doi: 10.1371/journal.pntd.0009556RašićG.FilipovićI.WeeksA. R.HoffmanA. A. (2014). Genome-wide SNPs lead to strong signals of geographic structure and relatedness patterns in the major arbovirus vector, Aedes aEgypti. BMC Genomics15, 275. doi: 10.1186/1471-2164-15-275RegilmeM. A. F.CarvajalT. M.HonnenA. C.AmalinD. M.WatanabeK. (2021). The influence of roads on the fine-scale population genetic structure of the dengue vector Aedes aEgypti (Linnaeus). PloS Negl. Trop. Dis.15 (2), e0009139. doi: 10.1371/journal.pntd.0009139RegilmeM. A. F.InukaiT.WatanabeK. (2022). Detection and phylogeny of Wolbachia in field-collected Aedes albopictus and Aedes aEgypti from Manila City, Philippines. Eur. J. Mol. Clin. Med.9 (03), 3060–3073.ReyesJ. I. L.SuzukiT.SuzukiY.WatanabeK. (2022). Discovery of prevalent natural Wolbachia in Aedes aEgypti in Metropolitan Manila, Philippines using locally designed primers: bacterial density is influenced by strain and host sex. bioRxiv. doi: 10.1101/2022.12.23.521724. [Preprint].RodriguesA. E. T.SilvaT. B. (2016). Molecular Diversity of Environmental Prokaryotes (United States: CRC Press).RuedaL. M. (2004). Zootaxa 589: Pictorial Keys for the Identification of Mosquitoes (Diptera: Culicinidae) Associated with Dengue Virus Transmission (Auckland: Magnolia Press), 60. doi: 10.11646/zootaxa.589.1.1RyanP. A.TurleyA. P.WilsonG.HurstT. P.RetzkiK.Brown-KenyonJ.. (2020). Establishment of wMel Wolbachia in Aedes aegypti mosquitoes and reduction of local dengue transmission in Cairns and surrounding locations in northern Queensland, Australia. Gates Open Res.3, 1547. doi: 10.12688/gatesopenres.13061.2SankarS. G.SundariT. W.AnandA. A. P. (2021). First report on the presence of natural Wolbachia population from major malarial vector mosquitoes Anopheles culicifacies s.l., and Anopheles stephensi from Tamil Nadu, India. bioRxiv. doi: 10.1101/2020.11.22.393652. [Preprint].SarmaD. K.KumarM.DhurveJ.PalN.SharmaP.JamesM. M.. (2022). Influence of host blood meal source on gut microbiota of wild caught aedes aEgypti, a dominant arboviral disease vector. Microorganisms10 (2), 332. doi: 10.3390/microorganisms10020332SchlöttererC.ToblerR.KoflerR.NolteV. (2014). Sequencing pools of individuals - mining genome-wide polymorphism data without big funding. Nat. Rev. Genet.11, 749–763. doi: 10.1038/nrg3803SchmidtT. L.BartonN. H.RašićG.TurleyA. P.MontgomeryB. L.Iturbe-OrmaetxeI.. (2017). Local introduction and heterogeneous spatial spread of dengue-suppressing Wolbachia through an urban population of Aedes aEgypti. PloS Biol.15 (5), e2001894. doi: 10.1371/journal.pbio.2001894SeemannT. (2014). Prokka: rapid prokaryotic genome annotation. Bioinformatics.30 (14), 2068–2069. doi: 10.1093/bioinformatics/btu153SherpaS.RiouxD.GoindinD.FouqueF.FrançoisO.DesprésL. (2018). At the origin of a worldwide invasion: unraveling the genetic makeup of the caribbean bridgehead populations of the dengue vector aedes aEgypti. Genome Biol. Evol.10 (1), 56–71. doi: 10.1093/gbe/evx267SieberK. B.BromleyR. E.Dunning HotoppJ. C. (2017). Lateral gene transfer between prokaryotes and eukaryotes. Exp. Cell Res.358 (2), 421–426. doi: 10.1016/j.yexcr.2017.02.009SilvaN. M.SantosN. C.MartinsI. C. (2020). Dengue and zika viruses: epidemiological history, potential therapies, and promising vaccines. Trop. Med. Infect. Dis.5 (4), 150. doi: 10.3390/tropicalmed5040150SimoesP. M.MialdeaG.ReissD.SagotM. F.CharlatS. (2011). Wolbachia detection: an assessment of standard PCR protocols. Mol. Ecol. Resources.11 (3), 567–572. doi: 10.1111/j.1755-0998.2010.02955.xSinhaA.LiZ.SunL.CarlowC. K. S. (2019). Complete Genome Sequence of the Wolbachia wAlbB Endosymbiont of Aedes albopictus. Genome Biol. Evol.11 (3), 706–720. doi: 10.1093/gbe/evz025SteineggerM.SödingJ. (2017). MMseqs2 enables sensitive protein sequence searching for the analysis of massive data sets. Nat. Biotechnol.11), 1026–1028. doi: 10.1038/nbt.3988StothardP.GrantJ. R.Van DomselaarG. (2019). Visualizing and comparing circular genomes using the CGView family of tools. Brief Bioinform.20 (4), 1576–1582. doi: 10.1093/bib/bbx081TeoC. H. J.LimP. K. C.VoonK.MakJ. W. (2017). Detection of dengue viruses and Wolbachia in Aedes aEgypti and Aedes albopictus larvae from four urban localities in Kuala Lumpur, Malaysia. Trop. Biomed.34 (3), 583–597.ThongsripongP.ChandlerJ. A.GreenA. B.KittayapongP.WilcoxB. A.KapanD. D.. (2017). Mosquito vector-associated microbiota: Metabarcoding bacteria and eukaryotic symbionts across habitat types in Thailand endemic for dengue and other arthropod-borne diseases. Ecol. Evol.8 (2), 1352–1368. doi: 10.1002/ece3.3676TorresR.HernandezE.FloresV.RamirezJ. L.JoyceA. L. (2020). Wolbachia in mosquitoes from the central valley of california, USA. Parasit Vectors.13 (1), 558. doi: 10.1186/s13071-020-04429-zTortosaP.CharlatS.LabbéP.DehecqJ. S.BarréH.WeillM. (2010). Wolbachia age-sex-specific density in Aedes albopictus: a host evolutionary response to cytoplasmic incompatibility? PloS One5 (3), e9700. doi: 10.1371/journal.pone.0009700UtariniA.IndrianiC.AhmadR. A.TantowijoyoW.ArguniE.AnsariM. R.. (2021). Efficacy of wolbachia-infected mosquito deployments for the control of dengue. N Engl. J. Med.384 (23), 2177–2186. doi: 10.1056/NEJMoa2030243WalkerT.JohnsonP. H.MoreiraL. A.Iturbe-OrmaetxeI.FrentiuF. D.McMenimanC. J.. (2011). The wMel Wolbachia strain blocks dengue and invades caged Aedes aEgypti populations. Nature.476 (7361), 450–453. doi: 10.1038/nature10355WongM. L.LiewJ. W. K.WongW. K.PramasivanS.Mohamed HassanN.Wan SulaimanW. Y.. (2020). Natural Wolbachia infection in field-collected Anopheles and other mosquito species from Malaysia. Parasit. Vectors13 (1), 414. doi: 10.1186/s13071-020-04277-xYangQ.ChungJ.RobinsonK. L.SchmidtT. L.RossP. A.LiangJ.. (2022). Sex-specific distribution and classification of Wolbachia infections and mitochondrial DNA haplogroups in Aedes albopictus from the Indo-Pacific. PloS Negl. Trop. Dis.16 (4), e0010139. doi: 10.1371/journal.pntd.0010139YangY.HeY.ZhuG.ZhangJ.GongZ.HuangS.. (2021). Prevalence and molecular characterization of Wolbachia in field-collected Aedes albopictus, Anopheles sinensis, Armigeres subalbatus, Culex pipiens and Cx. tritaeniorhynchus in China. PloS Negl. Trop. Dis.15 (10), e0009911. doi: 10.1371/journal.pntd.0009911ZhangH.GaoJ.MaZ.LiuY.WangG.LiuQ.. (2022). Wolbachia infection in field-collected Aedes aEgypti in Yunnan Province, southwestern China. Front. Cell Infect. Microbiol.12, 1082809. doi: 10.3389/fcimb.2022.1082809ZhengX.ZhangD.LiY.YangC.WuY.LiangX.. (2019). Incompatible and sterile insect techniques combined eliminate mosquitoes. Nature.572 (7767), 56–61. doi: 10.1038/s41586-019-1407-9ZhouW.RoussetF.O'NeilS. (1998). Phylogeny and PCR-based classification of Wolbachia strains using wsp gene sequences. Proc. Biol. Sci.265 (1395), 509–515. doi: 10.1098/rspb.1998.0324