Multiple introductions and overwintering shape the progressive invasion of Aedes albopictus beyond the Alps.

Aedes albopictus originates from Southeast Asia and is considered one of the most invasive species globally. This mosquito is a nuisance and a disease vector of significant public health relevance. In Europe, Ae. albopictus is firmly established and widespread south of the Alps, a mountain range that forms a formidable biogeographic barrier to many organisms. Recent reports of Ae. albopictus north of the Alps raise questions of (1) the origins of its recent invasion, and (2) if this mosquito has established overwintering populations north of the Alps. To answer these questions, we analyzed population genomic data from >4000 genome-wide SNPs obtained through double-digest restriction site-associated DNA sequencing. We collected SNP data from specimens from six sites in Switzerland, north and south of the Alps, and analyzed them together with specimens from other 33 European sites, five from the Americas, and five from its Asian native range. At a global level, we detected four genetic clusters with specimens from Indonesia, Brazil, and Japan as the most differentiated, whereas specimens from Europe, Hong Kong, and USA largely overlapped. Across the Alps, we detected a weak genetic structure and high levels of genetic admixture, supporting a scenario of rapid and human-aided dispersal along transportation routes. While the genetic pattern suggests frequent re-introductions into Switzerland from Italian sources, the recovery of a pair of full siblings in two consecutive years in Strasbourg, France, suggests the presence of an overwintering population north of the Alps. The suggestion of overwintering populations of Ae. albopictus north of the Alps and the expansion patterns identified points to an increased risk of further northward expansion and the need for increased surveillance of mosquito populations in Northern Europe.

INTRODUCTION

Reconstructing the history of biological invasions is fundamental to understand the evolutionary and ecological processes underlying successful invasions (Estoup & Guillemaud, 2010). The genetic structure of invasive populations reflects their introduction history, which includes their geographic origin, the number of introduction events (i.e., propagule pressure), and the number of individuals initiating the invasion (Garnas et al., 2016; Lockwood et al., 2005). A lack of genetic variation is expected in invading populations as the founder populations are often constituted by a limited number of individuals and experience pronounced genetic drift (Dlugosch & Parker, 2008). If the genetic variation of the founder population is too low, it may not be able to establish in a new environment and thus, it will disappear eventually (Facon et al., 2006). Indeed, previous findings suggest that successful biological invasions often originate from multiple rather than single introduction events (Dlugosch & Parker, 2008; Lockwood et al., 2005). Multiple introductions contribute to maintaining high genetic diversity and population size of the invading populations (Cristescu, 2015). This is especially true if introductions originate from geographically distant sources, as it increases the probability of introducing individuals with different genetic backgrounds (Rius & Darling, 2014). Once established, connectivity among introduced populations can additionally lead to admixture that further increases genetic variation, and this, in turn, may increase the probability of successful establishment and, ultimately, further spread (Slatkin, 1985).

For invasive vector species, knowledge of their dispersal dynamics, source populations, and introduction pathways is not only of academic interest but also of immediate relevance for public health. Understanding the invasion history allows better estimates of the risk of establishment of new populations, and thus provides important information for monitoring and control (Estoup & Guillemaud, 2010). A great example of a successful biological invader is Aedes albopictus (Skuse, 1894), the Asian tiger mosquito. It is considered one of the most invasive species worldwide (Global Invasive Species Database, 2020). Due to its vector competence for several arboviruses, including chikungunya, dengue, and Zika (Gratz, 2004; Wong et al., 2013), as well as dirofilarial worms (Cancrini et al., 2003), Ae. albopictus is of particular public health concern.

Aedes albopictus eggs can resist desiccation for long periods and overcome lower temperatures during winter in temperate zones through diapause (Hanson & Craig, 1994). These biological factors greatly facilitated the global expansion of this mosquito species together with human activities, which contributed to its expansion by enabling dispersal over long and shorter distances. Like other invasive Aedes species, Ae. albopictus is passively spread across continents primarily through the international trade of used tires into which mosquitoes had deposited eggs before shipment (Paupy et al., 2009). At the regional level, adult mosquitoes frequently hitch ride in vehicles and are subsequently displaced along roads (Egizi et al., 2016; Medlock et al., 2015).

Over the last four decades, Ae. albopictus has spread to every continent except Antarctica, while its native distribution range is in Southeast Asia, from tropical (e.g., Indian Ocean Islands and Indonesia; Bonizzoni et al., 2013) to temperate regions (Japan; Kobayashi et al., 2002). In mainland Europe, Ae. albopictus was first recorded in Albania in 1979 (Adhami & Reiter, 1998). In the Americas, it was first reported from Texas, USA, in 1985 (Sprenger & Wuithiranyagool, 1986) and 1 year later from the State of Rio de Janeiro, Brazil (Oswaldo, 1986). The populations in North America are considered to have served as bridgehead populations for secondary introductions into Europe (Garnas et al., 2016; Lombaert et al., 2010) at two sites in Northern Italy between 1990 and 1991 (Dalla Pozza & Majori, 1992; Sabatini et al., 1990). From there, the mosquito quickly spreads across Southern Europe (Sherpa et al., 2019).

To date, Ae. albopictus has firmly established across the Mediterranean region from Spain to Greece (ECDC, 2019), and from the sea to the foot of the Alps (Flacio et al., 2016). In addition, modeling studies, considering present and future climatic conditions, suggest that its range will be expanding even further north (Caminade et al., 2012; Kraemer et al., 2019). Indeed, isolated populations of Ae. albopictus have already been reported from north of the Alps in Southern Germany (Becker et al., 2013; Pluskota et al., 2008; Werner et al., 2012) and northern Switzerland (Biebinger, 2020) with mosquitoes frequently re‐introduced across the Alps along the highways from south to north (Fuehrer et al., 2020; Müller et al., 2020). Given the very patchy pattern of the reported Ae. albopictus populations and the uncertainties of the climatic models, the extent to which local populations north of the Alps are actually self‐sustainable, rather than temporarily established by re‐introduced individuals, is uncertain and their origins also remain largely unknown.

