# Multiple invasions, Wolbachia and human-aided transport drive the genetic variability of Aedes albopictus in the Iberian Peninsula

Federica Lucati, Sarah Delacour, John R.B. Palmer, Jenny Caner, Aitana

Oltra, Claudia Paredes-Esquivel, Simone Mariani, Santi Escartin, David A

Roiz, Francisco Collantes, et al.

► To cite this version:

Federica Lucati, Sarah Delacour, John R.B. Palmer, Jenny Caner, Aitana Oltra, et al. Multiple invasions, Wolbachia and human-aided transport drive the genetic variability of Aedes albopictus in the Iberian Peninsula. Scientific Reports, 2022, 12 (1), 20682 [15 p.]. 10.1038/s41598-022-24963-3. hal-04977320

HAL Id: hal-04977320

https://hal.science/hal-04977320v1

Submitted on 5 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

OPEN

# Multiple invasions, *Wolbachia* and human-aided transport drive the genetic variability of *Aedes albopictus* in the Iberian Peninsula

Federica Lucati<fnref n="1" />,²✉, Sarah Delacour³, John R.B. Palmer², Jenny Caner<fnref n="1" />, Aitana Oltra<fnref n="1" />, Claudia Paredes-Esquivel⁴, Simone Mariani<fnref n="1" />, Santi Escartin<fnref n="1" />,⁵, David Roiz⁶, Francisco Collantes⁷, Mikel Bengoa⁸, Tomàs Montalvo⁹,¹⁰, Juan Antonio Delgado⁷, Roger Eritja<fnref n="1" />,¹¹, Javier Lucientes³, Andreu Albó Timor<fnref n="1" />, Frederic Bartumeus<fnref n="1" />,¹¹,¹² & Marc Ventura<fnref n="1" />

The Asian tiger mosquito, *Aedes albopictus*, is one of the most invasive species in the world. Native to the tropical forests of Southeast Asia, over the past 30 years it has rapidly spread throughout tropical and temperate regions of the world. Its dramatic expansion has resulted in public health concerns as a consequence of its vector competence for at least 16 viruses. Previous studies showed that *Ae. albopictus* spread has been facilitated by human-mediated transportation, but much remains unknown about how this has affected its genetic attributes. Here we examined the factors that contributed to shaping the current genetic constitution of *Ae. albopictus* in the Iberian Peninsula, where the species was first found in 2004, by combining population genetics and Bayesian modelling. We found that both mitochondrial and nuclear DNA markers showed a lack of genetic structure and the presence of worldwide dominant haplotypes, suggesting regular introductions from abroad. Mitochondrial DNA showed little genetic diversity compared to nuclear DNA, likely explained by infection with maternally transmitted bacteria of the genus *Wolbachia*. Multilevel models revealed that greater mosquito fluxes (estimated from commuting patterns and tiger mosquito population distribution) and spatial proximity between sampling sites were associated with lower nuclear genetic distance, suggesting that rapid short- and medium-distance dispersal is facilitated by humans through vehicular traffic. This study highlights the significant role of human transportation in shaping the genetic attributes of *Ae. albopictus* and promoting regional gene flow, and underscores the need for a territorially integrated surveillance across scales of this disease-carrying mosquito.

The Asian tiger mosquito, *Aedes* (Stegomyia) *albopictus* (Skuse 1894) is a highly invasive species that originated in the tropical forests of Southeast Asia<fnref n="1" />. However, in the late 1970s it started a dramatic expansion throughout tropical and temperate regions of the world and it is now present in the five populated continents². This species is an ecological generalist capable of rapid evolution and, with the aid of man, speedy colonization of new habitats<fnref n="1" />. While in its native range *Ae. albopictus* inhabits forested areas, breeding in natural sites such as tree holes, bromeliads and bamboo stumps, it has now adapted to breed also in artificial man-made water containers from urban and suburban human settlements³. This species has an opportunistic feeding behaviour with a strong preference for mammals, especially humans⁴,⁵. Furthermore, in temperate regions, its eggs can survive the cold winters by entering diapause³. This ecological adaptability has important implications for the

<footnote n="1">¹Center for Advanced Studies of Blanes (CEAB-CSIC), Blanes, Catalonia, Spain. ²Universitat Pompeu Fabra (UPF), Barcelona, Catalonia, Spain. ³School of Veterinary Medicine, The AgriFood Institute of Aragon (IA2), University of Zaragoza, Zaragoza, Spain. ⁴Applied Zoology and Animal Conservation Group, Universitat de les Illes Balears, Palma, Mallorca, Spain. ⁵Associació Mediambiental Xatrac, Lloret de Mar, Catalonia, Spain. ⁶MIVEGEC, Institut de Recherche pour le Développement (IRD), Montpellier, France. ⁷Department of Zoology and Physical Anthropology, University of Murcia, Murcia, Spain. ⁸Anticimex, Sant Cugat del Valles, Spain. ⁹Agencia de Salut Pública de Barcelona, Barcelona, Catalonia, Spain. ¹⁰CIBER Epidemiología y Salud Pública (CIBERESP), Madrid, Spain. ¹¹CREAF, Cerdanyola del Vallès, Catalonia, Spain. ¹²Institut Català de Recerca i Estudis Avançats (ICREA), Barcelona, Catalonia, Spain. ¹³email: federalucati@hotmail.com</footnote>
epidemiology of several mosquito-borne diseases since the tiger mosquito is a competent laboratory vector of at least 16 viruses⁵,⁶. While *Ae. aegypti* is considered as the principal vector of dengue and Zika, *Ae. albopictus* is a less efficient epidemic vector (with local exceptions in some cases), having developed an enhanced transmission for chikungunya facilitated by genetic adaptation (E1-A226V substitution) of the ECSA strain⁶⁻⁹. *Ae. albopictus* has a role on sporadic autochthonous disease transmission of arboviruses (as chikungunya, dengue and Zika) in Europe, and therefore its surveillance and control are considered a regional priority¹⁰.

Invasive alien species have historically been spread by humans. However, advances in transportation logistics resulting in higher air traffic and sea-born trading are driving a more rapid dispersal of non-indigenous species in the world, including many vectors of human diseases¹¹. *Ae. albopictus* is one of top invasive alien species in the world, being considered, together with *Ae. aegypti*, the most costly invasive species¹²,¹³. Its expansion in Europe coincides with the wave of introductions of invasive species in the continent, which started almost 40 years ago as a result of globalization¹⁴. Considering its limited flight range (less than 200 m/day)¹⁵, a crucial element for its success is the longevity of its desiccation-resistant eggs¹⁶, which can be passively transported by humans through commercial shipping of used tires and aquatic plants³ and ground vehicles¹⁰,¹⁷. In Spain, *Ae. albopictus* was first detected in Sant Cugat del Vallès, Catalonia, in 2004¹⁸ as a result of the nuisances associated with its aggressive anthropophilic behaviour¹⁹. Almost twenty years later, this species is well established in the Mediterranean coast of Spain and it is now colonizing inner territories of the Iberian Peninsula¹⁷,²⁰.

A biological invasion is a three-step process, which involves: initial dispersal, establishment of self-sustaining populations and spread to neighbouring habitats²¹,²². It is, however, during the initial dispersal when proactive management efforts can be more cost-effective for preventing the establishment of invasive species²²,²³. For instance, successful surveillance in New Zealand in the 1990s intercepted the entrance of *Ae. albopictus*²⁴, which has not been reported so far in this country. On the other hand, once this has become established, it is highly difficult to eradicate¹⁶. In this respect, studying *Ae. albopictus* population genetic structure and identifying its dispersal routes, its main drivers and scales, is crucial to understand the tiger mosquito's spread and design integrated surveillance and preparedness strategies²⁵,²⁶. Previous population genetic studies pointed to a worldwide chaotic dispersion pattern in *Ae. albopictus*, with Europe harbouring several distinct genotypes that have been linked to multiple independent introductions (e.g. ²⁵,²⁷,²⁸). Furthermore, in Europe human transportation networks have been shown to have facilitated the spread of this mosquito and conditioned its genetic and demographic patterns²⁹⁻³².

As different marker types can tell different stories³³,³⁴, it is crucial to analyse both mitochondrial (mtDNA) and nuclear (nDNA) DNA data to understand what major factors contributed to shaping the genetic attributes of invasive populations. In this regard, the presence of maternally-transmitted endosymbiotic bacteria may be relevant, as they can distort mtDNA phylogenies and reduce mtDNA diversity, a process which may have little to no effects on nDNA variation³⁵. If a maternally-inherited symbiont confers a selective advantage to its hosts, the mtDNA variants originally associated with the symbiont can rapidly spread through the host population and go to fixation, thus resulting in a great increase in the frequency of a single or few mtDNA haplotypes³⁶. Such symbionts are widespread in many arthropod species, being *Wolbachia* the most common maternally-inherited symbiont on the planet³⁷,³⁸. For this reason, incorporation of nDNA data is essential to corroborate results derived from mtDNA.

