Microsatellites and mitochondrial evidence of multiple introductions of the invasive raccoon Procyon lotor in France

Raccoons (Procyon lotor) are worldwide invaders, due to deliberate or accidental releases, and their impacts exceed hundred of billions of Euros in Europe only. In France, raccoons have currently established three separate, expanding populations. Identifying the current spatial genetic structure, dispersal events and phylogeography of these populations is needed to infer the invasion history and identify management units. We used wild and captive individuals sampled in France and Belgium to characterize the genetic diversity and current population genetic structure of French raccoon populations and identify potential genetic connectivity with the Belgium population using both mitochondrial DNA and microsatellite loci. Results confirm that French populations are the result of at least three independent introductions. While the three populations display low genetic diversity and sign of recent bottleneck, they are still expanding, suggesting that in addition to their ecological plasticity, the remaining genetic diversity is sufficient to successfully adapt to their new environment and allow a quick colonization. Particular attention must be given to the North-Eastern population, which shows genetic admixture with the Belgium population, as admixed individuals may exhibit hybrid vigor facilitating their expansion. The comparison of captive and wild individuals did not allow to identify a potential captive origin of the wild populations. The current regulation in France allowing captivity in zoos without enforcement to tighten the biosecurity of detention facilities might dampen any management measure as few introduced founders might be enough to create new populations.


Introduction
Invasive species are considered among the most significant threats to global biodiversity (Bellard et al. 2016a;Pyšek et al. 2020), with invasive mammalian predators being the most damaging group of alien animal species (Bellard et al. 2016b). Through competition (Smith and Banks 2014), predation (Doherty et al. 2016), disease transmission (Chinchio et al. 2020), hybridization (Parker et al. 1999) and facilitation of other invasive species , they can cause the decline and extinction of native species with cascading effects throughout entire ecosystems (White et al. 2006). The economic consequences of invasive species are just as severe, with costs associated with mitigating their impacts exceeding €116 billion between 1960 and 2020 in Europe only (Haubrock et al. 2021).
The raccoon (Procyon lotor) is a medium-sized invasive carnivore species. Originating from Central and North America, it has now colonized many parts of the world (e.g., The Antilles, Helgen et al. 2008; Japan, Ikeda et al. 2004) after deliberate or accidental releases, including at least 27 European countries (Louppe et al. 2019;Salgado 2018). Their success is attributed to their high ecological plasticity in habitat selection and feeding habits (Louppe et al. 2021;Rulison et al. 2012), high reproductive potential (Asano et al. 2003), and the lack of natural predators (Salgado 2018). Despite a lack of comprehensive quantification of the damage by invasive raccoon populations, competition with native species (Kauhala 1996), predation on native species (e.g., Cichocki et al. 2021;Fiderer et al. 2019;Tricarico et al. 2021), and potential disease transmission (e.g., Beltrán-Beck et al. 2012;Frantz et al. 2021;Maas et al. 2022) have been suggested or evidenced to occur throughout Europe. As a consequence, the species is listed among the invasive species of European Union (EU) concern for which concerted action at EU level is required (EU Regulation 1143/2014).
The raccoon is considered as a non-native pest species in France, and is thus opportunistically trapped and hunted apart from any coherent national control strategy. However, to develop an effective management strategy such as a coordinated national plan, it is mandatory to gather basic ecological, genetic and phylogeographical information about the target species (Hohenlohe et al. 2021). Particularly, identifying the current spatial genetic structure and inferring dispersal pattern among populations, together with investigating the phylogeography, the historical processes responsible for the past to present geographic distributions of lineages, is needed to determine the invasion history (Le Roux 2022). It will also help to identify management units that, with appropriate care, could be eradicated with little chance of recolonization (Robertson and Gemmell 2004). Determining the level of genetic variability is also important to predict the future of established or newly-colonising populations (Abdelkrim et al. 2005;Miura 2007) as a high level of genetic variation is assumed to increase the adaptative potential of introduced species and their invasion success (Crawford and Whitney 2010;Le Roux 2022).
In France, a field investigation by Léger and Ruette (2014) showed that raccoons have currently established three separately expanding populations: (i) the oldest (1966) and biggest population in North-Eastern France (hereafter "North-Eastern France population"); (ii) the second population (1990) in the Massif Central mountains (hereafter "Massif Central population"); (iii) the third most recent (2007)  the South-West of France in Gironde region (hereafter "Gironde population"). These populations are now growing rapidly and are assumed to be the results of independent introductions (Léger and Ruette 2014). The first one may be the result of pet escapes from a former NATO military air settlement, and the last two are supposed to be the results of individuals escaped from zoos or private owners (Léger and Ruette 2014). Occasional reports are also made regularly throughout France (hereafter "Non-residential population"), probably corresponding to escapes from captivity that have not, a priori, founded permanent wild populations (Léger and Ruette 2014).
In this paper, we used wild and captive individuals sampled in France and Belgium (hereafter "Belgium population") to characterize the genetic diversity and the current population genetic structure of all French raccoon populations and to identify potential genetic connectivity with the Belgium population. We used mitochondrial DNA (mtDNA) control region (D-loop) and microsatellite loci, to determine genetic diversity and spatial genetic structure on different timeframes among studied populations. Indeed, while the variability of the more conserved mtDNA control region reflects maternal lineages present in the founder individuals, microsatellites, due to their higher variability, reflect more recent events that shaped the current genetic structure.
We first aimed at confirming the raccoon introduction history in France using mtDNA, and we hypothesized that each established population is the result of independent introductions, i.e., the main haplotype should be different in each population. We also compared wild and captive animals' haplotypes to determine whether zoos could have been the source of the established populations. We then assessed the spatial genetic structure and the level of genetic variation using microsatellites. We hypothesized that each established population is genetically different. Genetic diversity is supposed to be the highest in the North-Eastern France population due to its larger census population size, its oldest origin, and the potential admixture with the large Belgium population and German populations, while the most recently established Gironde population is supposed to have the lowest genetic diversity. Given the good dispersal capacity of the raccoon (Cullingham et al. 2008b;Rioux Paquette et al. 2014), we also expected the presence of migrants, i.e., individuals genetically assigned to one population but located in another population or between populations. Finally, we assessed the presence of genetic bottlenecks as invasive species are predicted to suffer from reductions in genetic diversity during founder events (Dlugosch and Parker 2008), and estimated the effective population size. We expected longer established populations to show less or no evidence of bottlenecks and larger effective population sizes.

