Introduction

Understanding the genetic processes associated with species introduction and establishment in new environments is fundamental to infer the evolutionary mechanisms underlying the dynamics and history of invasive processes. But, when an invasive species is a threat to human health, such as a disease vector, knowledge of population dynamics and demography is useful to interpret the causes and predict the risks of disease outbreaks. Aedes albopictus (Skuse 1894), the Asian tiger mosquito, is an emblematic example of an arbovirus vector that has rapidly and successfully invaded much of the globe supported by its life history traits and high propagule pressure driven by human activities. This mosquito took four decades to invade the world, while its counterpart Aedes aegypti did so in four centuries1. From its home range in tropical Southeast Asia, where it was originally a zoophilic forest species2, Ae. albopictus spread first to the islands in the Indian and Pacific Oceans3 and, during the 1980s, rapidly extended its range across temperate regions in Europe, the Americas and Africa1,4. This mosquito is able to tolerate climate/environment interactions that differ from its home range5,6,7. Its ability to diapause during unfavourable seasons and to lay desiccation-resistant eggs has facilitated its expansion in temperate regions8,9,10,11. It is listed among the top hundred most dangerous invasive species identified12. Its spread has created public health concerns as experimental infections have shown the species’ competence for at least 20 arboviruses and it is considered the main vector for chikungunya (CHIKV), and to a lesser degree for dengue (DENV) and Zika (ZIKV) viruses13,14,15,16,17,18,19,20,21. Anthropogenic activities, by creating new breeding and trophic niches of adaptation close to human habitations22,23, strengthened the association of Ae. albopictus with humans, with important public health implications. Indeed, beyond its original range, it has been increasingly involved in local autochthonous transmission of chikungunya in many places where it has become established, including La Réunion, continental Europe, Africa, the Americas and Japan, determining major CHIKV epidemics of severe, persistent, debilitating arthralgia24,25,26,27. Chikungunya virus (CHIKV, genus Alphavirus in the family Togaviridae), first discovered in 1952 in Tanzania28 and generally transmitted by Aedes mosquitoes, comprises three genetic lineages: Asian, West African and East/Central/South African (ECSA) which apparently evolved independently in different geographic areas29. The fact that Ae. albopictus is becoming the major vector of CHIKV is due to a series of adaptive mutations in the CHIKV ECSA lineage, the most known being an alanine to valine substitution at position 226 of the E1 glycoprotein (ECSA E1-A226V). This adaptive mutation, established in the Indian Ocean region in 2004–2005, dramatically increased the vector competence of Ae. albopictus for CHIKV, with respect to the previously more typical vector, Ae. aegypti. In association with the high invasive potential of Ae. albopictus, the ECSA E1-226V variant rapidly spread even in temperate areas such as Europe, causing CHIKV outbreaks.

The identification of the main dispersal routes of this mosquito vector can help in mitigating its spread and preventing/reducing the risks of arboviral outbreaks as once a vector is introduced, autochthonous transmission of its associated arboviruses may rapidly follow. But this implies knowledge of the demographic history of invading populations, their degree of connectivity and their vectorial capacities7. Several genetic, ecological and behavioural studies were performed to unravel different aspects of the invasion process and to infer the degree of genetic connectivity among populations at the micro- and macro-geographic levels23,30,31,32,33,34,35,36. However, vector competence for CHIKV has been found to be highly variable within and between populations, with transmission efficiencies that are strictly dependent on specific combinations of mosquito genome, viral genetic characteristics and temperature in a genotype-by-genotype-by-environment interaction37,38,39. On this basis, we previously proposed that admixture events and dynamics of different genetic backgrounds during the invasion process of this mosquito may have impacted competence for CHIKV in adventive populations32. In this study, we tested this hypothesis. We first provided a picture of the major events that determined, at the macrogeographic level, the demographic history of the adventive populations out of Southeast Asia starting from the 18th century. Secondly, we notably found that the population demographic history impacted their vector competence for CHIKV. We have demonstrated that the demographic history of populations was tightly related to CHIKV genotypes in a sort of adaptive coevolution, where genotype-by-genotype interactions influenced the corresponding Ae. albopictus vector competence phenotypes for this virus.

Results

Proximity does not impact population genetic connectivity

AMOVA analyses of molecular variance across mosquitoes from 25 localities spanning from Asia to Mediterranean area, the Americas and Central Africa (Table 1), indicated that 83% of variance occurred within populations, while only 11% and 6% were detected between populations and regions, respectively. The highest level of within-population variability was detected in the ancestral Southeast Asia region, where China displays the highest number of alleles and private alleles/individual. Outside this area, the adventive populations displayed different degrees of variability erosion (Supplementary Table 1). No signs of isolation by distance (IBD) were detected globally (y = 1e−06× + 0.0805; R2 = 0.015) nor within regions (R2 = 0.023 in South America, R2 = 0.15 in the Mediterranean area). Moreover, the population differentiation has been found to be heterogeneous not only between but also within geographical regions (Fig. 1, Supplementary Table 2). The most homogeneous area, in terms of non-significant or very low FST values, was the ancestral Southeast Asia (FST 0.014 – 0.069), which maintains low levels of differentiation with La Réunion and the Mediterranean Basin (FST 0.041 – 0.140), and with North America (FST 0.056 – 0.124). The Mediterranean area appears to be heterogeneous (FST 0.081 – 0.186). Within North America, Florida (VRB) was the least differentiated population (FST 0.065 − 0.070), sharing connectivity with the populations from Asia, La Réunion, and the Mediterranean area. It was the only population that maintained genetic connectivity with Central and South America. These latter two regions were the most heterogeneous based on FST values. The Brazilian and Argentinian populations had no genetic affinities with those from Southeast Asia, the Mediterranean area or North America, with the exception of Manaus (MAN) in the Brazilian State of Amazonas that, in turn, maintained low levels of genetic differentiation with the populations from Congo in Africa.