It is now widely accepted that a progress in the study of the ecology of *Ae. albopictus* requires an interdisciplinary approach¹¹,³⁹. For instance, population genetics and multilevel modelling analyses can be crucial to shed light on *Ae. albopictus* patterns of dispersal. Besides, it has been shown that genetic diversity in this species has been related to a higher adaptive potential²⁹ and a diversification of its interactions with the pathogens it carries³⁹. Understanding movement patterns in its populations in endemic areas is crucial to disentangle the dynamics of disease transmission between vectors and humans¹⁵. In this study we assess geographic patterns of mtDNA and nDNA variation in *Ae. albopictus* and use this information to investigate the species' dispersal routes across Spain. We also compare features of mtDNA and nDNA variation and assess whether they are incongruent, and test for *Wolbachia* presence to ascertain whether mtDNA diversity is affected by this maternally-transmitted parasite.

## Results

**Genetic variation and structure, and temporal trends of genetic diversity.** The nuclear DNA alignment (second internal transcribed spacer of ribosomal DNA—ITS2) included 470 sequences of 286 bp (424 sequences from Spain + 30 from France + 16 from Greece), while for the mitochondrial DNA alignment (cytochrome c oxidase gene subunit 1—COI) we obtained 471 sequences of 513 bp (424 from Spain + 40 from France + 7 from Greece) (Fig. 1, Table 1). Haplotype and nucleotide diversity estimates calculated at the province level ranged from 0.490 to 1 and from 0.004 to 0.064 for ITS2, respectively, and from 0 to 0.350 and from 0 to 0.0014 for COI (Table 1).

In the Iberian Peninsula, overall ITS2 genetic diversity was between three (haplotype diversity) and 10 (nucleotide diversity) times higher than for COI. Indeed, overall haplotype and nucleotide diversities were 0.7084 ± 0.025 and 0.00566 ± 0.00048 for ITS2, respectively, and 0.2196 ± 0.027 and 0.00054 ± 0.00008 for COI. We found 78 haplotypes defined by 83 polymorphic sites (21 parsimony informative) for ITS2, whereas only 18 haplotypes defined by 20 polymorphic sites (5 parsimony informative) were detected in the case of COI (Fig. 2).

No apparent association between haplotypes and geography was detected. Both ITS2 and COI haplotype networks showed the presence of worldwide dominant haplotypes and the absence of a clear genetic structure along the study area (Fig. 2). Specifically, in the ITS2 network, the most frequent and centrally placed haplotypes were observed in localities spanning the entire extension of the study area and were shared with other areas of the world. However, 89.7% of the haplotypes found were newly described (70/78), the majority being low-frequency haplotypes. In the case of COI, we obtained a much simpler star-like network defined by one
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Locations sampled for genetic analysis. Yellow circles are centered on sample locations with radius proportional to number of mosquitoes sampled at each location. Blue diamonds indicate locations of major commercial ports in Spain (ports with over 50,000 TEU -Twenty-Foot Equivalent Unit-, based on AAPA -American Association of Port Authorities- data from 2013). Grey lines indicate provinces of Spain. Land boundaries from Natural Earth. Province boundaries from GADM.</figcaption></figure>

<table><caption>**Table 1.** Basic genetic statistics of *Ae. albopictus* sampled provinces and countries for mtDNA (COI) and nuclear DNA (ITS2). N, sample size; S, number of polymorphic sites; HD, haplotype diversity; π, nucleotide diversity. Note that each province contains multiple sampling sites, as shown in Fig. 1.</caption><thead><tr><th>Province/Country</th><th>Year of first detection</th><th>Sampling period</th><th>N<sub>COI</sub></th><th>S<sub>COI</sub></th><th>HD<sub>COI</sub> (mean ± SE)</th><th>π<sub>COI</sub> (mean ± SE)</th><th>N<sub>ITS2</sub></th><th>S<sub>ITS2</sub></th><th>HD<sub>ITS2</sub> (mean ± SE)</th><th>π<sub>ITS2</sub> (mean ± SE)</th></tr></thead><tbody><tr><td>Barcelona</td><td>2004</td><td>2012–2015</td><td>104</td><td>10</td><td>0.262 ± 0.055</td><td>0.0008 ± 0.0003</td><td>104</td><td>32</td><td>0.693 ± 0.050</td><td>0.005 ± 0.001</td></tr><tr><td>Girona</td><td>2008</td><td>2011–2015</td><td>80</td><td>5</td><td>0.233 ± 0.062</td><td>0.0005 ± 0.0003</td><td>84</td><td>54</td><td>0.597 ± 0.059</td><td>0.006 ± 0.002</td></tr><tr><td>Tarragona</td><td>2005</td><td>2012–2015</td><td>7</td><td>1</td><td>0.286 ± 0.196</td><td>0.0006 ± 0.0005</td><td>4</td><td>2</td><td>0.833 ± 0.222</td><td>0.004 ± 0.002</td></tr><tr><td>Mallorca</td><td>2012</td><td>2014–2015</td><td>40</td><td>2</td><td>0.145 ± 0.074</td><td>0.0003 ± 0.0002</td><td>33</td><td>35</td><td>0.818 ± 0.067</td><td>0.012 ± 0.003</td></tr><tr><td>Huesca</td><td>2015</td><td>2015</td><td>3</td><td>0</td><td>0</td><td>0</td><td>5</td><td>10</td><td>1.000 ± 0.126</td><td>0.014 ± 0.005</td></tr><tr><td>Castelló</td><td>2010</td><td>2011–2015</td><td>7</td><td>1</td><td>0.286 ± 0.196</td><td>0.0006 ± 0.0005</td><td>7</td><td>5</td><td>0.524 ± 0.209</td><td>0.007 ± 0.003</td></tr><tr><td>València</td><td>2013</td><td>2011–2015</td><td>20</td><td>0</td><td>0</td><td>0</td><td>18</td><td>17</td><td>0.490 ± 0.142</td><td>0.006 ± 0.002</td></tr><tr><td>Alacant</td><td>2006</td><td>2011–2015</td><td>16</td><td>4</td><td>0.350 ± 0.148</td><td>0.0014 ± 0.0007</td><td>34</td><td>17</td><td>0.729 ± 0.078</td><td>0.005 ± 0.002</td></tr><tr><td>Murcia</td><td>2011</td><td>2012–2015</td><td>110</td><td>2</td><td>0.240 ± 0.049</td><td>0.0005 ± 0.0004</td><td>111</td><td>15</td><td>0.658 ± 0.031</td><td>0.004 ± 0.002</td></tr><tr><td>Málaga</td><td>2014</td><td>2014–2015</td><td>20</td><td>1</td><td>0.190 ± 0.108</td><td>0.0004 ± 0.0004</td><td>10</td><td>26</td><td>0.978 ± 0.059</td><td>0.020 ± 0.005</td></tr><tr><td>Granada</td><td>2014</td><td>2014–2015</td><td>4</td><td>0</td><td>0</td><td>0</td><td>6</td><td>14</td><td>0.800 ± 0.172</td><td>0.022 ± 0.007</td></tr><tr><td>Almería</td><td>2014</td><td>2014–2015</td><td>7</td><td>0</td><td>0</td><td>0</td><td>6</td><td>4</td><td>0.600 ± 0.215</td><td>0.005 ± 0.002</td></tr><tr><td>Gipuzkoa</td><td>2014</td><td>2014–2015</td><td>6</td><td>0</td><td>0</td><td>0</td><td>2</td><td>19</td><td>1.000 ± 0.500</td><td>0.064 ± 0.016</td></tr><tr><td>France</td><td>–</td><td>2014–2015</td><td>40</td><td>0</td><td>0</td><td>0</td><td>30</td><td>22</td><td>0.867 ± 0.043</td><td>0.012 ± 0.003</td></tr><tr><td>Greece</td><td>–</td><td>2013–2014</td><td>7</td><td>0</td><td>0</td><td>0</td><td>16</td><td>13</td><td>0.617 ± 0.135</td><td>0.006 ± 0.002</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Haplotype networks of COI mtDNA sequences (A) and ITS2 nuclear DNA sequences (B) analysed in *Ae. albopictus*. Each circle represents a unique haplotype and the circle area is proportional to the number of sequences of a given haplotype. Blue dots correspond to inferred unsampled haplotypes.</figcaption></figure>

main central haplotype connected to many low-frequency haplotypes with little genetic differentiation from the dominant sequence. Haplotypes were separated by no more than three nucleotide substitutions, although 83.3% of the haplotypes found were unique (15/18).

At the province level, we found a significant positive correlation between genetic diversity and colonisation time for COI, meaning that the provinces that were first colonised by *Ae. albopictus* bear higher mitochondrial genetic diversity ($R^2=0.497, p=0.007$ for haplotype diversity, $R^2=0.597, p=0.002$ for nucleotide diversity; Fig. 3). However, a deeper look showed that COI's variation in nucleotide diversity over time was almost non-existent ($\pi$ range: 0–0.0014), likely indicating a statistically significant but biologically irrelevant pattern. On the contrary, we found a significant negative correlation between ITS2 nucleotide diversity and colonisation time ($R^2=0.376, p=0.022$), which displayed higher variation ($\pi$ range: 0–0.022); the same relationship with haplotype diversity was not significant (Fig. 3).