High‐resolution population genetic markers are fundamental to accurately resolve invasion histories of target species, especially for species with a recent invasion on a fine geographical scale like Ae. albopictus (Cristescu, 2015). Previous studies attributed difficulties to reconstruct invasion histories to low resolution of genetic markers, such as mitochondrial DNA or microsatellites (Goubert et al., 2016; Manni et al., 2017). Genomic analysis based on the screening of thousands of genome‐wide single nucleotide polymorphisms (SNPs) using double‐digest restriction site‐associated DNA sequencing (ddRAD‐seq) allows for high‐resolution studies, enabling detection of patterns and levels of genetic differentiation for Ae. albopictus at different spatial resolutions ranging from global (Kotsakiozi et al., 2017), to continental (Pichler et al., 2019; Sherpa et al., 2019), and to city scales (Schmidt et al., 2017) studies. Here, we aimed at a higher resolution by using ddRADseq to identify a panel of 4000 SNPs to investigate the introduction of Ae. albopictus into Switzerland, to reconstruct the invasion history across the Alps and to evaluate if current populations are self‐sustained. The spatial scale of this study is about 300 km along the south–north axis across the Alps. To facilitate detection of both long‐ and short‐range dispersal events, we screened for genomic variations in specimens from six sites in Switzerland north and south of the Alps, 33 sites in Europe, 5 sites from the Americas, and 5 sites from its Asian native range. To evaluate temporal stability of the Ae. albopictus populations north of the Alps, we screened for variation in three population samples collected over two consecutive years.

MATERIALS AND METHODS

Sampling strategy

The sampling locations are reported in Table 1 and Figure 1, and all details on collection sites, time points, and methods are reported in Appendix S1 (Table S1). First, we investigated long‐range migration using a dataset consisting of 208 individuals from the native and invasive range (dataset named 1.native_invasive). Second, we assess dispersal at the European scale and genetic structuring across the Alps using a dataset consisting of a subset of 137 individuals, which included only European samples from 39 sites (dataset named 2.europe; Table 1 and Figure 1).

TABLE 1

Aedes albopictus specimens included in the two datasets of the present study

Dataset		Country	Collection sites (N)	N ind
1.native_invasive	2.europe	Switzerland (CH)	6	56
		Italy (IT)	15	34
		France (FR)	5	20
		Germany (DE)	3	9
		Liechtenstein (FL)	1	1
		Albania (AL)	4	8
		Greece (GR)	5	9
		Brazil (BR)	4	17
		USA (US)	6	6
		Japan (JP) a	1	11
		Indonesia (ID) a	12	14
		Hong Kong (HK) a	3	23

Note: N ind indicates the number of specimens included in the study prior any data filtering.

These specimens were included at the data analysis stage and are already published ddRAD data (Schmidt, Chung, Honnen, et al., 2020).

FIGURE 1

Aedes albopictus sampling sites. The pie charts represent collection sites, where the size of each pie represents how many individuals were collected in each location. The panels represent the sampling sites at the (a) global, (b) Europe, and (c) Swiss levels.

We used different sets of samples to address different questions. The 1.native_invasive dataset includes the core dataset (2.europe, see below), plus additional samples collected outside of the target study area to facilitate detection of potential long‐range introductions and their origins. The dataset comprises a total of 208 specimens, including 5 populations from the USA and Brazil, as they are considered to be a bridgehead for the European invasion (Battaglia et al., 2016), and 5 populations from the native range, Japan (Matsuyama), Indonesia (Bandung), and China (Hong Kong) (Schmidt, Chung, Honnen, et al., 2020), representing the three major genetic clusters previously detected in this species native range (Kotsakiozi et al., 2017; Sherpa et al., 2019). Table 1 details the new samples analyzed for this study and the ones with existing ddRAD data obtained from another already published study (Schmidt, Chung, Honnen, et al., 2020).

The samples collected specifically for this study constitute the core dataset (2.europe) and were collected during summer months in 2006, 2016, 2017, and 2018. This dataset includes samples from across the Alps in Switzerland, neighboring countries (Germany, France, Liechtenstein, and Italy), and from Albania and Greece. Samples from across the Alps included collections made along national highways around and in the city of Basel in Northern Switzerland (Figure 1c, CH‐North) and collections made in Ticino in southern Switzerland (Figure 1c, CH‐South). To allow assessment of overwintering ability and the presence of self‐sustainable populations, we also included collections made at multiple time points (2017 and 2018) from four locations across the Alps: Strasbourg in France, and Mendrisio, Luzern, and Basel in Switzerland. Samples from neighboring countries included Southern Germany (Baden‐Württemberg), eastern France (Haut‐Rhin), Liechtenstein, and 15 locations across Italy. We sampled Italy more comprehensively than the other neighboring countries because it is considered the most likely source of introduction into Switzerland (Sherpa et al., 2019, Figure 1b). We also included collections from Albania and Greece to provide larger geographic context and to include the region where Ae. albopictus was first reported in Europe (Adhami & Reiter, 1998).

The 160 specimens that did not have published ddRAD data available (see Table 1 for details on new vs. already published data) were collected as adults for DNA extraction of the full body or were collected as eggs or larvae and reared to the adult stage before DNA extraction. Larvae were caught by dipping into standing water and eggs were collected with ovitraps (for trap design, see Flacio et al., 2016). Larvae and eggs were reared to adults in a HPP110 constant climate chamber (Memmert GmbH + Co. KG, Schwabach, Germany), mimicking summer temperature and humidity regimes in southern Switzerland. To avoid sampling of siblings, one individual per dip/ovitrap was used. Adults were collected in Biogents Sentinel version 1 traps (Biogents AG, Regensburg, Germany), sent to us through citizen reports, or incidentally caught by the authors as human landing catches. Upon collection, samples were stored in 80% ethanol at 4°C until further processing.

In all analyses, mosquitoes collected from within the same city and in the same year were considered to be from one population. The minimum and maximum distances of samples within the same city are reported in the Appendix S1 (Table S1). This grouping of mosquitoes into populations for analysis is supported by previous estimates of Ae. albopictus dispersal that indicate highly localized and restricted active dispersal distances within urban areas (Vavassori et al., 2019). There is evidence that the sample sizes used are adequate because previous studies indicated that with >1000 SNPs, as few as two individuals per population provide adequate resolution to assess genetic differentiation and evolutionary relationships (Kotsakiozi et al., 2017; Nazareno et al., 2017; Willing et al., 2012). In our study, for some locations, only one individual was available (Table 1).

DNA extraction and ddRAD library construction

We extracted total genomic DNA from 160 individual mosquitoes—137 specimens in the core dataset (2.europe), as well as 23 specimens in the global dataset (1.native_invasive) that did not have published ddRAD data available (6 specimens from USA, 17 from Brazil, see Table 1 for details on new vs. already published data). DNA was extracted from adult mosquito specimens, using the High Pure PCR Template Preparation Kit (Roche, Rotkreuz, Switzerland), following the manufacturer's protocol with an additional step of RNAse treatment. DNA amounts were quantified with a Qubit 2.0 Fluorometer (Invitrogen, Carlsbad, CA, USA).

