The queenslandensis and the type Form of the Dengue Fever Mosquito (Aedes aegypti L.) Are Genomically Indistinguishable

Background The mosquito Aedes aegypti (L.) is a major vector of viral diseases like dengue fever, Zika and chikungunya. Aedes aegypti exhibits high morphological and behavioral variation, some of which is thought to be of epidemiological significance. Globally distributed domestic Ae. aegypti have often been grouped into (i) the very pale variety queenslandensis and (ii) the type form. Because the two color forms co-occur across most of their range, there is interest in understanding how freely they interbreed. This knowledge is particularly important for control strategies that rely on mating compatibilities between the release and target mosquitoes, such as Wolbachia releases and SIT. To address this question, we analyzed nuclear and mitochondrial genome-wide variation in the co-occurring pale and type Ae. aegypti from northern Queensland (Australia) and Singapore. Methods/Findings We typed 74 individuals at a 1170 bp-long mitochondrial sequence and at 16,569 nuclear SNPs using a customized double-digest RAD sequencing. 11/29 genotyped individuals from Singapore and 11/45 from Queensland were identified as var. queenslandensis based on the diagnostic scaling patterns. We found 24 different mitochondrial haplotypes, seven of which were shared between the two forms. Multivariate genetic clustering based on nuclear SNPs corresponded to individuals’ geographic location, not their color. Several family groups consisted of both forms and three queenslandensis individuals were Wolbachia infected, indicating previous breeding with the type form which has been used to introduce Wolbachia into Ae. aegypti populations. Conclusion Aedes aegypti queenslandensis are genomically indistinguishable from the type form, which points to these forms freely interbreeding at least in Australia and Singapore. Based on our findings, it is unlikely that the presence of very pale Ae. aegypti will affect the success of Aedes control programs based on Wolbachia-infected, sterile or RIDL mosquitoes.


Introduction
The mosquito Aedes aegypti (Linnaeus) is the most important arboviral vector in the tropics and subtropics [1]. Diseases transmitted by Ae. aegypti, like dengue fever and Zika, are on the rise [2], and some are reappearing. For instance, chikungunya has returned to the American tropics in 2013, after being absent for nearly 200 years [3]. Yellow fever was nearly eliminated thanks to an effective vaccine, but is now resurging in central and south Africa [4]. Such epidemiological trends highlight the need to persist with vector control efforts, which requires a thorough understanding of vector biology.
Nearly 60 years ago, Mattingly [5] noted that despite a vast body of literature, few mosquitoes have been"the subject of misconception. . ..in the minds of the general run of entomologists" like Aedes aegypti [5]. The species has a plethora of historical synonyms [6], mainly as a result of having extensive variation in body color and scaling patterns [7] which was also thought to correlate with behavioral differences (e.g. [8]). These issues motivated Mattingly [5] to revise the taxonomy of Ae. aegypti and create a foundation for the modern studies of this disease vector.
Mattingly [5] proposed the intraspecific classification of Ae. aegypti into three forms.
1. A very dark form that never has pale scales on the first abdominal tergite, avoids biting humans, prefers natural breeding habitats and is confined to sub-Saharan Africa. Mattingly gave this form a subspecies rank, Ae. aegypti spp. formosus (Walker).

