The fast invasion of Europe by the box tree moth: an additional example coupling multiple introduction events, bridgehead effects and admixture events

Identifying the invasion routes of non-native species is crucial to understanding invasions and customizing management strategies. The box tree moth, Cydalima perspectalis, is native to Asia and was recently accidentally introduced into Europe as a result of the ornamental plant trade. Over the last 15 years, it has spread across the continent and has reached the Caucasus and Iran. It is threatening Buxus trees in both urban areas and forests. To investigate the species’ invasion routes, native and invasive box tree moth populations were sampled, and moth’s genetic diversity and structure were compared using microsatellite markers. Our approximate Bayesian computation analyses strongly suggest that invasion pathways were complex. Primary introductions originating from eastern China probably occurred independently twice in Germany and once in the Netherlands. There were also possibly bridgehead effects, where at least three invasive populations may have served as sources for other invasive populations within Europe, with indication of admixture between the two primary invasive populations. The bridgehead populations were likely those in the countries that play a major role in the ornamental plant trade in Europe, notably Germany, the Netherlands, and Italy. All these invasion processes likely facilitated its fast expansion across Europe and illustrate the role played by the ornamental plant trade not only in the moth’s introduction from China but also in the species’ spread across Europe, leading to an invasion with a complex pattern.


Introduction
It is well known that human activities result in the transportation and introduction of non-native species across the world (Blackburn et al. 2011;Essl et al. 2015;Gippet et al. 2019). Insects represent a large share of these new introductions and the rate at which such insect species are being introduced is still climbing (Seebens et al. 2017;Bonnamour et al. 2021). In Europe, most of the non-native insects that have entered and established themselves since the 1990s are spreading faster than ever across the continent . The liberalization of commercial policies and travel in Europe, and the subsequent dismantling of border checkpoints within the European Union (EU), have facilitated trade, especially that of ornamental plants. These changes may have also simultaneously contributed to the dispersal of nonnative species across Europe . Understanding the invasion pathways of non-native insect species is a crucial part of identifying the factors driving invasions (Allendorf and Lundquist 2003;Bock et al. 2015). The knowledge gathered will help identifying trends on invasion pathways underlying faster spread observed for insects introduced accidentally (Sakai et al. 2001;Essl et al. 2015;Roques et al. 2016;Fahrner and Aukema 2018).
Human-mediated invasions often involve complex pathways (Garnas et al. 2016;Meurisse et al. 2019). To better understand invasion pathways, it is important to characterize the biology and ecology of nonnative species, to consider when they first appeared in their invasive ranges, and to examine any customs records of their interception , including the commercial products in which they occurred (Essl et al. 2015;Brockerhoff and Liebhold 2017). However, these data can be misleading, sparse, incomplete, or absent, making it difficult to trace invasion histories. In this context, molecular resources can help fill in the gaps (Estoup & Guillemaud, 2010;Lawson et al. 2011). Microsatellite markers can reveal whether the invasion process was characterized by certain common phenomena, such as bottlenecks, multiple introductions, or admixture events (Estoup and Guillemaud 2010;Lawson Handley et al. 2011;Bock et al. 2015;Dlugosch et al. 2015). Microsatellite markers can also highlight bridgehead effects, which occur when a primary invasive population serves as the source for secondary invasive populations (Lombaert et al. 2010;Lawson Handley et al. 2011). That said, it is not always straightforward to infer invasion routes using molecular data because appropriate analytical methods are often lacking to obtain detailed information (i.e. population source, number of independent introductions), as shown in previous studies Cristescu 2015). Luckily, approximate Bayesian computation (ABC) is an effective method for dealing with complex scenarios. It is a powerful tool that combines the use of simulated datasets and summary statistics (Beaumont et al. 2002;Csilléry et al. 2010;Bertorelle et al. 2010;Beaumont 2010), and it is a particularly good analytical fit for research on non-native species (e.g., Lombaert et al. 2010;van Boheemen et al. 2017;Fraimout et al. 2017;Lippens et al. 2017;Lesieur et al. 2019). Indeed, in the case of recently established invasive populations, geographical and historical information can be used to define and assess hypothetical invasion scenarios in order to infer introduction route(s) (Estoup and Guillemaud, 2010). Furthermore, ABC can be used to analyze large datasets and to estimate demographic, historical, and genetic parameters . This method was used to investigate the possible invasion routes of certain non-native insects (e.g., Ryan et al. 2019;Mutitu et al. 2020), including species that were observed spreading quickly across Europe ) like the Harlequin ladybird, Harmonia axyridis (Lombaert et al. 2010(Lombaert et al. , 2014a, the spotted wing drosophila, Drosophila suzukii (Fraimout et al. 2017), and the western conifer seed bug, Leptoglossus occidentalis (Lesieur et al. 2019). In these three cases, researchers observed that invasion dynamics and dispersal processes were complex, involving multiple introductions, admixture events, and bridgehead effects. However, such studies remain scarce, especially for species associated with the ornamental plant trade, which is one of the current major pathways for the introduction of non-native insects (Kenis et al. 2007; Roques 2010;Liebhold et al. 2016;Eschen et al. 2017).
The box tree moth (BTM; Cydalima perspectalis (Walker), Lepidoptera: Crambidae) is native to China, Japan, and Korea (Maruyama and Shinkaji 1987;Xiao et al. 2011;Kim and Park 2013). It is a multivoltine species that produces two to four generations per year in its invaded range (Nacambo et al. 2014;Göttig and Herz 2017). In the invaded range, damage has been recorded only on Buxus species Matošević et al. 2017;Ferracini et al. 2022). In its native range, the insect has been reported to develop on several Buxus species, and three other host plants (Ilex purpurea, Euonymus alatus and E. japonicus) ). However, the use of these three plants as host by BTM has recently been debated. Several performance experiments done with invasive populations showed that BTM larvae were only able to develop on Buxus species, suggesting that the moth is an herbivorous specialist of this plant genus (Brua 2014;Matošević et al. 2017;Ferracini et al. 2022). In Asia, approximately 40 native species of Buxus occur, among which at least 15 species are found in China (Kohler and Brückner 1989;Fang et al. 2011). In contrast, only up to 4 species, including the widespread B. sempervirens, are naturally present in Europe, the Caucasus and Iran (Di Domenico et al. 2012;Mitchell et al. 2018). These plants usually consist of shrubs that are naturally found in forest stands. Additionally, their slow growth pattern has developed gardeners' interest leading to a popular use as ornamental trees in urban environments.
BTM was accidentally introduced into Europe from Asia, likely as a result of the ornamental plant trade (Kenis et al. 2013). The species was first detected in 2007 in the German cities of Weil am Rhein and Kehl (state of Baden-Württemberg) (Krüger 2008). That same year, invasive populations also appeared in nearby cities in Switzerland (Leuthardt et al. 2010), as well as in plant nurseries in the Netherlands (Van der Straten and Muus 2010). Then, between 2007 and 2010, BTM was observed in several countries of western Europe, namely France (Feldtrauer et al. 2009), Belgium (Casteels et al. 2011), Italy (Bella 2013), and the United Kingdom (Salisbury et al. 2012). In the years that followed, the moth appeared in central and southeastern Europe (e.g., Croatia [Koren & Crne, 2012]; Hungary [Sáfián & Horváth, 2011]; and Greece [Strachinis et al. 2015]) as well as in Turkey [Hizal et al. 2012];southern Russia [Gninenko et al. 2014]; Georgia [Matsiakh et al. 2018] and Iran (Ahangaran 2016). At present, BTM has invaded more than 30 Eurasian countries (Bras et al. 2019), where it is threatening several Buxus species (Mitchell et al. 2018) both in urban and forest stands (Ferracini et al. 2022). It has been observed for the first time in 2018 in Canada (Frank 2019) and was also intercepted several time in the USA on nursery plants shipped from Canada in 2020 and 2021 (USDA-APHIS 2021).
In Europe, Buxus trees are widely used as ornamental plants, especially B. sempervirens, in public and private gardens, which results in their strong commercial trade (Kenis et al. 2018;Matošević, 2013;EPPO, 2012). It has been hypothesized that the Buxus tree trade between China and Europe led to the insect's arrival in Europe (Leuthardt et al. 2010;Casteels et al. 2011;Nacambo et al. 2014). Results from a study analyzing DNA sequences from native and invasive BTM populations suggest that the moth may have been introduced from eastern China on several occasions (Bras et al. 2019). Further evidence for the occurrence of multiple introductions is enhanced by the geographical disparity in the records of the moth's first appearance in 2007. Indeed, recently introduced populations were found simultaneously in Germany, the Netherlands, and Switzerland (Krüger 2008;Leuthardt et al. 2010;Van der Straten and Muus 2010). Moreover, the moth's population genetic structure in Europe (Bras et al. 2019) suggests that the species may have spread across the continent as a result of the ornamental plant trade among European countries (Kenis et al. 2013;Matošević, 2013). Indeed, several lines of evidence indicate that humans may have facilitated BTM's dispersal: (1) the species was intercepted in the Netherlands (EPPO, 2012), a major player in the horticulture industry (Eschen et al. 2017); (2) it was first observed in plant nurseries (Van der Straten and Muus 2010;Casteels et al. 2011;Salisbury et al. 2012); and (3) these initial records were associated with extremely disparate locations. However, because detailed data on the trade of specific plant genera are unavailable for many countries (Eschen et al. 2017), it has been challenging to identify the import/export routes of Buxus trees in Europe. In general, it has been difficult to trace both the pathways making up the ornamental plant trade into Europe (Dehnen-Schmutz et al. 2010) and the invasion routes of the pests introduced as a result, especially those of species not included in the EPPO alert list, such as BTM (Strachinis et al. 2015).
In this study, we aimed to characterize BTM's invasion routes into and across Europe. We assessed the evidence for the existence of multiple introductions, namely those hypothesized by Bras et al. (2019), as well as for bridgehead effects and admixture events, given that all these processes could have contributed 3868 A. Bras et al. to the species' successful introduction and fast expansion. Because historical records of documented introductions are scarce, we employed genetic methods (as recommended by Estoup & Guillemaud, 2010). First, we sampled and genotyped individuals from BTM's native and invaded ranges. Second, we used ABC to infer the species' invasion routes across Europe.