Table 1 Aedes albopictus populations sampled in areas in which the presence of this mosquito has been historically documented.
Fig. 1: Matrix of pairwise FST values among 25 Ae. albopictus population samples.
figure 1

FST values highlighted with a white circle are not significantly different from zero (P > 0.05) after Bonferroni correction.

Different lineages are present across the invaded areas

Six clusters/lineages represent the most parsimonious partitioning of the ancestry among individuals from the 25 populations (K1–K6 in Table 2) displaying a clear geographic distribution (Fig. 2, Supplementary Data 1). The genomes from the ancestral Asian populations, China and Japan, were distributed differently across clusters, but display the highest ancestry in K6. The other Asian population, Thailand, was fragmented mainly between K4 (55%) and K6 (24%). The ancestry heterogeneity of La Réunion and Mediterranean populations was obvious. The La Réunion mosquitoes from the southwest coast (RE) were mainly structured between K6 and K4 together with the majority of Mediterranean samples. By contrast, the two La Réunion populations from the north-eastern part (PROV and STD) were clustered with the French Mediterranean sample (BL) in K2. Hawaii was also fragmented between K6 and K4. North America appeared as a transition area in which populations such as Virginia (VA) still maintained membership in the two Asian clusters K4 (63%) and K6 (23%), while others such as Florida (VRB) and Missouri (TYS) exhibited ancestry profiles that transitioned to those typical of Central America (K3) and South America (K1). The Brazilian Amazonian (MAN) was fragmented across K1 (24%), K3 (17%) and notably the Southeast Asian K5 (45%). It is noteworthy that high percentages of K5 outside of Asia were recovered only in Manaus and in Congo, Africa (79%).

Table 2 Average coefficient of ancestry obtained from a STRUCTURE run with K = 6 for 702 individuals of Ae. albopictus from 25 samples collected in eight different geographical areas.
Fig. 2: Representation of the co-ancestry of Ae. albopictus mosquitoes from the ancestral and derived regions.
figure 2

Dates of invasion in the different regions are shown in the lower part of the figure.

Dynamics of invasion at the macrogeographic level

Five separate sequential Bayesian ABC analyses (Supplementary Table 3, Fig. 3) were compared to acquire a macro geographic picture of the demographic history of invasive populations out of the Southeast Asian home range to La Réunion, to Mediterranean Basin, to Americas and Africa. In the ancestral area, China has been confirmed as the ancestral population from which Thailand diverged, while Japan emerged from an admixture between China and Thailand.

Fig. 3: Graphical representation of the most-likely scenario of each set of scenarios describing the dynamics of samples within the native and derived areas using ABC methods.
figure 3

a Analysis 1: Southeast Asia, La Réunion and the Mediterranean invasion. b Analysis 2: The North America invasion. c Analysis 3: The South America Invasion. d The South America Invasion. e Analysis 4: Central American Invasion. f The Central Africa invasion.

Out of Southeast Asia to La Réunion and to Mediterranean basin (analysis 1b, Fig. 3a)

Mosquitoes from the southwest coast of La Réunion, St Pierre (RE), have previously been suggested as having an East Asian derivation32. We assessed whether the colonization of southern and northern parts of the island occurred independently as direct events from external South Asian source areas or via secondary internal derivation events (analyses 1a). The direct derivation of the northern PROV from the south-western RE was the most supported hypothesis. The STD population appeared to have a similar derivation. Manni and colleagues32 suggested that La Réunion/St Pierre (RE) contributed to the origin of the North Italian IT1 population (Cesena) through an admixture event with another North Italian population IT2, (previously derived from the North American Virginia (VA) population). On this basis, we tested whether La Réunion contributed also to the establishment of another Mediterranean population, Bar-sur-Loup (BL) in southern France. The hypothesis that La Réunion (PROV) in admixture events (74%) with the Italian IT1 (26%) contributed to the French Bar-sur-Loup population was strongly supported.

North America invasion

The South Atlantic states such as Florida and Virginia and the Midwestern state of Missouri are considered the first states to be invaded by Ae. albopictus according to historical records (Table 1). Thus, Manassas (VA) in Virginia, Vero Beach (VRB) in Florida and St Louis (TYS) in Missouri were considered as representative populations of these first invasions. After reconfirming that VA was derived from Japan (JP), the origin of Florida/VRB was determined. The most statistically supported scenario indicated that VRB in Florida was also a derivation of JP (analysis 2a, Fig. 3b). As for St Louis (TYS) in Missouri (analysis 2b), the strongest scenario supported its derivation from an admixture between JP and Florida (VRB) (66% and 34% respectively).

South America