2.
Ae. aegypti sensu stricto or the type form, distinctly paler and browner than spp. formosus, with pale scales restricted to the head and the first abdominal tergite. This form prefers to bite humans and to use artificial breeding containers, and is globally distributed. 3. A very pale form, Ae. aegypti queenslandensis (Theobold), with extension of the pale scaling on the thorax, tergites and legs, that co-occurs with the type form. Mattingly gave this form only a varietal rank (Ae. aegypti var. queenslandensis). This form was very common in the Mediterranean basin before it was eradicated from the region [9]. Aedes aegypti queenslandensis has also been considered the most domestic of the three forms, always breeding and resting very close to humans, both in and outside Africa [5,8].
A few years later, McClelland [7] reported a high level of variation in color and scaling within and among Ae. aegypti populations, suggesting that subdivision into forms seems oversimplistic and should be abandoned unless correlation between genetic and color variation can be demonstrated. In their latest review of the Ae. aegypti history, Powell and Tabachnick [9] pointed out that McClelland's recommendations have often been ignored for the past 45 years, despite the fact that multiple genetic marker systems (allozymes, microsatellites, nuclear and mitochondrial SNPs) have failed to find a clear differentiation between forms and markers [10][11][12][13].
Recently, Chan et al. [14] suggested that the DNA barcoding technique can be used to distinguish queenslandensis individuals from the type individuals in Singapore. The sequence divergence of 1.5%-1.9% between the two forms [14], although lower than a commonly adopted threshold of 3% for species delineation in insects [15], suggests that the two forms may not freely interbreed. Historical records indicate that the two forms have co-occurred in Singapore and other parts of south-east Asia and Australia for hundreds of generations [5,8]. In sympatry, genetic isolation can be maintained largely through pre-zygotic isolation mechanisms like incompatibilities in mating behavior [16]. For instance, molecular forms of the malarial mosquito, Anopheles gambiae, fly together in mating swarms but rarely hybridize due to flight-tone matching between males and females of the same form [17]. Similar incompatibilities in Ae. aegypti would have implications for control strategies that rely on successful mating between the release and target mosquitoes, like Wolbachia-based population replacement and suppression [18,19], releases of sterile males [20] or males with a RIDL construct [21].
To explore this further, we analyzed nuclear and mitochondrial genome-wide variation in the co-occurring pale and type Ae. aegypti from Singapore and northern Queensland (Australia). The RADseq approach we employed allows for detection of genetic structure and ancestry with power unparalleled by previous genetic studies of the Ae. aegypti forms [22]. Any association between genetic structuring (nuclear/mitochondrial) and the mosquito color/scaling would provide support for the hypothesis of restricted interbreeding between the type and the queenslandensis form, with implications for the implementation of biocontrol programs to suppress diseases transmitted by Ae. aegypti.

Ethics statement
The collection of wild mosquitoes in the study areas does not require specific field ethics approval. The sampling was not conducted on protected land, nor did it involve endangered or protected species. Consent was obtained from residents at each location where collections occurred on private property.

Sampling and identification
In Singapore, all samples were collected as larvae from the domestic breeding containers at nine locations during the second week of April 2015 (Fig 1, Table 1). These samples were collected during routine inspection by enforcement officers of the National Environment Agency (NEA), Singapore. Larvae were reared to the adult stage under standard laboratory conditions (25°± 1°C, 80 ± 10% relative humidity and 12 h light/dark cycle). In Townsville (northern Queensland), samples were collected as adults using Biogents Sentinel traps placed at 55 locations in January 2014 (Fig 1, Table 1). Adult mosquitoes were sexed and identified to form based on the key diagnostic color and scaling features, following Mattingly [5] and McClleland [7]. Eleven out of 44 mosquitoes (25%) from Singapore, and seven out of 99 mosquitoes (7%) from Townsville were identified as the queenslandensis form (Table 1). An additional four queenslandensis individuals collected in Cairns (northern Queensland) in December 2014 were included in the analyses (Table 1).

