Brazier Thomas

Cherif Emira
email: emira.cherif@ird.fr

Martin Jean-François

Gilles André

Blanchet Simon

ORCID Zhao Yahui

Combe Marine

ORCID Mccairns R J Scott

& Gozlan

Emira Cherif

Thomas Brazier

Jean-François Orcid

Simon Blanchet

Orcid

Rodolphe E Gozlan

R E G Wrote

The influence of native populations' genetic history on the reconstruction of invasion routes: The case of a highly invasive aquatic species

Keywords: Invasion pathways, aquatic biological invasion, demographic inference, population modeling, population genomics, Approximate Bayesian Computation

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

The influence of native populations' genetic history on the reconstruction of invasion routes: the case of a highly invasive aquatic species

Thomas Brazier, Emira Cherif, Jean-François Martin, André Gilles, Simon

Blanchet, Yahui Zhao, Marine Combe, R. J. S. Mccairns, Rodolphe E. Gozlan

To cite this version:

Thomas Brazier, Emira Cherif, Jean-François Martin, André Gilles, Simon Blanchet, et al.. The influence of native populations' genetic history on the reconstruction of invasion routes: the case of a highly invasive aquatic species. Biological Invasions, 2022Invasions, , 24, pp.2399Invasions, -2420. 10. 10.1007/s10530-022-02787-6. hal-03656986

INTRODUCTION

Invasive species are significant contributors to global change and often lead to biotic homogenization and biodiversity losses [START_REF] Iii | Consequences of changing biodiversity[END_REF][START_REF] Clavero | Invasive species are a leading cause of animal extinctions[END_REF][START_REF] Didham | Are invasive species the drivers of ecological change?[END_REF][START_REF] Villéger | Homogenization patterns of the world's freshwater fish faunas[END_REF][START_REF] Simberloff | Biological invasions: much progress plus several controversies[END_REF]. Globalization of trade and growth of worldwide transportation are the main drivers of non-native species introductions [START_REF] Hulme | Trade, transport and trouble: managing invasive species pathways in an era of globalization[END_REF]. Major ecological impacts from biological invasions include destabilization of trophic networks [START_REF] Stiers | Impact of three aquatic invasive species on native plants and macroinvertebrates in temperate ponds[END_REF][START_REF] Gallardo | Global ecological impacts of invasive species in aquatic ecosystems[END_REF], competition for resources and habitats [START_REF] Graebner | Invasive species grows faster, competes better, and shows greater evolution toward increased seed size and growth than exotic noninvasive congeners[END_REF][START_REF] Perdereau | Competition between invasive and indigenous species: an insular case study of subterranean termites[END_REF], predation on native species [START_REF] Salo | Alien predators are more dangerous than native predators to prey populations[END_REF]) and the transmission of novel infectious pathogens [START_REF] Crowl | The spread of invasive species and infectious disease as drivers of ecosystem change[END_REF]. Invasions can also have evolutionary consequences, such as genome introgression or altering selection pressures on native species [START_REF] Crispo | Broken barriers: Humaninduced changes to gene flow and introgression in animals: An examination of the ways in which humans increase genetic exchange among populations and species and the consequences for biodiversity[END_REF][START_REF] Mooney | The evolutionary impact of invasive species[END_REF][START_REF] Sinama | Nonhomogeneous combination of two porous genomes induces complex body shape trajectories in cyprinid hybrids[END_REF]Philips et al. 2006). Yet invasions offer an incredible framework to study adaptation to new environmental conditions and to understand how small introduced populations with a supposedly low level of genetic diversity manage to colonize large areas [START_REF] Bossdorf | Phenotypic and genetic differentiation between native and introduced plant populations[END_REF][START_REF] Peischl | Expansion load: recessive mutations and the role of standing genetic variation[END_REF][START_REF] Sax | Ecological and evolutionary insights from species invasions[END_REF][START_REF] Roman | Paradox lost: genetic diversity and the success of aquatic invasions[END_REF][START_REF] Facon | A general ecoevolutionary framework for understanding bioinvasions[END_REF], as well as to understand the spread and emergence of pathogens cointroduced together with their non-native hosts. Reconstructing the invasion history of a species is challenging. Typically, few individuals may constitute the initial colonizing group, but data on the amount of propagule pressure as well as its origin and time of first introduction are often lacking. Historical records of introduction events, census data or environmental monitoring projects often miss the first stages of invasion [START_REF] Holsbeek | A cryptic invasion within an invasion and widespread introgression in the European water frog complex: consequences of uncontrolled commercial trade and weak international legislation[END_REF][START_REF] Mergeay | Cryptic invasion and dispersal of an American Daphnia in East Africa[END_REF]. Molecular ecology allows the reconstruction of events that have not been directly observed but have left a genomic signature [START_REF] Cristescu | Genetic reconstructions of invasion history[END_REF]. Through the use of new demographic inference methods that permit deviations from the assumptions of the mutation-drift equilibrium model (which is often the case for recent colonization events), it is now possible to identify source populations or even to model the short-term demo-genetic processes underpinning invasion history [START_REF] Beichman | Using genomic data to infer historic population dynamics of nonmodel organisms[END_REF][START_REF] Cabrera | Inferring past demographic changes from contemporary genetic data: A simulation based evaluation of the ABC methods implemented in DIYABC[END_REF][START_REF] Shafer | Demographic inferences using short-read genomic data in an approximate Bayesian computation framework: In silico evaluation of power, biases and proof of concept in Atlantic walrus[END_REF][START_REF] Estoup | Reconstructing routes of invasion using genetic data : Why, how and so what?[END_REF]. Recent population genetic studies on invasion pathways have shown that invasion patterns may be more complex than previously thought [START_REF] Lombaert | Complementarity of statistical treatments to reconstruct worldwide routes of invasion: the case of the Asian ladybird Harmonia axyridis[END_REF]. A single introduction event is rare, and multiple introductions and/or admixture between source populations are often the cornerstone of successful invasions [START_REF] Estoup | Is there a genetic paradox of biological invasion?[END_REF][START_REF] Roman | Paradox lost: genetic diversity and the success of aquatic invasions[END_REF]. Since introduced populations inherit genetic variation primarily shaped by long-term evolutionary processes in the native area, incorporating this information in the interpretation of the results is important to avoid erroneous conclusions [START_REF] Rius | Phylogeography and the Description of Geographic Patterns in Invasion Genomics[END_REF].

The topmouth gudgeon, Pseudorasbora parva (Temminck & Schlegel), is a small freshwater cyprinid fish. More importantly, it is an interesting study system for biological invasion for several reasons: i) it is a highly successful invader that has spread and established very quickly in diverse environments throughout the world, ii) it has a well-documented history and iii) its invasion is associated with a pathogen co-introduction. It also has a very large native distribution in Eastern Asia (Eastern China, Taiwan, Korea and Japan) with broad environmental tolerance from continental climates to tropical ones [START_REF] Gozlan | Pseudorasbora parva temminck and schlegel (topmouth gudgeon). Handb Glob Freshw Invasive Species[END_REF][START_REF] Zhang | Species diversity and distribution of inland fishes in China[END_REF]. It was initially accidentally released outside its native range via aquaculture partnership exchanges of Chinese carp between China and former countries of the Soviet Bloc [START_REF] Gozlan | Pancontinental invasion of Pseudorasbora parva: towards a better understanding of freshwater fish invasions[END_REF]. During the 1960s, multiple introductions of P. parva took place all around the Black Sea area, followed by further introductions in the 80s in Eurasia and North Africa. However, introduction records lack precise information on the source populations.

After these initial phases of human-made introductions, natural local colonization of entire river networks occurred across major European rivers to the Middle East. Its life-history traits, which include early maturity (1 year) coupled with nest guarding behavior that ensures an increased likelihood of survival for juveniles and short longevity, have been identified as key traits to explain its rapid establishment and spread [START_REF] Gozlan | Pancontinental invasion of Pseudorasbora parva: towards a better understanding of freshwater fish invasions[END_REF][START_REF] Gozlan | Pseudorasbora parva temminck and schlegel (topmouth gudgeon). Handb Glob Freshw Invasive Species[END_REF], Gozlan et al. 2020), although its long-term invasion success depends on high genetic diversity to drive adaptation to the new environment [START_REF] Roman | Paradox lost: genetic diversity and the success of aquatic invasions[END_REF]. A major biodiversity problem associated with P. parva is the fungus-like pathogen, the rosette agent Sphaerothecum destruens, for which it acts as a vector, and that has subsequently spread to invaded areas causing mortality in a large number of native freshwater fish species in Europe [START_REF] Andreou | Associated disease risk from the introduced generalist pathogen Sphaerothecum destruens: management and policy implications[END_REF][START_REF] Combe | The rise of the rosette agent in Europe: An epidemiological enigma[END_REF].