The analysis of this area was considered before that of Central America, where the documented presence of this mosquito was more recent (Table 1). The tested hypothesis (analyses 3a–d, Fig. 3c, d) was whether the demographic histories of South American populations were related to those of North America, considering the concomitant historical presence of the species. Indeed, the most probable hypothesis indicated that the Brazilian coastal Parnamirim population (PNM) originated from Florida (VRB). Subsequent admixture events between Florida (VRB) and Parnamirim (PNM) (58% and 42% respectively) gave origin to the Rio de Janeiro state/Jurujuba (JRB) population. Among the Brazilian populations, the Amazonas State sample from Manaus (MAN) appeared to be the most differentiated. Its most likely origin (analysis 3d) was an admixture event between VRB-Florida and JP (52% and 48%, respectively). In Argentina, the mosquito has been recorded in Misiones Province near the border with Brazil. Thus, we tested whether the Misiones sample (MIA) was also a derivation from the Brazilian and/or Florida populations. The Argentinian mosquitoes appear to have originated as a split from the Brazilian JRB population (analysis 3c, Fig. 3c).

Central America

The relatively recent demographic histories of the populations Tapachula in Mexico (MXC) and Colòn in Panama (PAN) were considered under the hypothesis of a North American and/or a Brazilian derivation (analyses 4a, b, Fig. 3e). The most likely origin of MXC was a split from Florida-VRB. Regarding PAN, the most statistically supported scenario suggests that it originated from an admixture between VRB in Florida and Brazilian PNM.

Congo in Central Africa

The origin of the Brazzaville population in Congo was tested in analysis 5 (Fig. 3f). A very high posterior probability suggests that this central African population split from the Brazilian population of MAN (Manaus).

Population ancestry affects vector competence for CHIKV

Based on the above results, we selected 13 Ae. albopictus populations for which vector competence data obtained from experimental CHIKV infections were available (Table 1, Supplementary Table 4) to explore whether, and to what extent, the degree of ancestry detected among the populations impacted their dissemination efficiency (DE) and transmission efficiency (TE) for CHIKV. DE and TE respectively refer to the proportion of mosquitoes with infectious viral particles in the head and in the saliva among the number of mosquitoes tested.

A total of 1290 mosquitoes were considered in the dissemination analysis of which 1120 (86.8%) showed disseminated infection. CHIKV genotypes, dpi and blood titre were significantly associated with DE in univariate analysis (Table 3). When populations were grouped according to their major K-ancestry (as listed in Table 2), DE significantly differed according to K-ancestry in univariate analysis (P = 0.0001, Table 3). This remained true after adjusting on CHIKV genotypes, dpi and blood meal titre (Supplementary Table 5, Fig. 4a). As compared to mosquitoes with main K3-ancestry (North and Central America: OR = 1), those with major K6 (IT1, Italy), K2 (BL-France), K5 (Brazilian MAN and CONG1) and K1 (South America) ancestries showed significantly higher dissemination efficiencies (Fig. 4a). Dissemination of CHIKV strains from Asian genotype (CHIKV-ASIA) and, to a lesser degree, those of CHIKV ECSA 226 V strains significantly decreased across the K-ancestry lineages: from K6 (IT1, Italy) to the derived K2 (BL, France), to K5 (Manaus, Congo) to K1 (South America) and K3 (North and Central America) which displayed similar efficiencies. Another approach consisted in introducing all 6 K-ancestries when dichotomised. In univariate analysis, levels of K1, K5 and K6-ancestries >0.25 were significantly associated with DE. It appeared that these three dichotomised factors defined only four distinct population profiles (low K1, K5 and K6-ancestries, high K1-ancestry only, high K5-ancestry only and high K6-ancestry only); to simplify model interpretation, a new variable was defined allowing comparison of these 4 populations (Table 3). After adjusting for CHIKV genotype, dpi and blood meal titre, as compared to the populations with lower K5 and K6-ancestries, those with higher K5 or K6 presented higher DE (Table 3). Both multivariate models showed a tendency for DE to decrease from populations with major K6-ancestry to major K5-ancestry and finally to the other populations.

Table 3 Dissemination efficiency of CHIKV in Ae. albopictus lineages (Adjusted logistic regression).
Fig. 4: Graphical representation of the influence of K-ancestry on dissemination and transmission efficiencies.
figure 4

a Dissemination efficiency (n = 443, 13, 498, 0, 256 and 80 for K1–K6, respectively, overall P = 0.003); and b transmission efficiency (n = 401, 52, 424, 0, 203 and 79 for K1–K6, respectively, overall P = 0.0006). Circles represent the adjusted Odds Ratio estimate and the bar its 95% confidence interval obtained from the logistic regression presented in Tables 3 and 4, respectively.

Transmission efficiency (TE) for CHIKV was evaluated in 1159 mosquitoes of which 58.0% (672) could transmit the virus. CHIKV genotypes were significantly associated with TE, with ECSA E1-226V strains being more efficient (P = 0.0002) and ECSA (E1-226A) at a lower degree, while no association with TE was detected according to dpi and blood meal titres (Supplementary Table 6). Again, the K-ancestry lineages displayed different TEs (Supplementary Table 6). In univariate analysis, as compared to populations with major K3-ancestry, all other populations presented greater TEs. This remained true after adjusting for CHIKV genotype, as populations with major K2-ancestry (La Réunion, PROV and the derived French BL) displayed the greatest TEs (Table 4 and Fig. 4b). As for DE, another approach consisted in looking at the level of each six K-ancestry when dichotomised. In univariate analysis, the levels of K1, K2 and K3-ancestries were all significantly associated with TE (Supplementary Table 6). In multivariate analysis, a higher level of K3-ancestry was associated with lower TE, whereas higher levels of K2 were associated with a greater TE (Table 4). Both multivariate models showed a tendency for TE to reduce from K2 to K1 and finally K3.

