Tiger on the prowl: Invasion history and spatio-temporal genetic structure of the Asian tiger mosquito Aedes albopictus (Skuse 1894) in the Indo-Pacific

Background 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. Methodology/Principal findings 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. Conclusions/Significance 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.


Methodology/Principal findings
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 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111

Introduction
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][6][7][8][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 17 th and 18 th centuries, with international trade (particularly the used-tire and ornamental bamboo trades) further facilitating its global spread in the 20 th 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][21][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][24][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.

Ethics statement
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.

Mosquito samples
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]

Microsatellite analysis
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 H 2 O, 3μl 5X Mytaq buffer (Bioline, containing 5mM dNTPs and 15mM MgCl 2 ), 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 Table 1. Sample information for Aedes albopictus used in the microsatellite study, where n indicates the number of individuals per population (n total = 911, n pop = 50). Region/description shows broader geographic regions and descriptions referred to in text. Population indicates more specific collection sites and the year of collection in brackets.

Region/description
Population ( 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 F ST 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 [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 [ 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 (F ST , G" ST , Jost's D) on the full dataset and subsets. Final DAPC analyses were performed on both a full dataset (n ind = 911, n pop = 23; including all populations) and a reduced dataset (n ind = 458, n pop = 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

Invasion history-ABC
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][55][56][57][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 (N e ) ranged from 10 to 1,000,000 individuals (uniform distribution) [20][21][22] depending on if a population was modelled as going through a change in N e . For instance, each of the recently introduced populations was modelled to allow for a founder effect after its introduction using lower N e ranges (Fig 2, S3 Table). We additionally allowed for a change in N e 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 N e in order to see if the population had undergone any change in N e (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][21][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).

Fig 2. Invasion scenarios of Aedes albopictus in Australasia tested using approximate Bayesian computation (ABC).
One unsampled and six sampled populations were modelled, shown as colored lines in five different invasion scenarios. Time events (t1-t7) are not to scale, but their prior distributions are displayed as the year (except t7 which is shown in years before present (ybp, present = 2015)). Changes in effective population size (N e ) are represented as differently shaded lines, where db-db4 represent the duration; narrowing lines represent population bottlenecks that were given lower N e priors ranges; rate of admixture (ra) is also shown. All populations have samples at time = 0 (i.e. 2015) and asterisks represent additional temporal

COI analysis
The mitochondrial protein-coding gene COI was amplified using custom designed primers (albCOIF 5'-TTTCAACAAATCATAAAGATATTGG-3' and albCOIR 5'-TAAACTTCTGGA TGACCAAAAAATCA-3') for 259 random individuals across different populations. Each 25.3μL reaction consisted of 19μl H 2 O, 5μl 5X Mytaq buffer (Bioline, with pre-optimized concentrations of dNTPs and MgCl 2 ), 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].

Accession numbers
All new COI sequences generated in this study are available on Genbank: KY907195-KY907453.

Microsatellite genetic diversity and population structure
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 (R 2 = 0.02) between genetic (phiPT) and geographic distance (y = 0.0002x + 28.1). Pairwise estimates of F ST , 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 F ST estimates but the other measures are shown in S2A-S2C Table. Lowest F ST 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 (F ST = 0.011-0.103) (S2A Table). However, some comparisons between regions separated by vast geographical distances also showed low F ST scores, such as populations from the Solomon Islands with populations from mainland SE Asia (F ST = 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.

Invasion history-ABC
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 N e prior ranges for all founders (10-10,000)).
Each introduced population showed no relative change in N e due to large 95% CIs of posterior distributions, although median N e 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 N e since its introduction (due to overlapping 95% CIs), although there was a gradual increase in median N e 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).

COI haplotype networks and diversity
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. Displayed are the nine most prevalent COI haplotypes (of 52 in total) using data from ours and other studies, where each haplotype is represented as a different color and the size of the circle represents the number of individuals from a given region (which is plotted on the map). Note, the placement of circles does not correspond to the exact location of haplotypes, but represents the general region they are from; refer to S6 Table for  Spatio-temporal genetics of Ae. albopictus in the Indo-Pacific 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,  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).

Bottleneck tests
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).