Genetic structure and phylogeography, based upon classical mitochondrial and nuclear markers (microsatellites), have been extensively described in both native and invasive P.

parva populations, suggesting the existence of two genetic lineages within non-native populations in Europe [START_REF] Hardouin | Reconciling the biogeography of an invader through recent and historic genetic patterns: the case of topmouth gudgeon Pseudorasbora parva[END_REF][START_REF] Simon | Invasive cyprinid fish in Europe originate from the single introduction of an admixed source population followed by a complex pattern of spread[END_REF][START_REF] Simon | Human induced stepping-stone colonisation of an admixed founder population: the spread of topmouth gudgeon (Pseudorasbora parva) in Europe[END_REF]. This general pattern has been recently confirmed by a study using 13,785 single nucleotide polymorphisms in Slovakian and Turkish introduced populations (Baltazar- [START_REF] Soares | Genomic footprints of a biological invasion: Introduction from Asia and dispersal in Europe of the topmouth gudgeon (Pseudorasbora parva)[END_REF]). Yet, the finescale genetic structure of the native populations and the gene flow that has shaped it, as well as the demographic dynamics underpinning the invasion pathways, remain unresolved.

Assessing the source populations within the native range by extensive sampling, along with high-throughput genotyping, is crucial to detect high-resolution signatures of demographic history [START_REF] Muirhead | Identifying the source of species invasions: sampling intensity vs. genetic diversity[END_REF][START_REF] Beichman | Using genomic data to infer historic population dynamics of nonmodel organisms[END_REF][START_REF] Shafer | Demographic inferences using short-read genomic data in an approximate Bayesian computation framework: In silico evaluation of power, biases and proof of concept in Atlantic walrus[END_REF]. Although such large datasets require advanced statistical methods that are computationally demanding, Approximate Bayesian Computation (ABC) and associated Machine Learning algorithms take advantage of large sets of summary statistics, exploring huge parameter spaces with reduced computational effort [START_REF] Beaumont | Approximate Bayesian computation in population genetics[END_REF][START_REF] Cabrera | Inferring past demographic changes from contemporary genetic data: A simulation based evaluation of the ABC methods implemented in DIYABC[END_REF][START_REF] Pudlo | Reliable ABC model choice via random forests[END_REF]Raynal et al. 2017;Rey et al. 2015).

Here we took advantage of genomic signatures, ABC models and Machine Learning algorithms to i) characterize the genetic structure of native populations and test whether they aggregate on homogeneous and coherent demes and ii) retrace P. parva's introduction history from Asia to Europe. To do so, we first applied genetic clustering methods and explicit tests of admixture to describe the phylogeography of both native and invasive ranges. We then identified putative source populations of introduced demes with population assignment tests. We finally used ABC model-based procedures to infer the most probable demographic scenarios of P. parva invasion and reconstruct an integrative biological invasion history.

MATERIALS AND METHODS

Sampling material

We sampled P. parva from 21 discrete sites across its overall distribution in Asia, comprising sixteen different river catchments in the historical Chinese native range (Figure S1, Table 1, sites numbered from S1 to S18). Additional samples were obtained from Japan and neighboring Asian invasive populations in Tibet and South-East China (sites: S19, S20). We also sampled P. parva across thirteen locations within the invasive European and Middle-Eastern range (Figure S1, Table 1). The sampling methods included fish traps, electric fishing and micro-mesh seine netting where appropriate. In the field, the fish were euthanized with an overdose of anesthetic, initially preserved in ethanol and subsequently stored at -70 °C. Each fish was measured and a fin-clip taken and stored in ethanol.

GBS sequencing

In total, 858 DNA samples (746 individuals in total, with 75 individuals replicated across at least one sequencing lane) were genotyped for single nucleotide polymorphism (SNP) markers by first digesting genomic DNA with PstI, followed by genotyping-by-sequencing (GBS), yielding an average of 2,702,000 raw sequencing reads per sample. SNP calling was performed using programs comprising the Stacks (v1.46) bioinformatics de novo pipeline [START_REF] Catchen | Stacks: an analysis tool set for population genomics[END_REF]. Replicated individuals were used to estimating genotyping error and only markers with less than 1% errors were retained. A minimum read depth of 20 was required for each marker, ultimately yielding 3999 validated SNP markers. To reduce linkage disequilibrium amongst markers, only a single SNP (first position) was called for each locus stack. To prevent biases due to large proportions of missing data, populations with more than 70% missing data, individuals with more than 60% missing data and loci with more than 45% of missing data were removed from the dataset. After trimming, the final dataset contained 300 individuals from eighteen sites in the historic Asian range (including invasive populations in Tibet and South-West China) and 168 individuals from eleven invasive sites across Europe, Turkey and Iran, with on average 16 individuals per sampled site (Figure S1, Table 1).

Genetic diversity

Statistics of genetic diversity were estimated with the 'adegenet' and 'hierfstat' R packages [START_REF] Goudet | Hierfstat, a package for R to compute and test hierarchical F-statistics[END_REF][START_REF] Goudet | hierfstat: estimation and tests of hierarchical F-statistics[END_REF][START_REF] Jombart | adegenet: a R package for the multivariate analysis of genetic markers[END_REF][START_REF] Jombart | adegenet 1.3-1: new tools for the analysis of genome-wide SNP data[END_REF] for R version 3.5.3 (R Core Team, 2019) (Table S1). Sensitivity analyses on summary statistics showed that biases were minimized without loss of power when estimates were inferred from 2,000 to 3,000 loci with the lowest proportion of missing data (Figure S2). Consequently, the dataset for inferring population assignment consisted of 2,112 SNPs to reduce computational load (maximum 45% missing data per locus); however, up to 3000 SNPs (maximum ~50% missing data) were retained for other analyses to keep most of the genetic information. Simulations have shown that even loci with high levels of missing data (> 50%)

can retain meaningful information [START_REF] Huang | Unforeseen Consequences of Excluding Missing Data from Next-Generation Sequences: Simulation Study of RAD Sequences[END_REF], as it has been empirically assessed for population structure [START_REF] Chan | Species delimitation with gene flow: A methodological comparison and population genomics approach to elucidate cryptic species boundaries in Malaysian Torrent Frogs[END_REF].

Genetic clustering

Genetic clustering of sampled sites was assessed independently in the Asian range (native and non-native sites) and within the invasive European range. Putative demes were defined as groups of individuals sharing a gene pool. Results from an iterative K-means method and a model-based Bayesian clustering method were compared for cross-validation.

Discriminant Analysis of Principal Components (DAPC), implemented in the R package 'adegenet' [START_REF] Jombart | adegenet: a R package for the multivariate analysis of genetic markers[END_REF][START_REF] Jombart | adegenet 1.3-1: new tools for the analysis of genome-wide SNP data[END_REF] was first used to determine genetic clusters within native and invasive regions and then to predict membership of individuals within invasive sites to native clusters. The most probable number of genetic clusters, K, was searched within the distribution of K among 1,000 independent clustering iterations, based on the 'goodfit' and 'min' criteria of the Bayesian Information Criterion (BIC). Overfitting was prevented by a cross-validation procedure to define the optimal number of Principal Components required to discriminate amongst these K clusters.