Study area and sampling methods
We collected 322 raccoon samples between 2002 and 2020 ( Fig. 1). Most tissue samples were collected from hunted, trapped or road-killed raccoons in a collaborative partnership with local NGOs. Three main areas were sampled for wild animals in France, particularly in North-Eastern France (n = 120), Massif Central (n = 82) and Gironde (n = 31) as these three locations constitute the three established populations in France. Non-residential individuals (n = 10) which could not be geographically associated with one of the three established populations were also sampled, as well as individuals from the Belgium population (n = 33). The sample also included individuals from French Guiana (n = 4), USA (n = 1), Latvia (n = 1), UK (n = 1), from seven zoos in France (n = 35), and from one illegal private owner (n = 4) ( Table S1). All wild animal captures, and captive animal sampling were performed according to the appropriate French laws for wildlife management and animal welfare; the procedures were approved by the pertinent administration. All tissue samples were stored at − 20 °C or conditioned in pure alcohol (ethanol 99.8%) and hairs were stored in hermetic bags prior to DNA extraction.

DNA extraction
DNA extraction of the 317 tissue and five hair samples was conducted under sterile conditions and in a room free of DNA. Sample tubes were surrounded by both negative extraction controls (blanks) and positive tissue controls. Positive tissue controls consisted in tissue sample previously validated in terms of genotyping success on microsatellite markers. Sample tubes, as well as positive tissue controls and negative extraction controls, were lysed overnight at 56 °C, and their DNA was isolated and purified using purification columns and vacuum filtration according to the manufacturer's instructions (Nucleospin 96 Tissue Kit, Macherey-Nagel). DNA was eluted to obtain final concentrations between 20 and 100 ng/ µl. We also extracted five hair samples using the DNeasy Blood & Tissue kit (QIAGEN) with several modifications of the protocol: volumes from 1, 2 and 3 steps were doubled, and DNA was eluted in 50 µl AE buffer.