Whereas these province-level correlations may be limited, to some extent, by the different sample sizes in each province, we also analysed the genetic distances among the individual sampled mosquitoes. In agreement with the haplotype network, multidimensional scaling (MDS) analysis revealed no apparent evidence of distinct genetic groups or geographic consistency within the nuclear dataset, although it did show a certain level of drift over time (Fig. 4). Most of the samples were grouped into a unique cluster in the centre of the MDS, with the exception of a few samples collected in 2014 and 2015 at the geographic extremes of the study area that laid slightly outside the main cluster (i.e. samples from the southernmost provinces of Málaga, Granada and Murcia, from the northernmost provinces of Gipuzkoa and Girona, and from the Balearic Islands in the east). The two-dimensional MDS solution accounted for 30% of the variance in the data (goodness of fit = 0.302).

Again using the individual sampled mosquitoes, a Mantel test of correlation revealed a significant positive correlation between ITS2 pairwise genetic distances and spatial distance ($r=0.163, p=0.001$), and negative correlations with potential tiger mosquito flux ($r=-0.022, p=0.036$) and spatial proximity ($r=-0.042, p=0.001$). The Mantel test did not find evidence of correlation between ITS2 genetic distance and temporal distance.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Relationship between genetic diversity and year of first detection in the analysed provinces. The mitochondrial COI fragment is indicated by grey triangles and the nuclear ITS2 gene by black circles. Only provinces with four or more samples were included in the analysis. See Table 1 for province list.</figcaption></figure>

**Wolbachia infection.** From a total of 62 samples screened using the *wsp* and 16S markers, 49 (79%) from 25 localities tested positive for infection for both markers. Two additional individuals yielded positive amplifications for 16S only and were thus excluded in reporting *Wolbachia* prevalence.

**Modelling genetic distance from spatial distance, temporal distance and mosquito flux.** We use Bayesian multiple-membership multilevel zero-inflated beta regressions to model ITS2 genetic distance as a function of mosquito flux, spatial proximity, spatial distance, and temporal distance (see "Methods"). These models rely on each pair of sampled mosquitoes as the unit of analysis and they suggest that ITS2 genetic distance is better explained by the combination of potential tiger mosquito flux, spatial proximity, spatial distance and temporal distance than by any of these variables on their own or in smaller combinations. This is most clearly seen in our leave-one-out cross validation (LOO) comparison, in which the expected log pointwise predictive density for model 6 (M6) is 138 points lower than the next best model, with a standard error of only 20 (Supplementary Table S1). It is also seen in the Bayesian R-squared comparison, in which M6 also has the highest value, although in this case the differences are small compared to the standard errors (complete model comparisons are presented in Supplementary Fig. S1 and Supplementary Table S2). The slope estimates for the main effects in the equation for the mean ($\mu$) of the beta distribution are shown in Fig. 5A and Supplementary Table S2, and should be read in conjunction with the slope estimates for the main effects in the equation for the probability of zeros ($\pi$), shown in Fig. 5B and Supplementary Table S2. The coefficients on the spatial proximity and potential tiger mosquito flux variables are negative in all models of $\mu$ and positive in all models of $\pi$, even when both variables are included together (model 5 -M5- and M6), indicating that greater spatial proximity and mosquito fluxes are associated with lower ITS2 genetic distance between sampled mosquitoes, with a higher probability of zero distances. In contrast, the coefficient on the temporal distance variable is positive in all models of $\mu$ and negative in all models of $\pi$ (albeit with its posterior distribution overlapping zero in M5 and M6), indicating that greater passage of time between samples is associated with greater ITS2 genetic distance and with lower probability of zero distances (Supplementary Figs. S2, S3). In M6 we find that greater spatial distance is also associated with lower probabilities ($\pi$) of zero ITS2 genetic distance (Fig. 5B). Although greater spatial distance is also associated with lower ITS2 genetic distance in the model of $\mu$ (Fig. 5A), the combined effect of the two components is an overall positive relationship between spatial distance and ITS2 genetic distance (Fig. 6A).

Figure 6A shows the effects predicted by M6 of changes in inter-point distance (reflected simultaneously in the spatial distance and spatial proximity variables) on ITS2 genetic distance. The range of inter-point distance values used for these predictions is the same as that observed in the data. Potential mosquito flux is held at its
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Multidimensional scaling plots of ITS2 genetic distances among *Ae. albopictus* samples, coloured by province (A) and year (B) of sample collection. Plots show the 2-dimensional solution using classical scaling. Goodness of fit = 0.302.</figcaption></figure>

median, the sampling years are set to 2011 and 2015 (to give the widest range observed) and the sample pair (for purposes of the random intercepts) is arbitrarily selected. The overall pattern is non-linear, reflecting the combination of the spatial proximity and spatial distance variables. There is an initial steep increase in predicted ITS2 genetic distance for samples taken within a few meters of one another, which can be seen most clearly in the inset plot of Fig. 6A. Beyond these highly proximate samples, predicted ITS2 genetic distance then rises more gradually, driven by the spatial distance variable.

Figure 6B shows the effects predicted by M6 of changes in potential mosquito flux on ITS2 genetic distance. The range of fluxes used for these predictions is the same as that observed in the data. Inter-point distance is held at its median and used for calculating the spatial proximity and spatial distance variables, the sampling years are set to 2011 and 2015 (to give the widest range observed) and the sample pair (for purposes of the random intercepts) is arbitrarily selected. Here we see a steep drop in predicted ITS2 genetic distance as potential mosquito flux increases from 0 to several hundred per day, with effect of increased mosquito fluxes becoming weaker at higher values, reflecting the log-linear relationship used in the model. Although in this case the pattern is overwhelmed by the uncertainty of the predictions (the wide posterior predictive distributions shown in the lighter
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Estimated relationship between ITS2 pairwise genetic distance, spatial distance (sp_dist; geodesic distance between sample locations in meters), spatial proximity (sp_proxy; measured as the negative exponential of distance), potential tiger mosquito flux (mosqirim; estimated from commuting patterns and tiger mosquito population distribution) and temporal distance (yr_diff; measured as absolute difference between years in which samples were taken) on the beta mean (μ) parameter (A) and on the zeros (B) in the zero-inflated beta regression models. Parameters are estimated from a set of Bayesian multilevel zero-inflated Beta regressions with multiple-membership random intercepts for the samples and sampling years represented in each pair.</figcaption></figure>

shades of blue), this is the result of the combined uncertainty of the other variables in the model; the effect of potential mosquito flux, net of these other variables, is very clearly shown in the parameter estimates in Fig. 5.

Figure 7 shows the combined effects of inter-point distance and potential tiger mosquito flux changing together. We see, first, how the lowest values of inter-point distance (bottom edge of plot) correspond with the lowest predicted ITS2 genetic distances while the lowest values of potential tiger mosquito flux (left edge of plot) correspond with the highest. When both variables are at their lowest (bottom left corner), inter-point distance is determinative: predicted ITS2 genetic distance is low here regardless of potential tiger mosquito flux. Beyond these lowest values, we see how predicted tiger mosquito flux acts to hold predicted ITS2 genetic distance down even as distance increases. Although predicted ITS2 genetic distances never reach their lowest values at these longer distances, the potential tiger mosquito fluxes of around 26,000 mosquitoes per day (this is the potential flux, for example, between Barcelona municipality and El Prat de Llobregat) keep the predicted ITS2 genetic distance within the range of 0.0035–0.0036, even at distances of 100 km.

Note, finally, that although the ITS2 genetic distances and effect sizes shown in Figs. 5, 6 and 7 are small overall, this should be interpreted in light of the distribution of observed genetics distances, which ranged from 0 to only 0.07, with a standard deviation of 0.01.

## Discussion

*Aedes albopictus* has spread worldwide and particularly in Europe at a fast pace, making it one of the 100 most invasive species on Earth¹². In Spain, after *Ae. albopictus* was first found near Barcelona in 2004¹⁸, a continuous spread along the Mediterranean coast was observed¹⁷. All Mediterranean provinces are currently colonised, along with the Basque Country in the northern coast and several inland territories²⁰, where the species continues to expand. In light of this, understanding *Ae. albopictus* dispersal routes across scales is crucial for planning effective early warning surveillance in non-invaded areas and implementing surveillance and control activities in the areas already colonised. The same information can also be valuable in predicting the transmission risk of pathogens by this vector. Knowledge of the effect of vehicles and transport infrastructures in the genetic structure of vector populations and their overall spreading capacity can comprehensively show our potential as a natural selective force and the existing contradictions between globalization and our efforts to combat biological invasions and pests⁴⁰.