Additionally, we used STRUCTURE 2.3 to infer the number of genetic clusters and to estimate admixture between them [START_REF] Pritchard | Inference of population structure using multilocus genotype data[END_REF]. The admixture model with correlated allele frequencies was parameterized with a fixed Lambda value (parameter of the allele frequencies distribution) directly estimated from the data. Sampling location was set as a prior to improve inferences on weak genetic structure. STRUCTURE was performed for K=1-21 (number of sampled sites + 3) for Asian sites and K=1-14 for European and Middle-Eastern sites. Twenty replicates were computed for each K value, with 100,000 sampling iterations after a burn-in of 100,000. The most probable K was assessed after considering the smallest value of K minimizing differences of likelihood (i.e. the plateau method of [START_REF] Pritchard | Inference of population structure using multilocus genotype data[END_REF], the highest value of ΔK given K given using Evanno's method, computed with STRUCTURE HARVESTER [START_REF] Earl | STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method[END_REF][START_REF] Evanno | Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study[END_REF] Additionally, within-chain convergence for parameter estimates (i.e. ancestry coefficient Alpha and Ln Likelihood) was assessed with diagnostics implemented in the R package 'coda' [START_REF] Plummer | CODA: convergence diagnosis and output analysis for MCMC[END_REF].

Finally, based on the observed partitioning of genetic variance, sampled sites were pooled into demes (also referred to as populations) that made sense biologically (i.e. continuous gene flow among sites) and geographically. Admixture proportions inferred with STRUCTURE helped to delineate putative demes and were computed as the mean of the major ancestry coefficients (Q) in a given site. Sampled sites with a mean major ancestry coefficient (Q) under 70% were considered as admixed sites. Sampled sites that did not cluster well based on genetic markers were clustered into putative demes based on geographical and historical data. The consistency of our final clustering of sampled sites in genetic populations was assessed by a hierarchical AMOVA implemented in the 'poppr' R package [START_REF] Kamvar | Poppr: an R package for genetic analysis of populations with clonal, partially clonal, and/or sexual reproduction[END_REF][START_REF] Kamvar | Novel R tools for analysis of genome-wide population genetic data with emphasis on clonality[END_REF] and the significance of variance proportions was tested with 1,000 random permutations (Table S3). Maps were drawn from the R package 'maps' with the 'world' database [START_REF] Becker | maps: Draw Geographical Maps. R package version 3.3. 0[END_REF]. The main river network was drawn from 'RNaturalEarth' (South 2017).

Assignment to source populations with supervised machine learning

The R package 'AssignPOP' was used to assign invasive individuals to candidate source populations (i.e. native demes) with a machine learning classification algorithm (Chen et The training error rate was assessed as the reassignment success of known individuals to their source population. Source populations of invasive individuals were predicted from posterior assignment probabilities. Only individuals for whom the highest probability was at least twice that of the second were retained as confident assignments.

Inference of invasion history with ABC

Approximate Bayesian Computation (ABC) was used to infer past demographic events shaping contemporary genetic diversity (Supporting information). We simulated large datasets under various invasion and admixture demographic scenarios and estimated the probability that data were observed under a given demographic scenario [START_REF] Estoup | Estimation of demo-genetic model probabilities with Approximate Bayesian Computation using linear discriminant analysis on summary statistics[END_REF]. Scenarios were designed in a two-step hierarchical procedure of increasing complexity, with the second step derived from findings at step 1 (Figure S3, S4; see supplementary method for details). In the first step, we tested 3 independent sets of scenarios about the origins in the native range of 3 independent invasive demes (Western Europe, Eastern Europe and Iran; see Results). For each independent invasive deme, we tested if the source population was one of the three candidates or an admixture between the three candidates. To avoid the trap of infinite combinations of exhaustive scenarios, source populations predicted with population assignment served as initial candidate source populations. Choice of candidate source populations was also cross-validated with a Maximum-likelihood phylogenetic tree (Figure S8) inferred by the program TreeMix version 1.13 [START_REF] Pickrell | Inference of Population Splits and Mixtures from Genome-Wide Allele Frequency Data[END_REF]. As the Italian sample's origin was ambiguous in previous genetic clustering results, we performed the ABC analysis with and without Italian samples.

The comparison of both replicated scenarios allowed us to assess the sensitivity of the selected scenario to Italian individuals (Table 2). The second step was dedicated to resolving a more complex invasion pattern, comparing competing worldwide invasion pathways encompassing all invasive demes. Source populations of invasive demes were found to be admixed (results of step 1, described below). Therefore, plausible hypotheses to test at step 2 were: (1) three independent introductions from three independent admixed source populations;

(2) three independent introductions from a single admixed native population; or (3) a single continental introduction from an admixed native population.

Demo-genetic scenarios were simulated and summary statistics estimated, with DIYABC version 2.1 [START_REF] Cornuet | DIYABC v2.0: a software to make approximate Bayesian computation inferences about population history using single nucleotide polymorphism, DNA sequence and microsatellite data[END_REF]. Prior probabilities of scenarios were set to uniform.

Parameter prior distributions were first set to a biologically reasonable range of values, then confidence in priors checked (via DIYABC test of goodness of fit), with distributions refined iteratively until the simulated scenarios fit the data (final parameter spaces given in supplementary methods). For each scenario, 10,000 simulations were conducted, estimating all available summary statistics in DIYABC for each [START_REF] Cornuet | DIYABC v2.0: a software to make approximate Bayesian computation inferences about population history using single nucleotide polymorphism, DNA sequence and microsatellite data[END_REF]). The best model was selected using two machine-learning algorithms trained on the simulated datasets: a

Neural Network algorithm implemented in the R package 'abc' [START_REF] Csilléry | abc: an R package for approximate Bayesian computation (ABC)[END_REF]) and a Random Forest algorithm in the 'abcrf' package [START_REF] Pudlo | Reliable ABC model choice via random forests[END_REF]. One thousand Neural Networks were trained with 5 to 12 units of hidden neural layers. Parameter sets were weighted by an Epanechnikov kernel. The tolerance rate was set to 0.2. Other configuration values were set to the default value. Power of ABC model selection with Neural Network was evaluated by leave-one-out cross-validation repeated 100 times. Independently, 1,000 trees were grown in the Random Forest training set, with linear discriminant analysis scores added to summary statistics when it reduced the prior error rate. Power of Random Forest was evaluated by out-of-bag prior misclassification error rate. Lastly, the quality of the selected scenario was checked by comparing the marginal posterior predictive distributions to the observed values of the summary statistics. The marginal posterior predictive distribution was computed from 10,000 simulations under the selected scenario, parameterized with estimated posterior parameter distributions as priors.

The complexity of the selected scenario at step 2, combined with a restricted number of markers (3,000 SNPs), reduced the power to infer demographic parameters jointly with a reasonable confidence interval. Hence, demographic parameters of the selected scenarios were estimated at step 1 with a regression-based method adjusted by local linear regression [START_REF] Blum | Non-linear regression models for Approximate Bayesian Computation[END_REF][START_REF] Csilléry | abc: an R package for approximate Bayesian computation (ABC)[END_REF]. Parameters were weighted by an Epanechnikov kernel. One million simulations were produced under the selected scenario to explore parameter space, but only simulations closest to the observed dataset were retained for parameter estimation (i.e. tolerance rate). Confidence in parameter estimates was checked with leave-one-out cross-validation repeated 1,000 times, thus estimating prediction errors.

Genotype phasing and imputation

The haplotype phase and missing data of the 468 sequenced individuals were inferred using Beagle 5.1 [START_REF] Browning | Rapid and accurate haplotype phasing and missing-data inference for whole-genome association studies by use of localized haplotype clustering[END_REF][START_REF] Browning | A one-penny imputed genome from nextgeneration reference panels[END_REF]) to perform population migration modeling and admixture tests. Beagle uses the localized haplotype-cluster model and applies an iterative approach to infer the most likely haplotype pair for each individual.

At each iteration, phased input data are used to build a localized haplotype-cluster model.

Once the model is built, phased haplotypes for each individual are sampled from the induced diploid HMM conditional on the individual's genotypes. The sampled haplotypes are the input for the next iteration, and so forth. In the final iteration, the Viterbi algorithm selects the most likely haplotypes for each individual, conditional on the diploid HMM and the individual's genotype data. For each copy of each individual, missing alleles are randomly imputed according to allele frequencies, and the data for each individual are phased by randomly ordering the genotypes [START_REF] Browning | Rapid and accurate haplotype phasing and missing-data inference for whole-genome association studies by use of localized haplotype clustering[END_REF].

Native population migration modeling

The modeling of the native population splits and mixtures was performed using TreeMix, a statistical model inferring the patterns of population splits and mixtures in multiple populations [START_REF] Pickrell | Inference of Population Splits and Mixtures from Genome-Wide Allele Frequency Data[END_REF]. A Maximum-Likelihood tree was constructed using genome-wide allele frequency and genetic drift approximation. One thousand bootstraps were performed to assess the robustness of the inferred Maximum-Likelihood (ML) tree.

Migration edges were then added sequentially to connect pairs of populations when allele frequency covariance excess was detected. For this analysis, we tested mainland China populations. To strengthen the migration model, we added the genetically close-related Tibet sample to the north-central China deme and considered Japan's population as an outgroup.

We ran TreeMix with migration events ranging from zero to four and the TreeMix composite model incorporating known admixture with the '-cor_mig' option (Pickrell et al. 2012) based on STRUCTURE admixture results. The robustness of the tree and the migration edges were confirmed by 1,000 bootstraps using GNU Parallel [START_REF] Tange | Gnu parallel-the command-line power tool[END_REF]) and 'treemix.bootstrap' function implemented in the R package BITE [START_REF] Milanesi | BITE: an R package for biodiversity analyses[END_REF]).

Native population admixture test