Table 4 Transmission efficiency in Ae. albopictus lineages (Adjusted logistic regression).

Also, on the basis of population profiles, when compared in mosquitoes with low K2 and K3-ancestries, transmission was significantly higher in those with high K2-ancestry and lower in those with higher K3-ancestry (Supplementary Table 6).

Discussion

This study provides a picture of the major events which determined, at the macrogeographic level, the demographic history of the adventive populations established in the Indian Ocean island of La Réunion, in the Mediterranean area, in the Americas and in Central Africa, as a consequence of the rapid west-oriented invasive processes of Ae. albopictus out of South East Asia starting from the 18th century. We notably found that (1) admixture events had a major role in shaping the genetic makeup of the globally distributed adventive populations, impacting their demographic histories, and (2) the population demographic history impacted their ability to transmit CHIKV. It follows that knowledge of the relationship between demographic history and vector competence is crucial for assessing the risk of arbovirus outbreaks in areas where this mosquito has become established.

As far as it concerns the demographic history of this mosquito, six gene pools/lineages (clusters) were detected across the invaded regions at the macrogeographic level. The data clearly suggest that the invasion process of this mosquito fits a demographic model of historical lineage diversification through admixture events among different genetic compositions as summarized in Fig. 5. The lineage divergence times we inferred from genetic data (ABC, STRUCTURE) are congruent with the available historical documentation on the different steps of the invasion process, thus providing a picture of the history of this mosquito outside its home range over the last 200 years. In Southeast Asia, a high level of lineage-sharing is displayed by the ancestral Chinese and, to a lesser degree, by the derived Thai and Japanese populations, supporting high coancestry in this home range area (Fig. 2). The historical trade networks across East Asian maritime space, the ‘China seas’, may have promoted the diffusion of propagules to Thailand and Japan resulting in an increase in diversity and favouring expansion and adaptation. As illustrated in Fig. 5, we inferred that K6 (red) and K4 (green) were the first ancestral lineages that historically emerged from Southeast Asia. K6, highly represented in China, is the lineage that may correspond to the primary invasion route within East Asia (Thailand and Japan) and, out of this area, to Indian Ocean La Réunion in the 18th century, to Hawaii in 1895, to the Mediterranean area starting from 1979 and to the USA in 198540. K4 is highly present in Thailand and contributed to the ancient colonization of the west coast of La Réunion, to the USA, and to the recent outbreaks in the Mediterranean area, Greece and partially to Cesena (IT1)32. A spatial extension of K6 and K4 from the western side of La Réunion (St Pierre, RE) to the eastern part of the island (PROV and STD) gave rise to the K2 lineage diversification (orange). Habitat structuring and strong climatic differences between the east and western coasts of La Réunion may have supported this diversification41,42. K2 expanded its range out of La Réunion and became established in the French Mediterranean area of Bar Sur Loup (BL) in 2004. La Réunion is a major trading partner of Continental France and most of this international trade passes through the main port and the international airport at St Denis on the north-east coast where K2 was established (PROV and STD). K6 and K4 were almost lost in the invasion process of the Americas. In North America, they contributed to the diversification of the K3 lineage (grey in Fig. 5) which is the major component of the Ae. albopictus genomes in Florida and Missouri. Indeed, within a year of its introduction in Texas in 198543, it was found in Florida and then rapidly expanded throughout the south-east coast and part of north-central and north-east USA due to internal trade and multiple introductions from Japan/South Asia and Hawaii. According to our data and historical records, Florida played a key role in introducing this mosquito in South America in 1986 and later, in Central America, Mexico and Panama in 2002 and 2004, respectively44,45. Thus, Florida can be considered as a successful bridgehead invasive population which favoured the extension of K3 (grey in Fig. 5) in Central America and the divergence of K1 (cyan) in Brazil with its subsequent diffusion within that country and in Argentina. In South America, the Rio Grande do Norte (PNM) and Rio de Janeiro (JRB) states appear to have played a role in the spread of the K1 lineage. A particular case is provided by the genetic signature of Brazilian mosquitoes present in the Central Amazonian state (Manaus) since 2002. Unlike the other Brazilian mosquitoes, they display a higher proportion of the ancient Southeast Asian K5 lineage (0.450), with respect to the K1 Brazilian component which is unexpectedly low (0.240) and that of North Central America, K3, is even lower (0.167). The high K5 lineage component may be explained by the contribution that Southeast Asia/Japan provided to the origin of Manaus mosquitoes46. Manaus is one of the four free trade areas in Brazil and its port is a commercial centre for ocean-going shipping and the main transport hub for the Amazon basin. We found the K5 lineage well established in mosquitoes from Congo in Central Africa as a derivation from Amazonian Manaus. How this area could be implicated in the introduction of this mosquito to Congo is an open question although it is noteworthy that this inference was derived from a sample collected in 2011; the year of the first record of this mosquito in this area. Several authors have also proposed that propagules originating from South America could have contributed to an expansion in Africa.

