Out of the Black Sea: Phylogeography of the Invasive Killer Shrimp Dikerogammarus villosus across Europe

The amphipod Dikerogammarus villosus has colonized most of the European main inland water bodies in less than 20 years, having deteriorating effect on the local benthic communities. Our aim was to reveal the species phylogeography in the native Black Sea area, to define the source populations for the colonization routes in continental Europe and for the newly established UK populations. We tested for the loss of genetic diversity between source and invasive populations as well as along invasion route. We tested also for isolation by distance. Thirty three native and invasive populations were genotyped for mtDNA (COI, 16S) and seven polymorphic nuclear microsatellites to assess cryptic diversity (presence of deeply divergent lineages), historical demography, level of diversity within lineage (e.g., number of alleles), and population structure. A wide range of methods was used, including minimum spanning network, molecular clock, Bayesian clustering and Mantel test. Our results identified that sea level and salinity changes during Pleistocene impacted the species phylogeography in the Black Sea native region with four differentiated populations inhabiting, respectively, the Dnieper, Dniester, Danube deltas and Durungol liman. The invasion of continental Europe is associated with two sources, i.e., the Danube and Dnieper deltas, which gave origin to two independent invasion routes (Western and Eastern) for which no loss of diversity and no isolation by distance were observed. The UK population has originated in the Western Route and, despite very recent colonization, no drastic loss of diversity was observed. The results show that the invasion of the killer shrimp is not associated with the costs of loosing genetic diversity, which may contribute to the success of this invader in the newly colonized areas. Additionally, while it has not yet occurred, it might be expected that future interbreeding between the genetically diversified populations from two independent invasion routes will potentially even enhance this success.


Introduction
Biological invasions are the inherent symptom of global changes and a major threat to biodiversity [1][2][3]. Alien species may cause irreversible changes to invaded ecosystems, often resulting in reducing distribution or in extinction of native species through direct predation [4], food and shelters competition [5,6], transmission of parasites [7], or modifications of habitat [8].
Occurrence of the species in its native area is associated predominantly with brackish lagoons (limans) and lower reaches of large rivers draining to the Black Sea. Its phylogeographic history in the native area is unknown, although it could bring key information for understanding current invasion dynamic, as in the case of other Ponto-Caspian intruders [46,47]. The distribution in the invaded continental Europe was a subject of numerous studies and is well documented (summarised in Rewicz et al. [48]). Based on the distribution pattern, two major distinct routes for the invasion of D. villosus have been proposed [49]. The eastern route would encompass Dnieper, Prypiat, Bug and Vistula rivers (Fig. 1). The western one would be composed of the Danube, Rhine, main French rivers, but also some northern sites of central Europe, such as the Mittelland Canal in Germany, and Oder in Poland (Fig. 1). However, the existence of these two distinct routes has not been firmly tested, and numerous points are still subjects of a debate.
First, the origin and genetic diversity of populations found in northern and central Europe (i.e. Oder River and Mittelland Canal) is not clear. Colonization of these sites could be a westward expansion of the eastern route, from the Bug and Vistula rivers, based on the fact that other Ponto-Caspian invasive species followed such direction [49]. However, D. villosus has not been found in the waterways joining the Vistula and the Oder and the dates of first records of D. villosus from both river systems suggest that it was present in the Oder prior to colonization of the Vistula and Bug rivers [49]. Therefore, a secondary eastward extension of the western route was favored to explain D. villosus's presence in northern Germany and western Poland. If correct, there would be two fronts of invasion in Poland (eastern and western), presently 150 km apart and likely to get in contact in the near future [31].
These two fronts might be genetically distinct. First, they originated from different parts of the native area (Fig. 1). Second, genetic differentiation between source populations of Ponto-Caspian species is already known for other invaders, such as mysids [50] and gobies [47].
A second question is the global impact of the invasion process on the genetic variation of D. villosus. So far, only three studies dealt with its genetic diversity. They focused either on molecular identification of D. villosus versus two congeneric species present in the Danube [51,52] or on its invasion dynamics in south-western Europe [53]. The latter study suggested there was no loss of genetic diversity during the invasion process. However, this study was based on few molecular markers, and it is unknown if this pattern restricts to the western route of colonization, or if it is a general pattern for D. villosus invasion. In particular, genetic variation within the eastern route has not yet been explored. Additionally, the genetic variation of the populations recently established in UK [54], as well as their source, are unknown.
In this study, D. villosus populations from the western part of the Black Sea basin (native area), as well as invasive populations from the presumed western and eastern colonization routes and from UK were genotyped for mtDNA (COI, 16S) and seven polymorphic nuclear microsatellites, in order to answer the following questions: (1) What is the species phylogeography in the native Black Sea area, including the assessment of cryptic diversity (presence of deeply divergent lineages), historical demography, level of diversity and genetic differentiation between populations being potential sources for the two presumed colonization routes? (2) Is it possible to associate a distinctive genetic signature to the two presumed colonization routes in continental Europe? (3) Is there a loss of genetic diversity between source populations and colonized areas, and is there a loss of genetic diversity and an isolation by distance along the colonization routes? (4) What was the source(s) for the newly-established UK populations and was this colonization associated with a genetic bottleneck?