We selected one deme from each mainland China region (North, North East, Central, South

RESULTS

Complex genetic structure of the native range

We searched for the genetic structure of the native range in order to aggregate sampled sites into fewer homogeneous and consistent demes, hence simplifying further demographic inferences based on ABC modeling. The highest level of structure estimated by STRUCTURE was K=5-6 (Figure 1a, Table S3, Supplementary information). To select the most plausible number of genetic clusters K, multiple criteria were assessed and two clustering methods cross-validated (Table S3, Supplementary information). The modal value of the distribution of ΔK given K suggested K=3 in the native area, whereas the shape of the plateau of ln(Pr(X|K)) suggested a value of K between 5 and 7. At K=3, Northern China was well separated from Central China/Japan and Southern China. However, at K=6-7, a better resolution on Central China (S4, S6, S10, S11) and Tibet was obtained, with signals of a north-south admixture. The overall patterns of clustering and admixture were similar between K=6 and K=7 with a spurious cluster at K=7. Convergence tests assessed that most sampling chains were convergent, yet some of the sampling chains remained nonconvergent even after a burn-in of 100,000 and 100,000 sampling iterations. The highest values of H' was for K=3 (0.99), though H' reached 0.81 for K=6. Moreover, the DAPC approach yielded very similar results with a most plausible K=5. STRUCTURE and DAPC clustering were congruent in revealing the same separation between Northern and Southern

China. The K that we finally selected was that meeting statistical, geographical and historical compromises. Hence, K=6 was chosen as the better compromise between statistical parsimony and the highest level of genetic structure.

Definition of native putative demes

The definition of native putative demes based on genetic clustering was crucial to building invasion route scenarios (Supporting information). Most putative demes in the native region were consistent (South China, North China, Japan), despite high uncertainty in some sites (Figure 2a, Figure S6). S3 in particular had a small sample size and was strongly admixed, causing high uncertainty for its assignment to a putative deme. As a consequence, S3 was removed from the dataset for any further analysis. Two admixed demes, composed of sites with a major ancestry coefficient lower than 70%, were created in Central China. These groupings were also justified on the basis of previous studies reporting this region as a zone of secondary contact between Northern and Southern populations [START_REF] Hardouin | Reconciling the biogeography of an invader through recent and historic genetic patterns: the case of topmouth gudgeon Pseudorasbora parva[END_REF][START_REF] Simon | Invasive cyprinid fish in Europe originate from the single introduction of an admixed source population followed by a complex pattern of spread[END_REF]). S11 was assigned at 65% to the North East China deme, but its geographic proximity and connectivity with S10 supported the constitution of an admixed Central China deme encompassing S10 and S11. Lastly, the admixed S13 site was placed into the North China deme because of its location within the same river basin (Figure 2a).

While the Tibetan population is genetically representative of the Central China deme (ancestry coefficient > 75%), for historical reasons (i.e. recent Tibetan introduction of P.

parva and few commercial exchanges with Europe), this population was not considered part of the deme. On the other hand, S19 and S20 were genetically well clustered with S9 and S18 (ancestry coefficient >99%) and were considered part of the South China deme. Six native demes were finally defined: North China, North-East China, Central China, South-East China, South China and Japan (Figure 2a).

Two main gene flow directions shaped the native populations

In the obtained consensus tree (Figure 1b), almost all the nodes were well supported (75%-100%). The Central-North East China node was supported by 47% of the bootstraps. The model with four migration events showed stable migration edges after multiple runs of TreeMix and TreeMix composite models (Figure 1b). The migration edges showed two main S4).

Genetic structure of the invasive range

Genetic clustering in the invasive range was used to narrow the putative number of nonnative populations' origins to test (i.e. the number of scenarios). STRUCTURE results indicated a clear genetic structure in the European invasive range with K=3, supported by DAPC results, despite some uncertainty around the actual highest level of genetic clustering (Table S3). Indeed, for all values of K between 2 and 9, Turkey and Bulgaria clustered together without admixture; Iran formed a similarly distinct cluster (Figures S5a-S5b).

Conversely, Western Europe showed a pattern of admixture for all values of K and admixture persisted unless K was equal to the number of sites (Figures S5b-S5c). This indicated that Western Europe is most likely a single deme with a strong sub-structure.

Furthermore, hierarchical STRUCTURE analysis performed on the subset of West European sites confirmed a substructure at this scale. Within West European sites (Aus, Bel, Hun, Ita, Pol, Spa, UK), K=3-4 was the most probable grouping according to Evanno's method and K=5 according to the plateau method (Figure S5c). K=5 was displayed for the Western European sub-structure (Figure S5c) because spurious admixture indicated overfitting for higher values (Figure S5b). H' was high (0.98) for K=3-4-5 in the STRUCTURE analysis of the Western European deme. Under Geweke's diagnostic and trace plots, sampling chains seemed convergent or close to convergence, yet some seemed to show departure from stationarity. Heidelberg and Welsh's diagnostic supported that most of the Markov Chains were a stationary distribution. Additionally, Gelman and Rubin's diagnostic assessed that replicated chains converged on similar values. Finally, despite the strong sub-structuring observed, the Western European deme was considered as a single consistent genetic population in ABC scenarios. Hence, we retained three independent demes in the invasive range formed by Western Europe, Eastern Europe and Iran.

Population assignment to source populations

'AssignPOP' training was efficient, with an assignment accuracy greater than 90% (Figure S7) and all putative demes confidently discriminated by the training algorithm. Predictions with the SVM algorithm confidently assigned 100 individuals (59%) to a source population,

ABC inference of source populations and invasion pathways

ABC simulations were first used to infer the most probable source population in Asia of each one of the three main demes identified in the non-native area (step 1), followed by discrimination between different competing invasion pathways at a global scale (step 2).

Step

1: Source populations of non-native populations

Following prior calibration, all scenarios could be fitted to the data, with 112 summary statistics estimated to compare scenarios. Prior error rates were low for both Neural Network and Random Forest (Table 2), indicating a good predictive power. There was no confusion between scenarios, and the marginal posterior predictive distribution improved the goodness-of-fit of the selected model. Demographic inferences in Western Europe, Eastern

Europe and Iran (the three main demes identified from the clustering approaches) all predicted introductions issued from admixture events (scenario 4, Figure S4), supported by high posterior probabilities (Table 2). The choice of candidate source populations was congruent with phylogenies inferred with Maximum Likelihood in Treemix (Figure S8). Each invasive population formed a group with an Asian native population, with relatively short distances between them for European demes. The Iranian population was particular, exhibiting long drift from the most recent common ancestor shared with Japan.

Step 2: Invasion pathways As in step one, prior fitting to the data was achieved, although model selection was based on 256 estimated summary statistics. Prior error rates were also low (Table 2), giving confidence in subsequent inferences. Scenario 1, modelling three independent introductions from independent admixed populations leading to the three observed invasive demes, was selected with strong support from both Neural Network (99%) and Random Forest (86% of votes) algorithms (Table 2, Figure 2c). were an admixture between China (South East China, North East China) and Japan (Figure 2c).

Invasion process and founding populations

Parameters linked to invasion history (i.e. time of invasion and bottleneck severity of the founding population) were estimated with reasonable confidence intervals by local linear regression (tolerance rate of 0.05 for Eastern and Western Europe; 0.005 for Iran; Figure S9). The Iranian invasion was estimated to have occurred 66 generations ago (CI95% = 44; 90), with an effective population size of 217 (CI95% = 185; 253) for the founding group. The Eastern European invasion was estimated to be more recent at 39 generations ago (CI95% = 18; 62) and with a smaller effective population size 90 (CI95% = 32; 148). The effective population size of the Western European founding group required two million simulations to be estimated with a higher degree of confidence at 1,436 individuals (CI95% = 1,001; 1,825); furthermore, the time of divergence was estimated at 40 generations (CI95% = 20; 63).

DISCUSSION

Population genomics approaches used to inform biological invasions can be challenged by the discrepancy between the large timescale of genetic mutations and the smaller timescale of human-mediated invasions, leading to potentially flawed analyses [START_REF] Fitzpatrick | What can DNA tell us about biological invasions?[END_REF]).

Nevertheless, accurate assignments and estimates of geographic origin and number of introductions can be obtained through extensive sampling of native and invasive ranges, the use of high-resolution of molecular markers and a highly genetically structured source population [START_REF] Excoffier | Computer programs for population genetics data analysis: a survival guide[END_REF][START_REF] Dlugosch | Founding events in species invasions: Genetic variation, adaptive evolution, and the role of multiple introductions[END_REF]. Our extensive sampling across both native and invasive areas, the genome-wide markers used in this analysis, as well as previous knowledge on native genetic structuring [START_REF] Hardouin | Reconciling the biogeography of an invader through recent and historic genetic patterns: the case of topmouth gudgeon Pseudorasbora parva[END_REF], allowed us to highlight long-term phylogeography and recurrent gene flow. In turn, this knowledge of complex genetic structure in the native range was an essential element allowing us to successfully model probable invasion pathways. Our results suggest that the genomic diversity of invaders was shaped long before introduction by the presence of geographical barriers and by human-mediated gene flow in the native range. The combination of appropriate and complementary methods for populations departing from H-W equilibrium and appropriate clustering within demes allowed us to infer source populations and an invasion scenario despite high admixture within native demes, thus reconciling the genetic history of the native range with recent genomic patterns of invasion.

