Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Open Access

Peer-reviewed

Research Article

Glacial History of a Modern Invader: Phylogeography and Species Distribution Modelling of the Asian Tiger Mosquito Aedes albopictus

Abstract

Background

The tiger mosquito, Aedes albopictus, is one of the 100 most invasive species in the world and a vector of human diseases. In the last 30 years, it has spread from its native range in East Asia to Africa, Europe, and the Americas. Although this modern invasion has been the focus of many studies, the history of the species’ native populations remains poorly understood. Here, we aimed to assess the role of Pleistocene climatic changes in shaping the current distribution of the species in its native range.

Methodology/Principal Findings

We investigated the phylogeography, historical demography, and species distribution of Ae. albopictus native populations at the Last Glacial Maximum (LGM). Individuals from 16 localities from East Asia were analyzed for sequence variation at two mitochondrial genes. No phylogeographic structure was observed across the study area. Demographic analyses showed a signature of population expansion that started roughly 70,000 years BP. The occurrence of a continuous and climatically suitable area comprising Southeast China, Indochinese Peninsula, and Sundaland during LGM was indicated by species distribution modelling.

Conclusions/Significance

Our results suggest an evolutionary scenario in which, during the last glacial phase, Ae. albopictus did not experience a fragmentation phase but rather persisted in interconnected populations and experienced demographic growth. The wide ecological flexibility of the species probably played a crucial role in its response to glacial-induced environmental changes. Currently, there is little information on the impact of Pleistocene climatic changes on animal species in East Asia. Most of the studies focused on forest-associated species and suggested cycles of glacial fragmentation and post-glacial expansion. The case of Ae. albopictus, which exhibits a pattern not previously observed in the study area, adds an important piece to our understanding of the Pleistocene history of East Asian biota.

Citation: Porretta D, Mastrantonio V, Bellini R, Somboon P, Urbanelli S (2012) Glacial History of a Modern Invader: Phylogeography and Species Distribution Modelling of the Asian Tiger Mosquito Aedes albopictus. PLoS ONE 7(9): e44515. https://doi.org/10.1371/journal.pone.0044515

Editor: Luciano A. Moreira, Centro de Pesquisas René Rachou, Brazil

Received: April 13, 2012; Accepted: August 3, 2012; Published: September 6, 2012

Copyright: © Porretta et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Funding: This work was funded by the Ministero dell’Istruzione, dell’Università e della Ricerca (MIUR). The funders have no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

Introduction

The Pleistocene was a time of intensive climatic fluctuations, during which glacial–interglacial cycles occurred at irregular intervals with varying durations [1]. These climatic changes substantially influenced both the current distribution of species and their genetic diversity [1]–[5]. An understanding of how species responded to past environmental change has important implications for long-term conservation of biodiversity and could also be useful for understanding species’ responses to ongoing climatic change [2], [3], [6]–[10]. Understanding these responses may be of particular interest in the case of disease vectors and/or invasive species [11]–[13]. In recent decades, much light has been shed on the general features of how species coped with Pleistocene climatic changes [3], [4], [14]–[16]. However, our knowledge is uneven across the globe in consequence of the disparity in the numbers of studies conducted in different regions. The impact of these climatic changes has been well-documented for temperate taxa that inhabit North America and Europe, whereas it has received little attention for species of Asian regions [3], [17].

In Far East Asian regions, Pleistocene climatic changes caused dramatic environmental changes, although the glacial advances were not as extensive as those in Europe or the Americas [18]. During the Last Glacial Maximum (LGM; ∼21,000 years BP), at mid-latitude (30–40° N), steppe and desert biomes extended southward as far as the modern coast of Eastern China. The broadleaf evergreen warm mixed forest shifted southward into the lowlands as far as 24° N, and tropical forests were restricted to the northern part of the Indochinese Peninsula [18], [19]. Open habitats such as savannah extended from the Malay Peninsula southward [20].

Most phylogeographic studies of animal species in these regions have focused on forest-associated species to investigate whether and how they were affected by ice-age-driven changes in forest distribution, both at mid-latitude and in the Indochinese Peninsula [21]–[26]. Attempts to identify common patterns are being carried out for those areas where a greater number of studies have been conducted [23], [25], [26]. For instance, it has been suggested that in the Indochinese Peninsula, mosquito species associated with tropical forests experienced fragmentation phases into separate northeastern and northwestern refugia. Subsequent post-glacial expansion southward from glacial refugia has also been suggested, leading to the formation of a suture zone between intraspecific genetic lineages along the Thailand–Myanmar border [26].

Despite their abundance, animal species not strictly associated with forest habitats have been little studied; as a result, how they coped with Pleistocene climatic changes remains overlooked. To contribute towards filling this gap, we investigated the Pleistocene evolutionary history of the Asian tiger mosquito Aedes albopictus, a generalist species with wide ecological plasticity. It is adapted to both tropical and temperate zones and is distributed across East Asia, from Indonesia to Korea on the north–south axis and from Japan to India on the east–west axis. In temperate regions, Ae. albopictus overwinters by photoperiodic egg diapause ([27], [28] and references therein, [29]). It is common in forested areas as well as open habitats, in both suburban and rural areas. The two most typical microhabitats of Ae. albopictus appear to be tree holes and man-made containers of any material, from glass to stone, plastic, wood, or rubber [27]. The females need very little water for oviposition and feed on a wide range of vertebrate hosts [27], [30]. Evidence of the species’ extraordinary ecological flexibility and its ability to exploit different habitats comes from its recent history. Ae. albopictus is an invasive species that, in the last few decades, has spread from Asia to the Americas, Europe, Africa, and Australia, as a result of human activities [31]. During its worldwide spread, the species has successfully competed with co-occurring mosquito species and has adapted to a variety of different environmental conditions [28], [32]–[34].

Owing to its wide ecological plasticity, Ae. albopictus may not conform to the scenario of past fragmentation during glacial phases that has been suggested for other species [35]–[39]. Previous genetic studies on Ae. albopictus populations by using nuclear markers have shown a lack of genetic discontinuity across most of the geographic range of the species (pairwise Nei’s unbiased genetic distance [40] never higher than 0.08), which could support this hypothesis [41], [42]. Here, we used two sets of data to investigate how the species coped with Pleistocene climatic changes. First, we used sequence variation data for two mitochondrial DNA genes to investigate the phylogeography and historical demography of the species across East Asia. Second, we used species distribution modelling (SDM) to reconstruct the species distribution during LGM.