Fig. 5: Geographical representation of the global spread of Ae. albopictus highlighting the co-ancestry of the different derived populations and the dissemination and transmission efficiencies for Chikungunya virus (CHIKV).
figure 5

ECSA E1-226V: CHIKV strains from East-Central-South African genotype harbouring a valine at position 226 of E1 glycoprotein. ECSA: CHIKV strains from East-Central-South African genotype harbouring an alanine at position 226 of E1 glycoprotein. ASIA: CHIKV strains from Asian genotype.

Finally, as summarised in Fig. 5, the old colonized Mediterranean area appears to be a mixing pot where the introduction of different lineages or their components resulted in the mixing of different genomes with the consequent heterogeneity observed among the populations.

Geographic variation in vector competence for different CHIKV genotypes has been observed among Ae. albopictus populations, highlighting that distinct barriers to dissemination and to transmission may operate even at small geographic scales depending on the interactions between the vector’s genetic background and the viral genotype15,37,39,47. Our data demonstrated that the ability of an Ae. albopictus population to disseminate or transmit the different CHIKV genotypes is strongly associated with their degree of ancestry within a genetic lineage. It follows that if populations from different geographical areas are related by their demographic history notably in terms of ancestry, they will display similar levels of vector competence for a given CHIKV strain. This is true for the La Réunion PROV, STD and the French Mediterranean BL populations in K2; the Brazilian Amazonian Manaus and Congo populations in K5, and the North and Central American populations in K3 (Tables 3 and 4). On the other hand, this phenomenon generates variation in competence within regions in which different lineages or admixture of lineages coexist, like in the Mediterranean area.

If we consider the historical dynamics of Ae. albopictus expansion and the temporal and spatial history of CHIKV evolution and spread, it appears that they are bound by convergent adaptive evolution39,48,49. Indeed, the historical sequential divergence and admixture of Ae. albopictus lineages out of Southeast Asia determined changes in the level of populations’ competence for CHIKV, thus facilitating the establishment of CHIKV transmission cycles in new regions. All the mosquitoes from the genetic lineages examined were able to disseminate the Asian and ECSA E1-226V CHIKV genotypes, better than ECSA (E1-226A). However, there is a progressive and drastic loss of permissiveness of the mosquito barriers for CHIKV dissemination (i.e. antiviral immune responses at the midgut epithelium) regarding the Asian and ECSA E1-226V genotypes that parallels the historical lineage divergence from the ancestral Asian lineage K6 to K2, to K5 and finally, to the American K1 and K3 lineages. Indeed, the higher abilities to disseminate CHIKV Asian and ECSA E1-226V genotypes were observed in populations with higher levels of the ancestral K6 and K5 lineages (Fig. 5), suggesting that the loss or modification of ancestral components in the genome of invading populations resulted in enhanced tissue barriers to CHIKV dissemination in mosquitoes. In contrast to their dissemination efficiency, the lineages displayed reduced transmission efficiency for CHIKV-ASIA genotype with respect to both ECSA genotypes, especially ECSA E1-226V. These results can be explained by the constrained ability of Asian CHIKV strains to adapt to Ae. albopictus due to negative epistatic interactions between the E1 residue 98 (that harbours a threonine in Asian strains; E1-98T) and the E1-226 position49. The highest transmission efficiency was displayed by La Réunion K2 lineage followed by K6 which retains around two-thirds of transmission efficiency. Interestingly, the K2 lineage differentiated in La Réunion where the CHIKV E1-226V genotype was identified for the first time in 2005 and established as an adaptive mutation of ECSA lineage for this mosquito50,51. ECSA E1-226V, with a valine instead of an alanine at position 226 A of the E1 glycoprotein, significantly increased viral replication in midgut infection and transmission. The association of this mutation with Ae. albopictus invasiveness impacted CHIKV diffusion in a sort of convergent adaptive evolution between vector populations and viral genomes49. Moreover, the ECSA E1-226V mutation is thought to have arisen independently in other areas where Ae. albopictus is present, supporting this convergent evolution hypothesis49,52. We have found that the higher the level of the K2 lineage in a population is, the higher is its ability to transmit the ECSA E1-226V genotype. After 2007, the ECSA E1-226V genotype was identified in Southeast Asia and quickly spread, supported by the transmission efficiency of the Ae. albopictus ancestral lineages (K6 and K5), inducing major chikungunya epidemics53. Prior to 2007, Ae. albopictus had never been incriminated in outbreaks involving Asian CHIKV strains, despite the massive presence of the mosquito in this region. These outbreaks were vectored by Ae. aegypti which is also abundant in the area. One plausible reason would be that, as explained above, Asian strains of CHIKV were constrained in their ability to adapt to Ae. albopictus via acquisition of the E1-A226V mutation49. This inability could have facilitated the establishment of the Ae. albopictus-adapted CHIKV strains that harbour the E1-226V substitution in the region, leading to the ongoing shift from Asian to ECSA E1-A226V genotypes in Southeast Asia49.