The Asian history of P. parva shaped by paleogeography and anthropization

Paleogeographic influence

Large-scale genetic structure was consistent with the known phylogeography of the species [START_REF] Hardouin | Reconciling the biogeography of an invader through recent and historic genetic patterns: the case of topmouth gudgeon Pseudorasbora parva[END_REF][START_REF] Simon | Invasive cyprinid fish in Europe originate from the single introduction of an admixed source population followed by a complex pattern of spread[END_REF]. Likewise, the retained fine-scale structure in the native range was consistent with non-genetic data (historical, geographical and morphological), suggesting that the estimated value of K is representative of the actual genetic structure. Asian phylogeography revealed a complex genetic structure involving locally high gene flow and variable admixture. Results of STRUCTURE, TreeMix and Dstatistics taken together suggest that the current genetic structuring of the native P. parva range has been shaped by waves of gene flow originating from southern and northern populations (edges of the natural native distribution). Fish phylogeography in China has been greatly influenced by geological events [START_REF] Chiang | Diverse processes shape deep phylogeographical divergence in Cobitis sinensis (Teleostei: Cobitidae) in East Asia[END_REF][START_REF] Li | Studies on zoogeographical divisions for fresh water fishes of China[END_REF]) that can represent major barriers to gene flow [START_REF] Brandley | Does dispersal across an aquatic geographic barrier obscure phylogeographic structure in the diamond-backed watersnake (Nerodia rhombifer)?[END_REF]. The complex geological history of the South China landmass acted as a barrier to gene flows and induced vicariance in common cyprinid species [START_REF] Yang | The population history of Garra orientalis (Teleostei: Cyprinidae) using mitochondrial DNA and microsatellite data with approximate Bayesian computation[END_REF]. In China, notable genetic divergence occurred between the North and the South, with the Qinling Mountains acting as a strong biogeographic barrier between the temperate climate of the North and the subtropical climate of the South [START_REF] Yuan | Genetic structure of the tree peony (Paeonia rockii) and the Qinling Mountains as a geographic barrier driving the fragmentation of a large population[END_REF][START_REF] Dong | Tectonic evolution of the Qinling orogen, China: review and synthesis[END_REF] resulting in a wide range of local adaptations to various ecological conditions that might facilitate establishment in the invaded area.

Yet within both biogeographical regions, gene flow was restricted and populations were mainly differentiated per river basins. In China north to the Qinling Mountains, the two main demes (i.e. North China and North East China; Figure 2a) were influenced by two different major river networks, the Amur River and the Yellow River (S1, S2 and S3 are located on the Haihe and Huaihe Rivers but both rivers were influenced by the Yellow River in history).

South China (sites S9, S18; Figure 2a) was structured by the Yangtze River. Indeed, most native demes were congruent with the architecture of the main river networks. Historically, five Pseudorasbora species, based on morphological characteristics and their association with major river basins, were described as endemic to China [START_REF] Nichols | Chinese fresh-water fishes in the American Museum of Natural History's collections: a provisional check-list of the fresh-water fishes of China[END_REF]. These five species actually corresponded to the same P. parva species (Gozlan 2012), though genetic clustering was largely in agreement with the distribution of historical morphotypes [START_REF] Nichols | Chinese fresh-water fishes in the American Museum of Natural History's collections: a provisional check-list of the fresh-water fishes of China[END_REF]). Yet admixture was also observed in locations at the boundaries between river basins (e.g. site S16 at the frontier between the Yellow and Amur river basins), indicating that recurrent migration has happened between demes.

Wild populations in Japan split from Chinese continental populations 12.1 Mya [START_REF] Hardouin | Reconciling the biogeography of an invader through recent and historic genetic patterns: the case of topmouth gudgeon Pseudorasbora parva[END_REF], consistent with Japan's separation from the Eurasian continent 15 Mya [START_REF] Barnes | Origins of the Japanese Islands: The New" Big Picture[END_REF]). Yet our results demonstrated significant gene flow with the South East deme. Two mtDNA lineages co-exist in Japan, including one closely related to continental populations [START_REF] Watanabe | Phylogenetic relationships of the Japanese minnows, Pseudorasbora (Cyprinidae), as inferred from mitochondrial 16S rRNA gene sequences[END_REF], and recent hybridizations with Chinese populations have been described [START_REF] Hardouin | Reconciling the biogeography of an invader through recent and historic genetic patterns: the case of topmouth gudgeon Pseudorasbora parva[END_REF]). Across Honshū Island, P. parva rapidly expanded its distribution with translocations of commercial cyprinids into ponds, which may have facilitated the introduction and the spread of the Chinese lineage [START_REF] Konishi | Natural hybridization between endangered and introduced species of Pseudorasbora, with their genetic relationships and characteristics inferred from allozyme analyses[END_REF][START_REF] Konishi | Identifying conservation priority ponds of an endangered minnow, Pseudorasbora pumila, in the area invaded by Pseudorasbora parva[END_REF].

However, resolving the origins and contemporary phylogeography of the Japanese population within the native distribution of P. parva with sufficient power would require more extensive sampling in Japan.

Anthropogenic influence

Genetic structure with long-term admixture has been reported in South China for other fishes, explained by coastal land and tributaries between river basins, especially for the Yangtze and Pearl Rivers [START_REF] Yang | The population history of Garra orientalis (Teleostei: Cyprinidae) using mitochondrial DNA and microsatellite data with approximate Bayesian computation[END_REF]. Admixed populations may be the result of anthropogenic modifications to the hydrological landscape, especially new dispersal pathways such as canals. Recent secondary contacts between the North and South may have been facilitated by the increased structural connectivity between river basins. The construction of the Lingqu Canal 2,200 years ago connected the Yangtze River (sites S9, S18 in South China and to a lesser extent S6 and S4 in South East China; Figure 2a) to the Pearl River [START_REF] Fengshu | Preliminary study of geological geomorphological conditions of ling canal[END_REF]. Moreover, the Beijing-Hangzhou Grand Canal connected the Yangtze River to the Huai He, Hai He and Yellow rivers in the east of China (AD 581-618, Sui Dynasty), a region at the core of the two admixed regions (Central China, North East China and South East China: Figure 2a). Secondary contact due to human-induced corridors can increase local admixture between populations [START_REF] Crispo | Broken barriers: Humaninduced changes to gene flow and introgression in animals: An examination of the ways in which humans increase genetic exchange among populations and species and the consequences for biodiversity[END_REF], and canals serving as ecological corridors between two interbreeding cyprinid species have ultimately become hybrid zones [START_REF] Guivier | Canals as ecological corridors and hybridization zones for two cyprinid species[END_REF]). However, the most important vector for long-distance dispersal in the native range is aquaculture expansion throughout central China. The two most admixed demes correspond to areas of intensive aquaculture in China. In Central China, many ponds and reservoirs were intensively stocked in the 1950s-60s, with P. parva eventually becoming the dominant species [START_REF] Zhao | Current state of freshwater fisheries in China[END_REF][START_REF] Gong | Fishery in contemporary China[END_REF]. Massive human translocations for aquaculture have often led to biotic homogenization [START_REF] Olden | Ecological and evolutionary consequences of biotic homogenization[END_REF]) and increased pathogen dispersal into the pool of future invaders [START_REF] Price | Reconstructing the emergence of a lethal infectious disease of wildlife supports a key role for spread through translocations by humans[END_REF].