Sampling
Individual box tree moths were collected at 37 locations in the species' native and invaded ranges between 2012 and 2017 (Table 1). We preferentially sampled larvae and pupae, which were collected from different Buxus trees to minimize the likelihood of sampling closely related individuals (i.e. from the same mating pair). When this was not possible, we used pheromone traps to collect adults. We tried to obtain at least 20 specimens per location, but samples sizes were smaller at 9 sites (Table 1). All specimens were stored in 96% alcohol at -21 °C to preserve their DNA.
In BTM's native range, 13 populations were successfully sampled. The sampling effort covered most of BTM's putative range (Supplementary Material Figure S1). But no moths could be found in Japan or in China's southern provinces of Guangdong, Guangxi, and Fujian despite the species having been recorded there (Kawazu et al. 2007;Bras et al. 2019). In other parts of China and in South Korea, larvae and pupae were mostly collected in urban areas where Buxus trees (mostly B. microphylla) have been planted for ornamental purposes. One exception occurred in Fuyang (N-FU, Zhejiang province, China; see Table 1 for population codes), where moths were obtained from Buxus trees planted in a sentinel plant nursery established on the border between farmland and a natural forest (Kenis et al. 2018).
In BTM's invaded range, samples were obtained from 13 countries in Europe and Eastern Black Sea region. A total of 24 populations were sampled (Table 1; Supplementary Material Figure S2). We sampled at locations where the moth was first observed in Europe in 2007, namely Weil Am Rhein (I-WAR) and Kehl (I-KE) in Germany and Giessen (I-GI) in the Netherlands. Only larvae and pupae were collected, except at Legnaro (I-LE), Italy, where adults were caught using pheromone traps. All the samples were obtained from man-made habitats, with the exception of those for the Bzyb Valley (I-BV) and Mtirala Park (I-MNP) in Georgia and Solokhaul (I-SO) in Russia, where insects were collected in natural forests of Buxus colchica.