In the Mediterranean area, no CHIKV Asian genotype has been related to autochthonous transmission despite the intense influx of travellers returning from Southeast Asia and the Americas. It has been shown that exposure of European Ae. albopictus to low temperatures (20 °C) significantly increased the extrinsic incubation period and reduced the transmission efficiency of Asian CHIKV strains, while for strains from ECSA lineages these parameters remained unaffected15,54. Although other factors influencing vectorial capacity may interact and determine that an infectious bite is received by a susceptible human host7, these results firstly highlight the importance of three-way interactions between temperature, the mosquito population and the viral genotype. Secondly, they suggest a higher risk for ECSA strains to emerge in Europe. In addition, the presence in this old-colonized area of different Ae. albopictus lineages that are highly competent for ECSA genotypes (Supplementary Table 5; Supplementary Table 6), coupled to a growing number of imported cases55 led to autochthonous chikungunya transmission in France56,57,58 and Italy in 2007 and 201759,60,61,62. In all these autochthonous transmissions, the strains involved belonged exclusively to the ECSA lineage, which supports the CHIKV emergence hypothesis enounced above. Nevertheless, the heterogeneity of Ae. albopictus genomes in the Mediterranean region has resulted in a patchwork of populations displaying different levels of competence which makes risk assessment for CHIKV highly unpredictable at a finer scale in the area63. Beyond vector competence, context-specific differences in parameters with greater impact on vectorial capacity (i.e. longevity, density, host feeding preference)7 may challenge risk prediction for arboviruses and should be taken into account in further studies.

The most recent divergent American lineages, K1 (South America) and especially, the North and Central American lineage K3, displayed the lowest CHIKV dissemination and transmission efficiencies of all studied populations. Indeed, K3 is a heterogeneous lineage with USA populations such as Florida/Vero Beach that displayed lower CHIKV dissemination and transmission efficiencies (Supplementary Table 5; Supplementary Table 6). Vega-Rúa and colleagues37 and Honório and colleagues47 using the three CHIKV strains of this study (ASIA, ECSA E1-226A and ECSA E1-226V) performed comparative analyses of vector competence using Ae. albopictus and Ae. aegypti populations across the Americas. All the populations of both species were susceptible to infection and dissemination of the different CHIKV genotypes. However, the transmission rates were heterogeneous even among populations from the same species with Ae. albopictus from Brazil, and generally from South America, having higher transmission efficiency for ECSA genotypes when compared to their counterparts from Florida, USA. The first detection of CHIKV in the Americas was recorded in the Caribbean in late 2013, followed by a dramatic spread of the virus throughout the continent. It has been stated that CHIKV arrived in Brazil through two independent introductions: the Asian genotype from the Caribbean entered through the northern region, while the African ECSA genotype was imported through the north-east region64. Following their initial introduction, both genotypes established and expanded their range with CHIKV-ASIA being the prevalent genotype causing several outbreaks in relation to the presence of Ae. aegypti that was the main vector implicated in the transmission15,37. In 2014, local transmission of CHIKV-ASIA was detected in the Amazon region. A distinct lineage CHIKV-ECSA was detected simultaneously and it rapidly spread inducing replacement of the Asian genotype in the Brazilian Amazonian region65. Whether this replacement was due to the spread of Ae. albopictus in this region and its transmission capacity for these particular ECSA strains, which do not harbour the E1-226V mutation, is an open question that deserves further investigation.

The ancestral K5 lineage is differentiated in the African populations from Congo that became established in 2011 as a derivation of the K5 Amazonian component. Vazeille and colleagues in 201638 demonstrated that Ae. albopictus newly introduced in Congo were unable to transmit the ECSA E1-226V genotype with the same efficiency as La Réunion mosquitoes, highlighting that genetic differentiation among these two populations may have impacted their vector competence. Our results are in agreement with this hypothesis: Congo and La Réunion belong to two well differentiated lineages K5 and K2, respectively, with Congo having half the ECSA E1-226V transmission efficiency when compared to the La Réunion K2 lineage.

In conclusion, it is clear from these data that the demographic history of Ae. albopictus populations is tightly related to CHIKV genotypes in a sort of adaptive coevolution, where genotype-by-genotype interactions influenced the corresponding Ae. albopictus vector competence phenotypes for this virus. It follows that the knowledge of the vector demographic history dynamics combined with vector competence data can provide an accurate risk map for CHIKV. However, in addition to these genetic interactions, we must consider that the spread of this mosquito occurred across different environmental conditions. Thus, the variation in competence between the eco-geographic populations may be impacted by the interactions between unexplored genomic factors66,67 within genomes and lineages, and complex environmental landscapes.

Methods

Mosquito collections

Twenty-five mosquito populations were considered. Three were from the ancestral East Asian area and 22 were sampled across much of the expansion area: La Réunion in the Indian Ocean, the Mediterranean basin, Hawaii in the Pacific Ocean, North, Central and South America, and Central Africa (Table 1). The samples from Southeast Asia (Japan, China, Thailand), La Réunion (St Pierre), the Mediterranean basin (Albania, Greece, Italy1, Italy2), Hawaii and North America (Virginia) were previously described by our laboratory in terms of microsatellite allele frequencies, genetic structure and demographic history32. Their integration in this analysis aimed at obtaining a deeper genetic portrait of the demographic history of this species during its invasion process. For the 25 collection sites, samples were collected as eggs at the peak of density to avoid any seasonality effects due to the presence of diapause or pre-diapause periods10. As described in Manni and colleagues32, a standardized sampling protocol was adopted to minimize inbreeding effects at each breeding site. In each locality, 15–17 ovitraps were placed at a distance of least 500 m one to another, and at least 40 eggs/ovitrap were obtained.