Mitochondrial DNA amplification and sequencing
We amplified fragments of 566 base pairs from the mitochondrial control region which possesses the highest variation of the genes to confirm the species of the 322 tissue samples (i.e., Procyon cancrivorus or Procyon lotor) and for the phylogeographic analyses. We used the specific forward primer PLO-L15997; 5′-CCA TCA GCA CCC AAA GCT -3′ (Frantz et al. 2013) and the specific reverse primer PLO-CRL1; 5′ -CGC TTA AAC TTA TGT CCT GTA ACC -3′ (Cullingham et al. 2008a). PCRs were carried out in 25 μl volume containing 1.25 μl of each 10 µM primers, 2 μl of 2.5 mM dNTP, 2.5 μl of 10X Hifi buffer, 0.1 μl of 5 U/μl AccuStart Taq DNA polymerase (Quantabio), and approximately 20-25 ng of DNA. All PCR products were sequenced in both directions by Biofidal (Vaux-en-Velin, France). The sequences generated were visualized and analyzed using CLC Main Workbench (QIAGEN). Sequences were aligned in Seaview 4.7 (Gouy et al. 2010) using Muscle program with default parameters.

Microsatellite genotyping
A total of 271 individuals out of the 322 extracted were successfully amplified (Table S1). We amplified 236 wild animals: n = 23 individuals of the Gironde population, n = 72 from the Massif Central population, n = 110 from the North-Eastern population, n = 21 from the Belgium population and n = 10 from the Non-residential population (i.e., animals located in other parts of France). The remaining n = 35 come from zoos or private owners (i.e., captive animals).
For each DNA sample, 25 microsatellite markers and one marker for sex identification (ZFXY) ( Table S2) were amplified by two multiplex PCR (Polymerase Chain Reaction) A and B and analysed with an automated sequencer in two migrations. Three negative controls (blanks) and three positive DNA controls (raccoon DNA previously analysed and validated in terms of genotyping success and quality) were included per PCR reaction plate. PCR amplifications were performed in a 10 µl final volume containing 5 µl of mastermix Taq Polymerase (Multiplex PCR Kit, Qiagen), respectively 0.53 µL of a first pool of 13 pairs of primers or 0.53 µL of a second pool of 14 pairs of primers at a concentration of 0.06 to 0.60 (depending on the multiplex) and a mean of 30 ng of genomic DNA (Table S2). Our PCR thermal protocol consisted of 95 °C for 15 min, followed by seven touchdown cycles of 94 °C for 30 s, 61 °C to 55 °C for 90 s (decreasing 1 °C per cycle), and 72 °C for 60 s, then followed by 29 cycles of 94 °C for 30 s, 55 °C for 90 s, and 72 °C for 60 s, ending with an extension of 60 °C for 30 min.
PCR products were resolved on a calibrated ABI PRISM 3130 XL capillary sequencer (ThermoFisher Scientific) under denaturing conditions (HiDi Formamide, ThermoFisher Scientific) with an internal size marker that guarantees the same calibration for all samples. As all the samples were distributed on 3 plates and each plate contained the same positive reference controls (previously genotyped once), all positive controls were finally run 4 times on each marker and so guarantee both amplification and capillary resolution repeatability. As the positive and negative controls ensure the robustness of the process, each DNA sample was genotyped once as is usual for this type of samples. Electropherograms were analysed using GENEMAPPER 4.1 (ThermoFisher Scientific) independently by two analysts to determine the allele sizes for each marker of each individual. Reading errors were resolved, and in case of persistent disagreement, ambiguous results were considered as missing data. Genotype of each positive control (tissue controls and DNA controls) was compared to its known reference to ensure repeatability of extraction, PCR and analysis steps all along the process.
Comparison of the observed genotypes with the distribution of randomized genotypes generated with the program MICROCHECKER v.2.2.3 (Van Oosterhout et al. 2004) revealed that there were no null alleles in the data set.

Phylogeography: mitochondrial data analysis
A total of 322 D-loop sequences were obtained. The four samples of French Guiana were identified as Procyon cancrivorus and were thus removed from the analyses, and the remaining 318 samples were identified as Procyon lotor. We constructed two datasets to determine (1) the repartition and (2) the origins of the haplotypes present in populations of France and Belgium. The first one comprised only our sequences of P. lotor from France and Belgium. For the second dataset we collected all available data for P. lotor. We retrieved sequences from GenBank and haplotypes repartition from previous work (Biedrzycka et al. 2014;Cullingham et al. 2008a;Fischer et al. 2017;Frantz et al. 2013;Louppe et al. 2020). For all our analyzes we only kept the sequences as long as the sequences of the haplotypes determined in Culling- We identified the haplotypes present in our sequences, as well as the within population haplotypes diversity, and the differentiation among the three French and the Belgian populations, using DNAsp v6 (Rozas et al. 2017). Pairwise Fst based on haplotype frequency were estimated with Arlequin 3.5 (Excoffier and Lischer 2010).
We used the median-joining option (Bandelt et al. 1999) of PopArt (Leigh and Bryant 2015) to construct the networks for the two datasets. PopArt was also used to draw the distribution maps of the haplotypes present in France and Belgium.
Microsatellite data analysis and current genetic structure and diversity.

