article: Within the last century, increases in human movement and globalization of trade have facilitated the establishment of several highly invasive mosquito species in new geographic locations with concurrent major environmental, economic and health consequences. The Asian tiger mosquito, Aedes albopictus, is an extremely invasive and aggressive daytime-biting mosquito that is a major public health threat throughout its expanding range. We used 13 nuclear microsatellite loci (on 911 individuals) and mitochondrial COI sequences to gain a better understanding of the historical and contemporary movements of Ae. albopictus in the Indo-Pacific region and to characterize its population structure. Approximate Bayesian computation (ABC) was employed to test competing historical routes of invasion of Ae. albopictus within the Southeast (SE) Asian/Australasian region. Our ABC results show that Ae. albopictus was most likely introduced to New Guinea via mainland Southeast Asia, before colonizing the Solomon Islands via either Papua New Guinea or SE Asia. The analysis also supported that the recent incursion into northern Australia’s Torres Strait Islands was seeded chiefly from Indonesia. For the first time documented in this invasive species, we provide evidence of a recently colonized population (the Torres Strait Islands) that has undergone rapid temporal changes in its genetic makeup, which could be the result of genetic drift or represent a secondary invasion from an unknown source. There appears to be high spatial genetic structure and high gene flow between some geographically distant populations. The species' genetic structure in the region tends to favour a dispersal pattern driven mostly by human movements. Importantly, this study provides a more widespread sampling distribution of the species’ native range, revealing more spatial population structure than previously shown. Additionally, we present the most probable invasion history of this species in the Australasian region using ABC analysis.
  Many species of mosquitoes are amongst the most invasive pests in the world. They have a long history of human-mediated introductions [1,2] that have resulted in the spread of major epidemics (malaria, dengue, Zika, etc.) and the establishment of invading mosquitoes as a biting nuisance. The Asian tiger mosquito, Aedes albopictus (Skuse 1894) [3], is regarded as one of the most invasive mosquitoes in the world [4]. Native to tropical and subtropical Asia and multiple Western Pacific and Indian Ocean islands, Ae. albopictus now has a pan-global distribution [5–9]. Its initial movement from Southeast (SE) Asia toward the Indo-Malayan Peninsula and Indian Ocean islands may have resulted from the increase in human migration during the 17th and 18th centuries, with international trade (particularly the used-tire and ornamental bamboo trades) further facilitating its global spread in the 20th century [10]. It is among the primary vectors of several globally expanding and medically important arthropod borne viruses (arboviruses) –particularly dengue, chikungunya, yellow fever, and Zika—while also able to transmit at least 23 other arboviruses and canine heartworm [10]. Whilst at present, the yellow fever mosquito Aedes aegypti (Linnaeus 1762) [11] is responsible for most of the transmission of some of these important arboviruses, the increased cold tolerance of Ae. albopictus relative to Ae. aegypti suggests that it could extend the range of many of these diseases under the right circumstances [12,13]. Genetic techniques can provide critical information to infer the routes and sources of invasive species as well as informing on the demographic history and genetic composition of founding populations [8,14]. For mosquito vectors, this knowledge is not only useful for inferring invasion routes in order to focus biosecurity efforts, it can also inform us of the colonizing capacity, adaptability and behaviour of invading mosquito lineages [15,16]. However, the high dispersibility of Ae. albopictus mediated by human activities can make it challenging to detect genetic variation between populations due to high gene flow facilitated by these activities [17]. Microsatellite markers have been shown to be useful for exploring Ae. albopictus genetic patterns as they evolve rapidly and can often detect subtle population structure [8,18,19]. Recently, approximate Bayesian computation (ABC) analysis has proven a powerful tool in testing the probability of competing invasion scenarios. This can provide us with crucial information regarding the timing and origin of mosquito introductions, which has been used recently for Aedes mosquitoes [20–22]. Overall, the population genetics of Ae. albopictus through the SE Asian-Indo-Pacific region requires further exploration and samples from this region (particularly Australasia) are often lacking from global population genetic studies, despite the importance of this region for vector research [23–25]. A study by Beebe, Ambrose [18] explored part of the Indo-Australasian invasion by Ae. albopictus and provided an interesting scenario where there is high human connectivity (largely maritime) spanning both oceanic barriers and complex geographic landscapes. The current study expands on this work to include SE Asian native populations (Myanmar, Thailand, Malaysia, Singapore, Indonesia), as well as several younger populations that appear to have been introduced within the last six decades (Papua, Papua New Guinea (PNG), Solomon Islands, Fiji, Christmas Is. , Cocos (Keeling) Is. , Nauru, Torres Strait Islands (Australia) ). Additionally, we included populations from the United States (USA) and northern Asia (only for COI) to see how these populations fit into a broader geographic analysis. We use previously developed nuclear microsatellite markers [18] and mitochondrial cytochrome c oxidase subunit I (COI) sequences to investigate the population genetics of the Ae. albopictus within the Indo-Pacific region. Our primary aims were to uncover the most likely historical invasion route of Ae. albopictus into the Australasian region as well as to detail the genetic connectivity and population structure of Ae. albopictus throughout this broad geographic region that we refer to as the Indo-Pacific (the aforementioned populations). While there are some records of the progressive establishment of Ae. albopictus throughout this region, the origin/s and the precise timing of introductions require testing using genetic methods under a coalescent-based approach such as ABC analysis. Our secondary aim was to further investigate the 2005 colonization of the Torres Strait Islands, Australia [26]. Many of the Torres Strait Islands have undergone intense spraying efforts since the establishment of Ae. albopictus and the region experiences monsoon-dry seasons leading to regular population bottlenecks [27,28]. We hypothesized that the genetic changes in neutral alleles may be detectable over time in these newly invaded and small island populations.
  Mosquito collections involving HLC from the Solomon Islands (S1 Table) were approved by the Medical Research Ethics Committee in compliance with Australia’s National Statement on Ethical Conduct in Human Research (project no. 2011000603). Collectors involved in HLC took anti-malarial medication and wore long-sleeved, protective clothing. Both adult and larval samples were collected throughout Australasia, SE Asia, Indian Ocean and Pacific Ocean islands as well as in the United States (Table 1, Fig 1 (orange dots) ). Samples were stored in 70% ethanol or dried (adults) over silica beads. Samples were collected using human landing captures (HLC), human baited sweep netting, egg collections and sampling of aquatic habitats for larvae and pupae (S1 Table). For identification purposes some samples were reared to adults after field collection (S1 Table). The logistics of sourcing material across multiple international borders resulted in variability in collection methods and sample sizes (Table 1, S1 Table). Adult mosquitoes were identified morphologically [29] and for larval/pupal samples using either real-time PCR assays [30] or a PCR-restriction digest [31] to differentiate from Aedes scutellaris (Walker 1858) [32]. DNA was salt extracted [33] and diluted at 1: 10 with 1X TE buffer (Tris, EDTA). Thirteen nuclear microsatellite were used in this study. These markers were previously developed [18] and include two dinucleotide and 11 trinucleotide loci (see Beebe, Ambrose [18] for loci and primers). Some samples included from the previous study were amplified using a variation of the master mix (see Beebe, Ambrose [18]), other samples were amplified in a 15. 4μl reaction that consisted of 10. 8μl H2O, 3μl 5X Mytaq buffer (Bioline, containing 5mM dNTPs and 15mM MgCl2), 0. 1μl 10μM M13 tagged forward primer, 0. 2μl 10μM reverse primer, 0. 2μl M13 tagged fluorescent dye (VIC, NED, PET or FAM; S1 Table), 0. 01μl (1U) MyTaq polymerase and 1μl 1: 10 DNA template. PCR cycling used the same protocol as in Beebe, Ambrose [18]. Amplification was verified by running 1μl of PCR product on a 2% agarose gel stained with either GelRed (Biotium) or MidoriGreen (Bulldog Bio). Samples that amplified successfully were sent to Macrogen Inc. (Republic of Korea) for genotyping. GeneMarker v. 2. 4. 2 (SoftGenetics LLC [34]) was used to score alleles for each locus manually after passing the data through the standardization run wizard using the default animal fragment setting. Random selections of genotyped plates were rescored by a second person to assess consistency in scoring. In addition to the data collected in this study, we included microsatellite scores from samples in Beebe, Ambrose [18]. During the present study, it was uncovered that the Beebe, Ambrose [18] study used (in some cases) inconsistent fluorescent dyes for a given locus, which caused a dye-shift [35] resulting in inconsistently scored alleles. We regenotyped a random subset of individuals from each of the populations used in the study by Beebe, Ambrose [18] to ensure consistency with data collected from this study. The predictability of this dye-shift (based on the dyes used previously) enabled shifting of the allele scores from Beebe, Ambrose [18] for use in this study. Samples with fewer than nine scored loci of 13 total were removed before further study as we considered these poor quality samples; thus leaving 911 samples for final analyses (20% of samples were from Beebe, Ambrose [18]; Table 1). With the remaining dataset, missing values were replaced based on population allele frequencies using GenoDive v. 2. 0b27 [36]–based on preliminary analyses this did not significantly alter population structure and relationships between populations. Missing values were not replaced for the STRUCTURE analysis, calculation of HWE and for checking the presence of null alleles. Scored allele frequencies were checked for the presence of null alleles using MICRO-CHECKER [37] and for Hardy-Weinberg equilibrium using GenAlEx v. 6. 5 [38,39]. Additionally, we calculated fixation index (F), allelic richness (Na), number of effective alleles (Ne) and the observed (Ho) and expected (He; unbiased estimate: uHe) values of heterozygosity using GenAlEx v. 6. 5. Pairwise population indices of genetic variation for Jost’s D, G”ST and FST were also calculated between populations in GenAlEx v. 6. 5 (S2 Table). We used 9,999 permutations and an analysis of molecular variance (AMOVA) to assess significance. A Mantel test was also performed in GenAlEx v. 6. 5 on geographic and genetic distance (pairwise phiPT) using 9,999 permutation [38,39]. Population structure was investigated using the Bayesian program STRUCTURE v. 2. 3. 4 [40] to infer the most probable number of population clusters (K). Based on our preliminary runs (S1 File, Supplementary Methods: Preliminary STRUCTURE), final analyses were run at both a lower (K = 4) and upper (K = 9) K value. For both K values we used a burn-in of 100,000 and runtime of 2,000,000 generations per iteration (20 iterations). For K = 9, cluster membership probabilities were somewhat inconsistent across runs due to multimodality; 20 iterations helped to account for this [41]. We assessed whether the burn-in period was adequate by reviewing summary statistics in STRUCTURE [41]. CLUMPP v. 1. 1. 2 [42] was used to compile data from the 20 iterations for the independent values of K using the Greedy algorithm with 1,000 replicates. Final graphs were formed in DISTRUCT v. 1. 1 [43]. Discriminant analysis of principal components (DAPC) and correspondence analysis (CA) was used to further assess population structure. DAPC was implemented in R Studio v. 3. 2. 2 (RStudio Team 2015 [44]) using the adegenet 1. 4–1 package [45,46] using the whole microsatellite dataset, where group membership was defined by the populations outlined in Table 1 (see DAPC abbreviation). These populations were more broadly defined and differed slightly from those used in STRUCTURE, to allow for easier interpretation of the results presented here. Specifically, the Torres Strait populations were split into groups based on their genetic relationship to one another and geographic/temporal information to reduce clutter in plots (Table 1; DAPC abbreviation). Only populations that were genetically similar were grouped together, which was confirmed using STRUCTURE, DAPC, CA and pairwise tests for genetic distance (FST, G”ST, Jost’s D) on the full dataset and subsets. Final DAPC analyses were performed on both a full dataset (nind = 911, npop = 23; including all populations) and a reduced dataset (nind = 458, npop = 18; excluding Jakarta, Sumba, Timor-Leste, the Torres Strait Islands and Southern Fly Region) to help discriminate genetically similar populations. In adegenet, cross-validation was performed on each of our DAPCs independently, using a training dataset of 90% and a validation set of 10%, using 100 replicates. The number of PCs (n. pca) associated with the lowest root mean squared error (RMSE) was used as this was considered optimum [47]. Cross-validation suggested retaining 60 PCs for the full dataset and 40 PCs (n. pc) for the reduced dataset. We used five discriminant functions (n. da) for each of the analyses, but only the first three are plotted and discussed here as they explained the majority of variance (see Results). A correspondence analysis (CA) was also implemented in adegenet on the full dataset to investigate general trends and to complement DAPC, as visualization of the data is simplified in CA because within-population genetic diversity is not displayed. The Garza-Williamson index (M-ratio) was calculated for populations using our microsatellite dataset in Arlequin v. 3. 5. 2. 2 [48]. The M-ratio was used to investigate the demographic history of populations and to test for recent bottleneck events; wherein an index statistic closer to 1 suggests the population is in a stationary state whereas very low values suggests a population has gone through a genetic bottleneck in the past (with a critical value of 0. 68 indicating a bottleneck) [49,50]. A Wilcoxon test for heterozygosity excess was also conducted on populations to detect bottlenecks using BOTTLENECK v. 1. 2. 02 [51]. We used a two-phase model (TPM) of mutation with 10% infinite allele model and a 90% single step mutation model with 15% variance for 1000 iterations. A Wilcoxon signed rank test (two-tailed) was used to calculate significance (P < 0. 05). We tested the invasion history of Ae. albopictus in part of the study region (SE Asia/Australasia) —populations from the Indian Ocean, Fiji, Nauru and United States were not included due to insufficient sampling of these regions. Both COI sequences and the thirteen microsatellites loci were analysed together using ABC analysis in DIYABC v. 2. 1. 0 [52]. Due to the complexity of modelling each population separately in this region, we simplified our invasion scenarios by randomly subsampling individuals from distinct geographic and genetic groups (defined using our other analyses). These representative populations included: mainland SE Asia, Indonesia, Papua, PNG, the Solomon Islands and the Torres Strait Islands/Southern Fly Region. For both the PNG and the Torres Strait/Southern Fly Region populations, we included temporal sampling in our scenarios (asterisks, Fig 2). Each of these representative populations/sampling events was made up of 30 individuals except Papua which used all 20 samples from the only sampled population, Timika. In addition, we included an unsampled ancestral population (ANC, Fig 2) in our model that split into mainland SE Asian and Indonesian populations. Preliminary runs were carried out in accordance with Bertorelle, Benazzo [53] in order to optimise summary statistics, prior estimates and the scenarios tested. Final runs compared five invasion scenarios; the final summary statistics, prior estimates and parameter conditions used are outlined (S3 Table). Each scenario represents a plausible invasion route into the Australasian region. These were constructed using historical records regarding the timing and suspected sources of the different invasions [26,54–58]. Priors were sampled from a wide range of distributions based on these records—less certain time priors were given a wider distribution and standard deviation whereas more likely priors were assigned narrower estimates. The upper time bound for the divergence of mainland SE Asia and Indonesia from a common ancestor was based on Porretta, Mastrantonio [59] while the lower bounds allowed for the possibility of a more recent split associated with human migration [10] (S3 Table). Because DIYABC measures time in terms of the number of generations, we assumed 10 generations per year for Ae. albopictus (which typically ranges from 5–17 generations in the tropics). Estimates of effective population size (Ne) ranged from 10 to 1,000,000 individuals (uniform distribution) [20–22] depending on if a population was modelled as going through a change in Ne. For instance, each of the recently introduced populations was modelled to allow for a founder effect after its introduction using lower Ne ranges (Fig 2, S3 Table). We additionally allowed for a change in Ne in the Torres Strait Islands/Fly Region population due to the drastic temporal changes we observed in our other analyses—this would allow us to make a relative comparison of Ne in order to see if the population had undergone any change in Ne (indicative of bottleneck/expansion events) (Fig 2, S3 Table). For COI, we used the HKY mutation model [60] and sampled from a uniform distribution with mutation rates ranging between 7x10-10–1x10-7. For microsatellites, both di- and tri- nucleotide repeats were modelled separately due to the possibility of different repeat lengths having different mutation rates [61]. Both microsatellites used the default generalized stepwise mutation model and were assigned a loguniform distribution with mean mutation rates between 1x10-6–1x10-3 (S3 Table). All mutation rates were based on standard ranges for COI [62] and for Dipteran microsatellites [20–22,63]. We simulated 15,000,000 datasets and each of the five scenarios was given a uniform probability. The performance of the ABC approach was assessed using multiple methods in DIYABC (S1 File, Supplementary Methods: Performance of DIYABC). The mitochondrial protein-coding gene COI was amplified using custom designed primers (albCOIF 5’-TTTCAACAAATCATAAAGATATTGG-3’ and albCOIR 5’- TAAACTTCTGGATGACCAAAAAATCA-3’) for 259 random individuals across different populations. Each 25. 3μL reaction consisted of 19μl H2O, 5μl 5X Mytaq buffer (Bioline, with pre-optimized concentrations of dNTPs and MgCl2), 0. 1μl 100μM forward primer, 0. 1μl 100μM reverse primer, 0. 1μl MyTaq polymerase and 1μl 1: 10 DNA template. PCR used an initial denature of 94°C for 3 min, 35 cycles of denaturation at 95°C for 30 sec, primer annealing at 45°C for 40 sec, and primer extension at 72°C for 30 sec. Final elongation lasted 5 min at 72°C prior to cooling to 4°C. Amplification was confirmed using gel electrophoresis (as described previously) and PCR products were purified by adding 2μl per sample of a mixture containing equal amounts of Exonuclease I and Antarctic Phosphatase (New England Biolabs, Australia) before incubation at 37°C for 20 min and denaturation at 80°C for 10 min. Samples were sequenced by Macrogen Inc. (Republic of Korea) using Sanger sequencing. Additional COI sequences of Ae. albopictus were obtained from other studies and from Genbank (1044 sequences total, 259 produced in this study, S4 Table). Sequences were edited and aligned in Geneious v. 9. 0. 4 (http: //www. geneious. com, Kearse, Moir [64]) using the MAFFT alignment. The final alignment was trimmed to 445bp to incorporate the large number of COI sequences available from Genbank which were smaller than the ~700bp region sequenced in this study. All sequences were checked for stop codons in Geneious v. 9. 0. 4. TCS haplotype networks [65] were constructed using 1,000 iterations in PopArt v. 1. 7 (http: //popart. otago. ac. nz). In addition, we calculated Tajima’s D for populations with temporal data in PopArt v. 1. 7 to determine whether sequences were evolving randomly or non-randomly. Haplotype and nucleotide diversity was calculated using DnaSP v. 5. 10. 1 [66]. All new COI sequences generated in this study are available on Genbank: KY907195—KY907453.
  Allelic richness for microsatellites was highest in native populations of Ae. albopictus from Myanmar, Thailand, Malaysia, but also high in recently invaded areas such as some of the Torres Strait Islands and in PNG and the Solomon Islands (See Na in S5 Table). A Mantel test on the whole dataset showed a significant (P = 0. 0001) positive, but weak, correlation (R2 = 0. 02) between genetic (phiPT) and geographic distance (y = 0. 0002x + 28. 1). Pairwise estimates of FST, Jost’s D and G”ST all revealed similar results to each other and recovered mostly significant relationships between populations; here we discuss gene flow and genetic distance in regards to FST estimates but the other measures are shown in S2A–S2C Table. Lowest FST values were apparent between populations belonging to the same geographical region (for definition of regions see Table 1, Region/description), especially within mainland SE Asia (FST = 0. 011–0. 103) (S2A Table). However, some comparisons between regions separated by vast geographical distances also showed low FST scores, such as populations from the Solomon Islands with populations from mainland SE Asia (FST = 0. 050–0. 114) (S2A Table). The relationships between populations were mostly consistent with the results obtained in STRUCTURE and multivariate analyses, which are described in detail below. Within the study region, four to nine clusters were supported by the Evanno ΔK and log likelihood methods for inferring K. While K = 4 (Fig 1) represents the simplest summary of the genetic structure of Ae. albopictus in the region, we detected substantial substructure within these four main clusters which are apparent at K = 9 (S1 Fig). We discuss the data in the context of both values of K to avoid underestimating the degree of population structure within the study region. At K = 4, clusters mostly pertained to distinct but broad geographic boundaries although many populations and individuals show signs of admixture, despite the large geographic distances (Fig 1). The mainland SE Asian populations of Myanmar, Thailand, Malaysia and Singapore cluster with the USA (Hawaii and Atlanta; Fig 1B and 1C), La Réunion, as well as Fiji and Nauru (light purple; Fig 1). The second cluster (light green) contains populations from Indonesia (Jakarta and Sumba), Timor-Leste, the Southern Fly Region of PNG (Fig 1A) and several islands of the Torres Strait (especially collections following the first detection of Ae. albopictus in the straits in 2005 (collections between 2006–2014) ) (Fig 1A). An additional cluster (purple) is prominent within the Torres Strait region (Fig 1A) and represents populations on the islands collected more recently (2013–2015), suggesting temporal shifts in population structure have occurred on some islands (see Fig 1A: Ker, War). The fourth cluster (green) is composed of historically-established PNG populations, but note that some of these populations contain admixture with the SE Asian cluster (Fig 1). The island of Daru (Fig 1A: Dar), which is less than 5km from the Southern Fly Region incursion populations (light green cluster; Fig 1A; Sig, Kul, Mbd, Kat), is distinct and clusters with the historically established PNG populations. Timika and the Solomon Islands show genetic affinity to both PNG and SE Asian clusters (Fig 1). Indian Ocean islands (Fig 1: CK and CH) appear differentiated from each other and contain a notable degree of admixture, but Christmas Is. is more similar to SE Asia whereas Cocos (Keeling) Islands appear as an admixed population made up of the PNG and Indonesian clusters. When K = 9, the same broad population patterns are observed but some populations/regions become more distinct, including the USA and Hawaii, Solomon Islands, Timika, the Cocos (Keeling) Islands and Sumba/Torres Strait Island populations (S1 Fig). Relationships for this K value are described in S1 File, Supplementary Results, STRUCTURE (K = 9). For multivariate analyses, the DAPC on the full dataset (n. pc = 60, n. da = 5) explained 89% of variance, whereas the reduced dataset (n. pc = 40, n. da = 5) explained 79. 3% of the variance in the data. Eigenvalues for these first three PCs are 206. 37,110. 78 and 63. 21 for the full dataset (S2A Fig) and 99. 6,69. 2 and 44. 99 for the reduced dataset (S3 Fig); these values correspond to the ratio of between-group over within-group variance for each discriminant function. For the correspondence analysis (CA) we plotted the first three eigenvalues (0. 19,0. 15,0. 09), which indicate the proportion of variance explained by the first three PCs (S2B Fig). DAPC and CA results of the full dataset (S2 Fig) showed similar population differentiation as observed in STRUCTURE at K = 4. Due to the large number of populations, we describe population structure based on the broad clustering—populations are color coded with geographically close populations being more similar in color. Four major clusters of populations are noticeable when the first three PCs are plotted against each other (C1-C4, S2A Fig). However, there is considerable overlap between these clusters, particularly with C4 overlapping C2 and C3, suggesting that these individuals and populations are genetically similar and show signs of admixture (S2A Fig). Cluster 1 (C1) represents recent (2012–2015) collections from Torres Strait Islands and is the most distinct from the other clusters. It is most closely related to C2, which contains earlier collections (2007–2014) from the Torres Strait Islands and populations from the Southern Fly Region, Sumba, Timor-Leste and Jakarta (S2A Fig). The relationships uncovered in the DAPC of the full dataset were recovered in the CA and are more easily visualized, where the first three principal components (PCs) of the CA are plotted in three-dimensions (S2B Fig); however, note that the large amount of within-population variation (as displayed in the DAPC plots) is not shown. Because of the overlap of C3 and C4 clusters, we separately analysed these clusters by DAPC (referred to in the Materials & Methods as the reduced dataset) that excluded populations from the Torres Strait Islands, Jakarta, Sumba, Timor-Leste and the Southern Fly Region (S3A and S3B Fig) to explore substructure within these clusters (i. e. C3 & C4 in S2A Fig). Populations from PNG (Kiunga, Madang, Port Moresby and Daru) were similar to each other but distinct from the other populations (S3 Fig). The offshore PNG populations (Lihir Is. and Buka Is.) were somewhat differentiated from mainland PNG populations, although Buka Is. shares some overlap with both Port Moresby and Lihir Is. (S3 Fig). The Solomon Islands also appears similar to Lihir Is. and Buka Is. populations, but is somewhat distinct (S3 Fig). Mainland SE Asian populations appear genetically similar and tend to exhibit the most genetic overlap with other populations (S3 Fig). Nauru and Fiji are most similar to mainland SE Asian populations. In contrast, both USA populations (Hawaii and Atlanta) as well as La Réunion appear well differentiated from mainland SE Asian populations. Cocos (Keeling) Island and Timika share some overlap with each other, whereas Christmas Island shares overlap with both PNG and mainland SE Asia populations (S3 Fig). The scenario with the highest posterior probability using the logistic approach was scenario 4 (P = 0. 52 [95% CI: 0. 44,0. 59], Fig 2). In this scenario, Ae. albopictus colonized Papua and PNG in two separate events from mainland SE Asia, established in the Solomon Islands via PNG and more recently colonized the Torres Strait Islands/Southern Fly Region via Indonesia (Fig 2). The timing of each introduction event is shown in S3 Table and corresponds with historical records for the introduction of Ae. albopictus in the tested populations, although 95% confidence intervals (CI) suggest that some introduction dates could have been earlier than first observed (see Discussion). None of the other scenarios showed overlapping 95% CI with Scenario 4 (Fig 2), however, we detected moderate levels of type I (0. 45; probability that scenario 4 is rejected given that it is the ‘true’ scenario) and type II error (0. 48; probability of deciding scenario 4 is the ‘true’ scenario when it is not) that suggest Scenario 1 (P = 0. 34 [95% CI: 0. 31,0. 37]) could provide a plausible alternative invasion scenario for our data (S3 Table). Scenario 1 is identical to scenario 4, except that the Solomon Islands is modelled as originating from mainland SE Asia, rather than from PNG (Fig 2). Consequently, we discuss the Solomon Islands invasion based on both alternative origins and posterior estimates calculated under the combined scenarios are presented (in Discussion), although show no significant difference from scenario 4 alone (S3 Table). A preliminary analysis showed low support (P = 0. 001 [95% CI: 0. 00,0. 12]) for a scenario where the Solomon Islands introduction was modelled as an admixture event between mainland SE Asia and PNG compared to the five scenarios compared in our final analyses (but using slightly different Ne prior ranges for all founders (10–10,000) ). Each introduced population showed no relative change in Ne due to large 95% CIs of posterior distributions, although median Ne values were smaller for founding events (S3 Table). Likewise, the duration of the modelled bottleneck had large 95% CIs, but median values generally ranged from 16–30 generations (S3 Table). Overall, the Torres Strait Islands/Southern Fly Region population showed a stable Ne since its introduction (due to overlapping 95% CIs), although there was a gradual increase in median Ne over time, potentially suggesting growth of the population as a whole. (S3 Table). Our assessment of the performance of our ABC analysis was supported as fitting our observed data well (S4 and S5 Figs). A total of 52 COI haplotypes were identified from the 1044 individuals used for generating the TCS haplotype network, with 92% of individuals belonging to nine main haplotypes (H1-5,11,15,39,43) (Fig 3, Table 2, S6 Fig). The distribution of these haplotypes by specific population is shown in S6 Table. All new sequences generated from this study are available on Genbank (Accession no. KY907195—KY907453; see S4 Table for accession numbers of sequences from other studies). The COI haplotype network is less informative in regards to population structure than the microsatellite data, although it does highlight broader geographic relationships that are somewhat consistent with the microsatellite results. Of the nine main haplotypes, H1 has the most individuals, mostly from eastern Asia (China, Taiwan and Japan), USA (mainland USA and Hawaii), Madagascar and La Réunion (Fig 3, Table 2). Similarly, H39 is distributed in a similar temperate/subtropical region. Haplotype 3 consists primarily of individuals from the Torres Strait Islands, Fly Region, Indonesia, Timor-Leste, PNG and the Philippines. However, the majority of PNG sequences belonged to H5, which also includes individuals from the Solomon Islands, Indian Ocean islands (CK and CH), Singapore and Thailand. Another major haplotype, H4, includes the most diverse range of populations (in terms of geographic spread), although it mostly consists of mainland SE Asian populations and populations from the tropics. Of the additional COI haplotypes, many are exclusive or shared amongst close geographic regions (Fig 3, Table 2, S6 Table, S6 Fig), although others show no apparent geographic pattern. Haplotype diversity (Hd) for the total dataset was high (Hd = 0. 83) as was nucleotide diversity (π = 0. 0037); however, population measures of Hd and π varied considerably (Table 3). Neutrality tests (Tajima’s D) on all populations with temporal data were not significant (S7 Table). Past genetic bottleneck events were not consistently indicated using both Wilcoxon tests for heterozygosity excess and M-ratio, with the exception of the 2010 population from Waiben of the Torres Strait (M = 0. 68, P = 0. 032 (two-tailed Wilcoxon signed rank test for heterozygosity excess) (S8 Table). However, multiple populations showed signs of a bottleneck using a single method. The M-ratio indicated a bottleneck for some populations of the Torres Strait Islands (Mabuiag, Waiben, Ngurupai, Poruma, Iama), Port Moresby, Timika and Gizo, whereas the Wilcoxon test was significant for Waiben, Madang, Jakarta, Yangon and the Cocos (Keeling) Islands (S8 Table).
  Multiple population studies have attempted to capture the amount of genetic structure throughout the species’ native range. Allozyme studies have shown that Indonesian and Japanese populations of Ae. albopictus are likely distinct [67] and that SE Asia (Borneo, peninsula Malaysia) and southern Asian populations (India, Sri Lanka) can both be differentiated from northern Asian populations (China, Japan) [68]. While no studies have conducted a comprehensive analysis of the species’ full native range [8], the genetic differentiation of native Asian populations of Ae. albopictus may confer to both a north-south (Korea to Indonesia) and east-west (Japan to India) pattern of genetic differentiation. Our results partly support this pattern, with evidence for genetic differentiation separating northern Asia (COI data only, no microsatellite data available), SE Asia and Indonesia (both COI and microsatellite data). Within mainland SE Asia, our data revealed little to no population structure despite high genetic diversity and COI haplotype diversity, supporting the findings of other studies [59,69,70]. Using climatic modelling and two mitochondrial markers, Porretta, Mastrantonio [59] suggested that the low genetic structure across this mainland SE Asian region could be explained by the demographic growth between interconnected populations of Ae. albopictus preceding the last glacial maximum (LGM, occurring ~21,000 ybp), with the species’ ecological flexibility facilitating its success in the ecologically diverse Sundaland (exposed SE Asian landmass) during this period. Whilst their study lacked sampling from Indonesia, their data suggested climatically suitable habitat for Ae. albopictus existed across the southern range of Sundaland which later formed the Indonesian islands after a rise in sea levels. They hypothesized that the emergence of Sundaland during the LGM could have facilitated population connectivity across Indonesia and mainland Asia. In contrast, we found clear genetic differentiation between the Indonesian archipelago and mainland SE Asia with both our microsatellite and COI data, which could be explained by fragmentation and subsequent differentiation of Ae. albopictus populations following a rise in sea levels or potentially driven by human migration from the region [10]. Our ABC analysis supports a more historical split ~9,040 ybp (95% CI = 1,790–27,600 ybp) but further investigation with extensive native population sampling would be needed to clarify dating. The genetic homogeneity and high gene flow (FST = 0. 020–0. 103) in our microsatellite data across mainland SE Asia (including Malaysia, Singapore, Myanmar and Thailand) could also be explained by human-mediated gene flow associated with transportation infrastructure (land, air and sea); aircraft and road networks in particular may be major drivers in the connectivity of Ae. albopictus in this region given the inland location of many of these populations [71–73]. For our full dataset, the relationship between genetic and geographic distance was significant, however, the correlation was weak (R2 = 0. 02), highlighting the potential extant that human movements have had on Ae. albopictus population structure; although it does also highlight a minor trend of isolation by distance in our study region. A recent worldwide study that examined the mitogenome diversity of Ae. albopictus found three major haplogroups, two of which were implicated in the global spread of the species [9]. Some of these haplogroups appeared more prevalent in particular climatic and geographic regions, such as haplogroup A1a which chiefly characterized the tropics and A1a2 which is mostly distributed in temperate regions. Another haplogroup (A2) appeared important in the spread of Ae. albopictus from SE Asia toward Australasia, distinguishing many of the samples from the Philippines, PNG, Indonesia and the Torres Strait Islands [9]. A similar pattern was observed using our COI data when viewing the nine haplotypes that account for 92% of the sequences in our study. Haplotype three (H3) was found to be common in Philippines, Indonesia, the Torres Strait Islands and Southern Fly Region but also present in historically established PNG populations, Vietnam and an individual from Sepilok, Borneo (Malaysia). Consequently, the populations from Vietnam, Borneo, the Philippines and Sulawesi could represent important unsampled populations that may influence ABC results and should be considered in future studies, which could additionally explain why the posterior probability of our most likely invasion scenario was moderate (P = 0. 52). We also found that H1 was more prevalent in regions that experience temperate and subtropical climates, while H4 was widespread in the tropics (Fig 3). Due to the maternal inheritance of mtDNA, these results could suggest that the movement of females has been somewhat limited to their preadaptation to certain climatic regions. For example, females originating from a temperate region may be more being likely to successfully invade other temperate regions. It is possible that these genetic patterns for COI in Ae. albopictus are associated with the photoperiod response of different populations, which could contribute to higher gene flow and invasion success between climatically similar regions [17,74–76]. However, the photoperiodic response of Ae. albopictus in recently introduced regions within Australasia is lacking (especially for New Guinea, the Solomon Islands, Fiji and Nauru). Within the Torres Strait Islands it appears the population is of tropical, Indonesian origin and egg survival was lower in less humid conditions [77], supporting this conclusion. However, it is worth noting that there are multiple other COI haplotypes that do not conform to any obvious geographic patterns (S6 Fig), suggesting that there have been multiple introduction events into some locations from mainland Asia or potentially unsampled locations—a similar pattern which was also highlighted by other studies [9,74]. The possibility of insertion of mtDNA in the nuclear genome [78,79] was considered in this study and our COI sequences were assessed by examining chromatograms (no double peaks in chromatograms) and by checking for stop codons. It seems unlikely that there is nuclear insertion of mtDNA in samples from our study—but it remains a possibility and requires further research. We included several populations from the USA and Indian Ocean to explore how these populations fit into a broader geographic analysis with our samples, which are chiefly from SE Asia and Australasia. We did not include these in our ABC analysis because of the lack of temperate Asian populations, which have been shown as the source of USA introductions and because of insufficient sample sizes in Indian Ocean populations (Cocos (Keeling) Islands, Christmas Island and La Réunion). We included them in our other analyses to assist in future studies and have discussed them in S1 File, Supplementary Discussion. For the first time, we have used ABC analysis to compare various invasion scenarios of Ae. albopictus in the Australasian region and have additionally explored and characterised the genetic structure of a wide range of populations in the Indo-Pacific, which had not been compared under a single study. We uncovered notable temporal population structure in recently introduced Torres Strait Island populations. Importantly, this demonstrates some of the drastic changes that invading populations may undergo within a short time period (i. e. in less than a decade), which has substantial implications on the practicality and accuracy of using such genetic databases for estimating invasion sources. However, we also found that historically established populations of Ae. albopictus displayed stable population structure. Future studies that aim to address the global genetic structure of Ae. albopictus will need to consider the full native range of the species and the influence of temporal collections on population structure, especially newly established populations. Likewise, ABC analyses that account for complex scenarios (requiring thorough spatio-temporal sampling) and gene flow between populations will play a key role in better understanding the population dynamics of Ae. albopictus as well as other mosquito species that are highly associated with humans, such as Ae. aegypti. The standardisation of genotyping methods and sampling efforts will allow for more rigorous assessments of the global population structure of Ae. albopictus, given the scale of its current distribution which makes such population studies logistically challenging. This will prove essential in controlling the spread of Ae. albopictus and for assessing the health risk of different populations, given their variation in vector efficiency, physiology and behaviour [10,86–88].
summary: The Asian tiger mosquito, Aedes albopictus, is an aggressive mosquito that has expanded globally in the last century, chiefly due to the increase of human movements. It is a major public health concern due to its role in transmitting dengue, chikungunya and Zika viruses. Understanding how populations of Ae. albopictus are genetically related and how they have been introduced into new regions is important for controlling them and assessing their disease risk; few studies have explored this in the Indo-Pacific. In our study, we sampled a broader range of populations of Ae. albopictus in the Indo-Pacific to explore genetic patterns and to investigate the likely route of invasion through Australasia. We uncovered clear genetic groups throughout this region, but also found that some geographically distant populations are closely related, likely due to human-associated movements. We also found, that Ae. albopictus likely colonized New Guinea from mainland Southeast (SE) Asia before spreading to the Solomon Islands via either PNG or SE Asia. In contrast, the populations in Australia’s Torres Strait Islands were introduced from Indonesia. Interestingly, we found major genetic changes over time in some Torres Strait populations, less than a decade after its introduction, potentially reflective of a random reduction in genetic diversity (genetic drift) or a secondary invasion.
keywords: invertebrates
  solomon islands
  medicine and health sciences
  population genetics
  geographical locations
  animals
  indonesia
  genetic mapping
  population biology
  insect vectors
  fiji
  infectious diseases
  papua new guinea
  disease vectors
  insects
  arthropoda
  people and places
  mosquitoes
  haplotypes
  asia
  heredity
  oceania
  genetics
  biology and life sciences
  species interactions
  evolutionary biology
  organisms
title: Tiger on the prowl: Invasion history and spatio-temporal genetic structure of the Asian tiger mosquito Aedes albopictus (Skuse 1894) in the Indo-Pacific