The utility of integrating SDM into phylogeography has recently been highlighted by several researchers, who have shown that the inferences from one approach can be investigated and even potentially corroborated by another [15], [43]–[47]. Under a scenario of prolonged allopatric divergence in multiple refugia across East Asia during the glacial phase, we would expect to observe groups of genetically distinct populations, or highly diverging reciprocally monophyletic lineages or secondary contact zones resulting from admixture between populations previously restricted in separate refugia [2], [3]. Under this scenario, SDM would also show the occurrence of multiple geographically separated areas of suitable climatic conditions. Under an alternative scenario of a lack of population fragmentation during this period, we would expect a lack of the above genetic signatures, and SDM would show the occurrence of continuous suitable climatic conditions across the study area.

Materials and Methods

Ethics Statements

No specific permits were required for the described field studies. The localities, where the mosquitoes were collected, are no private or protected areas, and the study did not involve endangered or protected species.

Genetic Analyses

Adults of Ae. albopictus were collected from 16 localities through the East Asia to cover the range of distribution of the species in both tropical and temperate areas (Figure 1A, Table 1) and the refugia regions inferred for other animal species [21]–[26]. By using aspirators, mosquitoes resting near humans outdoors were captured in three-five places per site sampling. Then, in laboratory, mosquitoes were killed by exposure to ethyl acetate vapour, and were placed individually in eppendorf 1.5 ml tube with ethanol 80%.

thumbnail
Download:
Figure 1. Sampling sites and phylogenetic relationships among mitochondrial haplotypes.

(A) Map of East Asia showing the sampling sites of Aedes albopictus. For localities details see Table 1. (B) Statistical parsimony network, constructed using TCS software, showing phylogenetic relationships among the haplotypes found. Haplotypes are shown as circles with sizes corresponding to their frequencies in the total sample and colour corresponding to populations where they have been observed. Haplotypes are encoded as in Table 1. Dots indicate missing intermediate haplotypes.

https://doi.org/10.1371/journal.pone.0044515.g001

thumbnail
Download:
Table 1. Geographical origin, sample size and haplotypes’ distribution for the 16 sampled populations of Aedes albopictus (haplotypes in bold have been found in two or more localities; the number in brackets indicates how many time a haplotype has been observed at a particular locality).

https://doi.org/10.1371/journal.pone.0044515.t001

Total genomic DNA was extracted from single mosquitoes following the standard CTAB (cetyltrimethyl ammonium bromide) protocol [48]. Partial sequences of the mitochondrial DNA genes encoding for the cytochrome oxidase I (COI) and nicotinamide adenine dinucleotide dehydrogenase subunit 5 (ND5) were obtained through Polymerase Chain Reaction (PCR). The primers TW-J-1305 and TK-N-3782 [38] were used to amplify and sequence ∼2000 base pairs including the regions COI, tRNA-Leu and COII. The following primers were then designed to amplify and sequencing two regions of the COI gene: AealCOIa-f 5′- aaaaagatgtatttaaatttcggtctg-3′ and AealCOIa-r 5′-tgtaattgttactgctcatgctttt-3′; AealCOIb-f 5′-ttattacacaagaaagaggaaaaa-3′ and AealCOIb-r 5′-cattgcactaatctgccata- 3′. The ND5 gene region was amplified using the primers designed by Birungi & Munstermann [49]. Each PCR amplification was performed in 25 µl including 1x Buffer, 200 µM dNTPs, 2.5 mM MgCl2, primers at 0.2 µM, 0.5 units of Taq DNA Polymerase (GoTaq™ DNA Polymerase, Promega), and 20 ng of genomic DNA extracted from single mosquito. The PCR cycling procedure was: 95°C for 5 min followed by 34 cycles at 93°C for 1 min, annealing temperature: 55°C for COI and COIb; 50°C for ND5 for 1 min, 72°C for 1 min 30s, and a single final step at 72°C for 10 min. PCR sequences were obtained using an ABI PRISM 3700 DNA sequencer by Macrogen Inc. (www.macrogen.com). All individuals analysed were sequenced using both forward and reverse primers, and all haplotypes that were found unique (see Result section) were re-amplified by using high fidelity Taq DNA Polymerase (Phusion® High-Fidelity DNA Polymerase, Fermentas) and re-sequenced to check for consistency.

Sequences were edited and aligned using the software Chromas 2.33 (Technelysium Pty Ltd, Australia) and Clustalx 2.0 [50] respectively.

Nucleotide and amino-acidic polymorphisms of the COI and ND5 genes were assessed using the software Mega 5.0 [51]. For all the subsequent analyses the two genes were concatenated since the partition-homogeneity test – conducted with 1000 replicates using the software Paup* 4.0b10 [52] - did not reject the null hypothesis of homogeneity of their phylogenetic signals. The software jModelTest version 0.1.1 [53] was used to find the optimal model of sequence evolution for our data using the Akaike Information Criterion. The HKY+I+Γ substitution model with the alpha value of the gamma distribution = 0.812 and the proportion of invariable sites I = 0.862 was shown as the best-fit model for the data. Haplotype diversity (h) and nucleotide diversity (π) were estimated for each locality (with the exception of Nagasaki population where 4 individuals were sampled) and for the overall dataset, using the software Arlequin 3.1 [54].

The genealogical relationships between haplotypes were investigated by constructing a phylogenetic network [55]. We used the statistical parsimony algorithm described by Templeton et al. [56] and implemented in the TCS software [56], applying a 95% cutoff for the probability of a parsimonious connection. To check for consistency [57] we also used the median-joining (MJ) network algorithm as implemented in the Network 4.6.1.0 software (Fluxus Technology Ltd). The loops in the resulting phylogenetic network were resolved by applying the criteria described by Pfenninger & Posada [58].

The possible occurrence of groups of genetically distinct populations was investigated by Spatial Analysis of Molecular Variance as implemented in the software samova 1.0 [59]. This method defines groups of populations on the basis of genetic data, without any a priori hypothesis of the expected structure. Given a user-defined number K of groups of populations, samova uses a simulated annealing procedure to define groups of populations that are as genetically homogenous as possible (among population differentiation index, FSC, minimized) and maximally differentiated from one another (among groups differentiation index, FCT, maximized). We tested K values ranging from 2 to 15 and, to check for consistency, we conducted each analysis five times, with 10,000 independent annealing processes each.