Sample collection
Dikerogammarus villosus was collected from 33 sites, both in the native area (hereafter N) and the invaded part of Europe (Fig. 1, Table 1), during expeditions spanning 2002-2012. All the sampling sites were located in public and non-protected areas. No permissions were required for sampling. The study did not involve any endangered or protected species. In the native area, all suitable coastal habitats were surveyed along the western and northern coast of the Black Sea. In the invaded continental Europe, sampling covered both the putative western and eastern routes (hereafter WR and ER, respectively) and the recently invaded UK. One site in UK was sampled twice, in 2010 and 2012.

Testing for cryptic diversity
To visualize molecular divergence of mtDNA haplotypes, a Minimum Spanning Network was generated using ARLEQUIN 3.5.1.2 [60]. Pairwise Kimura 2 parameter (K2p) distances were estimated using MEGA 6.2 [61]. For analysis based on Bayesian inference, the AICM method of moments' estimator [62] was used to define best fitting model of evolution. The time calibrated phylogeny was reconstructed in BEAST, version 1.8.1 [63]. The Hasegawa, Kishino and Yano (HKY) model of evolution with proportion of invariable (I) and Yule speciation model were set for priors. The strict clock with rate 0.0142 proposed for the genus Gammarus was applied for the analyses [64]. Two runs of 20 M iterations of Markov chain Monte Carlo (MCMC) sampled each 1000 iterations were performed. Both runs were examined using Tracer v 1.6, all sampled parameters achieve sufficient sample sizes (ESS>200). Tree files were combined using Log-Combiner 1.8.1 [63], with removal of the non-stationary 10% burn-in phase. The maximum clade credibility tree was generated using TreeAnnotator 1.8.1 [63]. To add additional support for the tree topology, the same dataset was analyzed with Maximum Likelihood (ML) method based on the General Time Reversible (GTR) model [65] with 10000 bootstrap replicates. Model of evolution was selected using JMODELTEST2 [66]. ML analyses were performed in the MEGA 6 [61].

Historical demography within the native range based on mtDNA
To reveal historical demography in the Ponto-Caspian region we used 133 individuals from nine localities (Table 1). In order to assess the temporal changes of the effective population size in each of the three phylogeographic lineages (A-C, see results), a set of the Extended Bayesian Skyline Plot (EBSP) analyses [67] was performed in BEAST, version 1.8.1 [63]. The GTR model of evolution was used as the best fitting model. To ensure convergence, four runs of MCMC, 100M iterations long sampled each 1000 iterations, were performed. Both runs were examined using Tracer v 1.6, all sampled parameters achieved sufficient sample sizes (ESS>200).