Genotyping
All the specimens were dissected prior to DNA extraction. Using the DNeasy® Blood and Tissue Kit (Qiagen, Hilden, Germany), DNA was extracted from either the thoracic muscles, in the case of the larvae and pupae, or a leg, in the case of the adults. All vouchers are stored at − 21 °C at the URZF INRAE laboratory in Orléans, France. The DNA was diluted with ultrapure water when its concentration exceeded 50 ng/µL. Fifteen microsatellite markers developed by Bras et al. (2018) were used to genotype individuals from both the native and invaded ranges with the same multiplex PCR sets. PCR amplification and allele scoring were performed as described in Bras et al. (2018). In total, we genotyped 303 individuals from the native range and 593 individuals from the invaded range.
Genetic analyses: genetic diversity and population structure In a first step, we assessed the genetic diversity and the relationship between populations both from the native and invaded ranges. The software Arlequin v. 3.5 (Excoffier and Lischer 2010) was used to calculate mean allele number (An), observed heterozygosity (Ho), expected heterozygosity (He), and to test for significant deviation from Hardy-Weinberg proportions (HW), and linkage disequilibrium (LD). To deal with multiple pairwise comparisons between loci, a sequential Bonferroni correction (Rice 1989) was applied to the results for both HW and LD. Null alleles at each locus were estimated using FreeNA (Chapuis and Estoup 2007). Allelic richness (AR) for each population was calculated based on the minimum sample size of 14 individuals using FSTAT v. 2.9.3 (Goudet 1995). FreeNA software was also used to calculate pairwise F ST values using the excluding null allele correction method. We reconstructed a population-level neighbor-joining tree using  Population v. 1.2.32 (http:// bioin forma tics. org/ ~tryph on/ popul ations/) and Cavalli-Sforza and Edwards chord distances. A discriminant analysis of principal components (DACP) (Jombart et al. 2010) was employed to preliminarily explore the genetic structure of the dataset; the adegenet package (Jombart 2008) implemented in R software (R Development Core Team 2008) was used.
To further explore the dataset's genetic structure and to identify its potential number of genetic clusters, we used a Bayesian clustering approach implemented in STRU CTU RE (Pritchard et al. 2000). Three independent STRU CTU RE analyses were performed using (i) the data from all the specimens; (ii) the data from the native populations (eastern Asia); and (iii) the data from the invaded populations (Europe, Eastern Black Sea region). Parameter values were estimated using an admixture model where allele frequencies were correlated among populations. For each analysis, 15 replicate runs were performed for each value of cluster number (K). Every run comprised a burn-in period of 200,000 MCMC iterations that was followed by 1,000,000 iterations. For the full dataset and the dataset for the invasive populations, K was set to range between 1 and 15. For the dataset for the native populations, K was set to range between 1 and 10. The uppermost level of population structure was estimated for each analysis using the ΔK method (Evanno et al. 2005), which was implemented in STRU CTU RE Harvester (Earl and vonHoldt 2011). For each K, the mode among runs was checked using the Greedy algorithm in Clumpp (Jakobsson and Rosenberg 2007) via the online pipeline CLUMPAK (Kopelman et al. 2015). Genetic structure was visualized using Distruct v. 1.1 (Rosenberg 2004).

Inferring invasion scenarios with approximate Bayesian computation
We used ABC (Beaumont et al. 2002) to infer BTM's invasion routes into and across Europe. As ABC scenarios can be complex and thus difficult to assess (Estoup et al. 2012), we only used a selected number of native and invasive populations in our ABC analyses. More specifically, we focused on the populations that corresponded to the genetic clusters identified with STRU CTU RE, as per Lombaert et al. (2014a, b).
The STRU CTU RE analysis of the native populations revealed that each genetic cluster was represented by a single sampling region. Hence, there were three main genetic clusters in the moth's native range: one in South Korea (Seoul, N-SE), one in eastern China (Shanghai, N-SH), and one in southwestern China (Lijiang, N-LIJ). Because our sampling efforts were unsuccessful in several parts of BTM's native range (e.g., Japan, central China), we included a fourth, unsampled population (N-US) in the ABC analyses. No alternative datasets were available because we only obtained one population for South Korea and one population for southwestern China. We selected six populations for the STRU CTU RE analysis of the invaded range. In making our choices, we took into account geographical position and historical information (e.g., year of first introduction, country-specific data on the ornamental plant trade). First, we selected Weil Am Rhein (I-WAR), Kehl (I-KE), and Giessen (I-GI) (Van der Straten and Muus 2010) because they are where BTM was first observed in Europe. To obtain a better coverage of the invaded range, we also selected one population to represent each of the following sampling regions: western Europe, eastern Europe, and Eastern Black Sea region. Our choices were Lucca (I-LU), Budapest (I-BU), and Solokhaul (I-SO), respectively. Small sample sizes can negatively affect the rate at which the most probable scenario is identified . For Solokhaul (I-SO) for which we were only able to genotype 15 individuals, we thus adopted the methodology described in Lombaert et al. (2014a, b) and employed an alternative dataset, that associated with the next closest population, Bzyb Valley (I-BV).
The ABC analyses in themselves involved two steps. In the first step, we examined the oldest invasive populations recorded in Europe (Weil Am Rhein, I-WAR; Kehl, I-KE; and Giessen, I-GE) to better understand the pathways followed by BTM from eastern Asia. The goal was to (i) verify that eastern China was the source and (ii) determine whether the three first invasive populations in Europe, which were established in 2007 at distant locations from one another, resulted from a single or multiple introduction events (Table 2). Since these three populations appeared the same year, we conducted two separate ABC analyses. We initially assessed whether each invasive population had a Chinese origin. We then tested the relationship among these populations without accounting for admixture events since they appeared during the same year. In the second step, we used ABC analyses to clarify BTM's dispersal within Europe. We also wished to determine whether populations in the western, central, and eastern parts of the invasive range (Lucca, I-LU; Budapest, I-BU; and Solokhaul, I-SO, respectively) had resulted from direct introductions from China, another primary invasive population in Europe, or admixture events. The scenario topologies were based on the first dates of record for the moth. Based on the first set of ABC analyses, we defined scenarios where the targeted population was resulting of either an introduction from one population or an admixture between two defined populations. In total, we conducted eight successive ABC analyses comprising 4 to 36 scenarios. After each ABC analysis, the most likely scenario was used as the basis for subsequent analyses.
In our ABC analyses, the values for the historical parameters were drawn from prior distributions defined using historical data (Table 3). In summary, the years of introduction for BTM in Europe were drawn from uniform distributions bounded by 2000 and 2007, defined by the first three populations observed in Europe. Given that BTM produces two to four generations per year in the invaded range, we considered that one year was equal to three generations. The same range (i.e., 7 years) was used for the invasive populations observed after 2007. Following Lombaert et al. (2011), we simulated sub-structure by using of ghost populations in the native area to consider that the true source population might not have been sampled. Population divergence was immediately followed by a bottleneck period (DB i generations) during which there was a reduction in effective population size (effective population size at founding, Ninf i ). After the bottleneck period, effective population size climbed back up and stabilized (stable effective population size, N i ). The results were expressed using all the summary statistics available in DIYABC v. 2.1 (Cornuet et al. 2014). These statistics describe population-specific genetic variation (e.g., mean number of alleles), genetic variation between populations (e.g., F ST ), and admixture between populations. The total number of summary statistics ranged from 76 to 576 (Table 2).
The ABC Random Forest (ABC-RF) method described by Pudlo et al. (2016) was used to compare the tested scenarios for each ABC analysis to define the most likely. This method has the capacity to discriminate among scenarios more efficiently than traditional ABC methods, especially when complex invasion pathways are involved (see Fraimout et al. 2017). The ABC-RF method is a nonparametric, machinelearning tool which carries out the classification of hundreds of bootstrapped decision trees using summary statistics as predictor variables. It extracts the maximum of information from the entire set of summary statistics, avoiding choosing arbitrary a subset as required in other classical ABC methods. A "prior error rate" (PER) can also be estimated using the simulations that were not employed to build the trees. PER provides a direct means of cross validation. With DIYABC, we generated 10,000 microsatellite datasets per scenario, which we then used to grow a classification forest of 1000 trees. The scenario with the greatest number of votes was defined as being the most likely. When the number of votes was similar among several scenarios (i.e., scenarios differed by fewer than 50 votes), the two scenarios with the greatest numbers of votes were repeated using the ABC-RF method to verify the accuracy of scenario selection. We used the abcrf R package (Pudlo et al. 2016) to perform all the ABC-RF analyses. Finally, we inferred posterior probabilities (PPs) for the final scenario in the final analysis by performing Random Forest regression (Raynal et al. 2019) on classification forests of 1000 trees. Simultaneously, model checking was run with DIYABC to evaluate the fit between the most likely scenario and the observed dataset.

Genetic diversity of the box tree moth
Overall, we genotyped 896 individuals from 37 sampling sites. Genotyping success was high (98.85%), and there was a mean of 8.07 alleles per locus. In this full dataset, six cases of linkage disequilibrium were found after Bonferroni correction. A given pair of loci never displayed significant linkage disequilibrium in more than one population. The fifteen microsatellite markers were thus considered to be independent.
After Bonferroni correction, three populations in the native range (N-BE, N-XI, N-HU; Table 1) Table 1) had at least one locus that deviated significantly from Hardy-Weinberg proportions. When looking at Hardy-Weinberg proportions for each locus, only the locus BTM22 departed significantly from Hardy-Weinberg equilibrium in 11 populations (1 native population and 10 invasive populations). It also had a high frequency of null alleles (21%). For the other loci, mean frequencies of null alleles were less than or equal to 10%. We thus performed the analyses with and without BTM22 to determine whether inclusion of the locus introduced bias. Since the results were similar in both cases, we have reported the findings obtained with the full, fifteen-microsatellite dataset.
Mean allelic richness was higher in the native populations (4.23) than in the invasive populations (2.96) ( Table 1). For the native populations, the ranges of observed heterozygosity (Ho) and expected heterozygosity (He) were 0.43-0.56 and 0.50-0.57, respectively. Ho and He were lower in the invasive populations, and their ranges were 0.33-0.48 and 0.28-0.49, respectively. For the full dataset, the range for pairwise F ST was 0.01-0.38; it was smaller for native (0.01-0.14) than for invasive populations (0.01-0.38) Table 3 Historical, demographic, and genetic parameters used to assess box tree moth (Cydalima perspectalis) invasion scenarios via ABC analysis