Discussion
The population structure and genetic connectivity of Ae. albopictus within the Indo-Pacific region has been limited to a few studies that only examined regional structure or had restricted sampling within the species' range. In addition, the genetic characteristics of some populations examined in this study have been unexplored (e.g. Solomon Islands, many Indian and Pacific Ocean Islands, Papua-Indonesia). Using multiple lines of evidence (and both microsatellite and mitochondrial markers), we show high spatial genetic structure throughout this region. We used coalescent ABC analysis to test for the first time the likely invasion route of Ae. albopictus into the Australasian region and uncovered that the species  [18]. Furthermore, we provide evidence of rapid temporal shifts in population structure occurring less than a decade after the Asian tiger mosquito's introduction into Australia's Torres Strait Islands in 2005 [26]. In contrast, historically-introduced and native populations of Ae. albopictus showed less spatial population structure at a regional level, despite large geographic distances and international boundaries between some of these populations. Importantly, this study provides a widespread sampling distribution of the species' native range and revealed more spatial population structure than previously shown, as well as evidence for rapid temporal genetic change in newly established populations in the Torres Strait Islands.

Population structure in the Indo-Pacific
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 eastwest (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 (F ST = 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][72][73]. For our full dataset, the relationship between genetic and geographic distance was significant, however, the correlation was weak (R 2 = 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][75][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.

Invasion into Australasia
New Guinea. In New Guinea (PNG and Papua (Indonesia)), Ae. albopictus was first detected in Jayapura in the northeastern corner of Papua Province, Indonesia (formerly West New Guinea) in 1962 [55,57]. By the early 1970s it was reported in northern PNG near Madang and was established in the PNG capital Port Moresby by the 1980s [55,57,58]. There were reports of Ae. albopictus from New Guinea earlier but these have been considered doubtful and probably referred to Ae. scutellaris [56]. However, we considered these earlier dates in our ABC prior distributions for testing the timing of introduction events. We found support for a later arrival of Ae. albopictus into Papua (~1959 [95% CI = 1943, 1967]), rather than earlier dates from the 1920s, but likewise, the introduction could have been much earlier than the 1962 detection [55]. Additionally, our ABC analysis show that the Papuan population of Ae. albopictus probably originated from mainland SE Asia, rather than from our sampled Indonesian population. The subsequent spread into PNG was supported as originating from a similar mainland SE Asian source in a separate and later event, rather than from an introduction from Indonesia (including both Papua and non-New Guinea Indonesian populations in our study). The timing of the introduction into PNG in the ABC analysis (~1970 [95% CI = 1963, 1976]) corresponded with first detection dates from 1972 [57], but similar to Papua, could have been earlier than first observed which is plausible given the lack of surveys in PNG at the time.
Previous data [18] from PNG used inconsistent M13 dyes for a given microsatellite locus, causing dye shifts [35] resulting in some incorrectly scored alleles. This was uncovered and corrected in the present study (as outlined in the Materials & Methods) and revealed less substructuring of PNG populations than previously suggested [18]. We found that all historically established PNG populations were a homogenous genetic cluster (whereas in the previous study, Kiunga and Port Moresby both appeared as distinct genetic clusters using microsatellites), but this dye shift had no discernible influence on the general relationships between the other populations in their study [18]. Cooper, Waterson [58] suggested that the dispersal of Ae. albopictus in PNG is primarily driven by aircraft and coastal shipping as there is not an extensive road network in PNG. Surprisingly, the small island of Daru, offshore of southern PNG, showed no detectable genetic change using bottleneck statistics between samples analysed after 16 years (1992-2008), nor did Port Moresby over 4 years (1996)(1997)(1998)(1999); neutrality tests using COI revealed no contraction or expansion of these populations, but results were insignificant.
Solomon Islands, Fiji and Nauru. The introduction of Ae. albopictus into the Solomon Islands in 1979 and eventually into Fiji in 1989 was originally suspected to be from a progressive expansion, most likely from PNG via shipping [54,80]. Our ABC analysis suggested the Solomon Islands introduction originated from either PNG or mainland SE Asia in~1979 (95% CI = 1973CI = , 1983, corresponding with the first detection date-due to moderate type I and II error it is unclear which region was the primary source of the introduction (see Results). Major commercial shipping routes exist between the Solomon Islands and mainland SE Asia and currently the Solomon Islands has significant trade with Singapore and Malaysia [81] which could explain the ABC results and high between-region gene flow (F ST = 0.05-0.11). Large-scale logging operations have been prevalent in the Solomon Islands over the last 40 years, with a four-fold influx of logging industry multinationals, chiefly from SE Asia (especially Malaysia), during the early 1980s [82]. Shipping associated with this industry could have provided a prime opportunity for the introduction of Ae. albopictus from SE Asia. However, our STRUCTURE and DAPC plots also show the Solomon Islands has close genetic affinity to both mainland SE Asia and PNG clusters. Compared to SE Asia, PNG showed less, but relatively high gene flow with the Solomon Islands (F ST = 0.08-0.20). Papua New Guinea has strong historical and contemporary cultural ties with the Solomon Islands and other islands of the Melanesian region (spanning the islands from PNG to Fiji) and human movements between the regions have been ongoing, which may explain PNG's genetic ties to the Solomon Islands, particularly Lihir Island (S1 File, Supplementary Discussion: Lihir and Buka Islands). The majority of Solomon Islands COI sequences were H5, which is widespread throughout PNG; however, this haplotype is also observed in some mainland SE Asian populations (from Singapore and Thailand) as well as Indian Ocean Islands (Cocos (Keeling) Islands and Christmas Is.). Additionally, multiple private COI haplotypes were only found in the Solomon Islands, highlighting some unique structure in the Solomon Islands that was also shown at K = 9 in STRUCTURE (S1 Fig). For microsatellites we found there were no shared private alleles between the Solomon Islands and mainland SE Asia (i.e. alleles that were exclusively found in the two regions), but a single shared private allele was exclusive to the Solomon Islands and PNG. Overall, these results highlight that there have probably been multiple introductions into the Solomon Islands from both mainland SE Asia and PNG.
Due to insufficient sampling of Nauru and Fiji populations, we were unable to test these populations in our ABC analysis. Results are discussed in S1 File, Supplementary Discussion: Fiji and Nauru, but conclusions should be interpreted cautiously given small sample sizes. Aedes albopictus was only recently found on Nauru in April 2014 -first detected and morphologically identified by Michael Bangs, and genetically confirmed using PCR-restriction digest and COI in this study (see Materials and methods).

Torres Strait Islands and Southern Fly Region (Southern PNG).
In 1988, Ae. albopictus was detected in the coastal Southern Fly Region of PNG and on the nearby PNG island of Daru. However, a 1992 survey showed that while Ae. albopictus was still present on Daru, it was not found in/around the coastal villages of the Southern Fly Region. Consequently, Cooper, Waterson [58] highlighted that the species had failed to establish in these coastal locations after its initial introduction, despite an abundance of suitable larval habitats. However, a reintroduction occurred in 2005 at the same time as the Torres Strait incursion, which were suspected to be part of the same invasion wave into the region due to their genetic similarity [18]. Compared to Indonesian and historically established PNG populations, the Torres Strait/Fly Region invasive population was more genetically similar to Indonesian populations and was suspected to have been introduced by illegal, Indonesian fishing vessels [18]. Our ABC analysis supports the Indonesian origin of these introduced populations occurring~2005 (95% CI = 2003CI = , 2007. It is likely the Torres Strait and Fly Region were either seeded by a similar invasion source, or that local movements (mostly boat (no airstrips in the Fly Region)) after an initial invasion into the region have seeded new populations and homogenized population structure in the region. The presence of multiple shared COI haplotypes with PNG and mainland SE Asian populations also highlight that there have probably been multiple secondary introductions into the region, some potentially from unsampled populations. We found no support for initial founder effects after the Torres Strait/Fly Region introduction, but median N e showed a gradual increase after its introduction, potentially suggesting growth of the population over time. Note that in our ABC analysis this region is represented as a random subset of multiple populations (due to the complexity of modelling each population separately), so these results are more reflective of the stabilization of N e in the Torres Strait/Fly Region as a whole. In contrast, when we explored population structure on individual islands of the Torres Strait we found variable and sometimes drastic temporal changes in genetic structure (outlined below).
Our microsatellite results show high temporal genetic structure in the Torres Strait region over the past decade. Furthermore, our data indicate either a notable pattern of genetic drift on some islands (such as Keriri, Ngurupai and Warraber) or could represent a secondary introduction into the region from an unsampled population. We found evidence of genetic bottlenecks for several Torres Strait Islands, but only Waiben (2010 population) was supported by both bottleneck statistics used here (M-ratio and heterozygosity excess tests). Williamson-Natesan [83] highlighted that when a bottleneck is supported by the M-ratio, but not heterozygosity excess, it can indicate a more severe and older population bottleneck and vice versa for weaker and more recent bottleneck events. Consequently, we found evidence for older and more severe bottlenecks on Mabuiag, Waiben, Ngurupai from earlier collections after the initial invasion of Ae. albopictus into the region, as well as in more recent collections from Poruma and Iama. The dramatic wet-dry climate in the Torres Strait [77,84], intense spraying efforts and refocusing of control efforts (see S1 File, Supplementary Discussion: Control and surveying on the Torres Strait Islands; and [27,28]) could explain population bottlenecks on certain islands and the differing genetic trajectories of the islands. Additionally, initial founder effects (resulting from small invading populations) could have played a strong role in the pattern of genetic drift on the various islands. Laboratory experiments by Nicholson, Ritchie [77] showed that hotter, dryer conditions significantly reduced egg viability in Ae. albopictus from Masig Is. in the Torres Strait, so it is possible that seasonal changes have also caused genetic bottlenecks, but it is likely there are a multitude of factors at play. Both bottleneck tests used here have been shown to fail to detect recent bottlenecks (within 1-5 generations), particularly given small sample sizes and few loci [49,85]-to more accurately determine if a recent bottleneck has occurred more extensive temporal sampling will need to be conducted in the future. As previously mentioned, it remains a possibility that a secondary introduction from an unsampled source has caused the rapid temporal genetic changes observed on some islands. Indeed, Urbanelli, Bellini [67] showed that the Sulawesi region harbours highly distinct populations that can be distinguished from other Indonesian populations of Ae. albopictus and the region could be key to fully understanding the species' invasion of Australasia. Likewise, Southern PNG is highly connected to the Torres Strait Islands via local boating traffic [28] and the temporal genetic change detected could have originated from an unsampled population in this region.

Conclusion and future implications
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][87][88].
Supporting information S1   Table. Tajima's D and significance value (P value < 0.05) as well as nucleotide diversity (π) using COI sequences for populations of Aedes albopictus that had temporal collections from the present study. Sample sizes varied by year and are indicated. We found no significance (P > 0.05) in all calculations of Tajima's D. (XLSX) Populations from the Southern Fly Region, Nauru and Port Moresby (1996, 1997) were not included due to insufficient data (see Table 1 for n).  Fig. A) DAPC of the full dataset for 13 microsatellite loci for Aedes albopictus in the study region. Three-dimensional plots show the relationship between individuals belonging to 23 different populations (represented by colored dots, where the color of a dot corresponds to population) using the first 3 principal components (PC1-3). Each plot shows the same data, but is rotated along the horizontal plane. Four distinct clusters (C1-4) are indicated with dashed ellipses (not confidence intervals) and cluster membership of each population is denoted in the legend. This plot is chiefly to visualize the genetic relationships between the four main clusters and specific relationships are discussed in text. B) Correspondence analysis (CA) of the same data, but presenting population means rather than individual data points. Note that for both plots A and B, population definitions varied slightly from STRUCTURE analyses and are shown in Table 1