Population clustering
To explore the genetic structure across wild populations and to assess the source of putative immigrants, we focused on the 236 wild animals, and we searched for genetic groups using two approaches: multivariate analysis through discriminant analysis of principal components (DAPC) (Jombart et al. 2010), and population genetic model-based Bayesian clustering with the program STRU CTU RE (Hubisz et al. 2009;Pritchard et al. 2000). These two methods all aim to detect groups of genetically related individuals but differ in terms of population structure and the underlying evolutionary model.
We first used DAPC (Jombart et al. 2010) as implemented in the adegenet (Jombart 2008) package in R. DAPC maximizes differences among clusters while minimizing variation within and does not rely on a particular population genetics model, such as Hardy-Weinberg equilibrium. We applied the function find.clusters to determine the number of potential clusters K with 1 < K < 8. Minimization of a Bayesian Information Criterion (BIC) was used to identify the most probable number of clusters present in the data. We then used the DAPC itself to assign individuals to these clusters. In the data transformation step, 40 principal components of the PCA were retained based on the cross-validation procedure implemented with the function xvalDapc and discrimination was based on these 40 discriminant components. DAPC provides posterior membership probabilities (pp ik ) to  (Hubisz et al. 2009;Pritchard et al. 2000) was used assuming population admixture and correlated allele frequencies within populations. The Dirichlet parameter (λ), which describes how much the allele frequencies are correlated between populations, was inferred using the model of uncorrelated allele frequencies with K = 1 and its value was then used in the model of correlated allele frequencies (λ = 0.5687). Simulations were run with 10 6 Markov Chain Monte Carlo (MCMC) iterations after a burn-in of 10 5 iterations. K was varied from 1 to 8 and 20 independent simulations were run for each K value to check for consistency in the results. We used the ΔK criteria (Evanno et al. 2005) to identify the most probable number of clusters (K) present in the data. STRU CTU RE provides the relative contributions of subpopulations to the genome of the individual (q ik values) which sum to unity across clusters and were examined for geographic structure. We considered individual raccoons to be from a single cluster when their highest q ik value was > 0.8. Raccoons with the highest q ik value < 0.8 were considered as putative migrants or admixed individuals between two populations. Cluster assignments obtained from DAPC and STUC TUR E were compared for the sake of consistency of both methods.
Additionally, pairwise Fst values and their confidence interval between the genetic clusters identified with STRU CTU RE were calculated using the hierfstat (Goudet and Jombart 2020) package in R in order to quantify the level of genetic differentiation between the clusters. Only individuals with q ik value > 0.8 were considered.

Putative sources of wild population
To determine whether zoos or private owners could have been the source of established wild populations,  Jombart et al. 2010) in order to exemplify close genetic proximity between captive and wild individuals.

Genetic variation
All populations were expected to originate from an introduction of a small number of individuals and are thus expected to show a low level of genetic variation, traces of genetic bottleneck and low effective population size. In addition, the genetic variability of a population is correlated with its evolvability and adaptative potential. Therefore, for each genetic cluster identified with STRU CTU RE (only individuals with q ik value > 0.8), expected heterozygosity (He), observed heterozygosity (Ho), and rarefied allelic richness (Ar, El Mousadik and Petit 1996) were calculated using the hierfstat R package (Goudet and Jombart 2020). Inbreeding coefficient (Fis) values and significance were calculated using the diveRsity (Keenan et al. 2013) R package. The differences in He, Ho, Ar and Fis between the clusters were tested using linear mixed-effects models (LMMs). For each genetic diversity index, a model with the index as response, with the genetic clusters as fixed effect and the locus as a random effect was fitted using the lmer function of the lme4 (Bates et al. 2015) R package. In order to test the effect of the clusters, we ran a likelihood ratio test (LRT), comparing the model with the genetic clusters as fixed effect with a null model that excludes it. Differences among different genetic clusters were tested by Tukey HSD post-hoc comparisons using the function glht of the multcomp (Hothorn et al. 2008) R package. We also investigated the possibility of recent bottlenecks within each genetic cluster using the heterozygosity-excess approach (Cornuet and Luikart 1996) implemented in the BOTTLENECK v1.2.02 software (Piry et al. 1999). As some clusters had less than 20 polymorphic loci (see "Results" section), we used the one-step stepwise mutation model (SMM) and also the two-phase model (TPM) with a variance among multiple steps of 12, with 95% single step mutations and 5% multi-step mutations (Piry et al. 1999). Models were run with 10,000 replicates and significance of heterozygosity excess over all loci was determined with a Wilcoxon sign rank test.
Finally, we estimated effective population sizes (Ne) using the linkage disequilibrium method of Waples and Do (2008) as implemented in the NeEstimator v2.1 software (Do et al. 2014). We assumed random mating, screened out alleles with frequencies below 0.02, and used the parametric method for estimating confidence intervals (Jones et al. 2016).