Allelic/haplotypic diversity and differentiation
Diversity was assessed by calculating: (1) allelic-haplotypic (msat/mtDNA) diversity (k), (2) allelic richness (A r ) and private allelic richness (PA r ) corrected for a common sampling size using rarefaction approach [68]. Calculations were performed with HP-RARE 1.1 [69], differentiation in A r was tested using the non-parametric Mann-Whitney U-test in STATISTICA 10 [70], and put in brackets if significant. In addition, observed heterozygosity (H O ), expected heterozygosity (H E ) and fixation index (F IS ) were calculated, when appropriate, for microsatellite markers using FSTAT [71]. Pairwise differentiation was determined by two F ST estimators: Y ST with Tamura-Nei distance for mtDNA [72] and Y for microsatellites [73], both implemented in ARLEQUIN, statistical significance being measured using 10000 permutations. Genetic diversity and F ST were assessed either pooling sampling sites, or not, according the hypothesis tested, e.g. between fronts in Poland. Population structure was also analyzed using individual-based Bayesian clustering method implemented in STRUCTURE 2.3.4 [74]. Simulations were performed on the full data set including 29 populations and 876 individuals. Runs for each possible value of K (1 to 8) were repeated 20 times. Each run used a burn-in of 500000 iterations, a run length of 750000 iterations. All simulations were performed using the admixture and correlated allele frequencies models with no prior information. Selection of most probable value of K relied on the ΔK method developed by Evanno et al. [75].

Diversity and differentiation along Western Route (WR)
Based on 20 sites along the WR we tested if microsatellite differentiation increased positively with distance between sites (isolation-by-distance, hereafter IBD) but also if diversity (mean allelic richness) was associated with geographical distance from the source area (Danube delta). The distances were estimated using GOOGLE EARTH v.7.1.2. IBD was tested using Mantel test between F ST / (1-F ST ) and geographic distance as recommended by Rousset [76] for testing IBD in one-dimensional linear systems, with 100000 permutations, using the GENEPOP on the Web 4.2 [77] and ISOLDE software.

Phylogeography in the native Black Sea area
Out of 133 individuals from 9 sites in the native Black Sea region, a total of 17 haplotypes were identified based on concatenated (303+654 bp) 16S and COI mtDNA sequences (S1 Table). First, the difference observed between the most divergent haplotypes was only five nucleotides ( Fig. 2). Second, the mean overall K2p genetic distance between haplotypes was very low being 0.0009 (SD 0.0004). It showed clearly that there is no cryptic diversity involving highly divergent lineages. However, combination of the haplotype network (Fig. 2) and Bayesian phylogenetic reconstruction (Fig. 3A) revealed that the haplotypes may be grouped into three phylogenetic lineages. Their spatial distribution is partly structured geographically. Lineage A includes eight haplotypes (5)(6)(7)(8)(9)(10)(11)(12) specific to the Durungol liman in Turkey (2-N). Lineage B includes 7 haplotypes (1,4,13,(15)(16)(17) and lineage C includes the 3 remaining haplotypes (2, 3 and 14). The Bayesian chronogram showed that C diverged from A+B ca 280 kyr BP, while A and B split ca. 200 kyr BP (Fig. 3A). The results of the EBSP analysis indicated that population of D. villosus in the Durungol liman experienced steady growth for the last 20 ky, while populations of the remaining two lineages in the native area remained stable for most of the last 30k years, with accelerated growth starting less than 10 ky ago (Fig. 3B).
The highest diversity for mtDNA and msat was observed in Durungol (A r = 5.3 and 5.73 respectively); the locality harboring also the highest private allelic richness (PA r = 5.3 and 1.08) ( Table 2). The potential (i.e. a priori) source areas of invasion, i.e. Dnieper and Danube deltas,  Fig. 2). For msat these areas had similar diversities (A r = 3.94 and 3.57) and low private allelic richness (PA r = 0.05 and 0). Genetic differentiation differed from zero for all area pairwise comparisons, for both msat and mtDNA data (Table 3). However, the level of differentiation was heterogeneous, being the highest between Durungol and Dnieper (F ST = 0.180 and Y ST = 0.693) and the lowest between Danube and Dniester (F ST = 0.048 and Y ST = 0.103). The results of Bayesian clustering suggest that the four selected areas may represent four genetic clusters, although the division is not strict. The Durungol and Dnieper populations are the most homogeneous ones, while the Dniester and Danube populations show symptoms of migration or very recent common ancestry (Fig. 4A,B).

