Genomic data support multiple introductions and explosive demographic expansions in a highly invasive aquatic insect

The study of the genetic makeup and demographic fate of alien species is essential to understand their capacity to recover from founder effects, adapt to new environmental conditions and, ultimately, become invasive and potentially damaging. Here, we employ genomic data to gain insights into key demographic processes that might help to explain the extraordinarily successful invasion of the Western Mediterranean region by the North American boatman Trichocorixa verticalis (Hemiptera: Corixidae). Our analyses revealed the genetic distinctiveness of populations from the main areas comprising the invasive range and coalescent‐based simulations supported that they originated from independent introductions events probably involving different source populations. Testing of alternative demographic models indicated that all populations experienced a strong bottleneck followed by a recent and instantaneous demographic expansion that restored a large portion (>30%) of their ancestral effective population sizes shortly after introductions took place (<60 years ago). Considerable genetic admixture of some populations suggest that hypothetical barriers to dispersal (i.e., land and sea water) are permeable to gene flow and/or that they originated from introductions involving multiple lineages. This study demonstrates the repeated arrival of propagules with different origins and short time lags between arrival and establishment, emphasizing the extraordinary capacity of the species to recover from founder effects and genetically admix in invaded areas. This can explain the demonstrated capacity of this aquatic insect to spread and outcompete native species once it colonizes new suitable regions. Future genomic analyses of native range populations could help to infer the genetic makeup of introduced populations and track invasion routes.