Parameter prior distributions are indicated
The parameters Ninf i and DB i are associated with the invasive populations. Times (t) are expressed in numbers of past generations with a year being equal to three generations. The microsatellite loci were assumed to follow a generalized stepwise mutation model and were characterized using three parameters: the mutation rate (µ), the mean of the geometric distribution of mutation length (P), and the mean mutation rate for a single nucleotide insertion (SNI). Note that for [x; y], x and y are the bounds of the prior distributions. Assumptions: Nsup > N N-US ; t N-SHA > t GH-W ; t N-SHA > t GH-K ; t N-SHA > t GH-G ; t GH-W ≥ t I-WAR ; t GH-K ≥ t I-KE ; t GH-G ≥ t I-GI ; t I-GI > t I-LU ; t I-WAR > t I-LU ; t I-WAR > t I-BU ; t I-GI > t I-BU ; t I-GI > t I-SO ; t I-LU > t I-SO ; and t I-LU > t I-BV (Suppl. Mat. Table S1). Among the native populations, the south Korean population N-SE had the highest F ST values (0.06-0.14) followed by the southern Chinese populations N-GU (0.06-0.12) and N-LI (0.05-0.13). Compared to the invaded range, F ST values in Asia were increasing when geographic distances between populations were larger (Suppl. Mat. Figure S3).