Phylogeography and putative populations of origin
We obtained 322 D-loop sequences, 318 of Procyon lotor of and four of P. cancrivorus which were removed from the analyses. Almost all the D-loop sequences of P. lotor correspond to sequences already present in GenBank, except for the sample WB14008 from the USA. Our 318 sequences P. lotor correspond to 10 haplotypes, 9 of them present in France and Belgium (Fig. 2). For the wild samples, one haplotype is dominant in each of the three French areas: PLO13 in North-Eastern France (72%), PLO22 in Gironde (97%) and PL032 in Massif Central (95%). In North-Eastern France, close to the Belgian border, another haplotype (PLO2b) is also present in 23% of the samples. In Belgium there are two main haplotypes PLO2b (45%) and PLO13 (33%). Three haplotypes are present only in France (PLO32, PLO22, PLO17) and two only in Belgium (PLO2a, PLO110). The genetic diversity is higher in Belgium and North-Eastern France than in Gironde which is the most recently colonized French area (Table 1). Six haplotypes were found in captive animals (PLO2a, PLO2b, PLO13, PLO32, PLO110, PLO63) (Fig. 2). There are two main haplotypes in captive samples PLO2b (26%) and PLO13 (59%) which are also dominant in North-Eastern France and Belgium. PLO32 is dominant in Massif Central, PL02a and PLO110 were only found in Belgium and the last one PLO63 was only found in one zoo and not in the wild. Two zoos present an important haplotype diversity with 4 and 5 haplotypes for respectively 8 and 12 samples.
At a larger scale, the haplotypes present in France and Belgium belong to the three lineages described by Cullingham et al. (2008a) (Fig. 3). PLO22 (Gironde) belongs to the lineage I, PLO32 (Massif Central) to the lineage II, and the two main haplotypes present in North-Eastern France and Belgium (PLO13 and PLO2) respectively belong to the lineage II and III. Some of the haplotypes present in France (PLO02, PLO13, PLO32, PLO78 and PLO110) are also confirmed in the neighbouring European countries such as Germany or Spain whereas others are restricted to one country such as PLO22 (France), PLO66 (Spain) or PLO16 (Germany) (Fig. 3).
The Fst values estimated from D-loop haplotype frequencies (Table 2) are rather high and the lowest value (Fst = 0.476) was obtained between North-Eastern France and Belgium, two areas sharing several haplotypes.