RADseq genotyping
DNA was extracted from 29 individuals collected in Singapore (18 female type, 11 female queenslandensis) and 45 individuals from northern Queensland (17 male type, 17 female type, 11 female queenslandensis) ( Table 1). Qiagen Blood and Tissue DNA kit (Venlo, Limburg, NL) was used to extract DNA from a whole adult mosquito. 100 ng of DNA from each individual was used to construct the double-digest RAD library following a previously validated protocol [22]. In short, 100 units of the two frequently cutting enzymes (MluCI and NlaIII, New England Biolabs, Beverly MA, USA) were used to digest 100 ng of DNA during three hours of incubation at 37°C. 100 pM P1 and 300 pM P2 Illumina adapters with customized barcode sequences were ligated to the genomic fragments using 100 units of T4 ligase at 16°C overnight (New England Biolabs, Beverly, MA, USA). Pooled ligations were purified and size selected for fragments 300-450bp in length, using the 2% Pippin Prep cassette (Sage Sciences, Beverly, MA, USA). The final libraries (one for each geographic region) were enriched with 12 PCR cycles with standard Illumina primers and then sequenced in two HiSeq2500 lanes with the 100 bp paired-end chemistry.
Raw fastq sequences were processed within our customized pipeline. First, all reads were trimmed to the same length of 90 bp and removed if the base quality score was below 13 (FASTX Toolkit, http://hannonlab.cshl.edu/fastx_toolkit/index.html).High quality reads were then aligned to the reference mitochondrial genome [23] and the nuclear genome version AaegL1 [24] using the aligner Bowtie [25]. Uniquely aligned reads were passed to the refmap.pl program that runs the Stacks v.1.35 pipeline [26]. In addition to the samples from Singapore,  Table 1. Sample information. Sample ID, region (SNP-Singapore, QLD_T-Townsville, QLD_C-Cairns, Queensland), X, Y (longitude/latitude decimal degrees), collection (method/breeding container), sex (F-female, M-male), form (t-type, q-queenslandensis [5][7]), mitochondrial haplotype (mt hapl, Hap1-24), per individual proportion of heterozygous (het) nuclear loci, average (aver) locus depth, and proportion of missing (miss) loci. Townsville and Cairns, we included previously sequenced individuals: 15 from Rio de Janeiro (Brazil) [27], 15 from Gordonvale (northern Queensland), and 15 from Ho Chi Minh City (Vietnam) (S1 Table). This was done to compare the extent of genetic structuring within and among samples at a regional and global scale. Sexing of the larval samples from Brazil and Vietnam could not be done based on the external morphological features, so we employed a genetic sexing method based on the presence/absence of the male-specific RAD tags [28]. All 119 individuals were included in the creation of the RAD tag catalogues using the default Stacks parameters in the maximum likelihood model of SNP and genotype calling. The populations module was used to filter the catalogues and export data in the FASTA format (for the mitochondrial variation) and the variant calling format (VCF, for the nuclear variation).

Analyses of genetic diversity and structure
The mitochondrial haplotype richness within and among groups (Ae. aegypti forms and geographic regions) was calculated using the rarefaction method implemented in the program HPrare [29]. Phylogenetic relationship among mitochondrial haplotypes was estimated with the maximum likelihood approach in the program RAxML (GTRM + G, rapid bootstrap heuristic algorithm and thorough ML search) [30]. Haplotypes of three related Aedes species, for which  [32]. Haplotype sequence of the Ae. aegypti reference line (Liverpool, NC_010241.1) was also included in the analysis.
Parameters of data quality and diversity (RAD tag depth, percentage of missing data, heterozygosity averaged per individual) were compared between females of the two co-occurring forms using independent sample t-test. The level of nuclear genetic structuring was estimated using the non-spatial multivariate method DAPC [33] in the R package adegenet [34]. Rousset's genetic distance (â) and geographic distance between pairs of individuals were calculated in the program spagedi [35]. Color distance between pairs of individuals was treated as a binary value: 0 (same color/form) and 1 (different color/form).

Variation and phylogenetic relationship among mitochondrial haplotypes
From the mitochondrial RAD tag catalogue, we extracted 13 polymorphic tags that were shared between at least 80% of individuals (60/74, Table 1). Tags were distributed across eight different mitochondrial genes (COXI, Cytb, ATP6, ND1-2, ND4-6; S1 File). All 13 tags were concatenated into a final 1170 bp long sequence that was treated as a mitochondrial haplotype. We found 24 different haplotypes in samples from Singapore and Townsville. Haplotype richness did not differ between the two forms in either location (Singapore type = 5.13, queenslandensis = 5.07; Townsville type = 4.19, queenslandensis = 5.0). Moreover, seven haplotypes were shared between the two forms (Table 1, Fig 2).
There were 207 distinctive alignment patterns and 8.17% of undetermined characters in the dataset consisting of 24 haplotypes from Singapore and Queensland, one from the Liverpool  (Hap1-24) found in Aedes aegypti type and var. queenslandensis that co-occur in Singapore and northern Queensland, Australia. Sequences of the three outgroups (Ae. albopictus, Ae. vigilax, Ae. notoscriptus) and Ae. aegypti Liverpool strain were obtained from the NCBI nucleotide sequence/genome database, with the NCBI accession numbers listed in square brackets. The number of Ae. aegypti individuals with a given mitochondrial haplotype is listed in parentheses. A circle designates haplotypes found in Queensland, and a triangle those found in Singapore. Open symbols designate haplotypes found in the queenslandensis form, and filled symbols those found in the type form. Very Pale Ae. aegypti Are Not Genomically Separate strain and three from other Aedes species (outgroups). A phylogeny based on maximum likelihood revealed two highly statistically supported maternal lineages in Ae. aegypti: a basal clade (more similar to the outgroups) and a clade arising from it (a derived clade) (Fig 2). Nucleotide distance (p-distance) between the two clades ranged from 1.2% to 1.6% (S2 File). Importantly, haplotypes of the two Ae. aegypti forms were found in both clades, indicating no association between mitochondrial variation and color variation (Fig 2).
While our results do not support the tentative patterns suggested by Chan et al. [14], they match those from the most comprehensive mitochondrial phylogeny of the African and global Ae. aegypti generated to date [10]. Using the ND4 variation, Moore et al. [10] showed that Ae. aegypti populations outside Africa represent "mixtures" of mosquitoes from the basal clade and the derived clade, with the basal clade likely originating from West Africa and the derived clade mainly from East Africa. Our analyses of the mitochondrial genome-wide variation revealed the same matrilineage structure in populations from Singapore and northern Queensland (Fig 2). A lack of mitochondrial distinctiveness between the queenslandensis and the type form is also in line with the findings of Moore et al. [10], who could not separate the type and formosus forms into distinct mitochondrial clades despite their assumed subspecies rank.

Nuclear genetic structuring
We extracted nuclear RAD tags that were shared between at least 80% of individuals in the entire dataset (Singapore, Townsville, Gordonvale, Ho Chi Minh City and Rio de Janeiro). To avoid redundant information from the highly linked markers, we randomly selected one SNP per tag with a minor allele frequency greater than 5%, which gave a total of 16,569 markers for downstream analyses.
Parameters of data quality and diversity did not significantly differ between the co-occurring queenslandensis and type individuals, including the average: percentage of reads uniquely aligned to the reference genome (Singapore: t df, 27  Discriminant analysis of principal components (DAPC) showed a clear-cut differentiation of mosquitoes based on their geographic origin and not their color. When the entire dataset was considered, Ae. aegypti individuals formed genetic clusters that corresponded to their sampling region (i.e. Rio de Janeiro, Ho Chi Minh City, Singapore and northern Queensland) (Fig  3a). The only exceptions were three individuals in Singapore (K-M) that formed a distinct genetic group (Fig 3a). They were collected as larvae from the same breeding container, and two were identified as the type form and one as the queenslandensis form (Table 1, Fig 3b). Given their high relatedness (Supplemental file 4) and shared mitochondrial haplotype, as well as high nuclear differentiation from other mosquitoes in the region, it is likely that individuals K, L and M are offspring of an incursion female(s) not local to Australia and Vietnam. These individuals were found near the city port, suggesting a possible route of introduction.
Further analysis of genetic structuring within Singapore revealed that family groups were sampled within the breeding containers, some of which had both color forms (Fig 3b). Highly related queenslandensis and type pairs were found at four locations (Fig 3b), including the incursion family group (K,L,M). Most of the related individuals (24/28 pairs), however, had the same color (Fig 4). These results suggest that the color/scaling pattern is likely to represent a quantitative trait under some environmental influence (e.g. temperature, humidity, light, nutrient availability). The frequency of the color forms has been shown to vary between the dry and Very Pale Ae. aegypti Are Not Genomically Separate the wet season in Ae. aegypti populations from Surabaya, Indonesia [36]. Also, the dorsal abdominal scaling pattern responds to artificial selection [36,37] and multiple QTLs associated with this trait have been recently reported [37]. Individuals with the color/scalling patterns corresponding to the queenslandensis form have also been observed (albeit rarely) in our laboratory colonies which are maintained by occasional crossing to field-caught type Ae. aegypti (Jason Axford, personal communication).
Individuals from northern Queensland were grouped into three clusters corresponding to the three towns where the sampling occurred (Fig 3c). An exception was one queenslandensis individual from Cairns that was grouped with the type individuals from Gordonvale (Fig 3c). The two forms in Townsville could not be distinguished based on their nuclear genome-wide variation (Fig 3c). We found four pairs of closely related individuals: two queenslandensis and two type pairs (Fig 4, S3 File). In other words, all related pairs detected in Townsville were of the same form.
A lower frequency of related individuals in Townsville when compared to Singapore is not surprising given that different sampling methods were employed in these locations. Collection of multiple larvae from the same breeding container increases the chance of sampling family Pairs of Aedes aegypti collected in Singapore (upper graphs) and Townsville (lower graphs). A value of zero for Rousset's genetic distance (â) indicates a distance between a pair of individuals randomly drawn from a given sample, while a negative value indicates lower than average genetic distance between a pair (i.e. their higher relatedness). Color distance between pairs of individuals was treated as a binary value: 0 (same color/form) and 1 (different color/form). doi:10.1371/journal.pntd.0005096.g004 groups, as seen in Singapore and parts of Rio de Janeiro [22]. On the other hand, when BG-sentinel traps are used, the likelihood of related individuals being collected is low. In Townsville, 12.5% of pairs from the same trap were close relatives. Sampling effects are reflected in an elevated level of pairwise genetic distance over geographic distance for mosquitoes from Singapore when compared to Townsville (Fig 4). Such differences in genetic patterns could be erroneously interpreted as differences in the underlying processes (e.g. different dispersal rates), and highlight that sampling methods are crucial when inferring processes within and among Ae. aegypti populations.
In summary, we did not find any association between nuclear genetic variation and color/ scaling variation in Ae. aegypti from Singapore and northern Queensland. Our results are unlikely to be caused by a lack of power to detect genetic structure, given that more than 16,000 genome-wide SNPs allowed us to delineate family groups at a very fine spatial scale. In fact, several families had the queenslandensis and type members. A recent study of global Ae. aegypti populations at 12 microsatellite loci found that at least in one locality in Africa (Senegal) the two established forms (formosus and type) are interbreeding with no sign of genetic subdivision when brought into sympatry [11], so it is not surprising that the type and queenslandensis variety also form one genetic cluster. Our results also help explain the similar vectorial capacity for a dengue 2 viral strain of type and queenslandensis females originating from wild Ae. aegypti in Thailand [38].

Wolbachia infection
In addition to an absence of genetic structuring between the two Ae. aegypti aegypti forms, another line of evidence in support of ongoing interbreeding is the presence of Wolbachia in both forms. We detected this endosymbiotic bacterium in three (out of four) queenslandensis individuals from Cairns and 14 (out of 15) type individuals from Gordonvale, using a lightcycler assay for Wolbachia detection [39]. Wolbachia is not naturally found in Ae. aegypti, but was introduced into the populations in Gordonvale in 2011 and Cairns in 2013 in an effort to reduce dengue transmission [40,41]. This was done by releasing Wolbachia-infected females and males from a colony that originated from type Ae. aegypti [42]. Because the infection is transmitted from mother to offspring, the only way queenslandensis individuals could have become infected by Wolbachia is by mating with infected, type females. Given the high Wolbachia frequency (> 85%) in Cairns and Gordonvale at the time of our sampling [41], the presence of the infection in 3/4 individuals caught in Cairns, and 14/15 individuals caught in Gordonvale was expected.

Conclusion
Our analyses of mitochondrial and nuclear genome-wide variation and the Wolbachia infection indicate that Ae. aegypti queenslandensis and Ae. aegypti type mosquitoes interbreed freely, at least in Singapore and northern Queensland. These findings are of practical importance for control strategies that rely on successful mating between the released and target mosquitoes. Our results also re-enforce the recommendations by the early taxonomic authorities (Mattingly and McClelland) that the extant Ae. aegypti queenslandensis should not be ranked as a subspecies.