A comprehensive phylogeographic study of Arion vulgaris Moquin-Tandon, 1855 (Gastropoda: Pulmonata: Arionidae) in Europe

Arion vulgaris Moquin-Tandon, 1855 is regarded as one of the 100 most invasive species in Europe. The native distribution range of this species is uncertain, but for many years, the Iberian Peninsula has been considered as the area of origin. However, recent studies indicate that A. vulgaris probably originated from France. We have investigated the genetic structure of 33 European populations (Poland, Norway, Germany, France, Denmark, Switzerland) of this slug, based on two molecular markers, mitochondrial cytochrome c oxidase subunit I (COI, mtDNA) and nuclear zinc finger (ZF, nDNA). Our investigation included published data from two previous studies, giving a total of 95 populations of A. vulgaris from 26 countries. This comprehensive dataset shows comparable haplotype diversity in Central, North and Western Europe, and significantly lower haplotype diversity in the East. All haplotypes observed in the East can be found in the other regions, and haplotype diversity is highest in the Central and Western region. Moreover, there is strong isolation by distance in Central and Western Europe, and only very little in the East. Furthermore, the number of unique haplotypes was highest in France. This pattern strongly suggests that A. vulgaris has originated from a region spanning from France to Western Germany; hence, the slug is probably alien/invasive in other parts of Europe, where it occurs. Our results indicate the necessity to cover as much of the distribution range of a species as possible before making conclusive assumptions about its origin and alien status.