observe and study in real time ecological and evolutionary processes that, at best, can only indirectly be inferred from research conducted on natural populations (Lee, 2002;Moran & Alexander, 2014;Sakai et al., 2001;Sax et al., 2007). Understanding these processes can, in turn, be crucial to guide management practices aimed at preventing or mitigating the ecological and economic damage caused by exotic species (Estoup & Guillemaud, 2010;Facon et al., 2006;Yokomizo et al., 2017). Information on the ecoevolutionary dynamics, genetic makeup, and demographic fate of alien species is essential to determine their capacity to adapt to new environmental conditions, recover from founder effects, and disperse, key aspects to understand why some of them become invasive (Handley et al., 2011;Sax et al., 2007).
High-throughput sequencing now offers the possibility of inferring neutral and adaptive processes at an unprecedented resolution and testing refined hypotheses about the ecological and evolutionary pathways of biological invasions (Handley et al., 2011;Viard et al., 2016; for example, Ryan et al., 2019;Puckett & Munshi-South, 2019;McCulloch et al., 2021). Information on the genetic composition of invading populations can help to identify the number of source populations, determine invasion routes, and predict the intrinsic capacity of the species to spread locally (i.e., without humanassisted dispersal) once they arrive to non-native regions (Cristescu, 2015;Estoup & Guillemaud, 2010). For instance, the implications for management can be very different if the species spread from a single source population (i.e., one introduction event, high intrinsic capacity of expansion; for example, Cao et al., 2016;Trumbo et al., 2016;Zhang et al., 2019) or if the distribution in the non-native range is a consequence of multiple independent introduction events (i.e., this may suggest a more limited intrinsic capacity of expansion; for example, Hagenblad et al., 2015;Lejeusne et al., 2014;Miller et al., 2005;Xia et al., 2020). Determining the genetic background and origin of introduced populations is not only important to understand invasion routes and modes, but can also help to forecast admixture among previously isolated gene pools in the non-native range, and anticipate its potential negative impacts (Facon et al., 2008Kolbe et al., 2004). Interbreeding among introduced populations with different origins can quickly restore genetic diversity after strong founder effects, and increase the evolutionary potential of alien species through new allele combinations, which can ultimately boost their capacity to become invasive and outcompete native species with similar ecological requirements (Facon et al., 2008; but see Tsutsui et al., 2000). Genomic-based demographic reconstructions also have a high potential to understand the colonization history of invasive species, determine the mode of spread (e.g., gradual vs. explosive), and estimate the duration of lag phases between introduction and invasion if benchmark dates on the timing of introductions are available (Coutts et al., 2018;Kowarik, 1995;Puckett & Munshi-South, 2019;Rouget et al., 2016;Ryan et al., 2019). Direct estimates about the scale of gene flow are also essential to predict the capacity of the species to colonize new habitat patches and expand within the non-native range (Ciosi et al., 2011;Handley et al., 2011;Trumbo et al., 2016). When combined with spatially explicit information within a landscape genetic framework, such estimates can also help to identify barriers/corridors to dispersal and predict routes for future expansion Zalewski et al., 2009).
The boatman Trichocorixa verticalis verticalis (Fieber, 1851) (Hemiptera: Corixidae) is one of the few alien aquatic insects in the world and a paradigmatic example of a successful invasive insect (Carbonell et al., 2016;Fenoglio et al., 2016;Ricciardi, 2015). This corixid has an extended natural range in brackish water systems from Atlantic North America and the Caribbean region, but its distribution has recently expanded (<50 years ago) to similar habitats in southern Europe (Spain and Portugal), Africa (Morocco and South Africa), and Oceania (New Caledonia) (Carbonell et al., 2017;Fouzi et al., 2020;Rabitsch, 2008). This alien outcompetes native Corixidae (mainly Sigara spp.) in coastal wetlands (Carbonell et al. 2017(Carbonell et al. , 2020. The original dispersal pathway is unknown, but previous studies have suggested that the species probably arrived to non-native regions with maritime trade via ship ballast water or as accidental introductions accompanying the release of alien fish for mosquito control (Guareschi et al., 2013;Rabitsch, 2008). The fact that T. verticalis is one of few insects able to breed in hypersaline conditions and even survive in the open sea (Gunter & Christmas, 1959;Hutchinson, 1931) has also led to speculate that the species might even have reached the Western Palearctic via marine dispersion during severe storm events, although this possibility was considered to be highly unlikely (Sala & Boix, 2005 and references therein).
Different intrinsic factors have been suggested to make T. verticalis a successful invader, including high physiological plasticity (euryhalinity, broad thermal tolerance, etc.), short life cycle, extraordinary fecundity, and capacity to be passively transported as eggs (Carbonell et al., 2016(Carbonell et al., , 2020(Carbonell et al., , 2021Céspedes, Coccia, et al., 2019;Guareschi et al., 2013). However, no study has yet been performed on the demography and population genetics of the species in either the native or invasive ranges. As a result, nothing is known about the genetic makeup and demographic fate of introduced populations, patterns of genetic structure and dispersal at different spatial scales, or the potential role of different landscape features on promoting/limiting gene flow among populations.
In this study, we employ genomic data (ddRAD-seq) to shed light on key demographic processes that might underlie the extraordinarily successful invasion history of T. verticalis in the Western Mediterranean (Carbonell et al., 2017;Fouzi et al., 2020;Guareschi et al., 2013;Rodríguez-Pérez et al., 2009). Specifically, we first employed different approaches to quantify genetic structure (or its lack thereof) in populations within the non-native range in the Iberian Peninsula and Morocco. These analyses revealed significant genetic differentiation and the presence of three main genetic groups in the region, rejecting the hypothesis of genetic panmixia. Second, we performed coalescent-based simulations to infer the timing of divergence among the recovered genetic groups. This allowed us to test whether their divergence is compatible with in situ geographical diversification (e.g., due to strong genetic drift and/or genetic monopolization during serial colonization; De Meester et al., 2002;Trumbo et al., 2016) or if, on the contrary, the obtained inferences point to independent introduction events involving multiple source populations with different evolutionary histories (Cristescu, 2015;Handley et al., 2011). Third, we tested alternative demographic models to estimate the timing, magnitude, and rates of historical changes in effective population sizes, which can provide valuable insights about the capacity of the species to recover from founder effects (i.e., bottlenecks), the extent of time lags between introduction and invasion (i.e., invasion debt), and on whether demographic expansions during the invasion phase were gradual or explosive (Coutts et al., 2018;Essl et al., 2011;Roman & Darling, 2007;Rouget et al., 2016). Finally, we employed a spatially-explicit landscape genetic approach to infer potential corridors and barriers to dispersal (sea water, inland surface water, and land) that might explain genetic connectivity among populations and the spread of the species in the region (Handley et al., 2011;e.g., Zalewski et al., 2009;.

| Population sampling
In 2018, we used a D-framed pond net (0.5 mm mesh) to sample nine populations of Trichocorixa verticalis in the Iberian Peninsula (Portugal, n = 2 populations; Spain, n = 5 populations) and Morocco (n = 2 populations) (Table 1; Figure 1). We used occurrence records available in the literature to design sampling and collect specimens from populations representative of the known alien distribution range of the species in the Iberian Peninsula and northwestern Morocco (Carbonell et al., 2017;Guareschi et al., 2013). These include some previously unrecorded populations at wetlands where we suspected the species would be present, based on their coastal location and high salinities. We registered spatial coordinates using a global positioning system (GPS) and preserved whole specimens at -20℃ in 1500 μl ethanol 96% until needed for genomic analyses. Further details on sampling locations are provided in Table 1.

| Genomic library preparation and processing
We used NucleoSpin Tissue kits (Macherey-Nagel, Düren, Germany) to extract and purify DNA. We processed genomic DNA into one genomic library (72 individuals/library; 8 individuals/population;   (Catchen et al., 2013) to filter and assemble our sequences into de novo loci, call genotypes, calculate population genetic diversity (π, nucleotide diversity) and pairwise levels of genetic differentiation (F ST ), and export input files for all downstream analyses. Briefly, we used the "blacklist" option of stacks to remove outlier loci (i.e., those putatively under selection) identified using bayescan v. 2.1 (Foll & Gaggiotti, 2008) and the hierarchical and nonhierarchical island models (FDIST method; Excoffier et al., 2009) implemented in arlequin v. 3.5.2.2 (Excoffier & Lischer, 2010). We exported only the first SNP per RAD locus, and retained loci with a minimum stack depth ≥ 5 (m = 5), a minimum minor allele frequency (MAF) ≥ 0.01 (min_maf = 0.01) and that were represented in all populations (p = 9) and at least 50% of the individuals within each population (r = .5). For further details on data filtering and assembling, see Methods S1.

| Quantifying genetic structure
We analysed population genetic structure and admixture using the Bayesian Markov chain Monte Carlo clustering method implemented TA B L E 1 Locality, country, code, latitude, longitude, number of genotyped individuals (n), and genetic diversity (π) for each sampled population  (Pritchard et al., 2000). We ran structure analyses assuming correlated allele frequencies and admixture and without using prior population information (Hubisz et al., 2009). We conducted 15 independent runs for each value of K (from K = 1 to K = 10) to estimate the most likely number of genetic clusters with 200,000 MCMC cycles, following a burnin step of 100,000 iterations. We retained the 10 runs having the highest likelihood for each value of K and determined the number of genetic clusters that best describes our data according to log probabilities of the data (LnPr(X|K) (Pritchard et al., 2000) and the ΔK method (Evanno et al., 2005), as implemented in structure harvester (Earl & vonHoldt, 2012). We used clumpp v. 1.1.2 and the Greedy algorithm to align multiple runs of structure for the same K value (Jakobsson & Rosenberg, 2007) and distruct v. 1.1 (Rosenberg, 2004) to visualize the individuals´ probabilities of population membership in bar plots.
Complementary to Bayesian clustering analyses, we performed a principal component analysis (PCA) as implemented in the r v. 3.6.3 (R Core Team, 2020) package adegenet (Jombart, 2008). Before running PCAs, we replaced missing data by the mean frequency of the corresponding allele estimated across all samples (Jombart, 2008).

| Demographic analyses
We used the simulation-based approach implemented in fastsimcoal2 (Excoffier et al., 2013) to infer the evolutionary relationships among the three main genetic groups identified by structure and PCA analyses (see Results section) and estimate their divergence times under alternative scenarios of gene flow (e.g., Eaton et al., 2015;Lanier et al., 2015;Ortego et al., 2018). For these analyses, we selected the eight individuals from each genetic group that presented the lowest degree of admixed ancestry (q-value > 0.85; e.g., Zeng et al., 2018).
The specific individuals used for these analyses are indicated in the dataset describing the samples deposited in Figshare  . Colours indicate the main genetic cluster at which populations were assigned according to structure analyses (see Figure 1). Black dots label individuals selected for three-population analyses in fastsimcoal2. Population codes as described in Table 1 (Table S1; e.g., Figure 3).
We also used fastsimcoal2 to test alternative single-population models, which allowed us to investigate in detail the demography of each population while reducing the complexity and high amount of parameters that would require multiepoch models based on several populations (Liu et al., 2020;Zeng et al., 2018). We tested seven alternative single-population demographic models, including one-epoch We used the site frequency spectrum (SFS) and fastsimcoal2 to estimate the composite likelihood of the observed data given a specified model (Excoffier et al., 2013). We calculated the SFS for multiple-population (joint SFS) and single-population demographic analyses considering a single SNP per locus to avoid the effects of linkage disequilibrium. Because we did not include invariable sites in the SFS, we fixed the contemporary effective population size for one of the three demes (NWI) in multiple-population analyses and for each of the populations in single-population analyses to enable the estimation of other parameters in fastsimcoal2 (Excoffier et al., 2013). The effective population size fixed in the model was calculated from the level of nucleotide diversity (π) and estimates of mutation rate per site per generation (μ). Nucleotide diversity (π) was estimated from polymorphic and nonpolymorphic loci using stacks (Table 1). We used the mutation rate per site per generation (2.8 × 10 −9 ) estimated for Drosophila melanogaster (Keightley et al., 2014), which is similar to the spontaneous mutation rate estimated for the butterfly Heliconius melpomene (2.9 × 10 −9 ; Keightley et al., 2015). To remove all missing data for the calculation of the SFS, minimize errors with allele frequency estimates, and maximize the number of retained SNPs, each population group was downsampled to 63% of individuals in multiple-population analyses (i.e., five individuals) and 75% in single-population analyses (i.e., six individuals) using a custom Python script written by Qixin He and available on Dryad (Papadopoulou & Knowles, 2015).
Each model was run 100 independent times considering 100,000-250,000 simulations for the calculation of the composite likelihood, 10-40 expectation-conditional maximization (ECM) cycles, and a stopping criterion of 0.001 (Excoffier et al., 2013;Excoffier & Foll, 2011). We used an information-theoretic model selection approach based on the Akaike's information criterion (AIC) to determine the probability of each model given the observed data (Burnham & Anderson, 2002;e.g., Nater et al., 2015;Thome & Carstens, 2016). After the maximum likelihood was estimated for each model in every replicate, we calculated the AIC scores as detailed in Thome and Carstens (2016). AIC values for each model were rescaled (ΔAIC), calculating the difference between the AIC value of each model and the minimum AIC obtained among all competing models (i.e., the best model has ΔAIC = 0). Point estimates of the different demographic parameters for the best supported model were selected from the run with the highest maximum composite likelihood. Finally, we calculated confidence intervals (based on the percentile method; e.g., de Manuel et al., 2016) of parameter estimates from 100 non-parametric bootstrap data sets. We obtained the 100 pseudoreplicate data sets by sampling loci with replacement, calculated their respective SFS using the same downsampling scheme as for empirical data, and used each pseudoreplicate SFS to re-run the most supported model 100 independent times (e.g., Maier et al., 2019). We used parameter estimates from the highest likelihood run from each pseudoreplicate data set to obtain 95% confidence intervals.

| Exploratory landscape genetic analyses
We generated alternative spatially-explicit isolation-by-resistance (IBR) scenarios of population connectivity and tested which one F I G U R E 3 Demographic model most supported by threepopulation analyses in fastsimcoal2. Parameters include mutationscaled ancestral (θ ANC1 and θ ANC2 ) and contemporary (θ NWI , θ SWI and θ NWM ) effective population sizes, migration rates per generation (m 1 , m 2 and m 3 ), and timing of population divergence (T DIV1 and T DIV2 ) [Colour figure can be viewed at wileyonlinelibrary.com] is better supported by observed data of genetic differentiation (F ST ; Table S3; McRae, 2006). We applied circuit theory and used circuitscape v. 4.0.5 (McRae & Beier, 2007) to calculate resistance distance matrices between all pairs of populations under two main hypothetical scenarios of gene flow: (i) a "flat" landscape in which all cells have equal resistance (resistance = 1). This scenario is analogous to geographical distance but more appropriate for comparison with other competing models also generated with circuitscape (Noguerales et al., 2016;Ortego, Aguirre, et al., 2015); (ii) a spatially heterogeneous landscape defined according to the different resistance offered by three surface categories: sea water, inland surface water, and land/elevation. Elevation data was obtained from WorldClim (https://www.world clim.org/; Hijmans et al., 2005) and inland water was based on the "maximum water extent layer" available at the Global Surface Water Explorer (https://globa l-surfa cewater.appsp ot.com/), which provides detailed information about areas detected as water at least once by remote sensing during a 35- year period   (Pekel et al., 2016). All layers were transformed to 30 arc-sec (ca. 1 km) resolution for subsequent analyses.
In scenarios considering spatially heterogeneous landscapes, we assigned a minimum resistance value (=1) to inland water (i.e., the habitat occupied by the species) and a range of resistance values (1-100,000) to land and sea water (i.e., the two hypothetical barriers to gene flow). We tested all possible combinations of resistance values for land and sea water, which allowed us to identify the specific resistance values for these landscape features that provide the best fit to our data on genetic differentiation (e.g., Andrew et al., 2012;González-Serna et al., 2018;Ortego, Gugger, et al., 2015). In the scenario considering elevation rather than a homogeneous resistance value for land, the range of resistance values explored was only applied to sea water. This resulted in a total of 41 resistance scenarios plus the flat landscape (see Tables S4 and S5). Under each of these scenarios, we calculated resistance distance matrices between all pairs of populations considering an eight-neighbour cell connection scheme in circuitscape (McRae, 2006). All GIS calculations were performed using arcmap v. 10.5 (ESRI, Redlands, CA, USA).
We determined the fit of the different landscape resistance models to our observed data of genetic differentiation (pairwise F ST values; Table S3) using multiple matrix regressions with randomization (MMRR) as implemented in r v. 3.6.3 (Wang, 2013). The final model was selected following a backward procedure, initially fitting all explanatory terms and progressively eliminating nonsignificant variables until all retained variables were significant. The significance of the variables excluded from the model was tested again until no additional variable reached significance (e.g., Ortego, Gugger, et al., 2015). In order to explore the impact of spatial scale and the geographical separation of populations in different continental landmasses on the obtained inferences, we re-ran the analyses using four different data sets that excluded/included Moroccan populations (LARA and MERJ) and the highly isolated population from Aveiro (AVEI) (see Tables S4 and S5).

| Genomic data and basic genetic statistics
The average number of reads retained per individual after the different quality filtering steps was 1,639,507 (range = 666,878-3,347,506 reads; Figure S1). On average, this represented 85% (range = 73%-88%) of the total number of reads recovered for each individual ( Figure S1). After filtering loci (see Methods S1), the final data set contained 2490 SNPs with a proportion of missing data of 12.58%. Levels of genetic diversity (π) in populations of T. v. verticalis ranged between 0.0021 and 0.0028, with the highest estimates in populations from southwestern Iberia and northwestern Morocco and the lowest estimate in the highly isolated population from northwestern Iberia (AVEI; Table 1; Figure 1a). Accordingly, AVEI population presented significantly lower levels of genetic diversity (24% lower) than the rest of the study populations (one-sample t test: t = 24.32, p < .001; Figure 1a).

F I G U R E 4
Single-population demographic models tested using fastsimcoal2. Parameters include mutation-scaled historical (θ ANC , θ BOT1 and θ BOT2 ) and contemporary (θ CON ) effective population sizes, exponential growth rates (r, in models C, E and G), and timing of population size change (T BOT1 , T BOT2 , and T EXP ) 3.2 | Quantifying genetic structure structure analyses identified the most likely number of clusters as K = 2 according to the ΔK criterion, but LnPr(X|K) reached a plateau at K = 4 ( Figure S2). For K = 2, the two genetic clusters separated the highly isolated northwestern Iberian population (AVEI) from the rest of the populations (Figure 1b). structure analyses for K = 3 split

| Demographic analyses
Three-population demographic analyses in fastsimcoal2 supported the model considering a full-migration matrix and a sister relationship between the lineages present in northwestern and southwestern Iberia (Table S1; Figure 3). The second most supported model (ΔAIC = 2.12) was the one with the same topology, but excluding gene flow between populations from northwestern and southwestern Iberia (Table S1). Accordingly, estimates of gene flow between these two populations (m 1 ) were significantly smaller (i.e., 95% CI did not overlap) than those estimated among the two other population pairs (m 2 and m 3 ; Table 2; Figure 3). The initial split between the lineage of T. verticalis present in northwestern Morocco and the most recent common ancestor of the lineages introduced in the Iberian Peninsula (T DIV1 ) was estimated to happen ca. 4.7 M generations ago, whereas the split between the lineages introduced in northwestern and southwestern Iberia (T DIV2 ) took place ca. 95k generations ago (Table 2). Considering the generation time of the species (6 generations/year; Céspedes, Coccia, et al., 2019), these estimates correspond to ca. 787 and 16 ka, respectively (Table 2).
Single-population analyses identified the three-epoch model  (Table 3). This strong bottleneck was followed by very recent (T EXP < 60 years ago) instantaneous demographic expansion that led to a considerable recovery of contemporary effective population sizes (θ CON ) to >25% of θ ANC (Table 3).

| Exploratory landscape genetic analyses
Pairwise F ST values ranged between 0.043 (for the nearby populations of PEDR and CETI in southwestern Iberia) and 0.113 (for the geographically distant populations of AVEI in northwestern Iberia and LARA in northwestern Morocco) (Table S3). Most pairs of populations were significantly differentiated, with the only exceptions involving some pairs of populations from southwestern Iberia (Table   S3).
Landscape genetic analyses revealed contrasting effects of the different landscape features depending on the subset of populations considered (Tables S4 and S5; Figure 5). When the highly isolated population from northwestern Iberia (AVEI) was included in the analyses, genetic differentiation was exclusively explained by resistance offered by a flat landscape (all populations: R 2 = .735, t = 9.72, TA B L E 2 Parameters inferred from coalescent simulations with fastsimcoal2 under the most supported model (illustrated in Figure 3) (Céspedes, Coccia, et al., 2019). Effective population size of the NWI population (θ NWI ) was calculated from the level of nucleotide diversity (π; Table 1) and fixed in fastsimcoal2 analyses to enable the estimation of other parameters (see the Materials and Methods section for further details). The site frequency spectrum used in the analysis was based on 2345 variable SNPs. p = .006, Figure 5a; Iberian populations: R 2 = .754, t = 7.64, p = .048, Figure 5c). Note that this scenario is equivalent to geographical distance or an isolation-by-distance model (e.g., Noguerales et al., 2016). When AVEI was excluded from the analyses, the effects of landscape heterogeneity on spatial patterns of gene flow became evident. Analyses including Moroccan and Iberian populations (except AVEI) revealed that oceanic waters separating populations in the two continents represent a major barrier to dispersal relative to land and inland water (Figure 5b; Table S4). Accordingly, genetic differentiation was only significantly associated with resistance distances in those scenarios in which sea water offered ≥100 times more resistance than land (Table S4). The most supported model was the one considering that sea water offers 10 5 times more resistance to gene flow than land and inland water (R 2 = .418, t = 4.32, p = .027; Table S4; Figure 5b). When analyses focused on populations from southwestern Iberia, not separated by oceanic waters, the opposite pattern arose: genetic differentiation was only significantly associated with resistance distances in scenarios in which land offered more resistance than sea water and inland water (Table S5). In this case, the most supported model was the one considering that land offers 10 times more resistance to gene flow than sea and inland water (R 2 = .309, t = 2.41, p = .005; Table S5; Figure 5d).

| DISCUSS ION
Our  (Carbonell et al., 2017;Guareschi et al., 2013), which can explain the extensive natural range of the species along the eastern coast of North America and its extraordinary ability to colonize temporary wetlands on an annual basis . generations per year for the species (Céspedes, Coccia, et al., 2019). Contemporary effective population size for each population (θ CON ) was calculated from their respective levels of nucleotide diversity (π;

| Genetic makeup of introduced populations
Table 1) and fixed in fastsimcoal2 analyses to enable the estimation of other parameters (see the Materials and Methods section for further details). Population codes as described in Table 1. space, rejecting the hypothesis of genetic panmixia and supporting significant genetic structure within the invasive range of T. verticalis in the Iberian Peninsula and Morocco (Figures 1 and 2). These analyses showed the presence of three main genetic clusters separating northwestern Iberian, southwestern Iberian, and northwestern Moroccan populations. Coalescent-based estimation of divergence times for the three genetic groups indicate that these predate the onset of the Holocene (>12 ka; Table 2), which is incompatible with in situ geographical diversification and points to several introduction events involving different source populations with contrasting evolutionary histories (Roman & Darling, 2007). Previous studies have hypothesized that the species probably reached invaded ranges through ballast water or accompanying introductions of exotic fish species, specially Gambusia holbrooki (Guareschi et al., 2013;Sala & Boix, 2005;Vidal et al., 2010), these being dispersal pathways compatible with introductions from disparate source populations (Bailey, 2015). Our results indicate that the current distribution of T. verticalis in the region has been aided by multiple introductions from North America and suggest that the species might not spread so easily through other potentially suitable areas in the Western Palearctic (e.g., coastal wetlands in France, Italy, and Northern Africa) as suggested by species distribution modelling if anthropogenic dispersal vectors are disrupted (Guareschi et al., 2013).
Admixture among the three main genetic clusters revealed by structure analyses support secondary contact among gene pools that have probably remained isolated for extended periods of time (e.g., Stepien et al., 2005). Genetic admixture is expected to although other factors, such as differences in the initial size of introduced populations, availability of suitable habitats or extinctionrecolonization dynamics, cannot be discarded (e.g., Braasch et al., 2019;Erfmeier et al., 2013;Kaňuch et al., 2021). For instance, we sampled mainly permanent waters that may have held stable populations of T. verticalis for over a decade, but also temporary waters where the species experiences local extinctions and must recolonize when inundation occurs again . A good example of the latter is CARA, a site which dries out each summer, probably explaining its lower levels of genetic diversity compared to other nearby sites (Figure 1).  Table S2). These analyses F I G U R E 5 Heatmaps showing the goodness of fit for alternative models of landscape resistance. Colours show the coefficient of determination (R 2 ) for models analysing genetic differentiation in relation to isolation-by-resistance (IBR) distance matrices, calculated with circuitscape considering different combinations (from 1 to 100,000) of resistance values for sea water and land. also revealed that populations have restored a big portion (>30%) of their respective ancestral effective population sizes after demographic expansions, supporting the capacity of the species to quickly recover from dramatic demographic bottlenecks during the invasion process (Table 3). Considering the six generations per year estimated for the species in the invasive range (Céspedes, Coccia, et al., 2019), explosive demographic expansions were estimated to take place 24-53 years ago (95% CI = 19-85 years ago; Table 3).

| Demographic fate of introduced populations
Trichocorixa verticalis was first recorded in Andalucía region (Spain) in 2001 (Rodríguez-Pérez et al., 2009) and then identified in material collected from the Algarve region (Portugal) in 1997 (Sala & Boix, 2005), which establishes the first record of the species in the area 21 years before our genetic sampling (although the species is likely to have been overlooked for years before its discovery). Despite the expected time lag between introduction and detection of invasive species and the fact that T. verticalis is not a conspicuous species, and is unlikely to be identified by nonexperts in the group, estimates of demographic expansions are remarkably congruent with the first record for the species in the region. This congruence between the first detection of T. verticalis and genomic-based inference of the time of explosive demographic expansions suggests extraordinary short time lags between arrival and spread, which is in agreement with previous monitoring studies reporting a rapid local expansion of the species in the region, and points to its high capacity to become invasive once it disperses into non-native suitable areas (Fouzi et al., 2020;Rodríguez-Pérez et al., 2009). High physiological plasticity , extraordinary fecundity (Carbonell et al., 2016), short generation time (<60 days; Céspedes, Coccia, et al., 2019), all year round reproduction (Rodríguez-Pérez et al., 2009), and the capacity of the species to out-compete native corixids in some habitats (Carbonell et al., 2017(Carbonell et al., , 2020Rodríguez-Pérez et al., 2009) confer to T. verticalis all the attributes to be a successful invader, which might explain the explosive population expansions inferred by genomic-based demographic reconstructions (Guareschi et al., 2013). The simultaneous arrival of propagules from different source populations might have also contributed to a rapid increase in local levels of genetic diversity via admixture, a phenomenon that could potentially explain inferred demographic expansions shortly after colonization and the fast spread of the species documented in the Western Mediterranean (Guareschi et al., 2013;Rodríguez-Pérez et al., 2009;see also, Frankham, 2005;Kolbe et al., 2004;Roman & Darling, 2007).

| Insights from exploratory landscape genetic analyses
Our hierarchical landscape genetic analyses suggest that sea and land surfaces might act as barriers to dispersal, which is not surprising considering the impossibility of the species to form viable populations in these habitats. In this line, some pairs of coastal populations separated by >170 km (e.g., ODIE-PEDR) showed no significant genetic differentiation (Table S3), suggesting widespread dispersal across large geographical scales. This may reflect frequent stepping stone dispersal along the coast of Andalusia and Algarve regions, where there is little distance between sites occupied by the species owing to a series of low-lying deltas holding salt work complexes, inland lagoons and fish ponds. This contrasts with a near absence of suitable habitats between southern (PORT) and northern Portugal (AVEI) and between coastal (PORT, ODIE, CETI, and PEDR) and inland populations (CARA and ZARR) in southern Iberia, which might explain the significant genetic differentiation between these groups of populations (Table S3). It should be noted, however, that landscape genetic analyses and spatial patterns of genetic differentiation must be interpreted with extreme caution given the different evolutionary origins of populations (i.e., multiple introduction events) from the three main invaded areas in the Iberian Peninsula and Morocco. The study populations are most probably not yet at equilibrium and, thus, historical contingency rather than contemporary gene flow probably still dominate spatial patterns of genetic diversity and structure (Fitzpatrick et al., 2012). For instance, populations located in southern Iberia and northwestern Morocco are separated by sea water but their recent (probably < 50 years; Günther, 2004;Sala & Boix, 2005) and independent colonization histories probably  (Gunter & Christmas, 1959;Hutchinson, 1931) and, thus, passive dispersal through marine currents and possibly severe storms might play a role in the spread of the species at regional scales (Sailer, 1948;Sala & Boix, 2005); (ii) Low wing loading and high aspect ratios suggest that T. verticalis is probably a good flyer in comparison with native syntopic corixid species, which could facilitate dispersal among nearby water bodies and explain low levels of genetic differentiation among populations (F ST < 0.12; Table S3); (iii) Recent empirical evidence also suggests that T. verticalis could be passively dispersed by attachment of eggs to legs of waterbirds (i.e., epizoochory), a phenomenon that can greatly increase the dispersal distances that can be travelled by the species. Zoochory may potentially explain the recurrent colonization of inland temporary ponds located far away from their likely source populations in permanently flooded coastal wetlands (Carbonell et al., 2016(Carbonell et al., , 2021; (iv) Human-transportation (e.g., via ballast water) within the invasive range could still also facilitate gene flow between distant populations and contribute to genetic admixture among introduced lineages in different areas (Medley et al., 2015;Robinet et al., 2009;Valls et al., 2016); Finally, we cannot rule out the possibility that introduced populations have experienced local admixture after independent introduction events involving different proportions of multiple lineages, or even that some source populations from the native range already had admixed gene pools resulted from natural (i.e., contact/hybrid zones; Hewitt, 1988)

| CON CLUS IONS
Our analyses strongly support that multiple introduction events are responsible of the invasion of T. verticalis in the Western Mediterranean. This suggests that the expansion of the species at regional and global scales could be avoided or at least limited if current environmental regulations of most countries regarding the correct management of ballast water exchange and the prohibition of releasing alien aquatic organisms are implemented in practice and effectively enforced (e.g., Gollasch et al., 2007; see also Ormsby & Brenton-Rule, 2017). Genomic-based inferences of explosive demographic expansions immediately or shortly after arrival to non-native areas indicate that the lag phase between introduction and establishment is probably very short, which provides little opportunity for early detection, management or eradication (Coutts et al., 2018;Simberloff, 2003). Future genomic analyses of native range populations could help to infer the genetic makeup, sources and geographic origins of introduced populations, reconstruct the invasion history, and understand with more confidence the specific processes underlying the demographic dynamics of introduced populations (Cristescu, 2015;Estoup & Guillemaud, 2010). Detailed landscape genetic analyses in the native range could also contribute to confidently identify corridors/barriers to dispersal for the species and gain insights useful for understanding patterns of genetic connectivity and metapopulation dynamics in invaded areas (Handley et al., 2011). Likewise, continuous genetic monitoring of introduced populations might help to document the arrival of individuals from different source populations and provide a better understanding of the demographic and evolutionary dynamics (i.e., adaptation processes to water chemistry, colonization of temporary ponds, etc.) of recently established populations of the species (Carbonell et al., 2012;Céspedes, Valdecasas, et al., 2019;Coccia et al., 2016;Frisch et al., 2021).

ACK N OWLED G EM ENTS
We are grateful to Amparo Hidalgo-Galiana for preparing the

CO N FLI C T O F I NTE R E S T
The authors have no conflicts of interest to declare.