We constructed the ddRAD libraries following the protocol for Ae. albopictus described in Schmidt, Chung, Honnen, et al. (2020) and Schmidt, Chung, Van Rooyen, et al. (2020), an adaptation of the original protocol of Rašić et al. (2014). DNA was digested with the restriction enzymes NlaIII and MluCI (New England Biolabs, Beverly MA, USA). The size selection step targeted a fragment size between 350 and 450 bp. We allocated individuals from the same collection sites randomly across libraries (O'Leary et al., 2018).

We sequenced these same 160 individual mosquitoes—with pools of barcoded DNA of 56 specimens per library on an Illumina HiSeq 2500 system (Illumina Inc., San Diego, CA, USA) at the Department of Biosystems Science and Engineering (D‐BSSE), ETH Basel, Switzerland, using paired‐end, HiSeq Flow Cell v4, HiSeq SBS Kit v4, and with a 10% PHiX spike.

Data processing and SNP genotyping

We used the process_radtags function in STACKS v2.2 (Catchen et al., 2013) to de‐multiplex the raw reads and mapped them to the Ae. albopictus reference genome (Accession number: GCA_006516635.1) available on NCBI GenBank (Palatini et al., 2020) using the BWA‐MEM algorithm implemented in the Burrows‐Wheeler Aligner tool BWA v0.7.17 (Li & Durbin, 2009), allowing up to four mismatches. For SNP calling, we used the ref_map.pl wrapper in STACKS. The VCF file output was used to filter the data for sequencing and SNP call quality. Using VCFtools v1.9 (Danecek et al., 2011) and R version 4.0.3 (R Core Team, 2020), we excluded loci that mapped to repetitive regions of the genome, had more than 50% missing data, or did not exhibit allele balance. We included only bi‐allelic variants, with a maximum mean depth value of 30 and with a minimum allele count of three.

We used plink v1.9 (Chang et al., 2015) to include only individuals with less than 20% missing genotypes and a genotyping rate greater than 80% in iterative steps for the 1.native_invasive and 2.europe datasets, independently. We excluded tags with more than 10 SNPs and used the populations function in STACKS to obtain output files in VCF format. Since most of the downstream analyses require that SNPs are unlinked, we removed linked sites by excluding SNPs located within a window of 400 bp (i.e., option ‐‐thin 400) with VCFtools. The window size corresponded to our maximum fragment size, thus each SNP belongs to a single DNA fragment. After conducting a relatedness analysis, we excluded one individual per sibling pair from the analyses (see section below). The reduced dataset was split into two cleaned datasets: 1.native_invasive_cleaned, including 153 samples and 4714 loci and SNPs, and 2.europe_cleaned, including 93 samples and 6308 loci and SNPs (Table 3).

TABLE 3

Details on the datasets used in the study

Dataset	N ind	Loci (N)	SNPs (N)	Missing data (samples ± SD)	Missing data (locus ± SD)	Average read depth per individual ± SD	Average coverage per site ± SD
1.native_invasive	208	4930	23,240	11.1% ± 3.9	11.1% ± 5.8	12.9 ± 7.5	12.9 ± 6.6
1.native_invasive_cleaned	153	4714	4714	11.1% ± 4.0	11.1% ± 5.9	12.9 ± 7.3	13.1 ± 6.8
2.europe	137	9966	9960	12.18% ± 3.7	12.17% ± 6.0	13.1 ± 6.5	13.2 ± 6.7
2.europe_cleaned	93	6308	6308	10.6% ± 3.3	10.6% ± 6	12.2 ± 6.5	13.9 ± 6.6

Relatedness analysis

To exclude closely related individuals that could potentially bias the analysis of population structure, we calculated Loiselle's k (Loiselle et al., 1995), using the program SPAGeDi (Hardy & Vekemans, 2002) for the datasets 1.native_invasive and 2.europe. We identified putative full siblings based on pairwise k values of >0.1875, and putative half‐siblings with values ranging from 0.1875 > k > 0.0938, following Iacchei et al. (2013). The same cutoff values have also been used in a previous study on mosquitoes (Schmidt et al., 2018). In addition to SPAGeDi, we confirmed the putative relationships between individuals with two additional approaches. First, we confirmed relatedness analysis with the ‐‐relatedness2 flag of VCFtools (Danecek et al., 2011) based on the KING inference (Manichaikul et al., 2010) and selected only pairs of siblings identified by both SPAGeDI and VCFtools. Second, we used the software program ML‐Relate (Kalinowski et al., 2006) to confirm putative relationships as described in Schmidt et al. (2018). We run two specific hypotheses of putative relationships: we ran a first “standard” test assuming that the kinship category assigned using Loiselle's k was more likely than the next most likely kinship category. Second, we run a “conservative” test that assumed that the kinship category assigned using Loiselle's k was less likely to be correct. Thus, for pairs with k > 0.1875, statistical tests run with ML‐Relate would determine whether the identified pair was full siblings or half‐siblings, while for pairs with 0.1875 > k > 0.09375, tests would help determine whether the identified pair was full siblings, half‐siblings, or unrelated. Conservative and standard tests were run using 10,000 simulations of random genotype pairs.

Genetic structure

To assess population structure, we employed both model‐free and model‐based approaches. First, we employed the model‐free Principal Components Analysis (PCA) and Discriminant Analysis of Principal Components (DAPC) (Jombart et al., 2010) on the 1.native_invasive_cleaned and 2.europe_cleaned dataset, using the adegenet v2.0 package in R (Jombart & Ahmed, 2011). PCA is a multivariate analysis used to identify genetic clusters, without an assumption about the underlying population genetics model. The DAPC analysis maximizes the between group while minimizing the within‐group variance and computes a PCA, followed by a discriminant analysis to identify the number of genetic clusters (Jombart et al., 2010). We used the function find.clusters to estimate the number of clusters K and xvalDapc option to perform cross‐validation and to assess the most likely number of principal components to retain in the DAPC analysis. We used two model‐based approaches. The first approach was the maximum‐likelihood method implemented in the program ADMIXTURE v1.3.0 (Alexander et al., 2009). We used it to conduct ancestry analysis and to estimate the most likely number of evolutionary clusters K on the cleaned datasets (1.native_invasive_cleaned and 2.europe_cleaned). The second approach was to choose the most likely value for K using the ADMIXTURE's cross‐validation procedure. Genetic differentiation was further investigated with fineRADstructure v1.7.20 (Malinsky et al., 2018). This method enables fine‐scale population structure inference by using a Bayesian clustering approach and it has been shown to be especially informative in the case of recent gene flow between mosquito populations (Pichler et al., 2019). For this analysis, we only used the 1.native_invasive dataset because the algorithm takes into account haplotype information and uses all available SNPs allowing for a higher structural resolution (Malinsky et al., 2018).

Genetic differentiation, isolation by distance, and overwintering

To evaluate the degree of genetic differentiation, we estimated pairwise F st values (Weir & Cockerham, 1984) at country level on the dataset 1.native_invasive_cleaned with the R package HierFstat v0.5–10 (Goudet, 2005) and estimated the corresponding 95% confidence intervals by performing 1000 bootstraps over all loci. Next, we calculated individual inbreeding coefficients (F IS) by assessing their statistical significance with 1000 bootstrap samplings and estimated allelic richness (AR). We calculated the observed heterozygosity (H ) using VCFtools (−het). Individuals were grouped by country and differences in mean H between groups were tested with a non‐parametric Kruskal–Wallis (KW) test for statistical significances in R. Expected heterozygosity (H ) per country was computed using the R package adegenet v2.0.

To investigate patterns of genetic diversity, genotype frequency, and genetic differentiation across Europe (2.europe_cleaned dataset), we estimated F ST, F IS, AR, H , and H using the R package adegenet v2.0, F IS and pairwise F ST using the R package HierFstat (Goudet, 2005) and pairwise proportion of shared alleles (D ps) using the R package adegenet (propShared function). We visualized pairwise D ps with neighbor‐net networks with the software SplitsTree v5.0 (Huson & Bryant, 2006). To assess the impact of geographic distance on genetic differentiation, we performed a test for isolation by distance (IBD) with a Mantel test (Mantel, 1967) with 1000 permutations, using D ps and log‐transformed geographic distances as the input, r = 0 as the null hypothesis, and r > 0 as the alternative hypothesis. We used pairwise D ps rather than F ST in the test for IBD because this metric provides improved power to detect IBD at small geographical scales, with small genetic distances as expected due to the recent invasion history, high dispersal, and small sample sizes (Bowcock et al., 1994) (Shirk et al., 2017), which is characteristic for the 2.europe_cleaned dataset.

To assess overwintering ability and assess the presence of self‐sustainable populations, we estimated pairwise F ST values between individuals from Mendrisio, Luzern, and Basel that had been collected during two consecutive mosquito seasons in 2017 and 2018. To evaluate if Ae. albopictus overwinters in our study area, we compared the F ST values from the temporal comparison with the F ST values calculated among individuals from geographically distant locations. If overwintering does occur, the samples collected in the two different years would belong to the same population and, therefore, their F ST values should be considerably smaller than across geographically distant populations.

Genetic assignment test

In order to identify possible source population of the mosquitoes in Switzerland, we performed genetic assignment test with the R program assignPOP (Chen et al., 2018). We assigned individuals collected in France and in Northern Italy as source populations, considering their geographical proximity to Switzerland. We tested assignment accuracies via Monte Carlo cross‐validation based on the following parameters: proportion of individuals used in training set: 0.5, 0.7, and 0.9; proportion of loci used in training set: 0.25,0.5, and 1 and loci sample method FST; iterations: 30; and model: support vector machine.

RESULTS

SNP discovery

We sequenced 160 individuals obtaining a total of 828 million reads, with 5 million reads per sample on average, ranging from 7 thousand to 17 million. After filtering and removal of duplicate siblings, the dataset 1.native_invasive_cleaned included 153 individuals and 4714 SNPs and loci. The dataset 2.europe_cleaned included 93 samples and 6308 SNPs. An average of 3 million reads (73%) per individual aligned to the reference genome. Table 3 shows the details on the number of reads, individuals, SNPs, coverage, and level of missing data of each dataset.

We identified 15 full‐ and 10 half‐sibling pairs from the same collection sites (Table 4). One sibling from each pair was kept for all downstream analysis (i.e., 1.native_invasive_cleaned and 2.europe_cleaned). Although we had only 15 full siblings in our dataset, we found two siblings from the same location in two consecutive years. These were detected at a distance of 330 m in 2017 and 2018 in Strasbourg, France. The first sibling was collected from an egg sample in fall 2017, while the second one was caught as an adult during the second half of July 2018.

TABLE 4

List of individuals and their relative kinship identified in the 1.native_invasive and 2.europe dataset in the relatedness analysis

Individual ID 1	Individual ID 2	Loiselle k (SPAGeDi)	Kinship (VCFtools)	Kinship (ML‐Relate)
CHCAS2	CHCAS3	0.513879 (FS)	FS	FS
CHBE121	CHBE122	0.512615 (FS)	FS	FS
DEHD2	DEHD3	0.509926 (FS)	FS	FS
CHPR1	CHPR2	0.492123 (FS)	FS	FS
FRMA11	FRMA9	0.473704 (FS)	FS	FS
CHBE121	CHBE129	0.471598 (FS)	FS	FS
CHBE122	CHBE129	0.470286 (FS)	FS	FS
CHNE1	CHNE5	0.453771 (FS)	FS	FS
FRST10	FRST1	0.445381 (FS)	FS	FS
CHNE2	CHNE5	0.325529 (FS)	FS	FS
IPA4	IPA5	0.318006 (FS)	FS	FS
AAM6	AAM8	0.312713 (FS)	FS	FS
CHNE2	CHNE4	0.306556 (FS)	FS	FS
IBA2	IBA3	0.301730 (FS)	FS	FS
CHNE4	CHNE5	0.286318 (FS)	FS	FS
FRST1	FRST4	0.194805 (FS)	HS	HS
FRST4	FRST9	0.188862 (FS)	HS	HS
FRST10	FRST4	0.179891 (HS)	HS	HS
FRST1	FRST9	0.165271 (HS)	HS	HS
CHBA6	FRSL1	0.162914 (HS)	HS	HS
FRST4	FRST5	0.159820 (HS)	HS	HS
AAM3	AAM4	0.148871 (HS)	HS	HS
FRCO3	FRCO4	0.145682 (HS)	HS	HS
CHNE1	CHNE2	0.136978 (HS)	HS	HS
FRCO2	FRCO4	0.120723 (HS)	HS	HS
FRST10	FRST5	0.090178 (UR)	UR	UR
CHCO11	CHCO15	0.088116 (UR)	UR	UR
CHCAS1	CHCAS6	0.084865 (UR)	UR	UR
FRGR7	FRGR9	0.083812 (UR)	UR	UR
FRGR3	FRGR8	0.082819 (UR)	UR	UR
IBA3	IBA5	0.079113 (UR)	UR	UR
CHLU2	CHLU3	0.071667 (UR)	UR	UR
CHCAS3	CHCAS8	0.065919 (UR)	UR	UR
CHCO14	CHCO6	0.062289 (UR)	UR	UR
IC14	IC7	0.061200 (UR)	UR	UR
FRCO2	FRCO3	0.060119 (UR)	UR	UR
FRGR2	FRGR8	0.056716 (UR)	UR	UR
IVA2	IVA6	0.055244 (UR)	UR	UR

Note: Kinship was determined by three different methods, first using SPAGeDi where pairs of k > 0.1875 are identified as full siblings (FS), those of 0.0938 < k < 0.1875 as half‐siblings (HS) and those k < 0.0938 as unrelated individuals (UR), second using VCFtools (flag ‐relatedness2) based on a relationship inference algorithm and lastly ML‐Relate to run specific hypothesis tests of putative relationships assigned by SPAGeDi. We added to this table the values of 14 randomly selected individuals from the same population and their respective Loiselle K values. Sample information for each individual ID listed is reported in the Appendix S1.

At the global level, the DAPC analysis separated the specimens in four main clusters (1.native_invasive_cleaned dataset Figure 2b, Table 5): cluster 1 (gray) consists of the mosquitoes collected in Brazil and one specimen from Sicily (Sample ID: IMS3, Messina); cluster 2 (green) includes specimens collected in Europe, Hong Kong, and USA; clusters 3 (pink) and 4 (blue) comprise the specimens from Indonesia and Japan, respectively. A similar clustering was also recovered by the PCA analysis (Figure 2a): specimens from Indonesia, Brazil, and Japan were the most differentiated, whereas specimens from Europe, Hong Kong, and USA largely overlapped. The ADMIXTURE cross‐validation analysis supports the presence of two genetic clusters in the 1.native_invasive_cleaned dataset with individuals collected in Indonesia forming one cluster genetically distinct from all the other individuals (Figure 3 K = 2). However, as the cross‐validation errors were very similar for a number of K's (K = 1: 0.216, K = 2: 0.211, K = 3: 0.219, K = 4: 0.232, K = 5: 0.247), we further explored K = 2, 3, 4, and 5. K = 4 appears to be the most likely scenario with two clearly distinct genetic groups, including specimens collected in Indonesia and Hong Kong, and two genetic groups with a high level of genetic admixture, including specimens collected in Europe, Brazil, the USA, and Japan (Figure 3). The fineRADstructure analysis confirms the presence of four genetic clusters at the global scale (Figure 4).

FIGURE 2

Genetically distinct clusters of Ae. albopictus sampled populations. Four genetic clusters may be observed in the global dataset (1.native_invasive_cleaned) based on the PCA and DAPC analysis. (a) PCA with the percent variation explained by the first two principal components (b) Scatterplot showing the results from the DAPC (K = 4).

TABLE 5

Composition of the DAPC groups obtained for the 1.native_invasive_cleaned dataset (Figure 2b)

Country	Cluster 1	Cluster 2	Cluster 3	Cluster 4
Indonesia (ID)	0	0	14	0
Hong Kong (HK)	0	23	0	0
Japan (JA)	0	0	0	11
USA (US)	0	6	0	0
Brazil (BR)	6	0	0	0
Albania (AL)	0	6	0	0
Greece (GR)	0	8	0	0
Italy (IT)	1	22	0	0
Switzerland (CH)	0	39	0	0
France (FR)	0	13	0	0
Germany (DE)	0	3	0	0
Liechtenstein (FL)	0	1	0	0

FIGURE 3

Aedes albopictus ADMIXTURE barplot for all mosquito populations based on the results from the dataset 1.native_invasive_cleaned. Each bar represents one individual, while white vertical lines indicate separate countries. Country codes: AL: Albania, AU: Austria, BR: Brazil, CH: Switzerland, HK: Hong Kong, DE: Germany, FR: France, GR: Greece, ID: Indonesia, FL: Liechtenstein, IT: Italy, JP: Japan and US: USA.

FIGURE 4

Output of the fineRADstructure analysis of the 1.native_invasive dataset. The heat map indicates pairwise co‐ancestry between individuals, with black, blue, and purple representing the highest levels, red and orange indicating intermediate levels, and yellow representing the lowest levels of shared co‐ancestry. The tree on top of the heat map shows the inferred relationships between the specimens analyzed, with each tip corresponding to an individual. On the Y‐axis, country of origin with their sample collection site ID is reported if they create distinct clusters, otherwise are included in the Europe (rest) cluster. On the X‐axis, sample codes are encoded with their laboratory ID (see Appendix S1: Table S1). Siblings are depicted in black and blue colors.

At the European level, different analysis suggests a weak substructuring. In the ADMIXTURE analysis, the individuals collected in Europe, except for some individuals from Greece and Albania, shared high levels of genetic ancestry and showed the highest proportion of their genome being assigned to either cluster 2 or 4 (green and blue clusters in Figure 3). Two specimens collected in Albania and Greece, respectively, were also assigned to cluster 1, together with mosquitoes collected in Hong Kong. The results for other values of K (K = 2, 3, and 5) are reported in Figure 3. The fineRADstructure and the PCA analysis groups together specimens from Albania and Greece, separating them from the rest of the European samples (Figure 4). While the admixture analysis on the European dataset did not detect clear genetic structuring (i.e., optimal K = 1; see Figure 5), the DAPC separates the specimens into three clusters (Figure 6). Cluster 1 (Figure 6a – purple) includes all specimens collected in Albania with two additional specimens from Greece. Cluster 2 (Figure 6a – green) comprises all specimens collected from Northern Italy, except for two individuals from Como and Varese, samples from both southern and northern Switzerland, and one individual from Germany. Cluster 3 (Figure 6a – orange) contains specimens from South and North Italy, Greece, all the specimens from France, and 22 specimens from Switzerland (Table 6).

FIGURE 5

ADMIXTURE barplot obtained for the 2.europe_cleaned dataset for K = 2–3. Individuals represented by vertical bars along the plot grouped by country and collection site. The Y‐axis represents the probability of an individual to be assigned to a genetic cluster. Each cluster is given in a different color. Multi‐colored bars indicate admixed genetic ancestry in the respective individual. The white vertical lines indicate country limits. Country codes: AL: Albania, CH: Switzerland, DE: Germany, FR: France, GR: Greece, IT: Italy.

FIGURE 6

Genetic structure and differentiation of Ae. albopictus specimens collected in Europe. (a) Scatterplot showing the results of the DAPC (K = 3) on the 2.europe_cleaned dataset. Cluster 1‐purple includes mosquitoes collected in Albania with some specimens collected in Greece; cluster 2‐green includes mosquitoes collected in Northern Italy (with the exception of two specimens which clustered with cluster 3 (orange), mosquitoes collected in southern and northern Switzerland and one specimen from Germany. Cluster 3‐orange includes specimens collected in Italy‐Center‐South, Italy‐Sicily, Switzerland, and France. (b) Neighbor‐net network of D ps relative genetic distances among the specimens from Italy. The map shows the locations of the sampling in the region of Italy‐Center‐South and Italy‐Sicily. Specimens collected in Northern Italy are depicted with a green square (cluster 2 ‐ green) and the one collected in Central and Southern Italy with orange circles (cluster 3 ‐ orange). Specimens from the Italian island Sicily are not reported here. For the sample abbreviations, see Appendix S1, Table S1 Laboratory ID.

TABLE 6

Composition of the DAPC groups obtained for the 2.europe_cleaned dataset (Figure 6a)

Country	Cluster 1	Cluster 2	Cluster 3
Albania	5	0	1
Greece	2	1	3
Italy	0	6	16
Switzerland	0	12	27
France	0	0	13
Germany	0	2	1
Liechtenstein	0	0	1

The degree of differentiation between countries detected in the dataset 1.native_invasive_cleaned is low, with pairwise F ST values ranging from 0 to 0.21, with lower values between specimens from Italy, Switzerland, and France, and higher values between specimens from Indonesia and Switzerland (Table 7).

TABLE 7

Pairwise F ST values for the dataset 1.native_invasive_cleaned and the corresponding 95% confidence interval using 1000 bootstrapping in the top diagonal. Pairwise F ST was calculated only for populations with N ind ≥ 2

	AL	BR	FR	DE	GR	IT	CH	US	HK	JA	ID
AL	NA	0.10–0.14	0.03–0.05	−0.01–0.03	0.01–0.03	0.02–0.04	0.03–0.05	0.05–0.07	0.04–0.07	0.07–0.10	0.16–0.20
BR	0.12	NA	0.10–0.14	0.07–0.11	0.09–0.12	0.08–0.12	0.09–0.13	0.10–0.14	0.07–0.11	0.12–0.15	0.16–0.20
FR	0.04	0.12	NA	0.01–0.04	0.01–0.02	0.00–0.01	0.00–0.01	0.04–0.06	0.03–0.04	0.05–0.07	0.17–0.22
DE	0.01	0.09	0.02	NA	0.00–0.03	−0.01–0.02	0.00–0.03	0.02–0.05	0.02–0.05	0.03–0.06	0.12–0.16
GR	0.02	0.10	0.02	0.02	NA	0.01–0.02	0.01–0.02	0.04–0.06	0.02–0.04	0.06–0.08	0.15–0.19
IT	0.03	0.10	0.00	0.01	0.01	NA	0.00–0.00	0.03–0.05	0.03–0.04	0.05–0.07	0.18–0.22
CH	0.04	0.11	0.00	0.01	0.01	0.00	NA	0.03–0.05	0.03–0.04	0.05–0.07	0.19–0.24
US	0.06	0.12	0.05	0.04	0.05	0.04	0.04	NA	0.04–0.06	0.06–0.08	0.16–0.20
HK	0.05	0.09	0.03	0.04	0.03	0.03	0.04	0.05	NA	0.07–0.09	0.16–0.20
JA	0.08	0.13	0.06	0.04	0.07	0.06	0.06	0.07	0.08	NA	0.19–0.23
ID	0.18	0.18	0.19	0.14	0.17	0.20	0.21	0.18	0.18	0.20	NA

Abbreviations: AL, Albania; BR, Brazil; FR, France; DE, Germany; GR, Greece; IT, Italy; CH, Switzerland; US, USA; HK, Hong Kong; JA, Japan; ID, Indonesia.

Observed (H ) and expected (H ) heterozygosity ranged from 0.041 to 0.055 and from 0.048 to 0.075, respectively. H , H , and F IS within country for the dataset 1.native_invasive_cleaned are reported in Table 8. H differed among countries (K‐W H = 41.9, df = 10, p‐value <.001) and between the native and the invasive range (K‐W H = 6.3, df = 1, p‐value <.05). The highest heterozygosity measured within a country was among the Indonesian specimens (Table 8).

TABLE 8

Basic diversity statistics for the dataset 1.native_invasive_cleaned calculated at country level (N ind = 153, N snps = 4714)

Country	N ind	Private allele	Mean H_obs	Mean H_exp	F IS	F IS CI (0.025, 0.975)
Indonesia (ID)	14	643	0.055	0.075	0.288	0.227	0.339
Hong Kong (HK)	23	220	0.044	0.053	0.188	0.148	0.232
Japan (JA)	11	85	0.049	0.060	0.197	0.159	0.253
USA (US)	6	29	0.046	0.052	0.176	0.135	0.244
Brazil (BR)	6	87	0.041	0.051	0.212	0.169	0.293
Albania (AL)	6	15	0.046	0.051	0.116	0.116	0.220
Greece (GR)	8	25	0.046	0.050	0.125	0.009	0.179
Italy (IT)	23	74	0.045	0.054	1.166	0.127	0.205
Switzerland (CH)	39	117	0.046	0.055	0.156	0.132	0.190
France (FR)	13	30	0.042	0.048	0.131	0.096	0.181
Germany (DE)	3	11	0.049	0.049	0.138	0.113	0.237
Liechtenstein (FL)	1	1	–	–	–	–	–

Note: N ind, number of individuals; H_exp, expected heterozygosity; H_obs, Individual observed heterozygosity averaged per geographical regions; F IS, inbreeding coefficient with 95% confidence interval (F IS CI).

To further investigate the dispersion across the Alps and identify the presence of self‐sustaining populations, we grouped specimens from Italy and Switzerland according to the time since their first report of introduction into “long‐established” (i.e., established since 1990) and “recently‐established” (after 2003) populations. With this approach, we identified three groups of populations from Italy (Figure 1b i.e., IT‐North, IT‐Center‐South, and IT‐Sicily) and two groups of recent population from Switzerland (Figure 1c i.e., CH‐North and CH‐South) and compared their genetic diversity.

The recently established populations in Switzerland (CH‐North and CH‐South) did not show lower genetic diversity (i.e., H and H ) than the long‐established Italian populations (IT‐North, IT‐Center‐South, and IT‐Sicily) (K‐W H = 2.81, df = 4, p‐value = .59; Figure 7a). Pairwise F ST between collection sites in Switzerland ranged between 0 and 0.04 (Appendix S1, Table S3). For the Swiss locations, Mendrisio, Luzern, and Basel, we found a decrease of 0.004 in the heterozygosity between the specimens collected in 2018 versus the samples collected in 2017 (K‐W H = 4.59, df = 1, p‐value <.05; Figure 7b) and higher inbreeding coefficients within the three investigated sites (Figure 7c and d). Pairwise F ST values among the specimens collected in the three sites at multiple time points were similar (Mendrisio: 0.006, Luzern: 0, Basel: 0. Appendix S1, Table S3) and considerably smaller than values calculated between geographically distant populations (Appendix S1, Table S3). The specimens from Strasbourg, France, were excluded from this analysis because they were identified as full siblings (see relatedness analysis, Table 4).

FIGURE 7

Genomic diversity of Ae. albopictus collected in Italy and in Switzerland. (a) Individual observed heterozygosity (H_obs) estimated with VCFtools on the 2.europe_cleaned dataset. The individuals were grouped by geographical regions, including three regions in Italy and two regions in Switzerland, and difference in their mean heterozygosity (H O) was tested with the non‐parametric Kruskal–Wallis (KW) test. (b) H_obs in the samples collected in Switzerland from the same sites in 2017 and 2018. (c) F IS calculated between samples collected in the three sites in Switzerland in 2017 and in 2018. (d) Inbreeding coefficient F IS between specimens collected from the same sites in Switzerland in 2017 and 2018.

Pairwise D ps between individuals ranged between 0.88 and 0.93. We did not find any indication of isolation by distance among the samples collected in Italy and Switzerland in the 2.europe_cleaned dataset (Mantel R = −0.17, p‐value = .988; Figure 8) even if we only included the samples collected in mainland Italy, excluding samples from Sicily (Mantel R = −0.16, p‐value = .974; Figure 8). In contrast, the neighbor‐net tree based on D ps distances shows a separation between samples from Northern and Southern Italy (Figure 6b).

FIGURE 8

Isolation by distance (IBD) analysis using the 2.europe_cleaned dataset represented as scatterplots. (a) Correlation of genetic distances as proportion of shared alleles (D ps) and geographic distances on logarithmic scale for samples from Italy (excluding samples from Sicily). The correlation was assessed using a Mantel test, R = −0.16, p‐value = .974 based on 999 replicates. (b) Correlation between D ps genetic distances and logarithmic geographic distances for samples from Northern Italy and southern Switzerland. The correlation was assessed using a Mantel test based on 999 replicates, Mantel R = −0.17, p‐value = .988.

We performed genetic assignment tests on individuals from population collected in Switzerland, using the method implemented in assignPOP. Due to their geographical proximity, individuals collected in France and in Northern Italy were assigned as source populations. Assignment accuracies of individuals collected in Italy (pop_itnd) are relatively low, whereas those collected in France (pop_fr) are higher (Figure 9a). Simulations performed best when all loci and individuals were used. On average, 41% of the individuals collected in Switzerland were assigned to Northern Italy and 59% to France, but only 56% of individuals were assigned with a proportion of genetic constitution of >75%, which can be considered as effective assignment (Figure 9b).

FIGURE 9

Genetic assignment tests. (a) Assignment accuracy estimated by Monte Carlo cross‐validation based on the 2. europe_cleaned dataset. Assigned source populations were France (pop_fr) and Northern Italy (pop_itnd). Red horizontal lines indicate 0.33 null assignment rate, where the assignment accuracy is zero. (b) Membership probability of the individuals collected in Switzerland, organized from north to south. Individuals are sorted based on the probability of assignment to their original populations.

DISCUSSION

Our aim was to describe the invasion history of Ae. albopictus into Switzerland and across the Alps, and to estimate if the current populations are self‐sustaining. As a point of reference, we compared the genetic variability in mosquitoes from Switzerland to populations from Italy that have been established for over 25 years, to recently established populations in neighboring France, Germany, and Liechtenstein, and to populations from the mosquito's native range in Japan, Indonesia, and China. We found that populations from Switzerland had similar genetic variability to those from well‐established populations in Italy (Figure 7a, Table 2), and that there were no clear patterns of isolation by distance (Figure 8). We detected weak genetic structuring with a high level of genetic admixture, supporting a scenario of rapid expansion after introduction into Switzerland—both south and north of the Alps (Figures 2, 3, 4). These findings are in line with observations from the Swiss national monitoring program, suggesting human‐aided dispersal along main transportation routes (Müller et al., 2020). While the genetic pattern suggests frequent re‐introductions from Italian sources, the recovery of a pair of full siblings at a distance of 330 m in Strasbourg (France) in two consecutive years (Table 4) suggests the presence of an overwintering population north of the Alps. To our knowledge, this result is an indirect molecular evidence for establishment of a self‐sustaining population north of the Alps.

TABLE 2

Basic diversity statistics for the dataset 2.europe_cleaned

Geographical regions	N ind	Private alleles	Mean H O	Mean H E	AR	F IS	F is CI
							Upper	Lower
CH‐North	12	133	0.062	0.080	1.085	0.269	0.255	0.297
CH‐South	27	447	0.062	0.085	1.087	0.249	0.247	0.274
Italy‐North	9	122	0.063	0.082	1.088	0.241	0.257	0.298
Italy‐Center‐South	8	92	0.057	0.07	1.077	0.244	0.208	0.257
Italy‐Sicily	6	65	0.061	0.072	1.081	0.254	0.211	0.270
France	14	195	0.058	0.073	1.078	0.299	0.224	0.260
Germany	3	33	0.064	0.071	1.086	0.227	0.256	0.327
Albania	6	105	0.058	0.073	1.082	0.294	0.258	0.309
Greece	6	107	0.062	0.075	1.083	0.225	0.204	0.252

Abbreviations: N ind, number of individuals; Mean H , mean observed heterozygosity; Mean H , mean expected heterozygosity; AR, allelic richness; F IS, inbreeding coefficient with 95% confidence interval (F IS CI).

Across all of our specimens (within the 1.native_invasive_cleaned dataset), we detected the presence of four genetic clusters (Figures 2, 3, 4). High levels of shared ancestry were recorded between mosquitoes collected in France, Italy, Switzerland, Germany, and the USA, while the mosquitoes collected in Albania and Greece were genetically distinct from the rest of Europe (Figures 3 and 4). These results suggest that mainland Europe could have been invaded by mosquitoes originating via the USA to Italy as previously proposed (Battaglia et al., 2016; Sherpa et al., 2019; Zhong et al., 2013). While Albania was the first European country invaded by Ae. albopictus, our results suggest that samples collected in Albania are genetically closer to samples from Greece and the USA (Figure 3 K = 5 and Figure 4). Nevertheless, we may not completely rule out recent gene flow from Albania, while the genetic pattern as well as the geographical isolation of the country in the past rather supports the hypothesis that the invasion on mainland Europe goes back to an origin in the USA. Assigning the primary source with absolute rigor is very challenging considering the very recent colonization of this species in the study area. Our genetic assignment tests aiming to identify primary sources do not reveal the full picture (Figure 9) and, therefore, future studies should consider a denser sampling scheme across Italy, especially the northern regions. In addition to denser sampling, using whole genomes or a larger number of SNPs could help shedding more light on some of the recent invasion histories.

The approaches used to test genetic clustering in our European dataset did not yield entirely consistent results (Figures 5 and 6), suggesting that in Europe, there are at least three different clusters, with some genetic admixture between two of these clusters including specimens from Italy and Switzerland. This finding differs from previous studies (Pichler et al., 2019; Sherpa et al., 2019) that suggested two distinct genetic clusters in Italy, one comprising specimens from Northern Italy originating from the USA, and another one consisting of specimens from the central and southern areas that originated from admixture between the northern Italian genetic cluster and individuals from China. In our data, we also identified one mosquito from Sicily (Messina) that clustered together with mosquitoes collected in Brazil (Figure 2). Previous studies have suggested that the Brazilian populations are genetically distinct from the North American ones (Birungi & Munstermann, 2002) and are only partially related to the ones from the native range. An explanation for this apparent discrepancy might be undersampling of the native range (Kotsakiozi et al., 2017; Pichler et al., 2019), suggesting that the lineage that became invasive in Brazil has not yet been sampled in the native range. The Messina individual might have originated from the same yet unknown native lineage that also invaded Brazil, highlighting the need for future studies.

Overall, there is a genetic similarity of the mosquitoes collected in Switzerland to the one collected in Italy (Figure 3). This similarity, together with the close proximity of the two countries, the intense traffic of goods and people, and surveillance data and lack of isolation‐by‐distance, supports the hypothesis of the introduction of Ae. albopictus from Italy to Switzerland. This is in line with previously published studies indicating Italy as the main source for the European spread (Kraemer et al., 2019; Sherpa et al., 2019).

Interestingly, ddRAD‐seq SNP‐based studies with similar genomic resolution on a closely related species, Ae. aegypti, found strong spatial genetic structure at even small spatial scales (<4 km in Schmidt et al., 2018 and <200 m in Jasper et al., 2019), suggesting that the weak genetic structuring found in this study is not a result of low genomic resolution. This difference is more likely caused by different dispersal abilities and invasion histories of the two species. The global colonization of Ae. aegypti is older than in Ae. albopictus dating to 100 of years ago (Powell et al., 2018). In Ae. albopictus, reports of very high levels of differentiation among samples of recently invading populations at regional levels have been identified in Southern Russia, but heavily restricted gene flow or population exchange is reported between the different study sites (Konorov et al., 2021). The weak genetic structure, high levels of admixture, and lack of IBD found in this study for Ae. albopictus suggest rapid expansion most likely through human‐aided dispersal along transportation routes across the Alps. The human transportation network is known to have influenced and shaped the rapid spread of Ae. albopictus at regional levels (for a review see Medley et al., 2015 and Medlock et al., 2015). Switzerland is crossed by the European highways (E35 and E43). The E35 is a south–north European route that runs from Rome (Italy) to Amsterdam (the Netherlands), while the E43 connects Eastern Switzerland with Germany. Our results support the hypothesis that E35 has indeed acted as a key route of introduction of Ae. albopictus across the Alps, as previously suggested by surveillance data (Müller et al., 2020).

In Strasbourg, France, we collected a pair of full sibling in two consecutive years at 330 m of distance (Table 4). This finding is an indirect proof of overwintering of a mosquito population north of the Alps as well as the occurrence of skip oviposition behavior. Aedes albopictus is adapted to colder temperatures by producing dormant egg stages in fall that overwinter and hatch in the subsequent spring. Diapausing eggs have been described as the main mechanism enabling range expansion into regions at higher latitudes in North America and Northern Europe (for review see Armbruster, 2016; Batz et al., 2020). For overwintering populations, we would expect that mosquitoes caught in spring are closely related to mosquitoes caught in fall from the previous year. In Strasbourg, a full‐sibling pair was collected in two consecutive years supporting the hypothesis that the second pair hatched from eggs laid by the same mother and implying that they must have overwintered as diapausing egg. This result is supported by field observations from the local surveillance program as Ae. albopictus individuals were present in the same area in the preceding years, with the first detection dated in 2014 (Krupa et al., 2020), suggesting the presence of a self‐sustaining population. The proportion of the population actually overwintering and the proportion of individuals which are re‐introduced every year remain yet to be identified in this area. In Switzerland, we did not find such closely related siblings, but we observed a high genetic similarity in mosquitoes collected from the same sites in two consecutive years (Figure 7b–d). This, together with the decline of heterozygosity and the increase in the inbreeding coefficient between mosquitoes from the same sites, supports the presence of overwintering populations (Figure 7). A population that is continuously inbreeding locally is likely to have higher inbreeding coefficient as mating occurs between individuals related by descent and an overall decline in heterozygosity is expected (Rumball et al., 1994). The reporting by local surveillance activities of individuals early in the season in both 2018 and 2019 further supports the likely presence of overwintering eggs. In Germany, studies based on ecological (Kuhlisch et al., 2018; Pluskota et al., 2016) and genetic data (Lühken et al., 2020; Walther et al., 2017) also suggest the presence of overwintering populations across the country.

The recovery of one pair of full siblings between two consecutive years in Strasbourg also provides indirect evidence of skip oviposition (Table 4). Skip oviposition describes the behavior of a female mosquito depositing eggs in multiple breeding sites during a single gonotrophic cycle (Corbet & Chadee, 1993). Since the full siblings must have been from the same mother, we can conclude that the same female mosquito laid eggs of a single batch in two different breeding sites. This result confirms previous laboratory (Davis et al., 2015) and field (Davis et al., 2016) evidence, showing skip oviposition behavior in this species. This is especially relevant from a control perspective as this behavior could be potentially exploited to develop auto dissemination control measures (Caputo et al., 2012; Gaugler et al., 2012).

The high genetic variability of the mosquito populations (Figure 3, Table 2) across the Alps suggests multiple re‐introductions from different sources. The frequent re‐introductions of specimens from multiple sources are the likely cause of the high level of admixture found in our data (Figure 3), which is also contributing to maintain high genetic variation within local populations. This, in turn, might increase the probability of further spread. We found evidence of eggs going through diapause across the Alps, which suggests that the mosquito is potentially adapted to survive the colder winters. Taken together, the expansion patterns suggest that the Alps are not a barrier for Ae. albopictus and we may expect further spread in Central Europe. As a consequence, control measures should be designed to detect and target mosquitoes early in the season in order to prevent adults from hatching from diapausing eggs.

AUTHOR CONTRIBUTIONS

Laura Vavassori: Conceptualization (equal); data curation (lead); formal analysis (lead); investigation (lead); methodology (lead); project administration (lead); resources (lead); writing – original draft (lead); writing – review and editing (lead). Ann‐Christin Honnen: Data curation (equal); formal analysis (equal); methodology (equal); supervision (equal); validation (supporting); visualization (supporting); writing – review and editing (equal). Norah Saarman: Data curation (equal); formal analysis (equal); methodology (equal); resources (equal); software (lead); writing – review and editing (equal). Adalgisa Caccone: Formal analysis (equal); investigation (equal); methodology (equal); resources (equal); software (equal); supervision (equal); writing – review and editing (equal). Pie Müller: Conceptualization (lead); funding acquisition (lead); investigation (equal); methodology (equal); supervision (lead); visualization (equal); writing – review and editing (equal).

CONFLICT OF INTEREST

The authors declare no competing interests.

Appendix S1

Click here for additional data file.