To investigate signatures of past demographic changes we analysed the distribution of pairwise nucleotide differences between haplotypes (mismatch distribution) as implemented in the software Arlequin 3.1. In populations that have undergone a sudden demographic expansion, the mismatch distributions are expected to be smooth and bell-shaped [60]. The sum of square deviations (SSD) between estimated and observed mismatch distributions was used as goodness-of-fit statistics and its significance was tested through 1000 replicates. To further detect signatures of past demographic changes, we also used Fu’s FS neutrality test [61]. Large negative values of the parameter FS indicate an excess of singleton mutations, which is expected both under departures from selective neutrality and in populations which have experienced a recent demographic expansion. It was computed on Arlequin 3.1, and its significance was assessed through 1000 replicates. Finally, historical demographic trend was investigated using the Bayesian skyline plot (BSP) method by Drummond et al. [62], as implemented in the software beast 1.6.1 [63]. This method estimates the posterior distribution of the effective population size through time using a Markov chain Monte Carlo (MCMC) sampling procedure. We used the HKY+I+Γ substitution model of sequence evolution. To give an approximate time scale to the inferred demographic trend, we used a mutation rate of 0.0115 substitutions/site/lineage/Myr, based on the divergence rate of 2.3% per million years [64]. It may be considered as an appropriate rate for mosquitoes, as this divergence rate has been estimated from several arthropod datasets, including several Diptera. Furthermore, it has been recently used in different mosquito species and the demographic trend inferred was shown to accord with independent paleoenvironmental data [26]. A relaxed molecular clock model was used, with uncorrelated rates drawn from a log-normal distribution. Preliminary batch runs were carried out using the auto-optimization procedure implemented in beast in order to fine-tune MCMC parameters. The MCMCs were finally run for 20 million steps, and sampled every 1000 steps. The first two million steps were discarded as burn-in, after checking the likelihood history with tracer 1.4 [65]. The same software was used to analyse the results of the BSP analyses. Five replicates of the BSP analysis were run to check for consistency among different runs.

Species Distribution Modelling