Population clustering based on microsatellite variation
For both clustering methods, K = 4 was the most likely number of genetic clusters as evidenced by Evanno's index for STRU CTU RE (Fig. 4a) and the BIC plot for the DAPC (Fig. 4b). For the DAPC, all individuals had pp ik values > 0.8 (Fig. 5a) while 12 individuals had their q ik values < 0.8 with STRU CTU RE (Fig. 5b). Among the individuals with a membership probability > 0.8, the clustering solution found both by STRU CTU RE and DAPC were identical and had strong spatial consistency, with the Belgium, Gironde and Massif central populations belonging to three different clusters (cluster 1, 2 and 3 respectively), (Fig. 6, Fig. S1). The North-Eastern population was split in two, with individuals close to the Belgian border belonging to the Belgian cluster 1, and individuals further south belonging to a fourth genetic cluster (cluster 4). Among the 12 individuals which could not be assigned to a single cluster (i.e., with q ik values < 0.8), most of them showed admixture at most between two populations reflecting their spatial proximity (e.g., individuals of the North-Eastern population showing admixture between the North-Eastern and Belgium population, Fig. 6). However, admixture between distant populations was also identified with two individuals belonging to the clusters 2 and 4. Admixed individual showing mixture between the cluster 1 and 4 (Belgium and North-Eastern populations) was also found in the Gironde population, while an individual belonging to the Belgium cluster was found in the Massif central population (Fig. 6).
Because of its greater sensitivity, we decided to focus our analyzes on the STRU CTU RE results, keeping only individuals showing the strongest membership probability, i.e., q ik values > 0.8. Strong genetic differentiation was found among the four clusters with high significant pairwise Fst values (none of the confidence intervals overlapped 0, Table 3). The lowest value (Fst = 0.144) was found between the cluster 1 and 4, while the highest pairwise Fst values were found in each pair involving the cluster 3 (Gironde, Table 3).
The assignment of the 35 captive individuals to the four clusters identified by the DAPC revealed poor spatial constituency. Most of the individuals (30) sampled everywhere in France belonged to the cluster 1 (Belgium/North-Eastern populations, Fig. S1). Out

Fig. 5
Genetic clustering. Genetic clusters as inferred by a STRU CTU RE and b DAPC for the most likely number of genetic clusters K = 4. The relative ancestry proportion (q ik for STRU CTU RE) and membership probabilities (pp ik for DAPC) of each individual for each cluster are shown in different colors. Individuals are ordered first by geographic location (Belgium, North-Eastern France, Non-residential, Massif Central, and Gironde populations), and second by decreasing q ik and pp ik within each geographical population of the 5 remaining ones, three individuals belonging to the cluster 2 (Massif central population) were found close to Belgium population, and one individual belonging to the cluster 4 (North-Eastern populations) was found in a zoo of North-Western France. Only one individual showed some spatial constancy and was found close to the Massif central population to which it genetically belongs (Fig. S1).

Microsatellite genetic variation
The cluster 1 (Belgium/ North-Eastern) showed the highest genetic diversity with an average Ho of 0.61 (± 0.13), average He of 0.65 (± 0.12), with a mean Ar of 5.13 (± 1.37), and a total number of alleles of 150 (Table 4). The lowest genetic diversity was found for the cluster 3 (Gironde) with an average Ho of 0.39 (± 0.30), average He of 0.38 (± 0.27), with a mean Fig. 6 Spatial genetic structure. Assignment of the individuals by STRU CTU RE to the four genetic clusters. Individuals with a strong ancestry proportion (i.e., with q ik > 0.8) are displayed as a circle with a single color while individuals with a lower ancestry proportion (i.e., with q ik < 0.8) are displayed as pie charts partitioned into several colored parts proportionally to their ancestry proportion assigned to the four clusters Ar of 3.5 (± 0.87), and a total number of alleles of 57 (Table 4). Fis values were low, and only significant for the cluster 3 (Gironde , Table 4). For all the indexes, likelihood ratio test comparing the model including the clusters as a fixed factor with a model that excludes it were significant (all p < 10 -3 , Table 5). Post-hoc comparisons revealed that the cluster 1 (Belgium/North-Eastern) had the significantly highest genetic diversity while the cluster 3 (Gironde) had the lowest one (Fig. S2).
Contemporary effective population size estimations were low for all populations, with Ne ranging from 19.8 for the cluster 4 (North-Eastern), to 66.2 for the cluster 2 (Massif central, Table 4). Note that the sample size of cluster 3 (Gironde) was too low to reliably estimate the confidence interval. Contemporary bottleneck was detected within all clusters with both the TPM and SMM models (except for the TPM model for the cluster 1, p = 0.11, Table S3).