Moreover, consecutive introductions represent a type of punctuated gene flow that produces admixed genotypes with high genetic diversity but unpredictable evolutionary effects [START_REF] Crispo | Broken barriers: Humaninduced changes to gene flow and introgression in animals: An examination of the ways in which humans increase genetic exchange among populations and species and the consequences for biodiversity[END_REF][START_REF] Hasselman | Human disturbance causes the formation of a hybrid swarm between two naturally sympatric fish species[END_REF]. In native locations, admixture is known for negative effects such as outbreeding depression and loss of local adaptation [START_REF] Côte | Outbreeding depression in Atlantic salmon revealed by hypoxic stress during embryonic development[END_REF][START_REF] Huff | Mixed-source reintroductions lead to outbreeding depression in second-generation descendents of a native North American fish[END_REF][START_REF] Hufford | Inbreeding and outbreeding depression in Stylidium hispidum: implications for mixing seed sources for ecological restoration[END_REF]), but in novel environments, admixture may increase adaptive potential for translocated individuals [START_REF] Verhoeven | Population admixture, biological invasions and the balance between local adaptation and inbreeding depression[END_REF]. For example, admixture in the wild between divergent sculpin populations increased their genetic diversity, lineage differentiation, and facilitated colonization of new habitats [START_REF] Nolte | An invasive lineage of sculpins, Cottus sp.(Pisces, Teleostei) in the Rhine with new habitat adaptations has originated from hybridization between old phylogeographic groups[END_REF][START_REF] Nolte | Variable patterns of introgression in two sculpin hybrid zones suggest that genomic isolation differs among populations[END_REF][START_REF] Stemshorn | Rapid formation of distinct hybrid lineages after secondary contact of two fish species (Cottus sp.)[END_REF].

The invasion

Invasion starts in Asia

Invasive populations have been successfully established at the edges of the native distribution. In South-West China, sites S19 and S20 are recently introduced populations (1980s) from South-East China, and perfectly clustered with S9, S18 (Figure 2a; see also [START_REF] Hardouin | Reconciling the biogeography of an invader through recent and historic genetic patterns: the case of topmouth gudgeon Pseudorasbora parva[END_REF]. The Tibetan population was also introduced, commercial exchange being most probably the pathway for this invasion given that many living fishes are sold in Tibetan markets and are traditionally released into rivers [START_REF] Gozlan | Pseudorasbora parva temminck and schlegel (topmouth gudgeon). Handb Glob Freshw Invasive Species[END_REF]. Although Tibet shows similarities with Central China, a continuum of populations between them were not sampled, making it difficult to determine if Tibet should be clustered with this deme or another.

Out of Asia

Three invasive demes were consistently discriminated as (1) Iran, (2) a small Balkan-Anatolian deme (Eastern Europe) and (3) a large pan-European deme (Western Europe).

Population differentiation (FST) was overall higher between the three invasive demes, than between the invasive demes and their source populations, suggesting independent introductions with different genomic backgrounds (Figure 3, Table S5). In addition, clear genetic sub-structuring was observed in the Western European deme. More populations sampled in this region and more markers would then be essential to more accurately depict fine-scale genetic structure of the West European deme. The Iranian population was particular, forming a distinct deme completely isolated from all European populations, consistent with observations of [START_REF] Hardouin | Reconciling the biogeography of an invader through recent and historic genetic patterns: the case of topmouth gudgeon Pseudorasbora parva[END_REF], who revealed that Iran has a different mtDNA lineage. Differentiation between Eastern and Western Europe may be imputed to the former USSR era when on one side, the eastern introduction first opens and on the other commercial exchanges with Western Europe were limited [START_REF] Britton | Geo-politics and freshwater fish introductions: How the Cold War shaped Europe's fish allodiversity[END_REF]. The Danube river ensures structural connectivity which has favored colonization among a large central European region [START_REF] Weiss | Complex patterns of colonization and refugia revealed for European grayling Thymallus thymallus, based on complete sequencing of the mitochondrial DNA control region[END_REF], from Germany to Bulgaria through Austria and Hungary. However, regional biogeography of the Balkan region could also partially explain the observed population structure [START_REF] Oikonomou | Biogeography of freshwater fishes of the Balkan Peninsula[END_REF]). In addition, rapid range expansion of P. parva over long distances in the invasive range likely relied on humanmediated dispersal (e.g. unintentional releases, aquaculture exchanges) rather than natural dispersal, as P. parva favours lentic habitats. The observed genetic structure of invasive locations may suggest multiple independent introductions without subsequent gene flow.

This assumption is strengthened by an invasive population genetic diversity similar to the native populations (Table S1). A recent study showed a clear reduction in genetic diversity in the French invasive deme of P. parva compared to native demes, highlighting that the introduction was accompanied by a substantial loss of genetic diversity [START_REF] Combe | Towards unravelling the Rosette agent enigma: Spread and emergence of the co-invasive hostpathogen complex, Pseudorasbora parva-Sphaerothecum destruens[END_REF].

However, in a previous study, high genetic diversity among invasive P. parva has been reported, first by [START_REF] Simon | Human induced stepping-stone colonisation of an admixed founder population: the spread of topmouth gudgeon (Pseudorasbora parva) in Europe[END_REF], and later attributed by [START_REF] Hardouin | Reconciling the biogeography of an invader through recent and historic genetic patterns: the case of topmouth gudgeon Pseudorasbora parva[END_REF] to unbalanced sampling. Baltazar-Soares et al. ( 2020) have since highlighted a high genomic diversity attributed to a lack of admixture due to a recovery of the genetic bottleneck associated with the introduction. Here, with more balanced sampling between invasive and native populations, coupled with observed admixture in the invasive population (Figure 2b, Figure S5), our results suggest an invasive genetic diversity shaped by past invasions of admixed native populations.

It is common to observe strong genetic structure within demes for continental invasive fishes when using nuclear markers. This is partly explained by human driven long-distance dispersal and geographical barriers to gene flow [START_REF] Sanz | Genetic characterization of the invasive mosquitofish (Gambusia spp.) introduced to Europe: population structure and colonization routes[END_REF][START_REF] Simon | Invasive cyprinid fish in Europe originate from the single introduction of an admixed source population followed by a complex pattern of spread[END_REF]).

Although we detected signals of differentiation with nuclear SNP data, mtDNA has revealed that Western European locations were characterized by the same mitochondrial lineage and, thus, constitute a single genetic population sharing a common ancestor [START_REF] Simon | Invasive cyprinid fish in Europe originate from the single introduction of an admixed source population followed by a complex pattern of spread[END_REF]).

Low migration rates in a stepping-stone model induced by human-mediated secondary introductions might explain the relative differentiation between Western European locations, and might be a sign of an invasive bridgehead [START_REF] Van Boheemen | Multiple introductions, admixture and bridgehead invasion characterize the introduction history of Ambrosia artemisiifolia in Europe and Australia[END_REF]. Invasive populations could also act as source populations for a secondary invasion, even over longdistances [START_REF] Karsten | Deconstructing intercontinental invasion pathway hypotheses of the Mediterranean fruit fly (Ceratitis capitata) using a Bayesian inference approach: are port interceptions and quarantine protocols successfully preventing new invasions?[END_REF][START_REF] Lombaert | Complementarity of statistical treatments to reconstruct worldwide routes of invasion: the case of the Asian ladybird Harmonia axyridis[END_REF].

Origins of invasion: admixture and multiple introductions

Our models described major demographic events, and gave congruent results between the two independent Machine Learning algorithms, supporting three independent introductions at the origins of the three invasive demes. Multiple introductions would, thus, multiply the chances of successful establishment outside the native range. Results point towards a southern Chinese origin of the Eastern European deme, with alleles coming from and/or shared between the Central, South East and South China demes (Figure 2c). Similar patterns have been demonstrated via population assignment tests, revealing multiple origins in the same location indicative of multiple introductions or admixture within invasive demes [START_REF] Guillemaud | Biological invasions in agricultural settings: insights from evolutionary biology and population genetics[END_REF]. Studies have also highlighted that admixture is likely a common characteristic for invasive species [START_REF] Genton | High genetic diversity in French invasive populations of common ragweed, Ambrosia artemisiifolia, as a result of multiple sources of introduction[END_REF][START_REF] Kolbe | Genetic variation increases during biological invasion by a Cuban lizard[END_REF][START_REF] Rius | How important is intraspecific genetic admixture to the success of colonising populations?[END_REF]; however, it is often difficult to disentangle multiple introductions from admixture before introduction. The reconstruction of admixture history may depend on the genetic structure observed in the native area. Admixture after introduction may be inferred when source populations are clearly differentiated and admixture is demonstrated in the invasive area, as in the case of the mussel, Mytella charruana [START_REF] Gillis | Higher genetic diversity in introduced than in native populations of the mussel Mytella charruana: evidence of population admixture at introduction sites[END_REF]). However, this is not the case for P. parva, which shows high rates of admixture within the inferred source populations, suggesting that human-mediated gene flow within the native range shaped genetic variation in the pool of future invaders. The Central deme corresponded to the Wuhan region, at the core of Chinese aquaculture, with more than two million fishes produced annually, and frequent translocations from surrounding regions to restock ponds [START_REF] Zhao | Current state of freshwater fisheries in China[END_REF] [START_REF] Britton | Geo-politics and freshwater fish introductions: How the Cold War shaped Europe's fish allodiversity[END_REF]. On average, genetic divergence amongst source demes was relatively high (mean FST = 0.27; Figure 3, Table S5), suggesting limited gene flow and a lower likelihood of native admixture, compared to that for Eastern Europe source demes (mean FST = 0.13; Figure 3, Table S5). Moreover, individuals of the Hungarian location were systematically assigned to two genetic clusters (Figure S5c), making it the most admixed location of Western Europe, and suggested further that admixture most likely happened from multiple introductions in the same invasive area. Temporal dynamics of invasive mosquitoes in the US have shown the same admixture-like clustering pattern [START_REF] Fonseca | Fine-scale spatial and temporal population genetics of Aedes japonicus, a new US mosquito, reveal multiple introductions[END_REF]. Range expansion across European countries most probably began in the Hungarian admixed region, and its likely chronology has been reconstructed from census data [START_REF] Gozlan | Pancontinental invasion of Pseudorasbora parva: towards a better understanding of freshwater fish invasions[END_REF]. P. parva colonized Austria as early as 1982 [START_REF] Weber | Die ausbreitung der pseudokeilfleckbarben im donauraum[END_REF], then Germany from 1984 to 1987 [START_REF] Arnold | Pseudorasbora parva (Schlegel 1842) nun auch in der DDR[END_REF][START_REF] Lelek | Zustandsanalyse der Fischartengemeinschaften im Rhein[END_REF], and Belgium by 1992 [START_REF] Vandelannoote | Atlas van de Vlaamse beek-en riviervissen[END_REF]. Austrian, Belgian and Polish sites clustered together, suggesting a common ancestry. Their connectivity along the Danube and Rhine river systems may have facilitated natural dispersal [START_REF] Hegediš | Amur sleeper (Perccottus glenii Dubowski, 1877) spreading in the Danube River basin[END_REF][START_REF] Leuven | The river Rhine: a global highway for dispersal of aquatic invasive species[END_REF].