Our models indicate that at very small spatial scales (i.e. several meters) the genetic variation measured by ITS2 is sharply reduced, likely representing seasonal mosquito pools that come from a main source. Beyond
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Predicted ITS2 pairwise genetic distance as a function of the spatial distance and spatial proximity variables taken together (A) and potential tiger mosquito flux (B) in the zero-inflated beta regression Model 6. Panel A shows predictions for the range of inter-point distances in the modelled data (0–940 km), holding potential mosquito flux at its observed median, setting the sampling years to 2011 and 2015, and arbitrarily selecting a sample pair and its associated provinces for purposes of the model's random intercepts. The inset plot in this panel shows a close-up of the predictions at very small distances (0–10 m). Panel B shows predictions for the range of mosquito fluxes in the modelled data (0–672 km), holding inter-point distance at its median, setting the sampling years to 2011 and 2015, and arbitrarily selecting a sample pair and its associated provinces.</figcaption></figure>

these scales, genetic variability steadily increases with spatial distance, as a clear positive correlation between dispersal distance and genetic variation appears. This is reflected in the combined effects of the (small scale) spatial proximity variable and the linear spatial distance variable. More interestingly, our models suggest that human transportation has a role in shaping *Ae. albopictus* nuclear genetic structure by means of passive dispersal of adult tiger mosquitoes in cars: there is a clear negative relationship between mosquito flux and genetic variation. Previous studies have also highlighted that, although *Ae. albopictus* has low natural dispersal capabilities, human-aided transport (especially cars) has probably facilitated significantly the tiger mosquito's movement and invasion process<sup>29,30,32,41</sup>. Our findings are consistent with these but go a step further by suggesting that the “hitchhiking” of *Ae. albopictus* in cars observed in Eritja, Palmer, Roiz, Sanpera-Calbet and Bartumeus<sup>30</sup> actually helps to explain observed patterns of population genetics. While genetic variability increases with spatial distance, car transport can strongly reduce this effect. The high-resolution of our dataset makes it possible to show the significant role of human transportation in shaping the genetic constitution of *Ae. albopictus* and promoting regional gene flow. Mosquito movement can be affected by human activities like commuting and human-made structures like roads, which combined act as bridges for dispersal by favouring gene flow and promoting genetic mixing<sup>32</sup>.

At a much broader spatial scale, we find a general lack of genetic structure and geographic consistency among haplotypes across the large expanse of the Iberian Peninsula, with numerous haplotypes shared among several distant areas (e.g. the most abundant ITS2 haplotypes found in Spain were also detected in several other European and Asian countries). This suggests that human-mediated large-scale dispersal of *Ae. albopictus* is also common, and point to a pattern of regular introductions of the species from abroad, through e.g. transportation of used tires and aquatic plants, which allows for survival and establishment at long distances of whole batches of eggs. Previous genetic studies of *Ae. albopictus* have also showed little or no genetic structure according to geography, both in the native and introduced range of the species (reviewed in<sup>39</sup>). Such a genetic feature is likely the consequence of the high level of human-mediated spread from several genetically distinct source populations followed by global dispersal, and it is concordant to what has been found in other wide-ranging invasive insect species, especially those that are closely associated with humans, e.g. the German and American cockroaches (*Blattella germanica* and *Periplaneta americana*, respectively)<sup>42,43</sup> and the longhorn crazy ant *Paratrechina longicornis*<sup>44</sup>. Taken together, our results highlight the role of human activities in promoting unintentional mid- and long-distance dispersal and thus in shaping the current genetic structure of insect species commonly found in human-modified landscapes.

Our study has some limitations that should be considered. First, the modelled genetic diversity range, which was obtained from ITS2 pairwise genetic distances, is low. This likely has to do with the relatively low variability of the analysed genetic marker. Indeed, nuclear genes, as well as mitochondrial fragments, are expected to be less variable and bear lower resolving power than highly mutating markers, e.g. microsatellites. Second, although we rely on a relatively large number of sampling sites, additional sampling across larger areas of Spain could provide
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Predicted ITS2 pairwise genetic distance (indicated by fill colour) as a function of inter-point distance (the spatial distance and spatial proximity variables taken together) and potential tiger mosquito flux in the zero-inflated beta regression Model 6. Predictions are shown inter-point distances between 0 and 100 km and potential tiger mosquito fluxes between 0 and 30,000, setting the sampling years to 2011 and 2015, and arbitrarily selecting a sample pair and its associated provinces for purposes of the model's random intercepts.</figcaption></figure>

better balance and a wider range of values for the models. Third, intragenomic heterogeneity (i.e. presence of multiple haplotypes within the same individual) of ITS2 has been reported in several mosquito species including the genus *Aedes*45,46, and this can pose a challenge in DNA sequencing and analysis. In this study, all individuals presenting intragenomic heterogeneity were thus excluded from analysis. As for future directions, although ITS2 has proved to be a useful marker for studies on the spread of *Ae. albopictus*47, the reassessment of the species' genetic diversity and population structure through the use of molecular markers with greater variability and/or potential, such as microsatellites and single nucleotide polymorphisms (SNPs), could provide more detailed and deeper insights into the description of fine-scale dispersal patterns, gene flow and introduction routes.

Global patterns of mtDNA and nuclear variation were highly discordant, with mtDNA showing little genetic diversity and a single star-like haplotype network. This is in agreement with previous studies showing *Ae. albopictus* levels of nuclear variation within the range of most insects, but extremely low mtDNA variation both within and among populations48–51. Interestingly, we found a high infection rate (79%) of *Wolbachia* in the studied *Ae. albopictus* samples. *Wolbachia* is a genus of maternally-inherited endosymbiotic bacteria that is known to induce male killing, feminization, parthenogenesis and cytoplasmic incompatibility, which facilitate its spread within the arthropod population52. *Wolbachia* is capable of inducing selective sweeps in mtDNA, i.e. fixation of a single or few mtDNA haplotypes that may become widespread in the host population through cytoplasmic hitchhiking driven by *Wolbachia* invasions36. Selective sweeps on mtDNA have been shown to not only reduce haplotype diversity producing a characteristic single star-like network, but also to cause the remaining set of haplotypes to deviate from neutrality35. Within Culicidae, the natural presence of *Wolbachia* has been documented in more than 30 species (e.g.53–56), with *Ae. albopictus* harbouring significantly lower mtDNA diversity than the unfected species55. In light of this, our results point to *Wolbachia* as causative agent for the lack of mitochondrial
polymorphism here recovered, as suggested for other insect species, e.g. *Acraea* butterflies<sup>57</sup> and the cherry fruit fly *Rhagoletis cerasi*<sup>58</sup>. Nevertheless, it has to be noted that *Wolbachia* can show seasonal fluctuations in infection rates<sup>59</sup>, which were not addressed in this study. Hence, the infection rate reported here may not capture a representative overview of *Wolbachia* infection. Moreover, the low variability in mtDNA in the introduced ranges of *Ae. albopictus* could also be caused by demographic processes, such as genetic drift or population bottleneck during rapid colonization<sup>60</sup>. However, demographic processes cause changes in variation for both mitochondrial and nuclear markers, even though mtDNA is expected to show a stronger response<sup>61</sup>. Furthermore, the lack of mtDNA variation was also found in the native range of *Ae. albopictus*<sup>48</sup>, but see<sup>62</sup>, strengthening our hypothesis of a *Wolbachia*-induced selective sweep. Alternatively, we cannot rule out that the low COI diversity here detected may be at least partially due to the length of the analysed fragment, as higher diversity has been observed when targeting larger COI fragments (> 1300 bp)<sup>63,64</sup>. Further research is needed to test this hypothesis.

Mitochondrial DNA has been extensively used to shed light on the geographic origin of invasive *Ae. albopictus* populations and arthropods in general<sup>49,50,65</sup>, and disentangle their phylogeographic history<sup>62,66</sup>. Nevertheless, there is broad recognition that mtDNA can commonly be under selection, thus challenging the assumption of neutrality postulated by several population, phylogeographic and phylogenetic studies<sup>33,67,68</sup>. Selection can arise due to e.g. mito-nuclear co-evolution, adaptation of mtDNA to different climatic/environmental conditions, and selective sweeps of beneficial mtDNA haplotypes<sup>68</sup>. Mitochondrial DNA is a single, linked molecule with low or no recombination, meaning that selective processes such as selective sweeps can have a profound impact on the apparent rate of genetic drift<sup>33</sup>. Therefore, we suggest caution should be used in drawing conclusions from mtDNA alone and, whenever possible, aim for a multilocus approach to achieve a correct understanding of the genetic population structure and history of the study species.

## Methods