Population structure
In the unrooted neighbor-joining tree, the invasive and native populations formed distinct groups (Suppl. Mat. Figure S4). The structure of this tree was consistent with the geographical distribution of the native populations. For the invasive populations, populations at the sites where the moth was first observed in 2007 (Weil Am Rhein, I-WAR; Kehl, I-KE; Giessen, I-GI) occurred in separate groups, and most of the invasive populations were part of the Giessen (I-GI) group. However, the low bootstrap values at most nodes made it difficult to draw any robust conclusions about the relationships between groups. The DAPC analysis yielded the same pattern, except that Katerini (I-KA) did not group with the other invasive populations (Suppl. Mat. Figure S5). The STRU CTU RE analysis results were consistent with those from the neighbor-joining tree, the DAPC analysis, and the F-statistics analyses ( Fig. 1). At K = 2 (selected via the ΔK method), all the native populations as well as the invasive populations of Weil Am Rhein (I-WAR), Kehl (I-KE), and Delémont (I-DE) belonged to the same genetic cluster. The second genetic cluster included all the other invasive populations, with the exception of Budapest (I-BU) and Paris (I-PA). These latter two populations showed signs of having arisen from admixture events between the two clusters. At K = 3, the native populations formed a distinct group from the invasive populations.
The analysis focusing on the native populations found evidence of genetic structure (Fig. 1, Suppl. Mat. Figure S6). The uppermost level of structure was K = 2. South Korea and southwestern China were clearly differentiated from the other native populations. Populations in northern and eastern China were grouped together, and individuals from Xinyang (N-XI) and Lishui (N-LIS) were assigned to the first genetic cluster. At K = 3, the grouping pattern was consistent with population geography. Populations from South Korea and southwestern China formed two distinct genetic clusters, while populations from northern and eastern China grouped together and were mixed in an unclear pattern.
In the analysis focusing on the invasive populations, K = 2 was the optimal level of structure (Fig. 1, Suppl. Mat. Figure S6)