Adults from each field collection were identified using morphological keys68. The mosquito egg collections in the Americas, La Réunion, France, the North of Italy and in Congo/Central Africa were separated in two batches: one batch was used for population genetic studies and the remaining batch, that gave rise to several hundred individuals, was amplified in the laboratory under controlled conditions to obtain the offspring that was used in vector competence assessments for CHIKV (Table 1). Thus, the microsatellite genotyping analyses were based on an average of 30 G0 individuals for each population sample, while vector competence assays were performed on the corresponding G1 mosquito populations, with the exception of Brazaville (G2) and Cesena (G3).

Microsatellite analyses

We used previously validated SSRs as highly polymorphic markers to provide continuity with our published Ae. albopictus population data and because their use allowed accurate biogeographic data analyses especially at macrogeographic scale. Indeed, the availability of a population SSR genotype dataset was a reliable resource for ABC simulations as SSRs represented a good balance between being well characterized (i.e. in terms of mutation model) and having sufficient genetic diversity to provide statistical power, as we previously assessed31.

We used 11 highly informative SSR loci applied to 702 individuals from 25 localities to produce a SSR genotyping dataset with no major bias in statistical analyses69.

On this basis, genomic DNA extracted from each mosquito70 was genotyped at the following SSR loci: Aealbmic1, 2, 3, 5, 6, 9, 11, 14, 15, 16 and 1731. These loci were chosen due to their high polymorphism and because they are spread across the genome (they map to different scaffolds of the Ae. albopictus Foshan genome (ref. 71 and our unpublished data)). They have proven to be efficient markers to detect variability, even in relatively small samples32. PCR amplifications and fragment identifications were performed as previously described31. In order to account for genotyping errors, automated binning of allele lengths was performed with TANDEM v1.0972 followed by manual checking. When microsatellite amplification was not successful, or allele scoring was unclear, a new DNA extraction was performed. The generated SSR genotype dataset, in addition to offering continuity with previous studies, allows the formulation of reliable models of recent evolutionary history because of their well‐documented mutation rates73, high molecular diversity and high minor allele frequencies.

Variability and genetic diversity analyses

Linkage disequilibrium between pairs of loci in each sample (100 batches, 1000 iterations per batch) and deviation from within-population Hardy–Weinberg equilibrium (HWE) at each locus/sample combination were tested with GENEPOP V. 4.274. The statistical significance was assessed following Bonferroni corrections75. Variation within each locality was estimated in terms of average number of alleles (na), private alleles (np) and their respective frequencies (Ap) using GenAlEx 6.576. The average number of alleles (na/n) and private alleles (np/n) were also computed at the individual level. Gene diversity (Hs) and allelic richness (Rs) were computed with the program FSTAT V.2.9.3.277. Observed (Ho) and unbiased expected heterozygosity (uHe) values and Pairwise-FST were computed using Microsatellite Analyzer (MSA) V.4.0578. The statistical significance of each FST value was assessed by comparison of the assessed value with the value obtained using 10,000 matrix permutations and Bonferroni corrections.

Population structure

The structure and degree of ancestry in the 25 population samples were inferred using the Bayesian clustering analysis implemented in STRUCTURE V 2.3.2 under the admixture model and assuming independent allele frequencies. The “burn-in” phase was set to 500,000 iterations, then 1,000,000 Markov Chain Monte Carlo (MCMC) replications were considered. For each number of possible clusters (K), 20 independent runs were repeated. The cluster (K) interval was set between 1 and 25 (i.e. the number of considered populations).

The results were analysed using STRUCTURE HARVESTER79 and the most likely number of clusters was determined by plotting the log probability (L(K))80, as well as the ΔK over all the runs81. Once the best number of clusters, K, was inferred, CLUMPP software82 was used to merge the 20 runs and the results were plotted in DISTRUCT83.

Analyses of demographic History

The Approximate Bayesian Computation (ABC) method, as implemented in the software DIYABC v2.1, was used to disentangle the complex invasion patterns of Ae. albopictus at the macro-geographic level. This approach allows comparisons between competing hypotheses regarding the divergence of the populations on a global scale. DIYABC v2.1 implements a standard ABC method (ABC-LDA) which, although being more expensive computationally, was shown to perform comparably (for posterior probability) or slightly better (for prior error rate) than the more recent random forest method (ABC-RF)84 when large reference tables are used85, as in our case. Because of the huge number of possible evolutionary scenarios (25 factorial, i.e. 1.55e25) that could be tested with the 25 populations considered, we adopted a step-by-step approach considering subgroups of populations, starting from the ancestral Southeast Asian region to the invaded areas. Based on our published data, we considered the previously characterized samples from China, Thailand and Japan as representative populations from countries within the Southeast Asian native range.