The Species Distribution Model (SDM) for Ae. albopictus was developed using the “maximum entropy model” as implemented in Maxent 3.3.3e [66]. This software reconstructs the potential distribution of a species combining presence-only data with climatic layers, distinguishing presence from random [66]. Although our study is focused on East Asian regions, all SDMs were constructed including also the Indonesian Islands since, during the LGM, they were connected to the mainland Asia as a consequence of the Sundaland exposure (see Discussion). All models were developed using a total of 153 occurrence points of Ae. albopictus, including records from literature, records from the “Global Biodiversity Information Facility” (GBIF) online database (http://data.gbif.org) and our sampled localities. We selected all georeferenced localities, while among those without geographical coordinates, we used only localities for which exhaustive geographical data are furnished by the authors. These localities were then georeferenced using the software Google Earth.

The SDM was developed for Current, Last Glacial Maximum (LGM) and Last Inter-Glacial (LIG) conditions. To reconstruct the ecological niche of LGM conditions, we used two general circulation models (GCM): the Community Climate System Model 3 (CCSM 3) [67] and the Model for Interdisciplinary Research on Climate (MIROC 3.2) [68]. For all models, we used the bioclimatic variables available in WorldClim database (http://www.worldclim.org/). These 19 variables, derived from monthly temperature and rainfall values, were downloaded at spatial resolution of 2.5 arc-minutes for current and LGM conditions, and at spatial resolution of 30 arc-seconds for LIG conditions.

Pearson’s correlation matrix was generated for all 19 variables to quantify the correlation between them. For pairs that were highly correlated (correlation coefficient ≥0.75) we chose the variable more biologically meaningful for the species [27]. Finally, all SDMs were constructed using four variables, including: BIO6 (Min Temperature in Coldest Period); BIO8 (Mean Temperature of Wettest Quarter); BIO18 (Precipitation of Warmest quarter) and BIO19 (Precipitation of Coldest Quarter).

In Maxent, each run was conducted using the default convergence threshold (10−5) and maximum number of iterations (500) values. For the regularization parameter β, as suggested by Warren & Seifert [69], we tested the values 1, 3, 5, 7, 9, 11, 13, 15, 17 and 19 using the software EnmTools [70]. The best model, chosen using the corrected Akaike Information Criterion (AICc), was one for which β = 3. Twenty-five percent of record points were used for model testing, and the final SDMs performance was evaluated using the area under receiver operating characteristic (ROC) curve (AUC). An AUC score above 0.7 is considered good model performance [71]. All SDM predictions were visualized in DIVA GIS 5.2 [72]. To provide a consensus map of the species distribution at LGM we combined the CCSM and MIROC predictions (available under request to the authors) by including areas predicted as suitable by both models, while discarding areas where no agreement occurred [15]. We used the software ImageJ [73] to compute the area of suitable habitats of the species under the LIG and LGM conditions using the “minimum training presence” threshold for presence-absence.

Results

Genetic Data

A final sequence alignment including 1458 nucleotides (364 from the ND5 gene fragment, 405 and 689 base pairs respectively for AealCOIa and AealCOIb primer pairs) was obtained for 174 individuals. Sixty two haplotypes were found (GenBank Accession numbers: JQ436947-JQ437008) defined by 54 nucleotide substitutions (32 parsimony informative). Forty one nucleotide substitutions were located at the third codon position, five at the second codon position and eight at the first one. Thirteen non-synonymous nucleotide substitutions were found. Overall haplotype and nucleotide diversity estimates were 0.945 (±0.009) and 0.0021 (±0.0012), respectively.

No apparent association between haplotypes and geography was observed. The three more frequent and centrally placed haplotypes in the genealogical network (h1, h4, and h22), were found in populations spanning the entire east-west and north-south extension of the study area (Table 1, Figure 1).

Phylogenetic relationships among the 62 haplotypes found are shown in Figure 1B. Both statistical parsimony and median-joining networks showed the same topology, so only the statistical parsimony network is shown. The uncorrected mean sequence divergence among haplotypes was 0.30%. A large portion (39/62, 63%) of the haplotypes found were unique (see also Table 1). No haplotype-groups were apparent being all haplotypes connected to each other by no more than 3 nucleotide substitutions (Figure 1B).

By the spatial analysis of molecular variance (SAMOVA) no groups of genetically distinct populations were evidenced. All partitions of the sampling areas for each K value showed no significant increment in the FCT values. They showed a very narrow range (12/14 values were between 0.28 and 0.36) and the highest value (FCT = 0.40) was observed for K = 15, when all populations minus one were separated (Figure S1).

The haplotype diversity and nucleotide diversity estimates for each population are shown in Table 1. They ranged from 0.286 (sample 5, China) to 0.950 (sample 10, Vietnam) and from 0.0004 (sample 6, China) to 0.0023 (sample 9, Vietnam). The geographical distribution of haplotypes and the indices of population diversity showed the occurrence of an uneven distribution of genetic diversity throughout the study area (Table 1, Figure 1). The populations from Bhutan, Southern China, and Indochinese Peninsula showed the highest haplotype diversity. In these samples we also found all the most common haplotypes (h1, h4, h22) and those centrally placed in the network (h1, h4, h22, but also h25, h32).

Demographic analyses showed the occurrence of past demographic expansion using both separated (data not shown) and concatenated COI and ND5 genes (Figure 2). The mismatch distributions appeared smooth and bell-shaped. The SSD statistic supported a close fit of the observed distribution to the distribution expected under population expansion model (SSD = 0.002, P = 0.46) (Figure 2A). The past occurrence of a demographic expansion was also suggested by Fu’s Fs and the Bayesian skyline plot method (BSP). The former was large, negative and highly significant (Fs = −26.244, P<0.001), showing the occurrence of a larger than expected number of singleton mutations, as expected for populations that underwent a population expansion. The demographic trend inferred by the BSP also indicated the advent of a phase of demographic growth, which roughly started about 70,000 years BP until recently (Figure 2B).

thumbnail
Download:
Figure 2. Historical demographic analyses.

(A) Mismatch distribution computed using the software Arlequin 3.1. Histograms show the observed distribution; lines show the expected distribution under a model of sudden population expansion. (B) Bayesian Skyline Plot, constructed using the software beast 1.6.1. Population size (y-axis) is measured as the product of effective population size per generation length (Neτ). The solid line is the median estimate, and the grey areas show the 95% Highest Posterior Density (HPD) limits. Time (x axis) is expressed in years before present (BP).

https://doi.org/10.1371/journal.pone.0044515.g002

Species Distribution Modelling

The Maxent models predicted under LIG, LGM and current conditions of Ae. albopictus are shown on Figure 3. For the current condition model, the AUC for the training points was 0.935 and for test points was 0.899, indicating a good performance of the trained model. Our reconstruction under current conditions is concordant with those previously published by Benedict et al. [31] and Medley [34], which further supports the reliability of our predictions.

thumbnail
Download:
Figure 3. Bioclimatic models for Aedes albopictus.

(A) Last Interglacial (LIG; ∼140,000–120,000 years BP); (B) Last Glacial Maximum (LGM, ∼21,000 years BP); (C) current conditions. All models were developed using the “maximum entropy model” as implemented in the software Maxent 3.3.3e. In dark grey are shown the areas predicted as suitable by both the CCSM and MIROC models.

https://doi.org/10.1371/journal.pone.0044515.g003

The distribution inferred at LIG shows the occurrence of climatically suitable areas both in tropical and subtropical regions (i.e. Indonesia, Indochinese Peninsula) and with lesser degree in temperate regions (i.e. China and Japan), similarly to the current conditions (Figure 3A and C). The projection at the LGM predicted the occurrence of mostly continuous areas with climatically suitable conditions for the species comprising South-East China, Indochinese peninsula and the emerged Sundaland. The emerging of the Sunda shelf resulted at the LGM in an increase of the potential suitable area for Ae. albopictus with respect to the LIG of 42% (by using the minimum training presence threshold, 0.092). At the higher latitudes, in the northernmost part of its current range of distribution (Northern China and Japan), areas of predicted suitable habitats were mainly absent (Figure 3B).

Discussion

mtDNA Diversity in Native Populations of Aedes albopictus

Previous studies of Ae. albopictus mitochondrial diversity have shown that native populations of this species harbour little or no mitochondrial genetic diversity ([49], [74] and references therein). In particular, an analysis by Birungi & Munstermann of 405 base pairs of the ND5 mtDNA gene in samples derived from laboratory colonies revealed only 1 haplotype across an area spanning Japan, Malaysia, and Indonesia [49]. Armbruster et al. [74] suggested that the overall low mtDNA diversity observed with respect to nuclear markers in the native Asian populations was a result of a cytoplasmic ‘sweep’ caused by double infection of Ae. albopictus by the endosymbiont Wolbachia. An alternative explanation is that the mtDNA diversity of natural populations has been obscured by sampling bias: the Ae. albopictus individuals that were analysed were derived from laboratory colonies (F10–F23 generations) that did not represent the genetic diversity in natural populations. This explanation was supported by subsequent studies in which individuals from natural populations were used [75]–[79]. Our study lends further support to the hypothesis that the mtDNA diversity of natural populations has been overlooked. All the specimens used in this study showed positive results for infection with both wAlbA and wAlbB strains of Wolbachia ([80] and our unpublished data), as expected for Ae. albopictus populations [74]. We found not only a non-negligible genetic diversity–within the range of that observed in other mosquito species uninfected by Wolbachia [26], [48], [81]–but also an overall concordance between the genetic patterns observed at nuclear and mitochondrial markers; this would not be expected if mtDNA selective sweep had occurred [82]. Previous population genetics studies conducted using allozymes showed significant genetic differentiation between populations at a regional level, but a lack of discontinuity in the geographic distribution of genetic diversity at the level of mainland Asia as a whole [41], [42], [83]. Finally, the topology of the network showing the phylogenetic relationships between the mtDNA haplotypes (i.e., multiple starlike topology, Figure 1B) is quite different from that of the single star-like network expected under a scenario of recent selective sweep driven by maternally inherited symbionts such as Wolbachia [82]. Mitochondrial DNA has been widely used to infer the geographic origin of invasive populations of Ae. albopictus [49], [75], [77], [79]. The comparison between native and new colonizing populations is outside the scope of our study. However, our mtDNA survey adds new haplotypes that will supplement the poor genetic pool available to date and can be used in future studies.

Late Pleistocene History of Aedes albopictus

The genetic pattern observed in Ae. albopictus populations is hard to reconcile with a scenario of fragmentation in multiple refugia across East Asia during the last glacial phase, as has been proposed for forest-associated species. Indeed, we found neither phylogeographic structure nor highly diverging reciprocally monophyletic lineages, as would be expected if a prolonged allopatric divergence had occurred. It might be possible to explain the lack of phylogeographic structure by admixture between populations previously restricted in separate refugia. However, in this case, we would expect the existence of a genetic signature of prior allopatric fragmentation (i.e., the occurrence of diverging monophyletic lineages or secondary contact zones), which we did not observe.

In contrast, the genetic pattern observed best fits with a scenario in which the species never experienced a prolonged allopatric divergence but rather persisted in a wide area, most likely across Southeast China and the Indochinese Peninsula. In these regions, we found not only a lack of phylogeographic structure but also higher levels of genetic diversity (Figure 1), as well as the occurrence of the most common haplotypes and those centrally placed in the network (the signature of areas of the refugia) [14].

This scenario is also consistent with the reconstruction of the species distribution at the LGM (Figure 3B). Although we used a conservative criterion to combine the CCSM and MIROC predictions, the models predicted the occurrence of continuous climatically suitable areas not only in the southernmost regions of Sundaland but also through most of the Indochinese Peninsula northward as far as Southeast China and Bhutan (Figure 3B). Here, the occurrence of interconnected areas with suitable climatic conditions for the species was inferred. It is notable that paleoenvironmental reconstructions also showed the occurrence of different habitats (such as forests or open habitats) in these regions [18] that Ae. albopictus could have exploited and maintained large interconnected populations.

In the genetic pattern observed, we also found the signature for past demographic changes in Ae. albopictus populations. The demographic trend inferred from the BSP is that of a gradual expansion starting about 70,000 years BP, that is, during the last glacial phase including the LGM. Although the lack of a calibrated molecular clock demands caution, paleoenvironmental reconstructions of the region and SDM at the LGM support the hypothesis that populations of the species expanded during this period.

During the last glacial phase in Southeast Asian regions, the glacio-eustatic lowering of the sea level fully exposed the Sunda shelf (Sundaland), which connected the major islands of Sumatra, Java, and Borneo to mainland Asia by land bridges [20], [84]. At 70,000 years BP, the sea level was 20–50 m below its present level, and all these islands were joined to mainland Asia [85]. As shown by the SDM, the emergence of Sundaland led to a non-negligible increase (42%) in climatically suitable areas for the species compared to the LIG (Figure 3A–B). These emerged areas, according to paleoenvironmental reconstructions, were covered by tropical forests in areas of Sumatra, Java, and Eastern Borneo and by open habitats, such as savannah, from the Malay Peninsula southward [84]. The availability of large amounts of suitable habitat and the occurrence of suitable climatic conditions in the emerged Sundaland (Figure 3B) could have had a positive influence on Ae. albopictus populations and could account for the expansion inferred from genetic data. The effect of glacial-induced lowland emergence on population dynamics has recently been investigated in the Sardo-Corso region in a study by Bisconti et al. [86], which showed that the emergence of lowland plains during the last glacial phase counterbalanced the negative demographic consequences of climatic changes during the last glacial phase by providing new suitable habitats and led to net demographic expansions in a temperate and thermophilic species, the tree frog Hyla sarda.

The inferred demographic expansion of Ae. albopictus populations also encompasses the post-glacial phase (from 14,000 years BP until recently, Figure 2B). Following the LGM, the amelioration of climatic conditions at mid-latitude could have sustained the demographic expansion of the species into the northernmost areas of its current distribution (Central and Northern China and Japan) (see also Figure 3). The signature of this colonization event is still detectable in the geographic pattern of genetic diversity of the species: the northernmost populations (sample 1–6, Table 1) show lower levels of genetic diversity and a higher percentage of derived haplotypes than do the southernmost populations, as expected in recently originated populations [3]. It is possible to speculate that a factor beyond climatic conditions may have played a role in the post-glacial population dynamics of Ae. albopictus. Since approximately 15,000 years BP, human populations experienced a south-to-north expansion into East Asia [87]. Furthermore, several coastal settlements developed as a consequence of the post-glacial rise in sea level and the formation of coastal reef and lagoon/estuary systems, which allowed the exploitation of marine resources [88]. Finally, during the same period, in Southern China (at approximately 9,000 years BP), the transition from a hunting-gathering lifestyle to farming was initiated [89]. Ae. albopictus, unlike other tree-hole mosquitoes, is able to exploit human-associated habitats, so these human activities could have furnished the mosquito with an increased blood supply (humans themselves, along with their livestock) and suitable breeding sites; this could have favoured the inferred population expansion.

In this study, we focused on Ae. albopictus populations in the temperate regions of East Asia and in the Indochinese Peninsula. The results allow us to propose some hypotheses concerning the southernmost part of the range of the species, to be tested in future studies. SDM data showed the occurrence of climatically suitable areas across all of Sundaland. Therefore, it could be hypothesized that the emergence of Sundaland allowed population connectivity across the Indonesian islands and between these islands and mainland Asia. This scenario has recently been inferred for other mosquito species (such as Anopheles vagus s.l.) that can exploit open habitats [90] and for the fruit bat Cynopterus brachyotis [35]. Phylogeographic and historical demographic analyses based on extensive sampling throughout the Indonesian islands and New Guinea would allow us to test the possible occurrence of a similar scenario for Ae. albopictus. Our preliminary data show high genetic similarity in mitochondrial DNA among the populations from East Asia and 2 populations from Sulawesi Island (0.63%, mean sequence divergence) [91].

Conclusions

Previous genetic studies on Ae. albopictus populations that were conducted using nuclear markers have shown a lack of genetic discontinuity across most of the geographic range of the species [41], [42]. Here, by integrating phylogeographic and SDM data, we have provided strong support for an evolutionary scenario in which this generalist mosquito persisted in widely interconnected populations in mainland East Asia during the last glacial phase. A similar scenario has been found in species with wide ecological flexibility in Australia [36], in Southeast Asia [35] and in the Western Palaearctic [34], [39], [37], [92], most likely as a result of their ability to exploit different habitats or to adapt to different resources or different hosts.

Historical demographic analyses also showed that Ae. albopictus experienced a demographic expansion that started in the last glacial phase and lasted until recently. The population expansion following the LGM, as expected, was driven by amelioration of climatic conditions at mid-latitudes and perhaps also by environmental changes attributable to human activities. On the contrary, the population expansion that occurred before the LGM (from 70,000 years BP) is an interesting exception. SDM, along with paleoenvironmental reconstructions and life traits of the species, suggest a positive role of Sundaland in the dynamics of Ae. albopictus populations. Also in this case, the ecological flexibility and adaptability of the species most likely allowed it to exploit the wide range of habitats offered by the emergence of Sundaland. Therefore, the case study of Ae. albopictus contributes to our understanding of East Asian phylogeography and supports the findings of those studies that showed the positive role of the glacial-induced emergence of lowlands in maintaining populations during glacial phases [36], [38], [86].

Supporting Information

Figure S1.

Spatial analysis of molecular variance. Fixation indices, obtained using the software SAMOVA, for the best clustering option for each pre-defined values of K. FCT: variation among groups of populations; FSC: variation among populations within groups; FST: variation among populations among groups.

https://doi.org/10.1371/journal.pone.0044515.s001

(TIF)

Acknowledgments

We thanks Dr Yukiko Higa, Department of Vector Ecology and Environment, Institute of Tropical Medicine, Nagasaki University, Japan; Dr Hwa-Jen Teng, Center for Disease Control, Department of Health, Republic of China; Mr Tashi Tobgay, The Vector-Borne Disease Control Programme, Gelephu, Bhutan, for their support and collection of specimens; Dr Catherine Walton, Dr. Daniele Canestrelli for fruitful comments on the manuscript; Alessandra Spanò and Simone Sabatelli for technical assistance; Mark Eltenton and Frances Dorfman for the linguistic revision.

Author Contributions

Conceived and designed the experiments: DP SU. Performed the experiments: PS VM SU RB DP. Analyzed the data: DP VM. Contributed reagents/materials/analysis tools: PS VM SU RB DP. Wrote the paper: DP VM.

References

  1. 1. Quante M (2010) The Changing Climate: Pasr, Present, Future. In Relict Species: Phylogeography and Conservation Biology, Springer-Verlag Berlin Heidelberg. 9–56.
    • 2. Hewitt GM (2000) The genetic legacy of the Quaternary ice ages. Nature 405: 907–913.
    • 3. Hewitt GM (2004) Genetic consequences of climatic oscillations in the Quaternary. Phil Trans Roy Soc Lond B Biol Sci 359: 183–195.
    • 4. Schmitt T (2007) Molecular biogeography of Europe: Pleistocene cycles and postglacial trends. Front Zool 4: 11.
    • 5. Hofreiter M, Stewart J (2009) Ecological Change, Range Fluctuations and Review Population Dynamics during the Pleistocene. Curr Biol 19: R584–R594.
    • 6. Petit RJ, Aguinagalde I, de Beaulieu JL, Bittkau C, Brewer S, et al. (2003) Glacial refugia: hotspots but not melting pots of genetic diversity. Science 300: 1563–1565.
    • 7. Canestrelli D, Aloise G, Cecchetti S, Nascetti G (2010) Birth of a hotspot of intraspecific genetic diversity: notes from the underground. Mol Ecol 19: 5432–5451.
    • 8. Taberlet P, Cheddadi R (2002) Quaternary refugia and persistence of biodiversity. Science 297: 2009–2010.
    • 9. Jackson ST, Williams JW (2004) Modern Analogs in Quaternary Paleoecology: Here Today, Gone Yesterday, Gone Tomorrow? Annu Rev Earth Planet Sci 32: 495–537.
    • 10. Lawing AM, Polly PD (2011) Pleistocene Climate, Phylogeny, and Climate Envelope Models: An Integrative Approach to Better Understand Species’ Response to Climate Change. PLoS ONE 6(12): e28554 doi:10.1371/journal.pone.0028554.
    • 11. De La Rocque S, Rioux JA, Slingenbergh J (2008) Climate change: effects on animal disease systems and implications for surveillance and control. Rev Sci Tech 27(2): 339–54.
    • 12. Gray JS, Dautel H, Estrada-Peña A, Kahl O, Lindgren E (2009) Effects of climate change on ticks and tick-borne diseases in Europe. Interdiscip Perspec Infect Dis, 593232.
      • 13. Walther G-R, Roques A, Hulme PE, Sykes MT, Pysek P, et al. (2009) Alien species in a warmer world: risks and opportunities. Trends Ecol Evol 24(12): 686–693.
      • 14. Avise JC (2000) Phylogeography: the history and formation of species. Harvard University Press, Cambridge, MA.
        • 15. Waltari E, Hijmans RJ, Peterson AT, Nyári ÁS, Perkins SL, et al. (2007) Locating Pleistocene refugia: comparing phylogeographic and ecological niche model predictions. Plos ONE 2: 563–574.
        • 16. Svenning J-C, Normand S, Kageyama M (2008) Glacial refugia of temperate trees in Europe: insights from species distribution modelling. J of Ecol 96: 1117–1127.
        • 17. Beheregaray LB (2008) Twenty years of phylogeography: the state of the field and the challenges for the Southern Hemisphere. Mol Ecol 17: 3754–3774.
        • 18. Yu G, Gui F, Shi Y, Zheng Y (2007) Late marine isotope stage 3 palaeoclimate for East Asia: A data-model comparison. Paleogeog Palaeoclim Palaeoecol 250: 167–183.
        • 19. Rai N, Adams JM (2001) A GIS-based vegetation map of the world at the last glacial maximum (25,000–15,000 BP). Internet Archeology. Available: http://www.ncdc.noaa.gov/paleo/pubs/ray2001/ray2001.html.
          • 20. Bird MI, Taylor D, Hunt C (2005) Palaeoenvironments of insular Southeast Asia during the Last Glacial Period: a savanna corridor in Sundaland? Quat Sci Rev 24: 2228–2242.
          • 21. Iyengar A, Babu VN, Hedges S, Venkataraman AB, Maclean N, et al. (2005) Phylogeography, genetic structure, and diversity in the dhole (Cuon alpinus). Mol Ecol 14: 2281–2297.
          • 22. Yong-Chao S, Yung-Hau C, Sin-Che L, I-Min T (2007) Phylogeography of the giant wood spider (Nephila pilipes, Araneae) from Asian–Australian regions. J Biogeogr 34: 177–191.
          • 23. Fuchs J, Ericson PGP, Pasquet E (2008) Mitochondrial phylogeographic structure of the white-browed piculet Sasia ochracea: cryptic genetic differentiation and endemism in Indochina. J Biogeogr 35: 565–575.
          • 24. Vidya TNC, Sukumar R, Melnick DJ (2009) Range-wide mtDNA phylogeography yields insights into the origins of Asian elephants. Proc Roy Soc Lond B Biol Sci 276: 893–902.
          • 25. Sakka H, Quéré JP, Kartavtseva I, Pavlenko M, Chelomina, et al (2010) Comparative phylogeography of four Apodemus species (Mammalia: Rodentia) in the Asian Far East: evidence of Quaternary climatic changes in their genetic structure. Biol J Linn Soc 100: 797–821.
          • 26. Morgan K, O’loughlin SM, Chen B, Linton YM, Thongwat D, et al. (2011) Comparative phylogeography reveals a shared impact of pleistocene environmental change in shaping genetic diversity within nine Anopheles mosquito species across the Indo-Burma biodiversity hotspot. Mol Ecol 20: 4533–4549.
          • 27. Hawley WA (1988) The biology of Aedes albopictus. J Am Mosq Contr Ass 1: 1–39.
          • 28. Lounibos LP (2002) Invasions by insect vectors of human disease. Annu Rev Entom 47 233–266.
          • 29. Lounibos LP, Escher RL, Lourenço-de-Oliveira R (2003) Asymmetric evolution of photoperiodic diapause in temperate and tropical populations of Aedes albopictus (Diptera: Culicidae). Ann Entomol Soc Am 96: 512–518.
          • 30. Delatte H, Desvars A, Bouétard A, Bord S, Gimonneau G, et al. (2010) Blood-feeding behavior of Aedes albopictus, a vector of Chikungunya on La Réunion. Vector Borne Zoonotic Dis 10: 249–258.
          • 31. Benedict MQ, Levine R, Hawley W, Lounibos L (2007) Spread of the tiger: global risk of invasion by the mosquito Aedes albopictus. Vector Borne Zoon Dis 7: 76–85.
          • 32. Carrieri M, Bacchi M, Bellini R, Maini S (2003) On the competition occuring between Aedes albopictus and Culex pipiens (Diptera: Culicidae) in Italy. Environ Entomol 32(6): 1313–1321.
          • 33. Armistead JS, Arias JR, Nishimura N, Lounibos LP (2008) Interspecific larval competition between Aedes albopictus and Aedes japonicus (Diptera: Culicidae) in northern Virginia. J. Med Entomol 45: 629–637.
          • 34. Medley KA (2009) Niche shifts during the global invasion of the Asian tiger mosquito, Aedes albopictus Skuse (Culicidae), revealed by reciprocal distribution models. Global Ecol Biogeogr 19: 122–133.
          • 35. Campbell P, Schneider CJ, Adnan AM, Zubaid A, Kunz TH (2006) Comparative population structure of Cynopterus fruit bats in peninsular Malaysia and southern Thailand. Mol Ecol 15: 29–47.
          • 36. Burns EM, Eldridge MDB, Crayn DM, Houlden BA (2007) Low phylogeographic structure in a wide spread endangered Australian frog Litoria aurea (Anura: Hylidae). Conserv Genet 8: 17–32.
          • 37. Valdiosera CE, Nuria G, Anderung C, Dalén L, Crégut-Bonnoure E, et al. (2007) Staying out in the cold: glacial refugia and mitochondrial DNA phylogeography in ancient European brown bears. Mol Ecol 16: 5140–5148.
          • 38. Porretta D, Canestrelli D, Urbanelli S (2011a) Bellini R, Schaffner F, et al. (2011a) Southern crossroads of the Western Palaearctic during the Late Pleistocene and their imprints on current patterns of genetic diversity: insights from the mosquito Aedes caspius.. J Biogeogr 38: 20–30.
          • 39. Teacher AGF, Thomas JA, Barnes I (2011) Modern and ancient red fox (Vulpes vulpes) in Europe show an unusual lack of geographical and temporal structuring, and differing responses within the carnivores to historical climatic change. BMC Evol Biol 11: 214–223.
          • 40. Nei M (1978) Estimation of average heterozygosity and genetic distance from a small number of individuals. Genetics 89: 538–590.
          • 41. Black WC IV, Hawley WA, Rai KS, Craig GB Jr (1988) Breeding structure of a colonizing species: Aedes albopictus in peninsular Malaysia and Borneo. Heredity 61: 439–446.
          • 42. Kambhampati S, Black IV WC, Rai KS (1991) Geographic origin of United States and Brazilian populations of Aedes albopictus inferred from allozyme analysis. Heredity 67: 85–94.
          • 43. Svenning JC, Fløjgaard C, Marske KF, Nógues-Bravo D, Normand S (2011) Applications of species distribution modeling to paleobiology. Quat Sci Rev 30: 2930–2947.
          • 44. Hickerson MJ, Carstens BC, Cavender-Bares J, Crandall KA, Graham CH, et al. (2010) Review Phylogeography’s past, present, and future: 10 years after Avise, 2000. Mol Phylogenet Evol 54: 291–301.
          • 45. Kozak H, Graham CH, Wiens JJ (2008) Integrating GIS-based environmental data into evolutionary biology. Trends Ecol. Evol 23: 141–148.
          • 46. Knowles LL (2009) Statistical phylogeography. Annu Rev Ecol Evol Syst 40: 593–612.
          • 47. Richards CL, Bryan C, Carstens L, Knowles L (2007) Distribution modelling and statistical phylogeography: an integrative framework for generating and testing alternative biogeographical hypotheses. J Biogeogr 34: 1833–1845.
          • 48. Porretta D, Canestrelli D, Bellini R, Celli G, Urbanelli S (2007) Improving insect pest management through population genetic data: a case study of the mosquito Ochlerotatus caspius (Pallas). J Appl Ecol 44: 682–691.
          • 49. Birungi J, Munstermann LE (2002) Genetic structure of Aedes albopictus (Diptera: Culicidae) populations based on mitochondrial ND5 sequences: evidence for an independent invasion into Brazil and United States. Ann Entomol Soc Am 95: 125–132.
          • 50. Thompson JD, Gibson TJ, Plewniak F, Jeanmougin F, Higgins DG (1997) The Clustal_X windows interface: flexible strategies for multiple sequence alignment aided by quality analysis tools. Nucleic Acids Res 25: 4876–4882.
          • 51. Tamura K, Peterson D, Peterson N, Stecher G, Nei M, et al. (2011) MEGA5: molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Mol Biol Evol 28: 2731–2739.
          • 52. Swofford DL (2001) PAUP. Phylogenetic Analysis Using Parsimony (and Other Methods), Version 4. Sinauer Associates, Sunderland, MA.
            • 53. Posada D (2008) JMODELTEST: phylogenetic model averaging. Mol Biol Evol 25: 1253–1256.
            • 54. Excoffier L, Laval G, Schneider S (2005) ARLEQUIN (version 3.01): an integrated software package for population genetic data analysis. Evolutionary Bioinformatics Online 1: 47–50.
            • 55. Posada D, Crandall KA (2001) Intraspecific gene genealogies: trees grafting into networks. Trends Ecol Evol 16: 37–45.
            • 56. Templeton AR, Crandall KA, Sing CF (1992) A cladistic analysis of phenotypic associations with haplotypes inferred from restriction endonuclease mapping and DNA sequence data. III. Cladogram estimation. Genetics 132: 619–633.
            • 57. Salzburger W, Ewing GB, von Haeseler A (2011) The performance of phylogenetic algorithms in estimating haplotype genealogies. Mol Ecol 20: 1952–1963.
            • 58. Pfenninger M, Posada D (2002) Phylogeographic history of the land snail Candidula unifasciata (Poiret 1801) (Helicellinae, Stylommatophora): fragmentation, corridor migration and secondary contact. Evolution 56: 1776–1788.
            • 59. Dupanloup I, Schneider S, Excoffier L (2002) A simulated annealing approach to define the genetic structure of populations. Mol Ecol 11(12): 2571–81.
            • 60. Rogers AR, Harpending H (1992) Population growth makes waves in the distribution of pairwise genetic differences. Mol Biol Evol 9: 552–569.
            • 61. Fu YX (1997) Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics 147: 915–925.
            • 62. Drummond AJ, Rambaut A, Shapiro B, Pybus OG (2005) Bayesian coalescent inference of past population dynamics from molecular sequences. Mol Biol Evol 22: 1185–1192.
            • 63. Drummond A, Rambaut A (2007) Beast: Bayesian evolutionary analysis by sampling trees. BMC Evol Biol 7: 214.
            • 64. Brower AV (1994) Rapid morphological radiation and convergence among races of the butterfly Heliconius erato inferred from patterns of mitochondrial DNA evolution. Proc Ntl Acad Sci USA 91: 6491–6495.
            • 65. Rambaut A, Drummond AJ (2007) Tracer v1.4, Available: http://beast.bio.ed.ac.uk/Tracer.
              • 66. Phillips SJ, Anderson RP, Schapire RE (2006) Maximum entropy modelling of species geographic distributions. Ecol Modell 190: 231–259.
              • 67. Otto-Bliesner B, Brady E, Clauzet G, Tomas R, Levis S, et al. (2006) Last glacial maximum and Holocene climate in CCSM3. J Climate 19: 2526–2544.
              • 68. Hasumi H, Emori S (2004) K-1 Coupled Model (MIROC) Description. K-1 technical report, Center for Climate System Research, University of Tokyo, 34 pp.
                • 69. Warren DL, Seifert S (2011) Environmental niche modeling in Maxent: the importance of model complexity and the performance of model selection criteria. Ecol Applic 21: 335–342.
                • 70. Warren DL, Glor RE, Turelli M (2010) ENMTools: a toolbox for comparative studies of environmental niche models. Ecography 33(3): 607–611.
                • 71. Fielding AH, Bell JF (1997) A review of methods for the assessment of prediction errors in conservation presence/absence models. Environ Conserv 24: 38–49.
                • 72. Hijmans RJ, Guarino L, Jarvis A, O’Brien R, Mathur P, et al.. (2005). DIVA-GIS. Available: www.diva-gis.org/.
                  • 73. Rasband WS (1997–2008) ImageJ, U. S. National Institutes of Health, Bethesda, Maryland, USA. Available: http://rsb.info.nih.gov/ij/.
                    • 74. Armbruster P, Damsky Jr WE, Giordano R, Birungi J, Munstermann LE, et al. (2003) Infection of New and Old-World Aedes albopictus (Diptera: Culicidae) by the intracellular parasite Wolbachia: implications for host mitochondrial DNA evolution. J Med Entomol 40: 356–360.
                    • 75. Mousson L, Dauga C, Garrigues T, Schaffner F, Vazeille M, et al. (2005) Phylogeography of Aedes (Stegomyia) aegypti (L.) and Aedes (Stegomyia) albopictus (Skuse) (Diptera: Culicidae) based on mitochoindrial DNA variations. Genet Res 86: 1–11.
                    • 76. Chen H, Chen HB (2003) Sequences analysis of cytochrome C oxidase subunit I gene in Aedes albopictus from different geographic strains in China. Zhonghua Liu Xing Bing Xue Za Zhi 24: 491–493.
                    • 77. Maia RT, Scarpassa VM, Maciel-Litaiff LH, Tadei WP (2009) Reduced levels of genetic variation in Aedes albopictus (Diptera: Culicidae) from Manaus, Amazonas State, Brazil, based on analysis of the mitochondrial DNA ND5 gene. Genet Mol Res 8(3): 998–1007.
                    • 78. Usmani-Brown S, Cohnstaedt L, Munstermann LE (2009) Population genetics of Aedes albopictus (Diptera: Culicidae) invading populations, using mitochondrial nicotinamide adenine dinucleotide dehydrogenase subunit 5 sequences. Ann Entomol Soc Am 102(1): 144–150.
                    • 79. Kamgang B, Brengues C, Fontenille D, Njiokou F, Simard F, et al. (2011) Genetic Structure of the Tiger Mosquito, Aedes albopictus, in Cameroon (Central Africa). PLoS ONE 6(5): e20257.
                    • 80. Urbanelli S, Calvitti M, Porretta D (2005) Infezione del parassita intracellulare Wolbachia in Aedes albopictus: influenza sull’evoluzione del mtDNA dell’ospite. XIV Annual Meeting of the Italian Society of Ecology Torino 12–14 September.
                      • 81. Chen B, Pedro PM, Harbach RE, Somboon P, Walton C, et al. (2011) Mitochondrial DNA variation in malaria vector Anopheles minimus across China, Thailand and Vietnam: evolutionary hypothesis, population structure and population history. Heredity 106: 241–252.
                      • 82. Hurst GDD, Jiggins FM (2005) Problems with mitochondrial DNA as a marker in population, phylogeographic and phylogenetic studies: the effects of inherited symbionts. Proc Roy Soc Lond B Biol Sci 272: 1525–1534.
                      • 83. Urbanelli S, Bellini R, Carrieri M, Sallicandro P, Celli G (2000) Population structure of Aedes albopictus (Skuse): the mosquito which is colonizing Mediterranean countries. Heredity 84: 331–337.
                      • 84. Voris HK (2000) Maps of Pleistocene sea levels in Southeast Asia: shorelines, river systems and time durations. J Biogeogr 27: 1153–1167.
                      • 85. Bellwood P (1999) The changing nature of the Southeast Asian Environment. The Cambridge history of Southeast Asia vol. 1 (ed. by N. Tarling), 61–65. Cambridge University Press, UK.
                        • 86. Bisconti R, Canestrelli D, Colangelo P, Nascetti G (2011) Multiple lines of evidence for demographic and range expansion of a temperate species (Hyla sarda) during the last glaciation. Mol Ecol 20: 5313–5327.
                        • 87. Pope KO, Terrell JE (2008) Environmental setting of human migrations in the circum-Pacific region. J Biogeogr 35: 1–21.
                        • 88. Innes JB, Zong Y, Chen Z, Chen C, Wang Z, et al. (2009) Environmental history, palaeoecology and human activity at the early Neolithic forager/cultivator site at Kuahuqiao, Hangzhou, eastern China. Quat Sci Rev 28: 2277–2294.
                        • 89. Jones MK, Liu X (2009) Origins of agriculture in East Asia. Science 324: 730–731.
                        • 90. Zarowiecki M, Walton C, Torres E, McAlister E, Htun PT, et al. (2011) Pleistocene genetic connectivity in a widespread, open-habitat-adapted mosquito in the Indo-Oriental region. J Biogeogr 38: 1422–1432.
                        • 91. Porretta D, Bellini R, Sabatelli S, Samboon P, Urbanelli S (2011b) 59TH Annual meeting of the Entomological Society of America, Rino, NV, USA, November 13–16.
                          • 92. Porretta D, Mastrantonio V, Epis S, Bandi C, Sassera D, et al.. (2010) Connectivity among European Ixodes ricinus populations and its implication for pathogens diffusion: a genetic population approach. XX Annual Meeting of the Italian Society of Ecology, Rome, 27–30 September.