Assessing scenarios for the box tree moth's invasion of Europe
The first set of ABC analyses examined the populations at the locations where BTM was first observed in Europe: Weil Am Rhein (I-WAR) and Kehl (I-KE) in Germany and Giessen (I-GI) in the Netherlands. First, we found that, in each of the ABC analyses, the most likely scenario was that eastern China served as the source for these populations (PER: ~ 7.6% and PP ranging from 0.72 to 0.83; Table 2, Fig. 2a). Second, when we inferred the relationships between these three invasive populations, the most likely scenario was that there were three independent introductions of the moth from eastern China into Germany The fast invasion of Europe by the box tree moth: an additional example coupling multiple… Georgia (I-BV) (Fig. 2B). The best scenario for Lucca (I-LU) was that population arose from admixture between the Weil Am Rhein (I-WAR) and Giessen (I-GI) populations (PER: 29.46% and PP: 0.47; Table 2). That said, there was a difference of less than 50 votes between this scenario and a scenario in which only the Giessen (I-GI) population served as its source. We therefore carried out an ABC-RF analysis using just these two scenarios, and the results supported the admixture event between Fig. 1 Graphical representation of box tree moth (Cydalima perspectalis) population genetic structure estimated using the Bayesian clustering approach implemented in STRU CTU RE. Results for the full dataset, native range dataset, and invaded range dataset respectively are represented. The population codes appear below the plots (see Table 1). Each individual is represented by a vertical line, and each color represents a particular genetic cluster. In bold is the best K for each STRU CTU RE analysis (full dataset, native range alone and invaded range alone, respectively), as defined by STRU CTU RE Harvester Weil Am Rhein (I-WAR) and Giessen (I-GI) populations as the most likely scenario (PER: 9.13% and PP: 0.58). Consequently, this scenario was used in subsequent analyses. For Budapest (I-BU), the most likely scenario was that the population had arisen also from admixture between Weil Am Rhein (I-WAR) population and the Giessen (I-GI) population (PER: 33.56% and PP: 0.52). However, there was a difference of less than 50 votes between this latter scenario and two other scenarios: admixture occurred between the Shanghai (N-SHA) and Giessen (I-GI) populations or admixture occurred between the Shanghai (N-SHA) and Lucca (I-LU) populations. We used the same approach as with Lucca (I-LU), and the most likely scenario remained the occurrence of an admixture event between Weil Am Rhein (I-WAR) and Giessen (I-GI) (PER: 16.37% and PP: 0.62). Consequently, we used this scenario in subsequent analyses. Finally, we found an admixture event between Giessen (I-GI) and Lucca (I-LU) populations as the best scenario for Solokhaul (I-SO) population (PER: 34.98% and PP: 0.48%) whereas considering the alternative Bzyb Valley (I-BV) population, the scenario in which the Lucca (I-LU) population served as source was selected (PER: 27.21% and PP: 0.53; Fig. 2b). None of the ABC analyses resulted in highly probable scenarios where South Korea, southern China, or an unsampled population served as a source.
The posterior model checking performed on the final scenario of Solokhaul (I-SO) population and the alternative population, Bzyb Valley (I-BV), indicated that the selected model and posteriors fitted the observed genetic data ( Figure S7). As the PER and PP values were better for the alternative Bzyb Valley (I-BV) population, we chose it to infer the posterior The fast invasion of Europe by the box tree moth: an additional example coupling multiple… estimations of model parameters on the final scenario (Table S2).

Discussion
The observed genetic diversity and ABC results strongly support that the box tree moth (BTM) dispersed along complex invasion pathways as previously suggested by mtDNA diversity (Bras et al. 2019). We confirmed that C. perspectalis was introduced from Eastern China to Europe with at least three independent introductions in Germany and the Netherlands, two leader countries in ornamental plant trade. ABC simulations emphasize the likely role of both bridgehead, with at least three invasive populations serving as a source of colonists for remote new territories, and admixture events in the moth's spread across the European invaded range. Our findings emphasize the role played by the ornamental plant trade, which has been an important factor promoting the establishment and spread of non-native phytophagous insects (Eschen et al. 2017;Kenis et al. 2007;Roques, 2010). It certainly appears to have contributed to BTM's introduction and dispersal across Europe, The Caucasus and Iran.