Several competing invasion scenarios were considered taking into account the available historical information regarding the first records in invaded countries and population genetic data. The posterior probability for each scenario was computed. The competing scenarios were set using prior definitions and distribution of demographic parameters, as described in Supplementary Table 7. For some of these analyses, we referred to and reran the scenarios that we had previously validated32. We took into account the effective population size, the time range in which split or admixure events occurred, the number of founders that contributed to the establishment of adventive populations, the duration of the eventual bottlenecks occurring during colonization and the rate of admixture, if considered32. When no prior information was available, the prior parameters were kept deliberately broad. For the effective population sizes, broad priors [500–100,000] with uniform distributions were chosen for all the populations. For each event, the time range chosen was less broad if it was supported by historical data. The timing of events was expressed in numbers of generations back in time. Considering that the number of generations/year is dependent on bioclimatic conditions, we assumed around 4–7 generations/year in temperate regions, and 12–17 generations/year in tropical regions. In the cases where the climate situation is not so clear-cut, the broadest interval was considered (i.e. 4–17 generations/year).

As mentioned before32, although Ae. albopictus has great potential for rapid population growth86, uncertainty regarding the duration of the bottleneck was taken into account. Therefore, we assumed a bottleneck period with a uniform prior distribution bounded between 0 and 50 generations. During colonization, the number of founder individuals for each colonization event was described as NF and its size was drawn from a uniform distribution bounded between 1 and 100 individuals. If necessary, the rate of admixture between two populations was drawn in a uniform distribution and set from 0.001 to 0.999. For the genetic parameters, the microsatellite Generalized Stepwise Mutation model was considered, and a mean microsatellite mutation rate across loci was set between 10−5 and 10−3. The parameter of geometric distribution was specified with a uniform prior distribution bounded between 0.1–0.3. For each locus, the above considered mutation parameters were allowed to vary using Gamma distributions. Finally, the possibility of single nucleotide insertion/deletion mutations were considered with a mean frequency of 10−8–10−5 (Supplementary Table 7)32.

The genetic variation within and among the populations was summarized using sets of statistics conventionally used in ABC analyses: the mean number of alleles per locus, the mean expected heterozygosity, the mean allelic size variance, the overall pairwise FST values, and the Garza-Williamson index (number of alleles in a population divided by the range in allele size)87. This index helps to discriminate populations that have experienced recent losses in genetic variability from those that have been stable for a long time. One million simulated datasets were generated for each scenario. Competing scenarios were compared by calculating their posterior probabilities using a polychotomous logistic regression on the 1% of simulated data closest to the observed data88. Confidence in scenario choice was evaluated by computing type I and type II errors89. Once the most likely scenario was identified for each analysis, the posterior distributions of genetic and demographic parameters were estimated. This was achieved by computing a local linear regression on the 1% of the simulated data closest to our observed dataset, after applying a logit transformation to the parameters value. The goodness-of-fit of the estimation procedure was also evaluated by performing a model checking computation by generating 1000 pseudo-observed data sets with parameter values drawn from the posterior distribution given the most likely scenario32.

Vector competence data for CHIKV strains

Vector competence for CHIKV was assessed on adults derived from field-collected eggs from the Americas, La Réunion Island, France, North Italy and Congo as described above (Table 1). Batches of 5–7-day-old females were mono-infected using an artificial feeding system (Hemotek Ltd®, Blackburn, UK) with low passaged CHIKV strains (2–3 passages) belonging to either East-Central-South African genotype (ECSA) or the Asian genotype (see section Viral Strains for more details). Fully engorged females were kept for further estimation of DE, TE and their saliva was collected at the selected days post-infection (Supplementary Table 4) as described by Vega-Rúa and collegues37. The infectious status of head and saliva samples was determined by titration assays90. The blood meal titres, viral strains, days post-infection (dpi) as well as DE and TE data are summarized in Supplementary Table 4.

Viral strains

Six CHIKV strains isolated from human cases were used in infection assays. Five strains belonged to the ECSA genotype: CHIKV_0621 (Accession number AM258992) and CHIKV_05115 (Accession number AM258990) isolated in La Réunion, CHIKV_1909 (Accession number FR846305) from Southeast France, CHIKV_DRC and CHIKV_Congo from Congo38. All these ECSA strains were split into two different genotype groups according to the amino acid substitution at 226 position of E1 envelope protein: the group named ‘ECSA’ that contains the substitution E1-226A (Alanine) and the group named ‘ECSA E1-226V’ that contains a Valine at the same position. This mutation allows an increased transmission by more than 50-fold in Ae. albopictus50,51. The strain from Asian genotype used, CHIKV_NC, was isolated in New Caledonia in 2011.

Association between population ancestry and CHIKV competence

The effect of viral genotype, dpi, and population ancestry on DE and TE were investigated using logistic regression models. In case of complete prediction of the outcome in one category, a penalized likelihood estimation was considered91. For TE, only mosquitoes that had disseminated CHIKV infections were considered. The Odds Ratio (OR) has been considered as a measure of the risk of CHIKV dissemination and transmission. The main factor of interest was the degree of ancestry among populations (K-ancestry) and two approaches were used. First, each geographical population is defined by its major K-ancestry (Table 2). Second, all 6 K-ancestry components were considered but dichotomised as above or below the threshold of 0.25; this threshold was deliberately chosen so that some populations could be characterized by a mixture of several K-ancestors. For both dissemination and transmission, factors that presented P-values < 0.05 in univariate analysis were introduced in multivariate models. Two multivariate models were considered using either the major K-ancestry or the binary K-ancestry factors.

Statistics and reproducibility

The details regarding the statistical analyses, the tests, software used, sample sizes and number of replicates are described in the relevant sections of the “Methods”.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.