However, natural dispersal cannot explain the observed high genetic structure, nor the expansion over geographical barriers (e.g. the English Channel, the Alps), and so humanmediated dispersal is more likely [START_REF] Gozlan | Occurrence of the Asiatic cyprinid Pseudorasbora parva in England[END_REF][START_REF] Aparicio | Expansion of the invasive Pseudorasbora parva (Cyprinidae) in the Iberian Peninsula: first record in the Guadiana River basin[END_REF][START_REF] Caiola | First record of the Asiatic cyprinid Pseudorasbora parva in the Iberian Peninsula[END_REF]. Thus, Hungary could be the invasive bridgehead that initiated the colonization of Western Europe, followed by secondary bottlenecks associated with a series of founding events during westward range expansion. This invasive bridgehead effect scenario has also been suggested by [START_REF] Combe | Towards unravelling the Rosette agent enigma: Spread and emergence of the co-invasive hostpathogen complex, Pseudorasbora parva-Sphaerothecum destruens[END_REF] for a French deme, which probably originated from a successful invasive population. Surprisingly, Italy showed different origins than other Western European locations in population assignment tests (Central China instead of North-East China; Figure 2b), but was nonetheless considered in the Western European deme based on genetic clustering results. The origin of this small population was difficult to explain, as it may be part of an unsampled fourth invasive population linked with Slovakia that successfully spread across the Italian peninsula [START_REF] Carosi | Status of Pseudorasbora parva in the Tiber river basin (Umbria, central Italy) 20 years after its introduction[END_REF][START_REF] Záhorská | Reproductive parameters of invasive topmouth gudgeon Pseudorasbora parva (Temminck and Schlegel, 1846) from Slovakia[END_REF]. This uncertainty in the Italian origin underlies the necessity for the same extensive sampling in Europe as in Asia, especially in the Southeastern Europe that was believed to be the core of past introductions (e.g. Slovenia, Hungary, Romania, Bulgaria). This punctual sampling represents an inherent limitation to the definition of putative demes, and subsequently, to the number of invasion pathways inferred.

The Iranian population origin was complex to infer because of an ancient admixed origin between South East China and Japan. Iran showed a pattern of assignment similar to that of the Japanese admixture in nuclear markers between endemic and Chinese populations, potentially explaining why it was not possible to disentangle China and Japan as putative sources. In line with our results, a recent study [START_REF] Ganjali | West Asian colonisation of topmouth gudgeon, Pseudorasbora parva (Teleostei: Gobionidae): Genetic admixture at the crossroad of Europe and east Asia[END_REF], based on two mtDNA markers, identified three Iranian matrilineal haplotypes belonging to two distinct lineages, an older Japanese lineage and a Chinese lineage, corresponding to a recent natural dispersal from Azerbaijan. They also highlighted an admixture of highly divergent Japanese and Chinese lineages [START_REF] Hardouin | Reconciling the biogeography of an invader through recent and historic genetic patterns: the case of topmouth gudgeon Pseudorasbora parva[END_REF]). Many questions remain, mainly because the history of Japanese populations has not been fully resolved. The high number of private alleles in Iran (Table S1) and the long drift separating the Iranian population from the most recent common ancestor with Japan (Figure S8) might suggest that the true source population has not been sampled. Thus, additional sampling would be required to effectively determine the origins of Iranian P. parva.

A successful invasion

The surprisingly large effective population sizes in Western Europe likely enhanced the probability of introduced populations survival, and may explain its invasiveness.

Furthermore, the co-introduction of non-native host-pathogen systems such as P. parva-S.

destruens may create new host-pathogen interactions that can be detrimental to naïve native host species, hence providing a competitive advantage for the invaders [START_REF] Price | Parasite mediation in ecological interactions[END_REF]

CONCLUSION

Our results shed light on the importance of grasping the genetic history of the native range populations to better understand the effects of introductions and admixture on the success and adaptive potential of invasive populations. Our study also draws a picture of the complex demo-genetic history of an invasive species' source populations and its spread across recurrent global invasion pathways. Reconstruction of such invasion pathways is crucial for setting up conservation biology approaches and management to prevent further non-native species introductions and potential associated pathogens.