New insights into the box tree moth's introduction into Europe from China
Asia is one of the major sources of the non-native insects that have been introduced into Europe ( Roques 2010). To help prevent new introductions, it is important to pinpoint the source population of such species within their native ranges (Essl et al. 2015). Using mtDNA markers, Bras et al. (2019) arrived at the conclusion that eastern China could be the source of invasive BTM populations in Europe. Here, the Bayesian analysis of microsatellite markers yielded strong support for this hypothesis. Furthermore, the ABC analyses clearly excluded the possibility that South Korea, southwestern China, or an unsampled population served as a source for the invasive European populations that we sampled. These results are consistent with the large volumes of Buxus trees that have been imported from China to Europe between 2006(EPPO, 2012, and they fit with the fivefold increase of the European imports of plants from China over the period -2015(Eschen et al. 2015. Furthermore, the coastal economies in China has developed rapidly over the years with Eastern China's international trade accounted for example for 86.5% of the national total in 2013 (Wan and Yang 2016). It is also thought that eastern China has been the source of several other non-native insects introduced to Europe during the last decade (Roques et al. 2020), such as the yellow-legged hornet, Vespa velutina (Arca et al. 2015) and the spotted wing drosophila, Drosophila suzukii (Fraimout et al. 2017).
Given the moth's simultaneous appearance in 2007 in three different countries-Germany, the Netherlands, and Switzerland (Krüger 2008;Leuthardt et al. 2010;Van der Straten and Muus 2010)-and its patterns of mtDNA diversity (Bras et al. 2019), it seems likely that the insect was introduced multiple times from China. Our ABC analyses pointed out that BTM was likely independently introduced from eastern China into Germany on two occasions, corresponding to the two first records of its presence in Europe (Krüger 2008), and on one occasion in the Netherlands. Netherlands and Germany's importation from international countries represents the greatest volume of ornamental plants among EU countries (Eschen et al. 2017) which is supporting such invasive routes. More particularly, since the 2000s, imports from China to the Netherlands represent an important part of such commerce (van Valkenburg et al. 2014), with more than a million of Buxus trees being imported between 2006(EPPO, 2012. However, the prior error rate and posterior probabilities were moderate in our ABC analysis considering the three invasive population together. Statistical support for a scenario does not necessarily mean that it is true (Csilléry et al. 2010), which means other scenarios could also explain the data. In our case, the three initial invasive populations likely originated from the same source area in China and might have been simultaneously introduced within a short period of time. The estimated values of the time since introduction (t i ) for these populations were alike and their first observation records were no more than a year apart, thus supporting this hypothesis. As indicated by Benazzo et al. (2015), repeated introductions from the same source may lead to difficulties to discriminate correctly the number of independent introduction events. Thus, a scenario in which one of the initial invasive populations serving as a source for another invasive population-from Weil Am Rhein to Giessen or the other way around-cannot totally be excluded.
Non-native species can more easily become established when they are introduced on multiple occasions directly from their native range (Blackburn et al. 2015;Cristescu 2015), which could have favored BTM's invasion of Europe. Indeed, propagule pressure plays a crucial role in increasing the probability of establishment of non-native species (Allendorf and Lundquist 2003;Lockwood et al. 2005;Kaňuch et al. 2021). Our ABC analyses suggest that the bottleneck severity was moderate among the initial invasive populations, especially in Weil am Rhein (I-WAR), and Giessen (I-GI). A large number of individuals could have been introduced into Europe from China, which could explain the elevated haplotype diversity observed by Bras et al. (2019) across the moth's invaded range compared to the diversity in its native range. However, we cannot rule out the possibility that there have been recurrent introductions of the moth from China at the first European observation locations in 2007, given that the importation of Buxus trees has been continuous (EPPO, 2012). Unfortunately, ABC analyses cannot distinguish among such hypotheses because the genetic pattern associated with continual invasions may be genetically similar to that associated with a single introduction event with a large propagule size .
Admixture events and bridgehead effects were likely involved in the fast spread of the box tree moth across Europe Most of the non-native insects introduced into Europe over the past two decades-such as the spotted wing drosophila (Drosophila suzukii), the western conifer seed bug (Leptoglossus occidentalis), and the box tree moth (C. perspectalis)-appear to be spreading at unprecedented rates . The dispersal of both D. suzukii and L. occidentalis are the result of complex invasion dynamics, involving multiple introduction events, bridgehead effects and admixture events (Fraimout et al. 2017;Lesieur et al. 2019). Such genetic patterns may result from humanmediated dispersal, which facilitates the rapid spread of non-native species within the invaded range (Garnas et al. 2016;Bertelsmeier and Keller 2018). Thus, the ever-increasing trade of ornamental plants among European countries (Cadic and Widehem, 2006;Dehnen-Schmutz et al. 2010) may have played a crucial role in BTM's fast range expansion as its flight capacities cannot explain the species' spread across Europe, the Caucasus and up to Iran in less than ten years (Van der Straten and Muus 2010; Bras 2018).
Our ABC analyses suggest that bridgehead effects coupled with admixture events have been involved in BTM's dispersal in Europe. Indeed, they brought to light that the invasive populations in Weil Am Rhein (I-WAR, Germany), and in Giessen (I-GI, the Netherlands), have probably both served as source for the populations in Lucca (I-LU, Italy), and in Budapest (I-BU, Hungary), resulting in admixed populations. Our analyses also confirmed that the BTM population of southern Russia had an Italian origin, as suspected by Gninenko et al. (2014). Indeed, the moth was first observed on Buxus trees imported from Italy for the Olympic Village of the 2012 Games in Sochi. All the invasive populations identified as sources in our ABC analyses were present either in areas known as main centers for nurseries and trade in arboricultural products (i.e. Giessen, Lucca; Van der Straten and Muus 2010; Bella 2013) or port (i.e. Weil Am Rhein; Casteels et al. 2011) highlighting the role of human activities through the ornamental plant trade in BTM's dispersal in Europe. The three separate introductions of BTM into Europe combined with the expanding trade in ornamental plants among European countries (Cadic and Widehem 2006;Dehnen-Schmutz et al. 2010) could have favored admixture between invasive populations. Intraspecific genetic mixing is thought to enhance invasion success (Rius and Darling 2014;Bock et al. 2015), but the role that it played in Europe's invasion by BTM remains unclear. In fact, the STRU CTU RE analysis of the East Asian populations also indicated signs of admixture among eastern China populations. Because our samplings in the native range essentially proceeded from cities where planted Buxus trees could have been imported from various nurseries in China, we cannot exclude the coexistence in such man-made habitats of moth populations from different geographic origins, and their possible in situ hybridization. Moreover, the change of habitat from the forest to the city, and the subsequent selective pressures within these human-modified habitats may then have led to preselect populations adapted to human habitats. Such a two-step process with populations moving from natural habitats to human-disturbed habitats has already been shown to favor invasion success in some other invasive species such as the little fire ant, Wasmannia auropunctata (Foucaud et al. 2013). In the case of BTM, we can then hypothesis that both admixture and pre-adaptation to human-modified habitats in native populations occurring in man-made habitats might have also favored BTM's invasion success in similar habitats in Europe as the insect was firstly observed in cities (Krüger 2008;Leuthardt et al. 2010;Van der Straten and Muus 2010).
BTM's rapid dispersal from western Europe up to the Caucasus region seems likely to have been triggered by a combination of the three genetic processes we have highlighted through our ABC analyses: e.g., (i) multiple introductions, bridgehead effects, admixture events, (ii) the species' biological traits, such as its strong flight capabilities (Bras 2018), and (iii) the widespread availability of the host plant (Nacambo et al. 2014;Thammina et al. 2016;Mitchell et al. 2018). They may also favor the moth's establishment and spread in North America.
Tracing the dispersal routes of the invasive box tree moth Our ABC results revealed the complex invasion pathways followed by BTM. They also underscored the difficulties inherent in clearly defining such routes when invasive scenarios are based on sparse information, the invasion process is complex, and range expansion is rapid. Our analyses strongly suggest that BTM was independently introduced at least three times from the same source population in Asia. However, BTM invasive routes among Europe were harder to identify. Despite the robustness of the ABC-RF method, compared to more classical ABC analysis (Pudlo et al. 2016;Fraimout et al. 2017), determining whether several introductions resulted from the same source can be challenging, especially when the species in question is spreading rapidly . The short bottleneck period (identified by the posterior estimations) together with a late sampling of the invasive populations may have led to the formation of several few genetically differentiated outbreaks. This variability coupled with an identical population source in China could partly explain the moderate posterior probability and prior error rate values obtained for the most probable scenarios that we identified in Europe (Bermond et al. 2012;Benazzo et al. 2015;Hargrove et al. 2017). Moreover, DIYABC assumes that no recurrent migration has occurred between populations (Cornuet et al. 2014). However, some of the invasive populations that we examined were likely geographically connected, notably the Weil Am Rhein (I-WAR) and Kehl (I-KE) populations in Germany. Indeed, both were sampled several years after the moth's initial appearance in 2007. These factors, combined with BTM's strong flight capacities (Bras 2018), suggest that migration among populations could have contributed to population mixing which may have reduce the likelihood of accurate source assignment. Moreover, signs of admixture among eastern China populations could also have led to a first genetic blur of the boundaries between populations (Hamelin and Roe 2020).
Our study provides a new example of a complex invasion process involving multiple introduction events, bridgehead effects, and admixture events. The invasive scenarios that we identified for BTM are consistent with what is known about ornamental plant trade dynamics in Europe (Dehnen-Schmutz et al. 2010). Complex invasion processes mediated by human activities seem to be the rule rather than the exception for the rapidly dispersing of non-native insects introduced over recent decades into Europe (Fraimout et al. 2017;Lesieur et al. 2019;Lombaert et al. 2014a, b;Sherpa et al. 2019). The recent establishment of BTM in North America (Frank 2019; USDA-APHIS 2021) demonstrates the need of assessing invasion pathways of such species in order to develop appropriate management strategies to avoid further introductions.