**Sample collection and DNA extraction.** *Ae. albopictus* samples were collected during the period 2011–2015 at 140 locations encompassing most of the current species distribution in Spain (13 provinces mostly located along the Mediterranean coast; Fig. 1, Table 1). Additional samples were collected from two areas of France (Nice N 43° 38' 32" E 7° 5' 27" and Montpellier N 43° 40' 39" E 4° 2' 35") and two regions of Greece (the area of Pylaia, in Salonica N 40° 27' 22" E 23° 13' 20" and Sykia N 40° 2' 19" E 23°56' 22"), which were used as complementary data from these Mediterranean areas (Table 1). Sample collection was performed using different methods according to the life stage (adults, larvae or eggs) (see Supplementary Methods for details). All samples were stored in absolute ethanol and kept at -20 °C until genetic analyses.

DNA was extracted from the whole bodies of mosquitoes (adults or larvae) in a final volume of 250 µL using the HotShot protocol<sup>69</sup>.

**Ae. albopictus nuclear and mitochondrial gene sequencing.** We amplified two gene regions, including one nuclear ribosomal gene (ITS2) and one mitochondrial fragment (COI). The following primers were used for amplification and sequencing: for ITS2, primers ITS-CP-P1A (5'-GTGGATCCTGTGAACTGCAGG ACACATG-3') and ITS-CP-P1B (5'-GTGTCGACATGCTTAAATTTAGGGGTA-3')<sup>70</sup>, and for COI, primers LCOI490 (5'-GGTCAACAAATCATAAAGATATTGG-3') and HCO2198 (5'-TAAACTTCAGGGTGA CCAAAAAATCA-3')<sup>71</sup> or the degenerated primers ZplankF1_M13 (5'-TGTAAAACGACGGCCAGTTCTAS- WAATCATAARGATATTGG-3') and ZplankR1_M13 (5'-CAGGAAACAGCTATGACTTCAGGRTGRCCR AARAATCA-3')<sup>72</sup> that had the M13 primers attached following the suggestion of Ivanova, Zemlak, Hanner and Hebert<sup>73</sup>. Further details on amplification and sequencing are presented in the Supplementary Methods. In the case of ITS2, 22.7% (138) of the individuals could not be considered for analysis due to intragenomic heterogeneity (simultaneous presence of two or more haplotypes in the same individual). Resulting sequences were aligned with other relevant *Ae. albopictus* sequences from each gene (retrieved from GenBank) using the ClustalW algorithm in MEGA 7<sup>74</sup>.

**Wolbachia screening.** To test whether mtDNA variability could be affected by the presence of endosymbiotic bacteria of the genus *Wolbachia*, we analysed 62 randomly selected samples (13.2% of the overall samples). Individuals were selected aiming to cover all sampling years and all studied provinces. Two molecular markers were used for detecting *Wolbachia* infection, namely *wsp* and 16S rDNA. Primers used for amplification and sequencing were: for the *wsp* marker, primers 81F (5'-TGGTCCAATAAGTGATGAAGA-3') and 69R (5'-AAA AATTAAACGCTACTCCA-3')<sup>75,76</sup>, while for 16S, primers 16SF (5'-CGGGGGAAAAATTTATTGCT-3') and 16SR (5'-AGCTGTAATACAGAAAGTAAA-3')<sup>77,78</sup>. Amplification conditions followed Wiwatanaratanabutr<sup>53</sup> in the case of *wsp*, and Heddi, Grenier, Khatchadourian, Charles and Nardon<sup>78</sup> in the case of 16S. PCR products were visualized on a 1% agarose gel. To validate the results, two PCR replicates were run for each sample (following Carvajal, Hashimoto, Harnandika, Amalin and Watanabe<sup>54</sup>). A third replicate was run for samples that showed incongruent results based on the two prior replicates. *Wolbachia* infection was confirmed by two successful amplifications of both molecular markers. In light of the high infection rate detected in the selected samples (see Results), we decided to exclude COI from population structure and multilevel modelling analyses (see below), as mtDNA diversity was likely affected by the presence of *Wolbachia*.

**Genetic diversity and structure.** Genetic diversity indices including number of haplotypes, number of polymorphic sites, haplotype diversity and nucleotide diversity for mtDNA and nuclear DNA were estimated for the whole dataset and for each analysed province using DNASP 6.11.01<sup>79</sup>. To estimate gene genealogies, we used HAPLOVIEWER, which turns trees built from traditional phylogenetic methods into haplotype genealogies<sup>80</sup>. We estimated the phylogeny using a maximum-likelihood approach as implemented in RAxML 7.7.1<sup>81</sup>, with a
GTRCAT model rate of heterogeneity and no invariant sites for COI, and a gamma model rate of heterogeneity and invariant sites (GTR+G+I) for ITS2. The most appropriate model of nucleotide evolution was selected using jModelTest 2.1.3⁸² under the Akaike information criterion (AIC). Input data were COI or ITS2 sequences from each individual, subsequently collapsed into haplotypes. The best tree was selected for network construction in HAPLOVIEWER.

Population structure of nuclear DNA was characterized using classical multidimensional scaling (MDS), also known as principal coordinates analysis⁸³. Pairwise distances between ITS2 sequences were calculated in MEGA by estimation of evolutionary divergence over sequence pairs, using the Kimura 2-parameter substitution model. MDS analysis was then performed with the pairwise genetic distance matrix, using the cmdscale function of the stats package in R 4.1.2⁸⁴.

**Construction of covariates.** Our analysis of the drivers of genetic distance focused on four constructed covariates: spatial distance, spatial proximity, temporal distance, and potential tiger mosquito flux. We constructed the spatial distance variable as the geodesic distance between the locations of each sampled mosquito (inter-point distance) in meters, calculated using the Vincenty method on the WGS84 ellipsoid as implemented in the geosphere package for R⁸⁵. We constructed spatial proximity as the negative exponential of spatial distance in km, represented as e⁻⁶, where e is Euler's constant and d is spatial distance in km. This approach gives high values for spatial proximity (0.82–1) within the 200 m buffer often taken as the tiger mosquito's maximum flying distance, with values then quickly dropping and nearing zero by 5 km. The combination of these spatial distance and spatial proximity variables together is intended to capture the effects of inter-point distance at different scales. We constructed the temporal distance variable as the absolute number of years elapsed between the sampling times when each member of the sample pair was captured, rounded up to the nearest year (another way of thinking about this is as the number of mosquito seasons between each capture time). We constructed potential tiger mosquito flux as the potential daily bidirectional gross number of *Ae. albopictus* moving between each pair of municipalities. We estimated this by combining municipality-level *Ae. albopictus* risk estimates from the *Mosquito Alert* citizen science platform⁸⁶ with commuter flow estimates drawn from the Spanish Labour Force Survey (LFS) in a manner similar to that described in Eritja, Palmer, Roiz, Sanpera-Calbet and Bartumeus³⁰ (see Supplementary Methods for further details).

**Modelling the influence of spatial proximity, mosquito flux, and temporal distance on genetic distance.** We explored drivers of nuclear genetic structure by carrying out simple Mantel tests of correlations⁸⁷ between ITS2 pairwise genetic distance and potential tiger mosquito flux, spatial distance, spatial proximity, and temporal distance, using pairs of individual mosquitoes as the unit of analysis. We implemented all tests in the ecodist package for R⁸⁸, using 1000 permutations with bootstrap confidence limits estimated using 500 iterations.

We then used Bayesian multiple membership multilevel regressions to model ITS2 pairwise genetic distances between sampled mosquito pairs as a function of potential tiger mosquito flux between sampling sites, spatial distance and spatial proximity between sampling sites, and temporal distance.

We relied on zero-inflated beta regressions⁸⁹, in which the distribution of the stochastic component is a mixture of a beta distribution and a degenerate distribution in 0. Following Figueroa-Zúñiga, Arellano-Valle and Ferrari⁹⁰ and Branscum, Johnson and Thurmond⁹¹, we used the beta distribution for ITS2 genetic distance values in (0, 1) because it is a highly flexible distribution that is defined on that interval. As these authors have done, we parameterized the beta distribution in terms of a mean (μ) and a parameter (φ) that captures precision. The probability density function of a variable (y) in this parameterization is:

$$b(y|\mu, \phi) = \frac{\Gamma(\phi)}{\Gamma(\mu\phi)\Gamma(\mu\phi)\Gamma((1-\mu)\phi)} y^{\mu\phi-1} (1-y)^{(1-\mu)\phi-1}, 0 < y < 1 \quad (1)$$

Since the ITS2 genetic distances in our dataset include zeros, for which the beta distribution is undefined, we used a zero-inflation component in our models. In addition to making it possible to include the zeros in the analysis (without adding artificial noise to them), this approach also has the advantage of treating them separately, which should minimize any problems associated with the inadvertent sampling of siblings (see Supplementary Methods). As explained in Ospina and Ferrari⁹², the probability density function of a variable y for this zero-inflated beta (zib) mixture is:

$$\text{zib}(y|\pi, \mu, \phi) = \begin{cases} \pi, & \text{if ITS2 distance} = 0 \\ (1-\pi)b(y|\mu, \phi), & \text{if ITS2 distance} \in (0, 1) \end{cases} \quad (2)$$

where 0 < π < 1 is the probability of observing an ITS2 distance of 0, and μ and φ are defined as above.

We treated the observed ITS2 genetic distances between each sample pair i as independent random variables $y_i, ..., y_n$ drawn from this zero inflated beta distribution such that $y_i \sim \text{zib}(y|\pi_i, \mu_i, \phi_i)$. We used a logit link to model both π and μ, and a log link to model φ. We fit six models (M1–M6) with different combinations of mosquito flux, spatial proximity, temporal distance, and spatial distance as covariates. The models for μ are specified as following:
$$
\begin{align}
\text{M1: } \ln\left(\frac{\mu_i}{1-\mu_i}\right) &= \alpha_0 + \sum_{j \in \text{yr}(i)} \alpha_j + \sum_{k \in \text{ID}(i)} \alpha_k + \beta_1 \text{sp\_prox} \nonumber \\
\text{M2: } \ln\left(\frac{\mu_i}{1-\mu_i}\right) &= \alpha_0 + \sum_{j \in \text{yr}(i)} \alpha_j + \sum_{k \in \text{ID}(i)} \alpha_k + \beta_2 \ln_{\text{mosq\_flux}} \nonumber \\
\text{M3: } \ln\left(\frac{\mu_i}{1-\mu_i}\right) &= \alpha_0 + \sum_{j \in \text{yr}(i)} \alpha_j + \sum_{k \in \text{ID}(i)} \alpha_k + \beta_1 \text{sp\_prox} + \beta_3 \text{yr\_diff} \nonumber \\
\text{M4: } \ln\left(\frac{\mu_i}{1-\mu_i}\right) &= \alpha_0 + \sum_{j \in \text{yr}(i)} \alpha_j + \sum_{k \in \text{ID}(i)} \alpha_k + \beta_2 \ln_{\text{mosq\_flux}} + \beta_3 \text{yr\_diff} \nonumber \\
\text{M5: } \ln\left(\frac{\mu_i}{1-\mu_i}\right) &= \alpha_0 + \sum_{j \in \text{yr}(i)} \alpha_j + \sum_{k \in \text{ID}(i)} \alpha_k + \beta_1 \text{sp\_prox} + \beta_2 \ln_{\text{mosq\_flux}} + \beta_3 \text{yr\_diff} \nonumber \\
\text{M6: } \ln\left(\frac{\mu_i}{1-\mu_i}\right) &= \alpha_0 + \sum_{j \in \text{yr}(i)} \alpha_j + \sum_{k \in \text{ID}(i)} \alpha_k + \beta_1 \text{sp\_prox} + \beta_2 \ln_{\text{mosq\_flux}} + \beta_3 \text{yr\_diff} + \beta_4 \text{sp\_dist} \tag{3}
\end{align}
$$

where $\alpha_i$ represents a random intercept for each of the years in which the mosquitoes in pair $i$ were sampled, and $\alpha_k$ represents a random intercept for each of the sampled mosquitoes. This multiple membership approach allows us to account for the effects of sample year and sample itself, which could otherwise confound results given that each pair may be connected to other pairs by the sampling years and sampling that they share⁹³. $\beta_1, \beta_2, \beta_3$ and $\beta_4$ are the slopes on the spatial proximity (sp_prox), log mosquito flux (ln_mosq_flux), temporal distance (yr_diff) and spatial distance (sp_dist) variables, and $\alpha_0$ is the overall model intercept.

We explicitly model the $\phi$ parameter of our beta distributions as a function of the provinces in which each sample in the pair was taken. This allows us to explore how the precision of the beta distribution changes across geography. The choice of provinces as the areal units for this part of the model is based on these units representing patterns of human settlement and activity that should be of relevance to tiger mosquito spreading patterns, while also being large enough to avoid adding too many additional parameters to an already complicated model. The model for $\phi$ is the same in all six models:

$$ \ln(\phi_i) = \alpha_0 + \sum_{j \in prv(i)} \alpha_j \quad (4) $$

Finally, we modelled $\pi$ (the probability of ITS2 pairwise distance being equal to zero) as:

$$
\begin{align}
\text{M1: } \ln\left(\frac{\pi_i}{1-\pi_i}\right) &= \alpha_0 + \sum_{j \in \text{yr}(i)} \alpha_j + \beta_1 \text{sp\_prox} \nonumber \\
\text{M2: } \ln\left(\frac{\pi_i}{1-\pi_i}\right) &= \alpha_0 + \sum_{j \in \text{yr}(i)} \alpha_j + \beta_2 \ln_{\text{mosq\_flux}} \nonumber \\
\text{M3: } \ln\left(\frac{\pi_i}{1-\pi_i}\right) &= \alpha_0 + \sum_{j \in \text{yr}(i)} \alpha_j + \beta_1 \text{sp\_prox} + \beta_3 \text{yr\_diff} \nonumber \\
\text{M4: } \ln\left(\frac{\pi_i}{1-\pi_i}\right) &= \alpha_0 + \sum_{j \in \text{yr}(i)} \alpha_j + \beta_2 \ln_{\text{mosq\_flux}} + \beta_3 \text{yr\_diff} \nonumber \\
\text{M5: } \ln\left(\frac{\pi_i}{1-\pi_i}\right) &= \alpha_0 + \sum_{j \in \text{yr}(i)} \alpha_j + \beta_1 \text{sp\_prox} + \beta_2 \ln_{\text{mosq\_flux}} + \beta_3 \text{yr\_diff} \nonumber \\
\text{M6: } \ln\left(\frac{\pi_i}{1-\pi_i}\right) &= \alpha_0 + \sum_{j \in \text{yr}(i)} \alpha_j + \beta_1 \text{sp\_prox} + \beta_2 \ln_{\text{mosq\_flux}} + \beta_3 \text{yr\_diff} + \beta_4 \text{sp\_dist} \tag{5}
\end{align}
$$

The systematic part of the model for $\pi$ is almost identical to the systematic part of the model for $\mu$. The only difference is the absence of the sample random intercepts, which are excluded here for computational reasons. Details on model fitting and accuracy estimation are outlined in the Supplementary Methods.

## Data availability

Newly generated mitochondrial and nuclear DNA sequence data were deposited in GenBank under accession numbers OP060971-OP060988 and OP077008-OP077085.

Received: 21 July 2022; Accepted: 22 November 2022

Published online: 30 November 2022

## References

1. Hawley, W. A. The biology of *Aedes albopictus*. *J. Am. Mosq. Control Assoc.* **1**, 1–39 (1988).

2. Benedict, M. Q., Levine, R. S., Hawley, W. A. & Lounibos, L. P. Spread of the tiger: global risk of invasion by the mosquito *Aedes albopictus*. *Vector-Borne Zoonotic Dis.* **7**, 76–85 (2007).

3. Paupy, C., Delatte, H., Bagny, L., Corbel, V. & Fontenille, D. *Aedes albopictus*, an arbovirus vector: From the darkness to the light. *Microb. Infect.* **11**, 1177–1185 (2009).

4. Delatte, H. et al. Blood-feeding behavior of *Aedes albopictus*, a vector of Chikungunya on La Réunion. *Vector-Borne Zoonotic Dis.* **10**, 249–258 (2010).

5. Pereira-dos-Santos, T., Roiz, D., Lourenço-de-Oliveira, R. & Paupy, C. A systematic review: Is *Aedes albopictus* an efficient bridge vector for zoonotic arboviruses? *Pathogens* **9**, 266 (2020).

6. Gratz, N. Critical review of the vector status of *Aedes albopictus*. *Med. Vet. Entomol.* **18**, 215–227 (2004).

7. Grard, G. et al. Zika virus in Gabon (Central Africa)—2007: A new threat from *Aedes albopictus*? *PLoS Negl. Trop. Dis.* **8**, e6281 (2014).

8. Lambrechts, L., Scott, T. W. & Gubler, D. J. Consequences of the expanding global distribution of *Aedes albopictus* for dengue virus transmission. *PLoS Negl. Trop. Dis.* **4**, e646 (2010).

9. Lounibos, L. P. & Kramer, L. D. Invasiveness of *Aedes aegypti* and *Aedes albopictus* and vectorial capacity for chikungunya virus. *J. Infect. Dis.* **214**, S453–S458 (2016).

10. European Centre for Disease Prevention and Control (ECDC). *Vector Control with a Focus on Aedes aegypti and Aedes albopictus Mosquitoes: Literature Review and Analysis of Information* (ECDC, Stockholm, Sweden, 2017).

11. Tatem, A. J., Hay, S. I. & Rogers, D. J. Global traffic and disease vector dispersal. *PNAS* **103**, 6242–6247 (2006).
12. Lowe, S., Browne, M., Boudjelas, S. & De Poorter, M. *100 of the World's Worst Invasive Alien Species: A Selection From the Global Invasive Species Database*, Vol. 12 (Invasive Species Specialist Group, 2000).

13. Diagne, C. et al. High and rising economic costs of biological invasions worldwide. *Nature* **592**, 571–576 (2021).

14. Hulme, P. E. Trade, transport and trouble: Managing invasive species pathways in an era of globalization. *J. Appl. Ecol.* **46**, 10–18 (2009).

15. Marini, F., Caputo, B., Pombi, M., Tarsitani, G. & Della-Torre, A. Study of *Aedes albopictus* dispersal in Rome, Italy, using sticky traps in mark-release-recapture experiments. *Med. Vet. Entomol.* **24**, 361–368 (2010).

16. Bonizzoni, M., Gasperi, G., Chen, X. & James, A. A. The invasive mosquito species *Aedes albopictus*: current knowledge and future perspectives. *Trends Parasitol.* **29**, 460–468 (2013).

17. Collantes, F. et al. Review of ten-years presence of *Aedes albopictus* in Spain 2004–2014: Known distribution and public health concerns. *Parasit Vectors* **8**, 1–11 (2015).

18. Aranda, C., Eritja, R. & Roiz, D. First record and establishment of the mosquito *Aedes albopictus* in Spain. *Med. Vet. Entomol.* **20**, 150–152 (2006).

19. Giménez, N. et al. Introduction of *Aedes albopictus* in Spain: A new challenge for public health. *Gac. Sanit.* **21**, 25–28 (2007).

20. European Centre for Disease Prevention and Control and European Food Safety Authority. Mosquito maps [internet]. Stockholm: ECDC. https://ecdc.europa.eu/en/disease-vectors/surveillance-and-disease-data/mosquito-maps (2022).

21. Shigesada, N. & Kawasaki, K. *Biological Invasions: Theory and Practice* (Oxford University Press, 1997).

22. Puth, L. M. & Post, D. M. Studying invasion: Have we missed the boat? *Ecol. Lett.* **8**, 715–721 (2005).

23. Leung, B. et al. An ounce of prevention or a pound of cure: Bioeconomic risk analysis of invasive species. *Proc. R Soc. Lond. Ser. B Biol. Sci.* **269**, 2407–2413 (2002).

24. Lounibos, L. P. Invasions by insect vectors of human disease. *Annu. Rev. Entomol.* **47**, 233–266 (2002).

25. Manni, M. et al. Genetic evidence for a worldwide chaotic dispersion pattern of the arbovirus vector, *Aedes albopictus*. *PLoS Negl. Trop. Dis.* **11**, e0005332 (2017).

26. Roiz, D. et al. Integrated *Aedes* management for the control of *Aedes*-borne diseases. *PLoS Negl. Trop. Dis.* **12**, e0006845 (2018).

27. Lühken, R. et al. Microsatellite typing of *Aedes albopictus* (Diptera: Culicidae) populations from Germany suggests regular introductions. *Infect. Genet. Evol.* **8**, 104237 (2020).

28. Battaglia, V. et al. The worldwide spread of the tiger mosquito as revealed by mitogenome haplogroup diversity. *Front. Genet.* **7**, 208 (2016).

29. Medley, K. A., Jenkins, D. G. & Hoffman, E. A. Human-aided and natural dispersal drive gene flow across the range of an invasive mosquito. *Mol. Ecol.* **24**, 284–295 (2015).

30. Eritja, R., Palmer, J. R., Roiz, D., Sanpera-Calbet, I. & Bartumeus, F. Direct evidence of adult *Aedes albopictus* dispersal by car. *Sci. Rep.* **7**, 1–15 (2017).

31. Sherpa, S. et al. Unravelling the invasion history of the Asian tiger mosquito in Europe. *Mol. Ecol.* **28**, 2360–2377 (2019).

32. Swan, T. et al. A literature review of dispersal pathways of *Aedes albopictus* across different spatial scales: Implications for vector surveillance. *Parasit Vectors* **15**, 1–13 (2022).

33. Ballard, J. W. O. & Whitlock, M. C. The incomplete natural history of mitochondria. *Mol. Ecol.* **13**, 729–744. https://doi.org/10.1046/j.1365-294X.2003.02063.x (2004).

34. Toews, D. P. L. & Brelsford, A. The biogeography of mitochondrial and nuclear discordance in animals. *Mol. Ecol.* **21**, 3907–3930. https://doi.org/10.1111/j.1365-294X.2012.05664.x (2012).

35. Hurst, G. D. & Jiggins, F. M. Problems with mitochondrial DNA as a marker in population, phylogeographic and phylogenetic studies: The effects of inherited symbionts. *Proc. R. Soc. B: Biol. Sci.* **272**, 1525–1534 (2005).

36. Cariou, M., Duret, L. & Charlat, S. The global impact of *Wolbachia* on mitochondrial diversity and evolution. *J. Evol. Biol.* **30**, 2204–2210 (2017).

37. Zug, R. & Hammerstein, P. Still a host of hosts for *Wolbachia*: analysis of recent data suggests that 40% of terrestrial arthropod species are infected. *PLoS ONE* **7**, e38544 (2012).

38. Weinert, L. A., Araujo-Jnr, E. V., Ahmed, M. Z. & Welch, J. J. The incidence of bacterial endosymbionts in terrestrial arthropods. *Proc. R. Soc. B: Biol. Sci.* **282**, 20150249 (2015).

39. Goubert, C., Minard, G., Vieira, C. & Boulesteix, M. Population genetics of the Asian tiger mosquito *Aedes albopictus*, an invasive vector of human diseases. *Heredity* **117**, 125–134 (2016).

40. Western, D. Human-modified ecosystems and future evolution. *PNAS* **98**, 5458–5465 (2001).

41. Pech-May, A. et al. Population genetics and ecological niche of invasive *Aedes albopictus* in Mexico. *Acta Trop.* **157**, 30–41 (2016).

42. Vargo, E. L. et al. Hierarchical genetic analysis of German cockroach (*Blattella germanica*) populations from within buildings to across continents. *PLoS ONE* **9**, e102321 (2014).

43. von Beeren, C., Stoeckle, M. Y., Xia, J., Burke, G. & Kronauer, D. J. Interbreeding among deeply divergent mitochondrial lineages in the American cockroach (*Periplaneta americana*). *Sci. Rep.* **5**, 1–7 (2015).

44. Tseng, S.-P. et al. Genetic diversity and *Wolbachia* infection patterns in a globally distributed invasive ant. *Front. Genet.* **10**, 838 (2019).

45. Wesson, D. M., Porter, C. H. & Collins, F. H. Sequence and secondary structure comparisons of ITS rDNA in mosquitoes (Diptera: Culicidae). *Mol. Phylogen. Evol.* **1**, 253–269 (1992).

46. Mishra, S., Sharma, G., Das, M. K., Pande, V. & Singh, O. P. Intragenomic sequence variations in the second internal transcribed spacer (ITS2) ribosomal DNA of the malaria vector *Anopheles stephensi*. *PLoS ONE* **16**, e0253173 (2021).

47. Artigas, P. et al. *Aedes albopictus* diversity and relationships in south-western Europe and Brazil by rDNA/mtDNA and phenotypic analyses: ITS-2, a useful marker for spread studies. *Parasit Vectors* **14**, 1–23 (2021).

48. Armbruster, P. et al. 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 (2003).

49. Maia, R., Scarpassa, V. M., Maciel-Litaiff, L. & Tadei, W. P. Reduced levels of genetic variation in *Aedes albopictus* (Diptera: Culicidae) from Manaus, Amazonas State, Brazil, based on analysis of the mitochondrial DNA ND5 gene. *Gen. Mol. Res.* **2000**, 998–1007 (2009).

50. Birungi, J. & Munstermann, L. E. 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 (2002).

51. Kambhampati, S. & Rai, K. S. Mitochondrial DNA variation within and among populations of the mosquito *Aedes albopictus*. *Genome* **34**, 288–292 (1991).

52. Werren, J. H., Baldo, L. & Clark, M. E. *Wolbachia*: master manipulators of invertebrate biology. *Nat. Rev. Microbiol.* **6**, 741–751 (2008).

53. Wiwatanaratanabutr, I. Geographic distribution of wolbachial infections in mosquitoes from Thailand. *J. Invertebr. Pathol.* **114**, 337–340 (2013).

54. Carvajal, T. M., Hashimoto, K., Harnandika, R. K., Amalin, D. M. & Watanabe, K. Detection of *Wolbachia* in field-collected *Aedes aegypti* mosquitoes in metropolitan Manila, Philippines. *Parasit. Vectors* **12**, 1–9 (2019).

55. Atyame, C. M., Delsuc, F., Pasteur, N., Weill, M. & Duron, O. Diversification of *Wolbachia* endosymbiont in the *Culex pipiens* mosquito. *Mol. Biol. Evol.* **28**, 2761–2772 (2011).

56. Damiani, C. et al. *Wolbachia* in *Aedes koreicus*: Rare detections and possible implications. *Insects* **13**, 216 (2022).
57. Jiggins, F. M. Male-killing *Wolbachia* and mitochondrial DNA: Selective sweeps, hybrid introgression and parasite population dynamics. *Genetics* **164**, 5–12 (2003).

58. Schuler, H. et al. The hitchhiker's guide to Europe: The infection dynamics of an ongoing *Wolbachia* invasion and mitochondrial selective sweep in *Rhagoletis cerasi*. *Mol. Ecol*. **25**, 1595–1609 (2016).

59. Ross, P. A., Ritchie, S. A., Axford, J. K. & Hoffmann, A. A. Loss of cytoplasmic incompatibility in *Wolbachia*-infected *Aedes aegypti* under field conditions. *PLoS Negl. Trop. Dis*. **13**, e0007357 (2019).

60. Avise, J. C. *Phylogeography: The history and formation of species* (Harvard University Press, 2000).

61. Rokas, A., Atkinson, R. J., Brown, G. S., West, S. A. & Stone, G. N. Understanding patterns of genetic diversity in the oak gallwasp *Biorhiza pallida*: Demographic history or a *Wolbachia* selective sweep? *Heredity* **87**, 294–304 (2001).

62. Porretta, D., Mastrantonio, V., Bellini, R., Somboon, P. & Urbanelli, S. Glacial history of a modern invader: Phylogeography and species distribution modelling of the Asian tiger mosquito *Aedes albopictus*. *PLoS ONE* **7**, e44515. https://doi.org/10.1371/journal.pone.0044515 (2012).

63. Motoki, M. T. et al. Population genetics of *Aedes albopictus* (Diptera: Culicidae) in its native range in Lao People's Democratic Republic. *Parasit. Vectors* **12**, 1–12 (2019).

64. Zhong, D. et al. Genetic analysis of invasive *Aedes albopictus* populations in Los Angeles County, California and its potential public health impact. *PLoS ONE* **8**, e68586 (2013).

65. Usmani-Brown, S., Cohnstaedt, L. & Munstermann, L. E. Population genetics of *Aedes albopictus* (Diptera: Culicidae) invading populations, using mitochondrial nicotinamide adenine dinucleotide dehydrogenase subunit 5 sequences. *Ann. Entomol. Soc. Am.* **102**, 144–150 (2009).

66. Mousson, L. et al. Phylogeography of *Aedes* (*Stegomyia*) *aegypti* (L.) and *Aedes* (*Stegomyia*) *albopictus* (Skuse) (Diptera: Culicidae) based on mitochondrial DNA variations. *Genet. Res.* **86**, 1–11 (2005).

67. Bazin, E., Glémin, S. & Galtier, N. Population size does not influence mitochondrial genetic diversity in animals. *Science* **312**, 570–572. https://doi.org/10.1126/science.1122033 (2006).

68. Dowling, D. K., Friberg, U. & Lindell, J. Evolutionary implications of non-neutral mitochondrial genetic variation. *Ecol. Evol.* **23**, 546–554 (2008).

69. Montero-Pau, J., Gómez, A. & Muñoz, J. Application of an inexpensive and high-throughput genomic DNA extraction method for the molecular ecology of zooplanktonic diapausing eggs. *Limnol. Oceanogr. Methods* **6**, 218–222 (2008).

70. Porter, C. H. & Collins, F. H. Species-diagnostic differences in a ribosomal DNA internal transcribed spacer from the sibling species *Anopheles freeborni* and *Anopheles hermsi* (Diptera: Culicidae). *Am. J. Trop. Med. Cond.* **45**, 271–279 (1991).

71. Folmer, O., Black, M., Hoeh, W., Lutz, R. & Vrijenhoek, R. DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates. *Mol. Mar. Biol. Biotechnol.* **3**, 294–299 (1994).

72. Prosser, S., Martínez-Arce, A. & Elias-Gutiérrez, M. A new set of primers for COI amplification from freshwater microcrustaceans. *Mol. Ecol. Resour.* **13**, 1151–1155 (2013).

73. Ivanova, N. V., Zemlak, T. S., Hanner, R. H. & Hebert, P. D. Universal primer cocktails for fish DNA barcoding. *Mol. Ecol. Notes* **7**, 544–548 (2007).

74. Kumar, S., Stecher, G. & Tamura, K. MEGA7: Molecular evolutionary genetics analysis version 7.0 for bigger datasets. *Mol. Biol. Evol.* **33**, 1870–1874 (2016).

75. Zhou, W., Rousset, F. & O'Neill, S. Phylogeny and PCR-based classification of Wolbachia strains using wsp gene sequences. *Proc. R Soc. Lond. Ser. B Biol. Sci.* **265**, 509–515 (1998).

76. Braig, H. R., Zhou, W., Dobson, S. L. & O'Neill, S. L. Cloning and characterization of a gene encoding the major surface protein of the bacterial endosymbiont *Wolbachia pipiensis*. *J. Bacteriol.* **180**, 2373–2378 (1998).

77. Hu, Y. et al. Identification and molecular characterization of *Wolbachia* strains in natural populations of *Aedes albopictus* in China. *Parasit. Vectors* **13**, 1–14 (2020).

78. Heddi, A., Grenier, A.-M., Khatchadourian, C., Charles, H. & Nardon, P. Four intracellular genomes direct weevil biology: Nuclear, mitochondrial, principal endosymbiont, and *Wolbachia*. *PNAS* **96**, 6814–6819 (1999).

79. Rozas, J. et al. DnaSP 6: DNA sequence polymorphism analysis of large data sets. *Mol. Biol. Evol.* **34**, 3299–3302 (2017).

80. Salzburger, W., Ewing, G. B. & Von Haeseler, A. The performance of phylogenetic algorithms in estimating haplotype genealogies with migration. *Mol. Ecol.* **20**, 1952–1963 (2011).

81. Stamatakis, A. RAxML-VI-HPC: Maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models. *Bioinformatics* **22**, 2688–2690 (2006).

82. Darriba, D., Taboada, G. L., Doallo, R. & Posada, D. JModelTest 2: More models, new heuristics and parallel computing. *Nat. Methods* **9**, 772. https://doi.org/10.1038/nmeth.2109 (2012).

83. Gower, J. C. Some distance properties of latent root and vector methods used in multivariate analysis. *Biometrika* **53**, 325–338 (1966).

84. R Core Team. *R: A Language and Environment for Statistical Computing*. Vienna, Austria: R Foundation for Statistical Computing. http://www.R-project.org/. (2021).

85. Hijmans, R. J., Williams, E., Vennes, C. & Hijmans, M. R. J. Package 'geosphere'. Spher. Trigon. **1**, 5 (2017).

86. Palmer, J. R. et al. Citizen science provides a reliable and scalable tool to track disease-carrying mosquitoes. *Nat. Commun.* **8**, 1–13 (2017).

87. Mantel, N. & Valand, R. S. A technique of nonparametric multivariate analysis. *Biometrics* **1970**, 547–558 (1970).

88. Goslee, S. C. & Urban, D. L. The ecosist package for dissimilarity-based analysis of ecological data. *J. Stat. Softw.* **22**, 1–19 (2007).

89. Stewart, C. Zero-inflated beta distribution for modeling the proportions in quantitative fatty acid signature analysis. *J. Appl. Stat.* **40**, 985–992 (2013).

90. Figueroa-Zúñiga, J. I., Arellano-Valle, R. B. & Ferrari, S. L. Mixed beta regression: A Bayesian perspective. *Comput. Stat. Data Anal.* **61**, 137–147 (2013).

91. Branscum, A. J., Johnson, W. O. & Thurmond, M. C. Bayesian beta regression: Applications to household expenditure data and genetic distance between foot-and-mouth disease viruses. *Aust. N. Z. J. Stat.* **49**, 287–301 (2007).

92. Ospina, R. & Ferrari, S. L. Inflated beta distributions. *Stat. Pap.* **51**, 111–126 (2010).

93. Chung, H. & Beretvas, S. N. The impact of ignoring multiple membership data structures in multilevel models. *Br. J. Math. Stat. Psychol.* **65**, 185–200 (2012).

## Acknowledgements

We would like to thank those citizen scientists participating in the *Mosquito Alert* project who captured and sent us adult tiger mosquitoes for this study. We are grateful to ReNED (Red Nacional de Entomólogos Digitales) for supporting this work. We are also thankful to Victor Ojeda, Sandra Serra and Marc Pradell for their valuable help at the initial stage of the study. The research leading to these results has received funding from the Spanish Ministry of Economy and Competitiveness (MINECO, Plan Estatal I+D+I CGL2013-43139-R), "la Caixa" Foundation (ID 100010434) under agreement HR18-00336, and the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No. 853271).
**Author contributions**

J.R.B.P., F.B. and M.V. conceived and designed the study. S.D., J.R.B.P., A.O., C.P-E., S.M., S.E., D.R., F.C., M.B., T.M., J.A.D., R.E., J.L., and F.B. collected the samples. F.L., S.D., J.C., C.P-E. and A.A.T. analysed samples and data, under the supervision of J.R.B.P., F.B. and M.V. F.L., S.D., J.R.B.P., J.C. and C.P-E. wrote the first draft, S.M., S.E., D.R., M.B., T.M., J.A.D., R.E., J.L., A.A.T., F.B. and M.V. improved successive versions. All authors read and approved the final manuscript.

**Competing interests**

The authors declare no competing interests.

**Additional information**

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1038/s41598-022-24963-3.

**Correspondence** and requests for materials should be addressed to F.L.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_10.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2022