List of Tables and Figures

  East and South) to perform admixture tests. We used the four-population test of the Dstatistic[START_REF] Green | A draft sequence of the Neandertal genome[END_REF][START_REF] Patterson | Ancient admixture in human history[END_REF]) implemented in Popstats[START_REF] Skoglund | Genetic evidence for two founding populations of the Americas[END_REF] to test for admixture and gene flow directionality within the native populations. The notation used by Popstats for the D-statistic is D(O, P3; P1, P2), where O is the outgroup, P3 the test population and P1, P2 the sister populations. A significant negative D indicates that P3 exchanged genes with P1; conversely, a positive D indicates that P3 exchanged genes with P2[START_REF] Durand | Testing for ancient admixture between closely related populations[END_REF]). The D-statistic was estimated for each combination of demes. D-statistic significance was assessed by block jackknife of 5kb and the standard error (SE) was used to estimate the Z-score[START_REF] Skoglund | Genetic evidence for two founding populations of the Americas[END_REF].

  southern and northern origins of gene flow: gene flow from the South mainly to the Central and South East populations with 45% and 33% migration weights and gene flow from the North to the North East China population with 9% migration weight (Figure 1b). Admixture tests confirmed gene flow origins followed a directional pattern. The fourpopulation test assumed the population configuration D(O, P3; P1, P2), with Northern and Southern demes as our test populations (P3), the remaining demes as P1 and P2 and Japan as the outgroup O. D-statistics showed an excess of shared derived polymorphism, highlighting genes exchanged between i) the South population and the South East and Central ones and ii) between the North population and the Central and North East ones (Figure 1c, Table

  with a relative posterior probability >2. Assignment to source populations showed multiple origins in sampled sites, especially in Eastern Europe and Iran (Figure2b), congruent with DAPC. Posterior membership probabilities of invasive individuals assigned to native Asian clusters with DAPC showed a putative origin in North East China for the Western Europe deme (cluster) and South/South East China with admixture for Bulgaria and Turkey. The Italian site showed a different origin than the rest of Western Europe, closer to Eastern Europe. The Iranian population was separated from all others, linked to South East China or Japan.

  The three different invasive populations clearly formed three distinct groups with admixed origins. Western Europe origins were in a Northern part of China (North-East China admixed with Central China and North China). Removal of the ambiguous Italian sample did not change the selected scenario. Eastern Europe origins were in a Southern part of China (admixture between Central China, South East China and South China). Lastly, Iran origins

Fig. 1 Fig. 2

 12 Fig. 1 Genetic structure and population modeling in the native range. (a) Posterior admixture proportions of individuals estimated with STRUCTURE for 18 sampled sites (n=300) in the native range (Asia) and K=6 genetic clusters (proportions aggregated with CLUMPP from 20 independent replicates). (b) ML tree based on TreeMix with block size of 500 SNPs. The migration edges pointing from South China to Central China, South Central China, and North East China have migration weights of 45%, 33%, and 2,7%. The migration edge pointing from North China to North East China has a migration weight of 9%. (c) D-statistics for testing admixture in native populations. Populations are noted D(O, P3, P1, P2). Negative values indicate gene flow between P3 and P1 while positive values indicate gene flow between P3 and P2. Confidence interval at 95% of the D-statistic estimated by block jackknife of 5kb. Null expectation is a D-statistic of 0. Demes abbreviations are North=North China, South=South China, Central=Central China, NorthE=North East China, SouthE=South East China Fig. 2 Definition of demes in the native range, and source populations predicted for the invasive range. (a) (Native range, Asia) Mean admixture proportions estimated with STRUCTURE in the native range (Asia), for the chosen K=6. Native sampled sites are pooled into putative demes for subsequent analyses. Pie chart colors correspond to the genetic clusters inferred by STRUCTURE, as in Figure 1. (Invasive range, Europe & Middle-East) Assignment predictions of invasive individuals to native demes with AssignPOP's SVM algorithm with a relative posterior probability >2 (n=292 for training, n=100 for predictions). Pie chart colors correspond to the putative demes defined in the native range to which invasive individuals were assigned. (b) Posterior assignment probabilities to putative native demes estimated with AssignPop. (c) The demo-genetic scenario that was inferred with Approximate Bayesian Computation demonstrating three independent introductions from three independent admixed source populations. Branch lengths are not scaled. Bottleneck events are represented in thin red lines in branches. Colored branches correspond to invasive demes history.

Fig. 3

 3 Fig. 3 Population differentiation. Pairwise FST (Weir & Cockerham's FST ) between native and invasive demes. FST values were estimated from 3,000 loci.

Fig. 1 Fig. 2

 12 Fig. 1 Genetic structure and population modeling in the native range. (a) Posterior admixture proportions of individuals estimated with STRUCTURE for 18 sampled sites (n=300) in the native range (Asia) and K=6 genetic clusters (proportions aggregated with CLUMPP from 20 independent replicates). (b) ML tree based on TreeMix with block size of 500 SNPs. The migration edges pointing from South China to Central China, South Central China, and North East China have migration weights of 45%, 33%, and 2,7%. The migration edge pointing from North China to North East China has a migration weight of 9%. (c) D-statistics for testing admixture in native populations. Populations are noted D(O, P3, P1, P2). Negative values indicate gene flow between P3 and P1 while positive values indicate gene flow between P3 and P2. Confidence interval at 95% of the Dstatistic estimated by block jackknife of 5kb. Null expectation is a D-statistic of 0. Demes abbreviations are North=North China, South=South China, Central=Central China, NorthE=North East China, SouthE=South East China

Fig. 3

 3 Fig. 3 Population differentiation. Pairwise FST (Weir & Cockerham's FST ) between native and invasive demes. FST values were estimated from 3,000 loci.

  . Consequently, admixture in the native range before the translocation into Europe may be at the origin of the Eastern European deme. The original introduction probably took place in Romania, at the epicenter of carp aquaculture in the 1960s, in the Nucet Fisheries Research Centre in 1961[START_REF] Gozlan | Pancontinental invasion of Pseudorasbora parva: towards a better understanding of freshwater fish invasions[END_REF], although several other introductions took place at the same time throughout the Black Sea region. Introduction was most likely accidental, as P. parva is a common contaminant of carp stockings[START_REF] Wolter | Distribution history of non-native freshwater fish species in Germany: how invasive are they?[END_REF]. Assuming between one and two generations per year among eastern European populations[START_REF] Gozlan | Introduction of non-native freshwater fish: is it all bad?[END_REF], it complies with our estimates of the time of invasion (ca. 18 to 62 generations), and the time of Asian fish export increases. Likewise, development of aquaculture in this part of Europe may have subsequently facilitated the spread of P. parva in Turkey as early as 1982(Ekmekçi & Kırankaya 2006;[START_REF] Britton | Geo-politics and freshwater fish introductions: How the Cold War shaped Europe's fish allodiversity[END_REF].The Western Europe deme was probably colonized from a first introduction in Hungary, noticed in the Paks Fisheries Farm in 1963[START_REF] Gozlan | Pancontinental invasion of Pseudorasbora parva: towards a better understanding of freshwater fish invasions[END_REF], which then spread throughout Western Europe. The source populations mainly differed from those of the Eastern Europe deme, although both share a connection with Central China. Western Europe populations also include an admixed north Chinese origin to both the North East China and North China demes. The invasion time was approximately the same as for Eastern Europe, coincident with a trading cooperation period between countries of the Eastern Bloc and China

;

  [START_REF] Andreou | Associated disease risk from the introduced generalist pathogen Sphaerothecum destruens: management and policy implications[END_REF][START_REF] Combe | The rise of the rosette agent in Europe: An epidemiological enigma[END_REF] Vilcinskas & Knoll 2015). The gene flow resulting from the multiple introductions highlighted in Western Europe largely increased the effective population size as well as genetic diversity. However, effective population sizes were smaller in Iran and Eastern Europe, although the bottleneck they experienced was not severe, with 100 to 200 individuals in the founder populations. P. parva probably overcame the loss of diversity induced by bottlenecks due to admixture between native populations within its new environment and prior to range expansion. Our results demonstrated multiple source populations coming from a wide range of climatic conditions, the most flagrant being the admixture in the Wuhan region between subtropical and temperate populations, which represent different morphotypes[START_REF] Nichols | Chinese fresh-water fishes in the American Museum of Natural History's collections: a provisional check-list of the fresh-water fishes of China[END_REF][START_REF] Gozlan | Native drivers of fish life history traits are lost during the invasion process[END_REF]).

	Multiple introductions
	associated with admixture can either increase or decrease the adaptive potential of invasive
	populations (Barker et al. 2019; Verhoeven et al. 2011). However, if from genetically distinct
	sources, multiple introductions can mitigate the negative effects of bottlenecks associated
	with invasion, as genetically diverse populations are less affected by the deleterious effects
	of inbreeding depression (Verhoeven et al. 2011). Moreover, they can increase individual
	fitness through heterosis, contributing to the invasion's success (Rius & Darling 2014;
	Vallejo-Marín et al. 2021), which is most likely the case in P. parva. Hence, native
	populations with broad and diverse biogeographic distribution may act as genetic and
	ecological diversity reservoirs.

Table 1 .

 1 Population properties with site specific sample, geographic range, geographical coordinates, and (N) number of individuals.

Table 2 .

 2 Inference of the demographic model based on Neural Network and Random Forest. Selected models at each step are given for each algorithm (bold values are those of the selected model). 10,000 simulated datasets per scenario for training. Prior error rate estimated with leave-one-out cross-validation for Neural network and out-of-bag error rate for Random Forest. Random Forest does not provide posterior probability of each scenario, but rather the number of votes in favor to each one. Subsequently, posterior probability of the Random Forest corresponds to the posterior probability of the selected model.

Table 1 .

 1 Population properties with site specific sample, geographic range, geographical coordinates, and (N) number of individuals.

Table 2 .

 2 Inference of the demographic model based on Neural Network and Random Forest. Selected models at each step are given for each algorithm (bold values are those of the selected model). 10,000 simulated datasets per scenario for training. Prior error rate estimated with leave-one-out cross-validation for Neural network and out-of-bag error rate for Random Forest. Random Forest does not provide posterior probability of each scenario, but rather the number of votes in favor to each one. Subsequently, posterior probability of the Random Forest corresponds to the posterior probability of the selected model.

Acknowledgements

We would like to thank all co-authors of Gozlan et al. 2020 for their help in collecting P. parva samples across the native and invasive range as well as B. Haenfling. We would also like to thank Eric Lombaert for kindly agreeing to discuss the best practices of Approximate Bayesian Computation methods. This study has been conducted in accordance with the relevant animal or human ethics approvals.

This study is part of the project GENESIS ANR-AF 13-ADAP-0005-001 supported by the Agence Nationale de la Recherche. S.B. was co-funded by the project PROBIS (Biodiversa) and E.C. was funded by an IRD postdoctoral fellowship. R.E.G., A.G., J.F.M. and S.B. received funding for this research from ANR Bioadapt.

data and material

All sequencing data (GBS reads) are publicly-available under the BioProject ID PRJNA799666. The original genetic data yielded by genotyping (in vcf format), additional metadata (e.g. geographical coordinates) and all scripts necessary to reproduce the results presented in the paper are available in a publicly-available Github repository: https://github.com/ThomasBrazier/PseudorasboraGENESIS. Data are available from the authors upon reasonable request.