Colonization dynamics in Continental Europe
For mtDNA, ER and WR in Poland, i.e. fronts, are not differentiated from their respective putative sources in the native region, i.e. Dnieper (3+4+5) and Danube (6+7+8+9), respectively (Table 3), while differentiation between fronts was significant (Table 3). For msat, although all pairwise comparisons for differentiation were significant (Table 3), the level of differentiation between fronts in Poland and their respective putative areas of origin was low (F ST = 0.027 and = 0.019) compared to differentiation among sites belonging to different routes (0.157 < F ST < 0.215) ( Table 3). Bayesian clustering analysis showed clearly, that individuals from the western front (Poland), the western route and the source in Danube from a homogeneous genetic unit, while the eastern front with its putative source in Dnieper form another homogeneous genetic unit (Fig. 4A,B). Both fronts in Poland had the same level of diversity for msat compared to their putative source in the native region i.e. in Dnieper and Danube respectively ( Table 2). This is also true for mtDNA for beginning the end of the ER, but the WR in Poland seems to present less diversity than the Danube (6+7+8+9) ( Table 2). Along the WR, geographic distance from the source area did not explain msat diversity within sites (Fig. 5). Along the 4500 km long route, only one site located on the River Vah in Slovakia (site 16) had lower diversity. No isolation by distance (Fig. 6) was present as no significant correlation was detected between pairwise geographical distances and genetic differentiation (Mantel test, R 2 = 0.0047, P = 0.42).

Source population and diversity for the UK
We observed high genetic differentiation for both mtDNA (0.206 < Y ST < 0.298) and msat (0.138 < F ST < 0.172), for pairwise comparisons between ER in Poland and two pooled sites in the Netherlands (23+24) or any UK site. On the opposite, pairwise comparisons between two sites in the Netherlands and each UK site showed no significant differentiation for mtDNA (Table 4). For msat, lower level of differentiation was observed between the UK site 32 (A and B) and the Netherlands (0.026 < F ST < 0.035) than between the UK site and ER in Poland (0.138 < F ST < 0.172). The UK site 33 showed a less conclusive picture for F ST (Table 4). Bayesian clustering analysis (Fig. 4A,B) showed that the UK populations form a homogeneous genetic unit with the western route, and genetically different from the eastern route.
Only one mtDNA haplotype (haplotype 1), the most common in continental Europe, occurred in UK. Haplotype 4, while common in ER, was absent from the UK. In Grafham Water site (32A, 32B) we observed no loss of diversity for msat compared to the Netherlands with A r values being respectively 3.57 and 3.49, versus 3.8 (Table 2). Diversity in the more recent population (33) was 2.49 but did not differ in the statistical terms from the above A r values ( Table 2).