Introduction
Alien species can have negative effects on natural ecosystems including the displacement of endemic taxa (Sax et al. 2005;Davis 2009;Hatteland et al. 2013), which in turn can adversely affect food security by impacting ecosystem services. Therefore, alien and invasive species may be a threat to global biodiversity because of their impact on biological communities and economic costs (Walther et al. 2009;Kozłowski and Kozłowski 2011). Anthropogenic activity strongly increases the spread and impact of many aliens, especially gastropods and insects, which are easily spread with various trade goods. Among terrestrial invertebrates, Arion vulgaris Moquin-Tandon, 1855 (Gastropoda: Pulmonata: Arionidae), also known as the Iberian or Spanish slug, is regarded by DAISIE organization (Delivering Alien Invasive Species Inventories for Europe) as one of the 100 most invasive species in Europe (Rabitsch 2006). It is of huge economic significance because it causes serious damage to agricultural crops and orchard cultivation resulting in financial losses (Frank 1998;Gren et al. 2009;Kozłowski 2010Kozłowski , 2012. For many years A. vulgaris has been regarded as Arion lusitanicus (Mabille 1868), thought to have spread from the area of the Iberian Peninsula (Spain, Portugal and Azores) to other countries (Simroth 1891;Quick 1952Quick , 1960van Regteren Altena 1971;Chevallier 1972). However, Castillejo (1997) showed that A. lusitanicus sampled outside Iberia was very different from the A. lusitanicus (Mabille 1868) originally described from Portugal based on internal morphology. Phylogenetic studies by Quinteiro et al. (2005) and Hatteland et al. (2015) came to the same conclusion; thus, A. vulgaris has been recognized as distinct from A. lusitanicus species, as also suggested by Falkner et al. (2002). Due to the morphological similarities and confusing taxonomy of older studies, the name A. lusitanicus often refers to the so called pest slug A. vulgaris commonly occurring in Europe. Currently, A. vulgaris is considered to be common in many European countries, from south to north, and up to now has been reported from France (in 1855), Great Britain (1952), Germany (1970), Slovenia (1970), Italy (1971), Switzerland (1971), Austria (1972), Sweden (1975), Bulgaria (1983), Austria (1984), Norway (1988), Belgium (1989), the Netherlands (1989), Finland (1990), Denmark (1991), Poland (1993), Iceland (2003-2004, Greenland, Latvia, Lithuania, Faroe Islands and Romania (Weidema 2006;Zając et al. 2017).
A. vulgaris is regarded as an invasive species in most studies, whereas another two species of terrestrial slugs, which have similar external appearance-A. ater (L., 1758) and A. rufus (L., 1758)-are considered as native to North-Western Europe (Welter-Schultes 2012; Hatteland et al. 2015;Zając et al. 2017). These three species may co-exist in natural habitats and their differentiation is possible on the basis of molecular analyses, mating behaviour or genital morphology (Noble 1992;Quinteiro et al. 2005; Barr et al. 2009;Kałuski et al. 2011;Slotsbo et al. 2012;Hatteland et al. 2015). Moreover, in these three species, introgression has been observed and so far two types of hybrids have been discovered A. ater × A. rufus and A. vulgaris × A. rufus (Evans 1986;Roth et al. 2012;Dreijers et al. 2013;Hatteland et al. 2015). Introgression seems to be commonest between A. ater and A. rufus, and also occurs between A. ater and A. vulgaris (Hatteland et al. 2015).
Recently, two contradictory studies about the origin, migration routes and status (invasive or native) of A. vulgaris have been published (Pfenninger et al. 2014;Zemanova et al. 2016). Pfenninger et al. (2014) used two molecular markers, one mitochondrial, cytochrome c oxidase subunit I (COI), and one nuclear, zinc finger (ZF), together with species distribution modelling, to investigate the invasive status of this pest species. Based on these results, Pfenninger et al. (2014) suggested that A. vulgaris is probably native to Central Europe, because haplotypes from Spain or western France did not cooccur with those from Central Europe. On the other hand, a lack of data prevented testing for Eastern Europe as a place of origin for A. vulgaris. Zemanova et al. (2016) focused on A. vulgaris' status to verify whether it is invasive or native to Central Europe. Their analyses of mtDNA sequences (COI and ND1 markers) and microsatellites, obtained from European populations across the distribution range of this slug, indicated that classification of A. vulgaris as non-invasive would be premature.
The main aim of our study was to make a comprehensive comparison of genetic variability between and within European populations of A. vulgaris to shed light of the origin of this species. Moreover, this study aimed to explore and summarize the genetic structure of A. vulgaris populations in Europe.

Material collection
Adult specimens of Arion slugs were collected into plastic boxes in 2014 from 42 populations across eight European countries: Norway (13), Poland (12), France (4), Denmark (3), Germany (3), Switzerland (3) and Netherlands (1) (see Table 1, Supplementary Materials for details on sampling localities). Only mature slugs with genital morphology corresponding to the original description of A. vulgaris Moquin-Tandon, 1855 and A. lusitanicus (Welter-Schultes 2012) were collected (5-15 slugs from each locality). All individuals were sent to the Institute of Environmental Sciences (Jagiellonian University), or the Norwegian Institute of Bioeconomy Research (NIBIO), or University of Bergen, and subsequently prepared for further morphological and molecular analyses. Remaining soft body parts were stored in 70% ethanol and material was deposited in the Institute of Environmental Sciences (Jagiellonian University, Poland), in the Norwegian Institute of Bioeconomy Research (NIBIO, Norway) or in University of Bergen (Norway).

Species identification
Species identification was based on detailed examination of dissection of the distal genitalia under a stereomicroscope (Olympus SZX10). The main characters which allow A. vulgaris to be distinguished from other large arionids are as follows: a small atrium, almost symmetrical, one-partite, bursa copulatrix, a fallopian tube with a short, thin posterior end and a thick, rapidly expanding anterior end, a long oviduct with large, asymmetric ligula inside (Riedel and Wiktor 1974;Wiktor 2004;Hatteland et al. 2015). Slugs after dissection were preserved in 96% ethanol. During dissection, a piece of tissue of each specimen examined was preserved separately in 96% ethanol in − 80°C for DNA extraction.
DNA extraction, amplification and sequencing DNA extraction was performed according to Montero-Pau et al. (2008). The samples extracted following the latter procedure were initiated placing the tissue samples in 0.2-ml tubes with 50 μl of alkaline lysis buffer (pH = 12). Samples were placed at 95°C for 30 min., after which they were kept on ice for 4-5 min. A total of 50 μl of neutralizing solution (pH = 5) was added to each tube and samples vortexed and centrifuged. PCR generated fragments of the cytochrome oxidase subunit I (COI, mtDNA), and of zinc finger (ZF, nDNA) gene using the following primers, LCO1490 (5′-GGTC AACAAATCATAAAGATATTGG-3′), HCO2198 (5′-TAAA CTTCAGGGTGACCAAAAAATCA-3′) (Folmer et al. 1994), and ZF forward (5′-TGTTTGGAATCTAGAGCCTG A-3′), ZF reverse (5′-TTTGTTGCTGTCCTGCTTGT-3′), respectively (Pfenninger et al. 2014). These two markers have a considerable mutation rate making them appropriate for intraspecific comparisons, as well as for comparisons between closely related species. PCR of each sample was performed in a 10 μl reaction mixture, consisting of 1 μl of template DNA, 0.2 μl of each primer, 1 μl of 10× buffer, 1 μl of 25 mM MgCl 2 , 6.5 μl of ddH 2 O, 0.1 μl of 20 mM dNTP and 0.1 μl of Taq-Polymerase. PCR conditions for COI consisted of 5-min denaturation at 95°C, followed by 30-s denaturation at 95°C, 1 min annealing at 48°C, and 90-s elongation at 72°C for 35 cycles, followed by a final elongation step for 10 min. PCR conditions for ZF consisted of 1min denaturation at 95°C and then 30-s denaturation at 95°C, 45 s annealing at 52°C and 1-min elongation at 72°C for 38 cycles, followed by a final elongation step for 3 min. A 1-μl sample of PCR product was run on a 1.4% agarose gel for 30 min at 100 V to check DNA quality. PCR products were diluted depending on this quality control and then sequenced in one direction. The sequencing reactions were performed in the BiK-F Laboratory Centre, Frankfurt am Main, Germany, in Molecular Ecology Lab, Institute of Environmental Sciences, Jagiellonian University, Kraków, Poland, and at the Sequencing Facility, University of Bergen, Bergen, Norway. A total of 307 COI sequences and 285 ZF sequences were obtained from tissue samples. Unique haplotype sequences which were not detected in previously published studies of both molecular markers were deposited in GenBank with the following accession numbers: COI: MK598594-607, MK617770, ZF: MK679547-MK679552.  (Katoh et al. 2002;Katoh and Toh 2008), using default settings. All obtained alignments were checked manually against non-conservative alignments in BioEdit 5.0.0. (Hall 1999). Aligned sequences were trimmed, and in case of COI translated into protein sequences in MEGA version 7, to check against pseudogenes (Kumar et al. 2016). Uncorrected pairwise distances were calculated in MEGA version 7 (Kumar et al. 2016).

Comparative and phylogenetic analyses
Phylogenetic trees were constructed based on single COI sequences representing each discovered haplotype within A. vulgaris, and published COI sequences of other species from the genus Arion; Arianta arbustorum (Accession number AF296945) was used as the outgroup (Table 4, Supplementary Materials). In total, 156 sequences were used for construction of the phylogenetic tree. The obtained alignment (482 bp) was divided into three data blocks, constituting three separate codon positions using PartitionFinder version 2.1.1 under the Akaike Information Criterion (AIC) (Lanfear et al. 2016). The best scheme of partitioning and substitution models was chosen for posterior phylogenetic analysis, and the analysis run to test all possible models implemented in the program. As best-fit partitioning scheme, PartitionFinder suggested retaining three predefined partitions separately. The best fit-models for these partitions were as follows: GTR + I + G for the first codon position, GTR + I + G for the second codon position and TVM + G for the third codon position. Models GTR, GTR + I, GTR + G and GTR + I + G were also tested, because RAxML allows for only a single model of rate heterogeneity (of the GTR family) in partitioned analyses (Stamatakis 2014). The best fit-model for all partitions in this analysis was GTR + I + G.
Maximum-likelihood (ML) topology for COI was constructed using RAxML v8.0.19 (Stamatakis 2014). Strength of support for internal nodes of ML construction was measured using 1000 rapid bootstrap replicates. Bootstrap (BS) support values ≥ 70% on the final tree were regarded as significant statistical support. Bayesian inference (BI) marginal posterior probabilities were calculated using MrBayes v3.2  with 1 cold and 3 heated Markov chains for 10,000,000 generations, and trees were sampled every 1000th generation. Each simulation was run twice. Convergence of Bayesian analyses was estimated using Tracer v. 1.3 (Rambaut et al. 2014). All trees sampled before the likelihood values stabilized were discarded as burn-in. All final consensus trees were viewed and visualized by FigTree v.1.4.3, available at: http:// tree.bio.ed.ac.uk/software/figtree.
Median joining haplotype networks of COI (482 bp length of alignment) and ZF (404 bp length of alignment) for A. vulgaris were performed in PopART (Bandelt et al. 1999). The haplotype networks were constructed from all haplotype sequences presented in each studied population (see Tables 5 and 6, Supplementary Materials for population list). Haplotypes are identified by numbers; circles without a number indicate a hypothetical intermediate haplotype necessary to link observed haplotypes. Hatch marks in the network represent single mutations. In order to estimate the population genetic parameters, number of haplotypes (NH), polymorphic sites (S), nucleotide diversity (π) and haplotype diversity (H) were calculated in DnaSP v.5.10 (Librado and Rozas 2009).
For all localities with sample size larger or equal to two individuals of A. vulgaris, haplotype diversity was plotted as bubbles against latitude and longitude for both loci. The larger the bubble, the higher haplotype diversity was observed at a particular locality. After this step, all populations were divided into four regional groups: East, West, Central and North based on their geographical coordinates. A total of 81 (for COI marker) and 54 (for ZF marker) populations of A. vulgaris were partitioned into a Western (longitude < 8°, latitude < 58°), Central (longitude > 8°and < 12°, latitude < 58°), Northern (comprising all populations in Scandinavia north of the Oresund) (latitude > 58°) and an Eastern (longitude > 8°and latitude < 58°) groups ( Fig. 1; Tables 1, 2 and 3 in Supplementary Materials). Division into regional groups is contextual and refers to the partitioning of sampling areas, not to the geographical division of Europe. Haplotype diversity was calculated once again for populations grouped in larger parts (West, Central, East, North groups) and visualized as separate figures for COI and ZF markers. For microsatellites, allelic richness was calculated, and for these analyses, data was adopted from the study of Zemanova et al. (2016). Similarly, allelic richness per region was visualized as figure. In order to obtain the statistical significance of the differences between regions, Mann-Whitney U test was performed because the data was not normally distributed.
Population differentiation was measured as population pairwise Φ ST values, estimated in ARLEQUIN v. 3.5.1.2 (Excoffier and Lischer 2010). To account for the effect of geography on population differentiation (isolation-by-distance), we regressed the population pairwise Φ ST s against the pairwise geographical distances. All statistical analyses were performed with the software package PAST v. 3 (Hammer et al. 2001).

Species identification
Identification of collected Arion species was done on the basis of morphology of the reproductive system together with molecular confirmation (sequences of cytochrome c oxidase subunit I (COI)). We investigated 39 European populations, 33 of them belong to A. vulgaris (Norway-11, Poland-9, Germany-3, France-4, Denmark-3, Switzerland-3). Additionally, we have detected three populations of A. rufus (Poland: Wronki, Nakielno, Zbąszyń) and three populations of hybrids: two of them were A. ater × A. rufus (Norway) and one A. vulgaris × A. rufus (Norway) (  Fig. 2 for haplotype network). The 482-bp COI fragment has 46 polymorphic sites (9.54%). Haplotypes 1 and 3 were the most frequent and present in most countries, occurring in 25.26% (24 localities) and 55.79% (53 localities) populations of its total number, respectively (Table 5, Supplementary Materials). Haplotype 1 was found in Central, Northern and Western Europe (Austria, Denmark, Faroe Islands, Germany, Hungary, Netherlands, Norway, Poland), whereas haplotype 3 was found in all countries except Spain. However, haplotype 3 has been found at several localities near the Pyrenees, in close proximity to the border with Spain (Boussenac, Selonnet localities).
The second molecular marker used in this study was ZF, the 404-bp fragment has 25 polymorphic sites (6.19%).We have analysed 220 of our original sequences of A. vulgaris together with 154 sequences deposited in GenBank as Arion sp. from Pfenninger et al. (2014). Some individuals identified as A. vulgaris based on mtDNA and genital morphology had nDNA typical for A. rufus (Bodø and Neauphle-le-Vieux populations), possibly indicating introgression. These latter sequences have been removed from the analyses. A total of 24 distinct haplotypes were identified (Fig. 3). Haplotype 1 was the most common, being present in 44 of 54 localities, representing all countries except Slovenia (Table 6,   Supplementary Materials). The second most common haplotype (haplotype 2, one mutational step away from haplotype 1) was present in 23 populations from seven countries (Austria, Denmark, France, Germany, Norway, Poland and Switzerland). Most haplotypes are separated by only one or two nucleotide substitutions (Fig. 3). Some haplotypes were limited to particular countries, for example, haplotypes 4 and 6 occur only in Denmark, whereas others have been more widely geographically distributed. In total, 15 of 24 haplotypes occurred in a single locality, mainly in Central Europe (France and Germany) ( Table 6, Supplementary Materials). The highest genetic diversity, defined by number of haplotypes in a population, was observed in Germany (Ottenhöfen-10 haplotypes, Munich 3-4, Nieder-Olm-4, Munich 2-3, Schlüchtern-3), Denmark (Aarhus, Silkeborg, Lisbjerg-4), France (Colmar-4, Boissy-Fressnoy, Neauphle-le-Vieux, Ornans-3), Austria (Bichlbach-3), a n d N o r w a y ( H o r t e n -3 ) ( Ta b l e 1 , Table 6 i n Supplementary Materials). The highest genetic distance between haplotypes within populations was observed in German localities: Ottenhöfen (range from 0.25 to 2.97%), Munich 3, Nieder-Olm (0.25-1.24%) and French: Ornans (0.25-1.24%). The genetic pairwise distances between all ZF haplotypes ranged from 0.2 to 3.2% (1.1% on average). The other population genetic parameters were π = 0.00502, H = 0.808. The average genetic distance (calculated on the basis of COI sequences available at GenBank) between species of the Arion genus is 18.2%. For all localities with sample size larger or equal to two individuals of A. vulgaris against latitude and longitude, haplotype diversity was calculated for COI and ZF markers, separately (Fig. 4). This resulted in a clear spatial pattern of diversity at both loci. There are populations with zero diversity all over the range (probably due to a pronounced metapopulation structure, where new suitable habitat patches are colonised by a low number of individuals), whereas there are many more localities devoid of haplotype diversity in the East part of Europe. However, this pattern is not so visible in the Northern localities.
All studied populations of A. vulgaris were partitioned into four, approximately equally strong regional groups (West, Central, East and North) based on their geographical coordinates ( Fig. 1). At all markers, there was no notable differences among the Western and the Central groups, while the Eastern group showed definitely the lowest mean values (Table 1, Fig. 5a, b). Statistically significant differences in haplotype diversity per regions were observed between Eastern and Western region (Z = − 2.67228; p = 0.007534) in the case of COI and between Central and Eastern regions (Z = 3.069854; p = 0.002142) for the ZF marker. Haplotype diversity for COI (a) and ZF (b) markers calculated for all localities of A. vulgaris with sample size larger or equal to two individuals against latitude and longitude. The larger bubble size, the higher haplotype diversity was observed at a particular locality Table 1 Haplotype diversity and allelic richness for A. vulgaris populations divided into four regional groups for COI, ZF markers and microsatellites, respectively. Data of allelic richness were adopted from the study of Zemanova et al. (2016)  This overall pattern is confirmed by microsatellite data from the study of Zemanova et al. (2016). Partitioning populations into the same regions showed that there is no notable difference in allelic richness between Western and Central regions at the microsatellite loci (Fig. 5c). The Northern populations showed the lowest values of allelic richness in comparison with other groups (Fig. 5c, Table 1). All obtained results for allelic richness comparisons between regions were not statistically significant, probably due to insufficient sample size (Table 2).
For COI marker, population differentiation (F ST ) is equally strong among the Western (F ST = 0.42265) and the Central (F ST = 0.38347) groups, but it is much higher in the East (F ST = 0.50369) and much lower in the North (F ST = 0.28689). In the East, this is compatible with a scenario of many, independent anthropogenic introductions from diverse sources.
Looking at the spatial population structure, it is remarkable that the overall relation between population pairwise F ST and the geographic distance is relatively low (Mantel's r = 0.130). The picture changed, however, when looking at the presumed ancestral regions (West and Central groups, Mantel's r = 0.237) and the North (Mantel's r = 0.322) separately. In the East, this relation breaks down (Mantel's r = 0.05). The more sophisticated model selection between a discreet and a spatially continuous population model as performed in Pfenninger et al. (2014) was tried, but could not be successfully repeated because of too little information content in the data in conjunction with too many sampling sites and too many individuals. The Monte-Carlo-Markov Chain did not converge for several important parameters despite runs of more than 50 million generations and several successive prior tuning attempts. It should be noted, however, that the root location parameters did converge in all runs and indicated a root location for the COI marker at 49.11°N and 7.56°E (the Saar region).

Phylogenetic analyses
Both phylogenetic trees constructed for COI by maximumlikelihood (ML) and Bayesian inference (BI) produced similar topologies. In the case of Bayesian analysis, phylogenetic trees were prepared in two ways: the first presents clades for particular Arion species which were collapsed regardless of value of sequence divergence (Fig. 6). The Arion species investigated here were represented by 21 strongly supported (1.00 posterior probabilities) clades. All obtained haplotypes for A. vulgaris clustered into one clade. This clade clustered together with sequences defined in GenBank as A. lusitanicus, which probably are in fact sequences of A. vulgaris. Differences in naming are caused by changes in nomenclature proposed by Falkner et al. (2002) of the common pest slug in Europe. The phylogenetic analysis confirm that A. vulgaris is the closest relative of A. rufus and A. ater. Some sequences of A. rufus (AY987900, AY987903) clustered with the A. ater clade; the same pattern was observed for one A. ater sequence (JF950513) which clustered with the A. rufus clade. The A. subfuscus clade shows the highest intraspecific variability within the Arion genus. Additionally, the Bayesian phylogenetic tree has been shown to represent collapsed clades with    . 6 A COI-based Bayesian phylogenetic tree for Arion species constructed from currently available sequences with Arianta arbustorum sequence as an outgroup (see also up to 3% sequence divergence, with the exception of the A. vulgaris clade (Fig. 1, Supplementary Materials). For some Arion species, more than one clade has been specified (A. ater, A. subfuscus, A. flagellus), with the genetic distance between sequences of these species higher than 3%. In total, 28 clades, strongly supported in most cases (1.00 posterior probabilities), have been distinguished for these Arion species. The maximum likelihood phylogenetic tree showed a similar pattern of clade distribution, but was poorly supported by bootstrap values (Fig. 2, Supplementary Materials).

Discussion
Our study comprehensively summarizes the genetic structure and variability within and between European populations of A. vulgaris. Our results confirm a possible origin of A. vulgaris in Western and Central Europe (i.e. France, Benelux and Western Germany). Moreover, the data suggest that A. vulgaris' expansion eastwards may be as a consequence of anthropogenic activity after 1989, when the iron curtain in Europe limiting the transport of goods ceased to exist. The haplotype distribution of A. vulgaris shows clear spatial patterning of diversity at both loci. There are populations with zero diversity all over the range, probably caused by a pronounced meta-population structure, where new suitable habitat patches are colonised by a low number of individuals. Among all COI haplotypes described in our study, haplotype 3 was the most common, being present in 19 out of 26 countries.
Arion vulgaris has been detected from the area of the Iberian Peninsula by Zemanova et al. (2016) but from only three sites located near the border between France and Spain where single individuals were collected. The Iberian Peninsula is unlikely to be the place of the origin of A. vulgaris due to the low densities recorded in this area (Pfenninger et al. 2014;Zemanova et al. 2016). In fact, in this area, the endemic A. lusitanicus occurred. According to other studies, A. vulgaris was observed along the Pyrenees (Quinteiro et al. 2005;Hatteland et al. 2015). Quinteiro et al. (2005) discovered that specimens designated as A. lusitanicus clustered into two different clades based on mitochondrial DNA; the first one closely related to A. ater and A. rufus, and the second to A. nobrei and A. fuligineus. Slugs from the first clade which seem to be A. vulgaris have been collected from Italy (Alpi Carniche, Rivolato), France (Mont. Noire, Central Massif) and Spain (Girona), while A. lusitanicus was recorded from Serra da Arràbida in Portugal. These findings suggest that A. vulgaris and A. lusitanicus are two separate species, both present on the Iberian Peninsula (Quinteiro et al. 2005).
Our analysis revealed the presence of 11 distinct COI haplotypes (45.83% of all haplotypes found in the study) in Polish populations. None of them was unique to Poland, and were also found in other parts of Europe (Table 3). Arion vulgaris was recorded in Poland for the first time in 1993 in Podkarpacie, in areas near Albigowa and Markowa to the east of Rzeszów, but it was probably present previously, in the early 1980s Kornobis 1994, 1995). In Poland, the most common haplotypes were 1 and 3, and none of them appeared in Rzeszów, which is located in close proximity to the first area where A. vulgaris was observed. Analyses of haplotype distribution in Poland suggest that A. vulgaris has been introduced multiple times from different parts of Europe, as also suggested by Soroka et al. (2009) based on mtDNA sequences.
In Norway, A. vulgaris was observed for the first time in 1988 in three localities: Langesund (south of Norway), Krakerøy (southeast) and Molde (southwest) (von Proschwitz and Winge 1994). It was probably introduced as a consequence of anthropogenic activity, with shipments from the Netherlands (Hatteland et al. 2013). A. vulgaris became a serious pest in many parts of southern Norway during the 1990s. Observations in Sweden showed that by increasing in population density, introduced A. vulgaris may outcompete native A. ater (von Proschwitz 1997). The same observation has been reported from Surrey in the UK, where a localized population of A. vulgaris has become established (Davies 1987). Hatteland et al. (2013) studied the distribution of this pest slug in Norway and discovered that A. vulgaris is common in large parts of southern Norway and currently occurs as far north as Troms county, north of the Arctic Circle (Hatteland et al. 2013). In Norway, we have detected eight distinct COI haplotypes, the most frequent being haplotype 1, which is also present in southern Europe. In Molde, one of the localities returning an early Norwegian record of A. vulgaris, haplotypes 1 and 3 were present.
Based on the comparison of haplotype diversity among species, it was suggested that A. vulgaris, with the lowest diversity was invasive, as compared with A. ater and A. rufus (Zemanova et al. 2016). The degree of genetic variation at a locus within a species depends on the complex interplay of mutation rate, migration rate, effective population size, population structure, age of the species, its demographic history, and acting selective forces. To deduce one particular cause from an interspecific comparison of intraspecific genetic variability is therefore not sufficient. Despite the fact that a high number of haplotypes have been detected in A. vulgaris populations, genetic differentiation is relatively low and can be considered as intraspecific genetic variation (Watanabe and Chiba 2001;Sánchez et al. 2011).
Former studies investigating the invasive status and putative origin were either limited by the population samples (Pfenninger et al. 2014) or by the number of loci (ND1 and COI) (Zemanova et al. 2016), which are both mitochondrial and as such segregate as a single locus as well as the choice of analysis (mtDNA mismatch distribution analysis) (Zemanova et al. 2016), which does not seem to be suitable for inference of population historic events (Schenekar and Weiss 2011;Grant 2015), and putative boom-bust cycles and metapopulation structure which may account for the low genetic diversity in some German and French populations have not been taken into account in previous analyses (ABC by Zemanova et al. 2016).
In conclusion, our comprehensive data set plus the two previous studies (Pfenninger et al. 2014;Zemanova et al. 2016) conformed well with the previous assumption of a long-lasting presence in the regions of West and Central Europe, i.e. France and Western Germany (Pfenninger et al. 2014), while in Eastern Europe, the spatial pattern of genetic diversity is compatible with the welldocumented recent, anthropogenically driven invasion (or perhaps better expansion, as a region immediately adjacent to the native range was colonized). Unfortunately, we had no data from the former GDR available. It is tempting to assume that the expansion into the East is related to the fall of the iron curtain in the early 1990s with a subsequent enormous increase in trade between Western and Eastern Europe.
Phylogenies showed that all detected haplotypes of A. vulgaris clustered in one clade, together with sequences designated as A. lusitanicus deposited in GenBank. Both Bayesian and ML phylogenetic trees confirm that A. vulgaris is the closest relative of A. rufus and A. ater, a finding reminiscent of earlier studies (Quinteiro et al. 2005;Hatteland et al. 2015). Distinction of separate clades among A. subfuscus and A. flagellus indicate that observed high genetic variability within particular taxa may result from incorrect identification or presence of cryptic species complexes (Pinceel et al. 2004;Pfenninger and Schwenk 2007;Jordaens et al. 2010).
Our work suggests that A. vulgaris may be native to France and Western Germany, based on the genetic diversity patterns of both mitochondrial and nuclear loci. Thus, A. vulgaris is probably alien/invasive in other parts of Europe where it occurs. Lower genetic diversity in A. vulgaris was observed in Eastern Europe, which definitely justifies the classification as an alien species in this region. However, if the species has a wider native distribution in Central Europe, its biological traits have changed drastically in the last decades to explain why it has been recorded as a major pest in the last 40-50 years, as suggested by Cameron (2016). Future studies need to explore further the understanding of this and other invasive species bridging population genetic studies with other species traits such as reproductive output, food preferences and interactions with natural enemies.
Acknowledgements The first author of this work thanks Daniel Stec (Jagiellonian University) for his support and helpful comments during writing of this manuscript.
Funding information This work was supported by Senckenberg Biodiversity and Climate Research Centre (BIK-F) in Frankfurt am Main (Germany), grants from the Polish-Norwegian Scientific Cooperation (Pol-Nor/201888/77) and Jagiellonian University (K/ZDS/ 006320).
Data availability The molecular genetic datasets generated during and/or analysed during the current study are available in the GenBank repository, [https://www.ncbi.nlm.nih.gov/genbank/].

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.