Discussion
In this study we used both mitochondrial and microsatellite data to analyze the genetic diversity and current spatial genetic structure of all the French raccoon populations to determine their introduction history and their invasiveness potential. We also investigated whether French and Belgium populations were connected.   Table 5 For each genetic diversity index, values of the models including the clusters as a fixed factor and the model that excludes it The number of parameters (K), Akaike's information criterion corrected for small sample size (AICc), Akaike differences (ΔAICc), and log-likelihood (LL) are indicated for each model. Results of the likelihood ratio test comparing the model including the clusters as a fixed factor with a model that excludes it are also shown, as well as the marginal (R 2 m ) and the conditional (R 2 c ) R 2 Spatial genetic structure and genetic diversity Both the DAPC and STRU CTU RE converged to the same clustering solution, evidencing the presence of three genetic clusters in France. These three genetic clusters have a strong spatial consistency with the geographical distribution of the French populations and are genetically different from the Belgium population which led to a fourth cluster (but see North-Eastern France population). This tends to confirm that, without evident geographic barriers, these populations are the result of at least three independent introductions (Léger and Ruette 2014): they are not the result of a natural colonization from the Belgium population nor from the oldest North-Eastern one. The same pattern was observed using mtDNA: populations from France and Belgium are well differentiated with different haplotype panels, even if the difference tends to be less important between Belgium and neighboring French populations. This tends also to plead for multiple introductions. This pattern seems to be typical of raccoon invasions and has been evidenced in several countries such as Spain (Alda et al. 2013), Germany (Fischer et al. 2015) and central Europe (Biedrzycka et al. 2014;Fischer et al. 2017). Even though it is not directly comparable due to differences in the microsatellites used in this study, the genetic diversity of the three detected French populations was relatively low (He: 0.38-0.50; Ar: 2.28-3.5) compared to that of the native range of the raccoon (e.g., He: 0.41-0.90;Ar: 4.06-11.58, Santonastaso et al. 2012), of the Belgium population (He: 0.65; Ar: 5.13), or even to that of larger introduced populations in Europe (e.g., He: 0.54-0.0.68; Ar: 4.41-6.25, Biedrzycka et al. 2014). Mitochondrial genetic diversity in France and Belgium (8 haplotypes for 280 wild samples) was also rather low compared to the diversity observed in the native range of the raccoon (76 haplotypes in 311 samples, Cullingham et al. 2008a) but a bit higher than the ones observed in other European populations such as in Spain (4 haplotypes) or Central Europe (6 haplotypes) (Fischer et al. 2017). These low estimations are in line with the detection of contemporary bottleneck in each of these populations and with the estimated effective population size which are 5-43 fold lower (Ne: 19.8-66.2) than the estimations in the native range (e.g., Ne: 370.1-854.0, Alda et al. 2013). These results likely correspond to a founder effect (i.e., population founded by a small number of individuals) exacerbating the effects of stochastic processes like genetic drift on genetic variation (Dlugosch and Parker 2008) and are still detectable 60 years after first introduction in North-Eastern France.
Despite these low levels of genetic diversity, the three populations are still expanding (Léger and Ruette 2014). Small founding populations can lose genetic variation via drift, and inbreeding among members of small populations can lead to low fitness. However, many genetically depauperate introduced populations are able to continue their expansion, to adapt, and are incredibly successful, a phenomenon called the "genetic paradox of invasions" (Allendorf and Lundquist 2003). Considering the genetic bottlenecks detected and the succesful invasion, the French raccoon populations seem to be a good example of the invasion paradox (Estoup et al. 2016). However, despite being lower than in the species native range, the observed genetic diversity is in the range of terrestrial vertebrates (DeWoody and Avise 2000). Furthermore, it has also been shown that variation at key parts of the genome permits rapid adaptation even in a population with low overall genetic diversity (Tepolt et al. 2021). It is thus likely that, in addition to their generalist dietary and habitat requirements and their great behavioral flexibility (Daniels et al. 2019), the remaining genetic diversity is enough to successfully adapt to a new environment and allow a quick colonization. A landscape genomics study (Balkenhol et al. 2019) would be needed to detect regions of the genome associated with adaptation to environmental variables and help to understand the mechanisms supporting quick adaptation to new environment of invasive species.
The geographical origins in the raccoon native range could not be precisely determined. The haplotypes present in France belong to three lineages described in North America, but their spatial distributions overlap in several states of the US (Cullingham et al. 2008a). Two haplotypes PLO2 and PLO13 dominant in North-Eastern France and Central Europe (Fischer et al. 2017) are also widespread in North America (Cullingham et al. 2008a). The other ones have a more restricted distribution (SE USA (PLO17), NE USA (PLO22, PLO32), Florida (PLO78), Cullingham et al. 2008a) but are still present in several other states. The haplotype PLO110 was only described from Germany and not sampled in the native range (Frantz et al. 2013).