Discussion
The Ponto-Caspian region has been recognized as the most prominent donor of non-indigenous hydrobionts to Europe and to the North American Great Lakes system. Their taxonomic spectrum is wide including amphipods, mysids, cladocerans, gastropods and fishes [79][80][81][82]. Phylogeography and population genetics patterns of these invaders may help in understanding colonization dynamics and in controlling their further spread [16]. Our results confirm that the invasion scheme for one of these species, D. villosus, is complex, with multiple routes, and that the loss of genetic diversity during the course of colonization is weak. We evidenced that   Phylogeography and contemporary genetic structure in the native region Level of cryptic diversity for invertebrates in the Ponto-Caspian region is highly variable. Deep level of divergence, but below the species threshold, was observed for e.g. cladocerans [78] and mysids [79]. Cryptic diversity was revealed e.g. in monkey goby [80]. Although the phenomenon is known to occur in several amphipods [81,82], no cryptic diversity was detected in D. villosus in its native region. The 17 mtDNA haplotypes showed shallow divergence with an overall K2p genetic distance of 0.0009, far below the threshold of 0.03-0.055 identified between crustacean species [83]. Such low divergence may be related to very recent history of the species within the Black Sea region (see below).
Although shallow, the divergence between mtDNA haplotypes is geographically structured with lineage A including a set of 8 haplotypes (out of 17) and being restricted to the Durungol liman. This pool of haplotypes separated from others ca 200 kyr BP. In addition, the other sampled areas are also characterized by high private haplotypic richness. Overall, the level of differentiation (Y ST ) between populations is high in the Black Sea area. The turbulent Pleistocene hydro-geological history of the region with recurrent changes of sea level and salinity may be among the most powerful driving forces explaining this pattern [46,50,78,79,[84][85][86]. During the last 670 kyr, there were at least 12 significant saline water intrusions from the Mediterranean Sea, and eight intrusions from the Caspian Lake to the Black Sea [87]. These events caused water level fluctuations and substantial salinity shifts from nearly fresh to full marine conditions that could cause shift ranges and population fragmentations in oligohaline hydrobionts inhabiting this basin [88]. During fully marine salinity stages, slightly brackish estuaries and limans may have become isolated refugia and differentiation centers for local aquatic fauna. The dating of divergence between Durungol liman and others sites coincides with one of the most prominent salinity raises [87]. In other sites presence of shared haplotypes reflects probably both recent and historical migration events among various areas in the native region. However, the overall presence of private haplotypes and high differentiation level indicate possible founder effects at the time of colonization. It could be followed either by restricted gene flow (which is confirmed by msat results) or even by allopatric divergence during stages of raised salinity. The results of EBSP analyses support very recent post-Pleistocene demographic  [50,79,84], cladocerans [78] and gammarids [46,78,86].
The msat analyses pointed out high level of differentiation (F ST ) between the four native sampled areas, that were divided in four genetic clusters in Bayesian analyses. If the Durungol liman is clearly isolated from the other sampled areas, the latter show connectivity as pointed out by unclear Bayesian assignment of some individuals. One of the explanations can be a shorter geographic distance between these populations. The north-western Black Sea is also the shallowest and least saline, due to massive sedimentation and inflow of riverine waters from e. g. Dnieper and Danube rivers [89,90]. Chances for migrations between these populations are high, including anthropogenic transport due to the high ship traffic between local ports.

Invasion routes and dynamics in continental Europe
In our study, combination of mtDNA and msat analyses clearly identified Danube and Dnieper deltas as differentiated sources for the two invasion routes we named "Western" and "Eastern", respectively. Dikerogammarus villosus has been highly monitored throughout Europe due to its detrimental impact on the ecosystem. Therefore, accurate map of invasion progress can easily be drawn and converted into the most likely scenario for colonization routes [33,[89][90][91][92]. Agreement of our results with putative routes might seem trivial at first sight. However, few studies upon other species pointed out that molecular data identified routes that were different from the most likely, census-based, scenarios [9,93]. Based on geographic invasion patterns of several aquatic species, Bij de Vaate et al. [91] defined three invasion corridors from the Ponto-Caspian region into continental Europe i.e. the northern (Volga River, Beloye, Onega and Ladoga lakes, Neva River to the Baltic Sea), the central (Dnieper, Pripyat, Pripyat-Bug channel, Vistula, Oder, Mittelland canal) and the southern (Danube, Rhine) one. Numerous species invasions fitted this pattern [49,94]. Contrary to other species, D. villosus used only the eastern part of the central corridor and has not passed the Bydgoski channel in Poland which is connecting the Vistula and the Oder rivers (Fig. 1). On the other hand, the western part of the central corridor was colonized eastwards by population which came up the entire southern corridor westward. Possibly, the Bydgoski channel, with its prominently soft bottom, slow current and abundant vegetation is not prone to be colonized by D. villosus [95]. However, the closely related invader, D. haemobaphes, along with some other invasive amphipods, such as Echinogammarus ischnus and Chelicorophium curvispinum, managed to pass the Bydgoski channel and use the entire central corridor westwards. This channel was an important shipping route until mid-20th century, but the traffic now is heavily limited [96]. Yet, we cannot exclude possibility of future contact between these two distinct populations of D. villosus. The two fronts in Poland are differentiated and characterized by level of diversity analogous to their source regions. If the two fronts meet, hybridization will occur as the two populations are not phylogenetically and ecologically divergent which implies the absence of reproduction barrier. This may result in producing a potential "super-hybrid"an even more effective invader, as it was observed in other cases [25,26]. Thus, the situation deserves particular surveillance and management to avoid contact between these two fronts.
The WR has a length of about 4500 km from the source population in the Black Sea region to the invasion front in Poland. All mtDNA haplotypes found in the native range were observed in the invaded area, 3 out of 8 non-frequent haplotypes being present in the last 1000 km. In addition, in the Mittelland Canal we found one haplotype (haplotype 18) not even encountered in the native area, probably due to its very low frequency. For microsatellites, we found no loss of diversity along the route and no isolation by distance. Globally, it suggests that no bottlenecks occurred along WR. Similar conclusions were made by Müller et al. [51] and Wattier et al. [53], who conducted research on a smaller scale or with fewer genetic markers. Even if reduction of diversity in the invaded areas was often expected in the literature in 20th century, numerous studies since then have shown it might be far from being the rule [97][98][99]. Lack of diversity loss may result from a very large propagule size i.e. large founding population, and/or propagule frequency i.e. recurrent waves of invaders. We suspect the latter to play an important role for D. villosus. Indeed, no loss of diversity is observed while genetic differentiation between sites is present. This suggests that recurrent waves are both maintaining allelic diversity at a high level, and reshuffling allelic frequencies what generates differentiation.
Only one site (site 16) in the middle section of River Vah, was characterized by both low diversity and high level of differentiation from other sites. Strong founder effect is likely to explain this pattern as, probably, the site was sampled very recently after D. villosus first colonization [48,100].

Recent UK overseas conquest
The overseas introduction of D. villosus into UK in 2010 was noticed in popular media [54]. It proved clearly that large funds spent on biosecurity programs to prevent the spread of invasive species (i.e. the procedure "check, clean and dry") have been insufficient to stop the killer shrimp [101,102]. Both F ST and Bayesian clustering of msat allowed us to exclude the populations from ER as donors for the UK sites. In addition, the haplotype 4, both frequent and private to ER, was not detected in UK, although our sampling size was limited in this case.
Diversity in any UK population was not different from the continent and no bottleneck effect was observed. Apparently, the propagule pressure was high enough to alleviate diversity loss. We are not able to conclude whether the introduction to UK was a single event followed by secondary colonization or multiple introductions. Anyway, the killer shrimp is spreading very efficiently throughout the UK. Furthermore, another congeneric invader, the "demon shrimp", D. haemobaphes, has already been recorded in UK [103]. Based on several possible expansion models of both D. villosus and D. haemobaphes it has been estimated that more than 60% of the UK waterbodies is suitable and vulnerable to colonisation by these two invaders [101,[104][105][106][107]. Moreover, high popularity of water sports may further accelerate the invasion [108], due to high ability of the killer shrimp to spread via boating and diving equipment [32].

Conclusions
Our results identified impact of the Pleistocene sea level and salinity fluctuations on the phylogeographic structure of Dikerogammarus villosus in the Black Sea native region and presence of two differentiated source populations, i.e. the Danube and Dnieper deltas. These sources are associated with two independent invasion routes (Western and Eastern) in continental Europe for which no loss of diversity is observed. We can expect further spread of the killer shrimp in continental Europe, even in smaller tributaries. MacNeil & Platvoet [95] pointed out that solid objects, like concrete fish passages, could be used by D. villosus as mainstay in smaller tributaries. This may pose a threat to native gammarids occupying such refuges [109].
The UK population has probably originated in the Western Route and despite very recent colonization, no drastic loss of diversity was observed. This recent overseas conquest provides rather non-optimistic message, accounting that the UK authorities implemented preventive biosecurity protocols and risk assessments. The Great Lakes of North America are likely to be the next step, since other Ponto-Caspian invertebrates already managed to reach them [110,111].
Finally, the Dniester native area is characterized by high msat allelic diversity (including private alleles). Thus, even if at the moment it is not a source population for the colonization of Europe, it may act as a donor for other source areas if the anthropogenic transport increases, enhancing the local genetic diversity.