The North-Eastern France population case
Both the mtDNA and microsatellites clustering results showed that North-Eastern France and Belgium populations are both spatially and genetically connected as evidenced by the presence of several admixed individuals. The ongoing gene flow between these two populations, also evidenced by the lowest pairwise Fst value, calls for the control of the raccoon in the region. Previous studies have shown that mixing different lineages following multiple introductions and their spatial expansion is an important factor of many successful biological invasions (reviewed in Forsman 2014). For example, Wagner et al. (2017) experimentally demonstrated that invasions of the beetle Callosobruchus maculatus resulting from more than one source spread farther than single-source invasions. Admixture may increase the genetic diversity of the merged population and thus, may increase the evolutionary potential to respond to novel selection pressures, thus facilitating their expansion (Dlugosch and Parker 2008;Keller and Taylor 2010). Additionally, admixture may cause heterosis (fitness advantage associated with heterozygosity) which in turn may have a 'catapult effect', i.e., a fitness benefit during the early stages of the introduction that increases the odds of long-term establishment (Drake 2006). Control programs should particularly target these two populations and aim to reduce the opportunity for further admixture events.

Dispersal event and captive sources of the wild populations
We evidenced genetic admixture between North-Eastern France and Massif Central populations for two individuals, and between Massif Central and Gironde populations for two other individuals (Fig. 6) suggesting that these populations have been recently connected through some long-distance dispersal events (> 200 km). This is coherent with raccoon's movement capacities as they have been observed travelling distances from 20 to 200 km (Cullingham et al. 2009;Dharmarajan et al. 2009). Despite being recent, these movements and the resulting admixture raised concerns about a potential increase in the invasiveness of the raccoon in these regions (see previous paragraph). In its native range, raccoon dispersal has been shown to depend on the landscape structure and composition (Moncrief et al. 2017;Rioux Paquette et al. 2014). Thus, a landscape genetics study (Manel et al. 2003) is particularly recommended to identify the environmental variables that could favor or impede raccoon dispersal and identify high connectivity areas that should be intensively monitored and where the species should be controlled to prevent gene flow among established populations.
Another individual showing admixture between the Belgium and North-Eastern France populations has been caught in the Gironde population. Given the large geographical distance between these populations, this admixture is unlikely to be the result of natural dispersal but rather of a translocation by humans or an escape from a zoo. While neither the assignment of the 35 captive individuals using the DAPC nor their mitochondrial lineages allowed to identify a potential captive origin of the wild populations, individuals sampled in a nearby zoo showed the same haplotype and belonged to the same genetic cluster as this admixed individual. The way in which the raccoons are kept in this zoo should be inspected. Most of the captive individuals were genetically related to the Belgian population. It should also be noted that one individual has been sampled in a zoo close to the Massif central population to which it belongs genetically. The history of captive individuals in this zoo should be investigated more thoroughly to determine whether former individuals could be at the origin of this population or if, conversely, the zoo has collected individuals from wild populations. If occasional introductions from private breeding or keeping of animals continue, the control of the established populations is likely to be of limited success. Indeed, given the detected bottleneck and the low effective population size, it is clear that populations have been founded by a very small number of individuals, showing that the release of few individuals can results in large established populations if nothing is done, and even a low number of individuals will allow the populations to rapidly recover (Alda et al. 2013). It is thus essential to track illegal or accidental raccoon releases and potential range expansion.

Conclusion
We demonstrated that the French established raccoon populations are the result of independent introductions, during which they have lost an important part of their genetic diversity. Despite this, they continue to expand, suggesting that the remaining diversity is enough to successfully invade to their new environment. Additionally, admixture areas have been identified with the largest population in Belgium (Wallonia) and among the populations in France, suggesting the risk of an increase in raccoon invasiveness in these areas. Intensive management measures should target these areas and try to limit gene flow among populations. Genetics can be used to identify priority areas for raccoon control, and to genetically inform control campaigns on progress (e.g., effective population sizes, bottlenecks). To set up an effective control campaign will require setting common management objectives, information flow and coordination of control actions across international borders (Adriaens et al. 2019). However, such management program can only be successful if it is combined with a strict control of the captive raccoon population in Belgium and France. The EU Invasive Alien Species regulation 1143/2014 prohibited its sale or exchange, reproduction, transport, to private owners. The keeping of raccoons in captivity is also prohibited unless the animal was legally kept before the regulation came into force. In France, despite the EU Regulation, captive breeding of raccoons in zoos is still possible, and there is no enforcement of increased biosecurity in detention facilities. Such biosecurity measures might not be sufficient to prevent new introductions as few introduced individuals might be enough to create new populations.