Ponto-Caspian amphipods have high invasion potential, and as many as 13 morphospecies originating from that region have been recorded as alien or invasive elsewhere (Holdich and Pöckl 2007). This number includes three species of the genus Dikerogammarus (Ricciardi and MacIsaac 2000; Cristescu et al. 2004; Berezina 2007). One of them, D. haemobaphes (nicknamed the demon shrimp), is native to the lower courses of large rivers in the Black and the Caspian drainage areas, to their brackish water lagoons as well as to the Caspian Sea (Sars 1894; Carauşu 1943) (Fig. 1). Already Mordukhai-Boltovskoi (1964) forecasted possible quick range expansion of the demon shrimp in Western Europe, but until 1980 the species was reported only from the lower and middle Dnieper, Dniester, Don, and Volga rivers (Jażdżewski 1980). In 1976, D. haemobaphes was found in the upper Danube in Germany (Bij de Vaate et al. 2002). In the 1990s, its rapid expansion in Germany had begun. In 1993, it was observed in the Main-Danube canal (Schleuter et al. 1994) and already in 1994 in the upper Rhine River (Schöll et al. 1995). At the same time, in Russia, it spread up the Volga River until Rybinsky Reservoir and invaded the upper Moskva River (Berezina 2007). Dikerogammarus haemobaphes was first recorded in the Vistula River in Poland in 1997, and, during the monitoring of its middle and lower course in 1998 and 1999, it appeared to have established populations there (Konopacka 1998; Bij de Vaate et al. 2002). At the end of the 20^th century, the species was commonly found also in the lower and middle Oder River, Vistula Lagoon, and in the Masurian Lake District in Poland (Jażdżewski 2003; Grabowski et al. 2007b). In 2008, the species was recorded for the first time in the Meuse, Moselle, and Seine rivers in France (Labat et al. 2011). In 2011, it was found in two isolated Alpine lakes in Austria and Germany (present study). In Switzerland the species was reported from the Murtensee in 2011 and from the Neuenburgersee in 2016 (Altermatt et al. 2019). The latest record of this species outside its native range came from 2012 when it was found in the River Severn in Great Britain and, subsequently, it has established sustainable populations in the central part of the country (Aldrige 2013; Etxabe et al. 2015; Constable and Birkby 2016; Johns et al. 2018). As a result, the species is presently distributed in the majority of large European rivers, and the history of its range extension suggests that it has used all three invasion corridors proposed by Bij de Vaate et al. (2002) (Fig. 1).

The demon shrimp has a high potential for invasion (Grabowski et al. 2007a; Bacela et al. 2009; Bacela-Spychalska and Van der Velde 2013). However, another invasive species of the genus Dikerogammarus, the so-called killer shrimp D. villosus, was often considered a more successful invader (Dick et al. 2002; Rewicz et al. 2015; Kobak et al. 2016), attracting more scientific attention than D. haemobaphes. Only two molecular studies on D. haemobaphes have been reported so far by Müller et al. (2002) and Cristescu and Hebert (2005). However, the number of localities and individuals sampled in both studies was too low to track the invasion history of the demon shrimp. Recent studies of the congeneric D. villosus revealed the absence of cryptic species and no genetic diversity loss along invasion corridors (Wattier et al. 2007; Rewicz et al. 2015).

Based on the thorough sampling across Europe and using two mitochondrial and one nuclear marker, our study aims to reveal the phylogeographic structure and historical population dynamics of D. haemobaphes in its native Black Sea basin and in the invaded range. Taking into account what is currently known about the recently studied and closely related killer shrimp, the history of the demon shrimp invasion and the geological history of the native area, we hypothesized that (i) it is plausible that in the native region the demon shrimp encompasses several weakly divergent lineages dating back to Pleistocene, (ii) there are possibly two sources, the Dnieper and the Danube deltas and, hence, two independent routes for the species invasion to Central and Western Europe – the so-called Central Corridor and the Southern Corridor of invasion, respectively, (iii) learning from the invasion history of the closely related killer shrimp, we can expect that the populations of the demon shrimp will also show signs of demographic and spatial expansion and no loss of genetic diversity, except the UK population, in comparison to the source population.

Dikerogammarus haemobaphes was collected at 55 sites in both its native (23 sites) and invaded (32 sites) ranges between 2000 and 2014 (Table 1). Some localities were geographically very close (ca. 20 km distance) to each other (same site number in Table 1, but one of the two or three indicated with letter A or B), so we decided to combine them with the closest sampling site to simplify further analysis (Table 1). As a result, our dataset comprises 41 localities, of which 19 belonged to the native and 22 to the invaded range (Fig. 2). All the sampling sites were located in public and non-protected areas. In the native area, sampling was done along the western and northern coast of the Black Sea either in limans of large coastal rivers or lower courses of large rivers (e.g., Danube, Dniester, Dnieper). The area around the Azov Sea was included in the sampling campaign, but the demon shrimp was not found in any of the 15 sampling sites visited. The Caspian Sea was not included in our study. The sampling in the invaded areas covered all three invasion corridors used by aquatic species to spread in Central, Western, and Northern Europe (see Bij de Vaate et al. 2002) including the recently colonized waters in the UK.

Figure 2.

A geographic distribution of D. haemobaphes haplotypes in the native and invaded area. Numbers in triangles (native range) or circles (invaded range) represent sampling localities coded as in Table 1. Sites with lineage A indicated with purple. n = number of sequences analyzed in the site B the minimum spanning network of mtDNA haplotypes identified in D. haemobaphes. The size of circle denotes the frequency of each haplotype. The numbers correspond to the haplotype numbers. GDU – geo-demographic unit.

The total DNA was extracted from 434 individuals according to the standard phenol-chloroform method (for details see Hillis et al. 1996). Air-dried DNA pellets were eluted in 100 μl of TE buffer, pH 8.0, stored at 4 °C until amplification, and subsequently at −20 °C for long-term storage. Two mtDNA markers: a gene for the 16S ribosomal RNA (16S rRNA; ca. 320 bp fragment) and the standard barcoding fragment of the cytochrome c oxidase subunit I gene (COI; 658 bp fragment) were amplified. LR-J-GAM/LR-N-GAM primers (Müller et al. 2002) and reaction conditions after Grabowski et al. (2012) were used for 16S rRNA amplification. COI gene was amplified using LCO1490/HCO2198 (Folmer et al. 1994) and UCOIR/UCOIF (Costa et al. 2009) primers and reaction conditions following Hou et al. (2007). Sequences were obtained using BigDye sequencing protocol on the Applied Biosystems 3730×l capillary sequencer by Macrogen Inc., Korea. Sequencing of the COI gene was performed unidirectionally (with the forward primer). In cases when the quality of the obtained sequence was not sufficient, a reverse sequencing was applied as well. Sequencing of the 16S gene was bidirectional. Sequences were edited, aligned with the ClustalW algorithm (Chenna et al. 2003) using BioEdit 7.2.5, and trimmed to the length of the shortest one. The resulting 433 sequences of 16S (302 bp) and COI (598 bp) were subsequently concatenated for the purpose of the analyses.

Haplotypes were retrieved using DNA SP v5 both for individual markers and for the concatenated sequences (Librado and Rozas 2009). Then, at least two individuals of lineage A as well as each detected geo-demographic unit (GDU, defined based on the spatial distribution of well-supported lineages on the chronogram, see Results for details) were amplified for the additional nuclear marker, 28S rRNA, for phylogeny reconstruction. The nuclear marker was amplified with 28F and 28R primers and reaction conditions published by Hou et al. (2007).

Relevant voucher information, taxonomic classifications, and the COI barcode sequences are publicly accessible through the public data set “DHAEMOBA” (https://doi.org/10.5883/DS-DHAEMOBA) on the Barcode of Life Data Systems (BOLD; www.boldsystems.org) (Ratnasingham and Hebert 2007). Newly generated sequences were also deposited in GenBank: COI: MN342874–MN343144, MN343146–MN343307; 16S: MN343308–MN343578, MN343580–MN343741; 28S: MN343743–MN343753.

Bayesian time-calibrated phylogeny was reconstructed in BEAST, version 2.5.2 (Bouckaert et al. 2014) to infer the time frame of the D. haemobaphes diversification. The COI and 16S mitochondrial markers were used as separate partitions in the Bayesian inference analysis (BI). The priors to evolutionary models were set using bModelTest (Bouckaert and Drummond 2017). The strict clock was calibrated with the COI mutation rate set at 0.01773 substitutions/site/Ma as proposed for amphipods by Copilaş-Ciocianu et al. (2019), given that very similar rates were reported for the family Gammaridae in other studies (Mamos et al. 2016; Grabowski et al. 2017). Four runs of Markov chain Monte Carlo (MCMC), each 20 million iterations long and sampled every 1000 iterations, were performed. Runs were examined using Tracer v 1.6, and all the sampled parameters achieved a sufficient sample size (ESS > 200). Tree files were combined using Log-Combiner 1.8.1 (Drummond et al. 2012), with the removal of the non-stationary 25% burn-in phase. The maximum clade credibility chronogram was generated using TreeAnnotator 2.5.2 (Bouckaert et al. 2014). For additional support of tree topology, the concatenated COI + 16S dataset was analyzed with the Maximum Likelihood (ML) method, using the Tamura 3-parameter model (Tamura 1992) selected through the Bayesian Information Criterion (BIC) with 10000 bootstrap replicates. ML analysis was performed in MEGA X (Kumar et al. 2018).

To visualize molecular divergence of mtDNA haplotypes, a Minimum Spanning Network was generated using Arlequin 3.5.1.2 (Excoffier and Lischer 2010). Pairwise Kimura 2-parameter (K2P) distances were estimated using MEGA X (Kumar et al. 2018).

Four molecular species delimitation methods were applied to reveal the potential Molecular Operational Taxonomic Units (MOTUs) that could represent putative cryptic species within the studied demon shrimp populations. Two methods were distance-based: Barcode Index Number (BIN) System (Ratnasingham and Hebert 2013), and the barcode-gap approach using the Automatic Barcode Gap Discovery (ABGD) (Puillandre et al. 2012). The following two were tree-based, a phylogenetic approach using Generalized Mixed Yule Coalescent (GMYC) model-based method (Pons et al. 2006), according to Monaghan et al. (2009) and the Bayesian implementation of the Poisson Tree Processes (bPTP) (Zhang et al. 2013). For both tree-based methods, we have used the already generated Bayesian time-calibrated phylogeny.

The BIN method is implemented as part of the Barcode of Life Data system (BOLD; Ratnasingham and Hebert 2007). Newly submitted sequences are compared together with sequences already available in BOLD. Sequences are clustered according to their molecular divergence using algorithms that aim at finding discontinuities between clusters. Each cluster is ascribed a globally unique and specific identifier (aka Barcode Index Number or BIN), already available or newly created if the submitted sequences do not cluster with previously known BINs. Each BIN is registered in BOLD.

The ABGD method is based upon pairwise distance measures. With this method, the sequences are partitioned into groups (MOTUs), such that the distance between two sequences from two different groups will always be larger than a given threshold distance (i.e., barcode gap). We used primary partitions as a principle for group definition, as they are typically stable on a broader range of prior values, minimize the number of false-positives (over-split species), and are usually close to the number of taxa described by taxonomists (Puillandre et al. 2012). The default value of 0.001 was used as the minimum intraspecific distance. As there is currently no consensus about which maximum intraspecific distance is reflecting delimitation of species, neither based on morphology (Costa et al. 2007; Weiss et al. 2014; Katouzian et al. 2016) nor on reproductive barrier (Lagrue et al. 2014) we explored a set of values up to 0.1. The standard Kimura 2-parameter (K2P) model correction was applied (Hebert et al. 2003).

The GMYC method defines MOTUs through identification of the switch from intraspecific branching patterns (coalescent) to interspecific species branching patterns (Yule process) on a phylogenetic tree. First, a log-likelihood ratio test is performed to assess if the GMYC model fits the observed data significantly better than the null model of a single coalescent species. If there is evidence for overlooked species inside the phylogenetic tree, the threshold model is tested for the observed data to estimate the boundary between intra- and interspecific branching patterns. The Bayesian tree was uploaded into the R (R Core Team, 2013) software package ‘SPLITS’ (Species Limits by Threshold Statistics) (Ezard et al. 2009) and analyzed using the single threshold model.

The bPTP is another phylogeny-based method. The bPTP incorporates the number of substitutions in the model of speciation and assumes that the probability that a substitution gives rise to a speciation event follows a Poisson distribution. The branch lengths of the input tree are supposed to be generated by two independent classes of the Poisson events, one corresponding to speciation and the other to coalescence. Additionally, the bPTP adds Bayesian support values (BS) for the delimited species (Zhang et al. 2013). The analysis was performed on the bPTP webserver (available at http://www.species.h-its.org/ptp/) with 500 000 iterations of MCMC and 10% burn-in.

The 28S rRNA nuclear phylogeny was reconstructed for Dikerogammarus haemobaphes, D. bispinosus Martynov, 1925, and D. villosus using the Maximum Likelihood method and the Tamura 3-parameter model selected through BIC as a best fitting model (Tamura 1992). All positions with less than 95% site coverage were eliminated. The bootstrap test was done with 500 replicates (Saitou and Nei 1987). These analyses were done in MEGA X (Kumar et al. 2018). Sequences of D. villosus (KF478495, KF478496) and Pontogammarus robustoides (Sars, 1894) (KF478447) used as an outgroup were retrieved from GenBank.

To reveal the historical demography of D. haemobaphes, we used the Extended Bayesian Skyline Plot (eBSP) (Heled and Drummond 2008) constructed in BEAST, version 2.5.2 (Bouckaert et al. 2014), using COI and 16S as separate partitions. The eBSP was performed for the lineage A and three out of four geo-demographic units (GDUs) identified in the lineage B (see results). Due to a small number of individuals representing the clade from the Moskva River in Russia (GDU B4), this unit was not explored. The clock rate and model selection were performed the same way as in the case of time-calibrated phylogeny reconstruction. The population scaling factor was set to 0.5. To ensure convergence, four runs of MCMC, each 100 million iterations long and sampled every 5000 iterations, were performed. Runs were examined using Tracer v 1.6; all the sampled parameters achieved sufficient sample sizes (ESS > 200) and presented congruent results. The final figures were generated in R (R Core Team, 2013).

We have also examined the current demographic status of the lineage A and three GDUs of lineage B of D. haemobaphes in Arlequin 3.5 (Excoffier and Lischer 2010) with mismatch distribution, supplemented by the selective neutrality tests, i.e., Tajima’s D (Tajima 1989) and Fu’s Fs (Fu 1997) as indicators of population expansion. We verified the validity of both models (sudden demographic and spatial expansion) with the sum of squared deviations (SSD) between observed and expected mismatches and by the Harpending’s Raggedness statistics (Harp).

The haplotype diversity was assessed by calculating the allelic richness (A_r), and the private allelic richness (PA_r), where haplotypes equaled to alleles, corrected for a common sampling size using a rarefaction approach (Leberg 2002) in case of localities with at least three individuals. Calculations were done using Hp-Rare 1.1 (Kalinowski 2005).

We used 19 sites along the Central Corridor to test for a positive correlation between genetic differentiation and the distance between sites (isolation-by-distance, IBD). Besides, the data from these sites were used to check whether diversity (haplotype diversity, nucleotide diversity, A_r, PA_r, see Table 1, Suppl. material 1: Table S1) was associated with geographical distance from the source area (Dnieprovski Liman). We evaluated the trend in haplotype diversity, nucleotide diversity, allelic richness, and private allelic richness along the Central Corridor of the invasion by calculating Pearson Correlation Coefficient. The distances (Suppl. material 1: Table S1) were estimated using Google Earth v.7.1.2 as a measurement of shortest distance along the waterway (Google Earth, option path, zoomed enough to fit the line in the riverbed) of Central Corridor starting from our sampling locality 11A (Dnieprovski Liman) to locality 32 (Oder). IBD was tested using the Mantel test between F_st/(1-F_st) and geographic distance as recommended by Rousset (1997) for testing IBD in one-dimensional linear systems, with 100000 permutations, using the ISOLDE software embedded in the GenePop on the Web 4.2 package (Raymond and Rousset 1995).

In the dataset composed of 433 sequences, we identified 15 haplotypes of 16S (302 bp) and 27 haplotypes of COI (598 bp). The concatenation of both fragments (900 bp) resulted in the recognition of 39 mtDNA haplotypes derived from different combinations of the 16S and COI haplotypes (Suppl. material 2: Table S2, Suppl. material 3: Table S3). The pairwise K2P distance between the haplotypes ranged from 0.002 to 0.183 for COI and from 0.003 to 0.099 for 16S; the highest values suggesting the presence of cryptic diversity within the studied morpho-species.

All the COI sequences, represented by 27 haplotypes, were uploaded into the Barcode of Life Data System (BOLD), resulting in two Barcode Index Numbers (BINs) obtained (Fig. 3). The BIN AAX9262 was attributed to 406 sequences (23 COI haplotypes) and was already reported in BOLD and GenBank associated with D. haemobaphes. The mean K2P distance within this BIN was 0.0023, and the maximal K2P distance was 0.0263. The BIN ADB9467, associated with 27 sequences (four COI haplotypes), has not been reported before. The mean K2P distance within this BIN is 0.0013, and the maximal K2P distance is 0.0051. The closest BIN to ADB9467 is the AAX9262, with an average K2P distance of 0.176.

Figure 3.

Time-calibrated phylogeny of Dikerogammarus haemobaphes and the results of different species delimitation methods. Maximum clade credibility chronogram was inferred from a strict molecular clock model based on the COI + 16S data set of D. haemobaphes. The numbers given next to the respective main nodes indicate Bayesian posterior probabilities (> 0.5) and ML bootstrap values (> 50%). B1-4 represent the geo-demographic units within the lineage B. BIN – Barcode Index Number (BIN) System, ABGD haplotypes – Automatic Barcode Gap Discovery based on haplotypes, ABGD sequences – Automatic Barcode Gap Discovery based on all sequences, GMYC – Generalized Mixed Yule Coalescent, bPTP – Bayesian Poisson Tree Processes.

Automatic Barcode Gap Discovery (ABGD) was performed on two data sets (i) 27 sequences, each representing one of the 27 COI haplotypes and (ii) 433 sequences, representing all individuals sequenced. The first analysis indicated the existence of two MOTUs (lineages A, B), divergent by the 0.183 maximum K2P distance, that may be interpreted as two potential cryptic species (Fig. 3). The individuals belonging to the lineage A included only four haplotypes (27 sequences) and were found in five localities, exclusively in the native range of Dikerogammarus species (two sites in the lower Danube, two in the lower Dniester and one on the Kerch Peninsula). The lineage B consisted of 23 haplotypes (406 sequences) widely distributed both in the native and the invaded range of D. haemobaphes. The second analysis resulted in partitioning all COI sequences into three groups. One corresponded to the lineage A, while within the lineage B, the sequences from individuals found in the River Moskva in Russia were separated as yet another MOTU.

The GMYC analysis also rejected the null model of the D. haemobaphes sequences representing a single species. The single threshold method of lineage identification allowed to recognize five entities: lineage A and lineage B, the latter further divided into four sub-lineages (Fig. 3). The estimated number of species in bPTP ranged from three to five, the best Bayesian support (>0.987) being for three (Fig. 3).

All the five methods used for species delimitation are congruent at recognizing the lineage A as distinct from any other MOTU. Given the extremely high K2P-distance (max 0.183) between lineage A and all the others, we can conclude the presence of cryptic diversity within D. haemobaphes. In addition, since this lineage is restricted only to the native range, we do not use it in tracking the demon shrimp invasion in Europe. Lineage B was recognized by BINs and, partially, by ABGD as one MOTU, while other methods subdivided it into smaller groups. Given the weak phylogenetic support for these sublineages and lack of congruent delimitation results (Fig. 3) we decided not to treat them as cryptic species, yet taking into account their allopatry (see below), we refer to them further as geo-demographic units (GDU B1–4).

The phylogeny reconstruction based on the 28S rRNA confirmed the separation of D. villosus and D. bispinosus from D. haemobaphes, as well as the clear separation of lineages A and B, further supporting the presence of two cryptic species within D. haemobaphes (Fig. 4). Within the morphospecies D. haemobaphes, the lineage A considerably differed from lineage B and, although with quite a low bootstrap support, showed more affinities to D. bispinosus than to lineage B.

Figure 4.

Maximum Likelihood tree of Dikerogammarus spp. based on the 28S rRNA gene. Sequences representing D. haemobaphes lineage A, each of the geo-demographic units (GDU B1-4) within D. haemobaphes lineage B as well as D. bispinosus and D. villosus were used. The evolutionary history was inferred by using the Maximum Likelihood method and the Tamura 3-parameter model. All positions with less than 95% site coverage were eliminated. Bootstrap test was applied with 500 replicates. Information about the D. bispinosus sequence may be found in public BOLD dataset: DHAEMOBA. Dikerogammarus villosus and Pontogammarus robustoides sequences were retrieved from GenBank.

The Bayesian phylogenetic reconstruction indicated the existence of two main lineages, A and B (Fig. 3). Lineage A is constituted by four haplotypes (36–39). These haplotypes were found in the Dniester river and its liman, the lower course of the Danube river, and in one locality at the Crimean peninsula. Lineage B can be, taking into account the results of MOTU delimitation together with the allopatric distribution of delimited entities, further divided into four geo-demographic units. The geo-demographic unit B1 (GDU B1) includes 12 haplotypes (haplotype numbers 1–12) endemic to Durugöl (Duruşu) liman in Turkey; GDU B2 includes 14 haplotypes (haplotype numbers 13–17 and 24–32) found both in native as well as recently invaded range. GDU B3 is represented by six haplotypes (haplotype numbers 18–23) found only in the Dniester river. The remaining three haplotypes (haplotype numbers 33–35) made up GDU B4 and were found exclusively in the Moskva River in Russia. The Bayesian chronogram showed that lineage A diverged from the rest ca. 5.1 Ma, whereas GDU B4 diverged from GDU B1-3 ca. 750 ka, GDU B1 diverged from GDU B2-3 ca. 300 ka, while the split between GDU B2 and GDU B3 was ca. 250 ka.

The results of eBSP analysis indicated that the size of populations within the lineage A remained stable between 80 and 20 ka, followed by a slight decrease from 20 to 15 ka and a large decrease during the last 10 ka (Fig. 5).

Figure 5.

Multilocus extended Bayesian Skyline Plots for four lineages / geo-demographic units (GDU) of Dikerogammarus haemobaphes. Solid lines indicate the median posterior effective populations size through time; dotted lines indicate the 95% highest posterior density interval for each estimate.

Within the lineage B, the GDU B1 from the Durugöl (Duruşu) liman in Turkey experienced slight growth from 80 ka till 50 ka. Later, the population remained stable. On the contrary, the widespread and invasive clade GDU B2 retained stable size till ca. 2.5 ka, and after that experienced accelerated growth. The Dniester population (GDU B3) was stable until ca. 10 ka; later steady but slow growth can be observed. There were not enough sequences of GDU B4 to perform the eBSP analysis.

Mismatch analysis showed that neither demographic nor spatial expansion model could be rejected for the lineage A or for any GDU of the lineage B (Table 2). Both neutrality tests suggested the recent expansion only in the case of the GDU B1 from the Durugöl (Duruşu) liman, while the Fu’s Fs suggests a recent expansion of the GDU B3 from the Dniester. Interestingly, neither of the neutrality tests indicated recent population size expansion for the widespread and invasive GDU B2.

The analysis of mtDNA haplotype distribution across the distribution range of the demon shrimp clearly showed the two invasion routes to Central and Western Europe (Fig. 2). The most frequently observed in the invaded area were haplotypes 14 and 17. However, while haplotype 14 was found in several localities both along the Central and Southern Corridor, haplotype 17 and haplotype 15 were present only in the North-Western Black Sea drainage (excluding the Danube and its delta) and in the Central Corridor. Only the mtDNA haplotype 17 was found in the UK, suggesting that the local population derives from the Central Corridor. We found four private haplotypes (haplotypes 28–31) in the invaded range in Belarus and Poland. They were not encountered in the native area, probably due to their generally very low frequency. Additionally, it is worth noting that haplotype 32 was found exclusively along the Southern Corridor (Danube drainage). Single haplotypes were recorded in the two Alpine lakes. In the case of Starnbergersee in Germany (site 35), it was haplotype 14, while in the Traunsee in Austria (site 36), only haplotype 32 was found. The population of Balaton lake consisted predominantly of individuals bearing haplotype 14 with a small share of individuals with haplotype 32.

The highest diversity within the D. haemobaphes lineage B expressed by the number of mtDNA haplotypes (12) as well as by the highest private allelic richness was observed in the Durugöl (Duruşu) liman (A_r = 4.03 and PA_r = 4.03, respectively) (Table 1). In all the other localities within the native range of the species, the number of detected haplotypes did not exceed five. The allelic richness was higher than three in two cases, in the Dniester (station 5/5A) and in the Dnieper at Kerch Peninsula (station 13) with the values 3.30 and 3.27, respectively. The private allelic richness was high in the Akkarzhanka River (station 8, PA_r = 1.06) and in the Dnieper at station 16/16A (PA_r = 0.87).

In the invaded range, the number of haplotypes was the highest in the Oder River (station 31/31A, six haplotypes) followed by the station 24 in the Bug River with four haplotypes recorded (Table 1). In all the other stations, the number of detected haplotypes did not exceed three. Allelic richness ranged from one to 3.31, with the highest values observed in the Bug River at station 24 (A_r = 3.31) followed by the station 22/22A in the same river and by the station 27 in the Vistula River (in both cases A_r equaled 2.86). The highest private allelic richness was observed in the Vistula River at station 26 (PA_r = 0.59). Slightly lower values (PA_r = 0.52) characterized the channel in Dubay village (Dnieper River system, station 21) and in the Oder River at the station 31/31A.

No loss of allelic richness was observed along the Central Corridor, ca. 2500 km from the source population to the furthest locality studied, as the highest value of this parameter was observed in the Bug River (station 24, A_r = 3.31) in Poland, ca. 1960 km from the source area (Table 1, Suppl. material 1: Table S1). We observed a weak isolation-by-distance (IBD) effect (Mantel test, R² = 0.0369, p = 0.037) (Fig. 6A).

A slight decrease, although not statistically significant, in haplotype diversity (R² = 0.006, r (17) = −0.07; p = 0.76), nucleotide diversity (R² = 0.018, r (17) = −0.13; p = 0.59), allelic richness (R² = 0.001, r (17) = −0.02; p = 0.92), private allelic richness (R² = 0.131, r (17) = −0.36; p = 0.13) were observed along the Central Corridor of invasion (starting in the source of the corridor at the Dnieprovski Liman locality 11A, belonging to the native range of the species) (Fig. 6B–E, Suppl. material 1: Table S1).

Figure 6.

A plot of pairwise F_st/(1-F_st) versus pairwise linear distance of 19 populations of D. haemobaphes encompassing source populations for the Central Corridor and Central Corridor itself. Plot of: Haplotype diversity (B), nucleotide diversity (C), allelic richness (D) and private allelic richness (E) within 19 populations of D. haemobaphes along source populations for the Central Corridor and Central Corridor itself (see Suppl. material 1: Table S1 for details).

We observed the presence of two highly divergent lineages among the individuals of the demon shrimp (Fig. 3). The COI K2P distance between these two lineages was 0.183, which exceeded the threshold limit widely used for amphipod species delineation and observed even between many morphologically well-defined species (Costa et al. 2009; Raupach et al. 2015; Lobo et al. 2017). The distinctness of these two divergent lineages was also supported by the nuclear (28S) gene (Fig. 4). The divergence between these two lineages possibly dates back to 5.1 Ma (Fig. 3). This coincides with the eustatic marine regression and the formation of separate Dacian and Kimmerian basins in the present Ponto-Caspian area (Cristescu et al. 2003; Stoica et al. 2013). The considerable extension of freshwater habitats in this area during the Pliocene may have induced the evolution of species with freshwater preferences. Among them could be the lineage A, as our findings suggest that it occurs in fresh and not in brackish waters of the north-west Black Sea basin.

Dikerogammarus haemobaphes was originally described from the brackish waters of the Caspian Sea, while Martynov (1919), based on the material from the lower Don River, described the riverine form Dikerogammarus haemobaphes morpha fluviatilis to distinguish it from the brackish water population. Some authors followed this separation (e.g., Cărăuşu 1943; Cărăuşu et al. 1955) and since then the above-mentioned form has often been treated as a valid species (Straškraba 1962; Barnard and Barnard 1983; Jażdżewski and Konopacka 1988; Özbek and Özkan 2011). However, the morphological features separating D. haemobaphes and Dikerogammarus fluviatilis were extremely poorly defined by Martynov (1919), so the existence of the latter species remains unclear, and a thorough taxonomic revision was recommended (Jażdżewski and Konopacka 1988). Our study of the material molecularly assigned to the lineage A did not reveal any morphological evidence that would allow to classify it as D. fluviatilis or as any other species than D. haemobaphes. Given the wide distribution of the lineage B in the native region and taking into account that it includes the presumably Caspian lineage from the Moskva River, we believe that this lineage most probably represents the “real” D. haemobaphes. Thus the lineage A is an undescribed cryptic lineage known so far only from the lower Danube and other rivers in the north-western Black Sea region.

The results of eBSP analyses indicate a recent (less than 20 ka) substantial decrease in the effective population size within the lineage A that may explain its rarity and lack of invasion potential. Such demographic decline is hard to explain, but it may be associated with the late Pleistocene/early Holocene marine transgression and subsequent salinity rise in consequence of reconnection of the former Pontic Lake to the Mediterranean Sea that made up the present Black Sea (Svitoch et al. 2000; Badertscher et al. 2011). This might have forced the numerous freshwater species to move up the rivers emptying to the Black Sea, disconnecting their populations and causing demographic declines.

The divergence pattern within the D. haemobaphes lineage B in its native range is shallow but reflected, to some extent, by the geographic distribution of GDUs. For example, GDU B4 found exclusively in the Moskva River (Caspian Basin) apparently diverged from the others at about 750 ka (Fig. 3), already after the Black and Caspian basins disconnected (Cristescu et al. 2003). A similar separation of genetic lineages was observed for several other local crustacean taxa, including amphipods, mysids, and some cladocerans (Cristescu et al. 2003, 2004; Cristescu and Hebert 2005; Audzijonyte et al. 2009). Geo-demographic unit B1, represented by 12 haplotypes that emerged ca. 200–300 ka was, in our study, restricted to a single locality in the Durugöl (Duruşu) liman in Turkey. A similar pattern was observed for D. villosus (Rewicz et al. 2015), and our findings support the hypothesis that the populations of freshwater and oligohaline crustaceans inhabiting the Ponto-Caspian region have started to diverge during the significant salinity rise that caused a barrier for gene flow between biota “imprisoned” in the nearly freshwater lagoons and limans (Badertscher et al. 2011; Rewicz et al. 2015).

The eBSP analysis showed different historical demographic patterns in the remaining geo-demographic units. In the case of GDU B1, restricted to the Durugöl Liman, the effective population size remained relatively stable. The two GDUs found in the North Black Sea basin (B2 and B3) seem to experience demographic growth in the last couple of thousands of years. These results stay in contrast with the ones obtained for D. villosus; in this case, the post-Pleistocene demographic boost was observed in all the studied lineages of the species (Rewicz et al. 2015).

The distribution of mtDNA haplotypes along the expansion routes clearly points out to the Dnieper River as a donor of invaders for the Central Corridor and, subsequently, for Western Europe. It is supported by the haplotype composition of the studied populations, with haplotypes, such as the most common haplotype 17, unique for the Dniester, Dnieper, Bug, Vistula, Oder and all along the Mittellandkanal in Germany but absent from the Danube and from the Alpine lakes. In the UK we found only haplotype 17, which suggests that the Central Corridor was also the most probable source for the introduction of D. haemobaphes to the UK.

In contrast, the Danubian expansion route is characterized by only two haplotypes, haplotype 14 that is shared with the Central Corridor and the private haplotype 32. Such poor haplotype diversity in the Southern Corridor is hard to explain. Already in the 1950s, the demon shrimp was reported from the middle Danube above Györ in Hungary, and in 1994, it was found as high as the Danube-Main canal in Germany (Straškraba 1962; Schleuter et al. 1994; Borza et al. 2015). In the early 1990s it was abundant in the German part of the Danube, but soon after the appearance of the killer shrimp, the population of demon shrimp in that area declined (Weinzierl et al. 1996; Kley and Maier 2006; Borza et al. 2018). Similar observations were made in other rivers: Rhine, Lahn, Drava, Oder, and Vistula (Weinzierl et al. 1996; Bernauer and Jansen 2006; Kley and Maier 2006; Grabowski et al. 2007b; Chen et al. 2012; Ćuk et al. 2019; own unpublished data), suggesting effective replacement of this species by D. villosus, even if in some conditions co-existence of D. haemobaphes and D. villosus was observed (Borza et al. 2017; Hellmann et al. 2017). Taking into account all the above, we cannot exclude a possibility that the haplotype diversity in the Southern Corridor had once been richer but pauperized after the populations of the demon shrimp declined following lengthy interaction with the killer shrimp.

Interestingly, while the demon shrimp was a successful invader reaching Western Europe via the Central Corridor, the killer shrimp expanded to this area via the Southern Corridor. The population of D. villosus from the lower Dnieper expanded its range along the Central Corridor only to central Poland (Rewicz et al. 2015). The interspecific competition between both congeneric invasive species may be the reason for these differences. The killer shrimp is known to be a very active and robust competitor, which may force less competitive species (such as D. haemobaphes) to leave their optimal ecological niches and extend their ranges into habitats yet unoccupied by competitors (Kobak et al. 2016 and references therein). Additionally, the presence of stronger competitors may lead to niche partitioning as it was observed for Dikerogammarus spp. in the Danube (Borza et al. 2017). Such factors have possibly promoted range extension of the less competitive species and may help to explain why D. haemobaphes conquered rivers along the Central Corridor (Kobak et al. 2016).

The Southern Corridor seems to be the most probable source from which demon shrimp colonized the three studied lakes: the large lowland Balaton Lake in Hungary, as well as the two submontane Alpine lakes – the Traunsee in Austria and the Starnbergersee in Germany. The lake Balaton is connected to the Danube by the Sió canal, allowing different organisms to spread, and the presence of alien species was reported from the lake already since the 1930s (Muskó et al. 2007). The spontaneous range extension through the Danube-Balaton waterway is most probable, but an unintentional introduction of D. haemobaphes together with another crustacean species, Limnomysis benedeni, used as a food source for fish, could also be possible (Muskó 1992; Muskó et al. 2007). The two submontane Alpine lakes also belong to the Danube drainage area, but they are connected to the Danube by relatively small tributaries, namely fast-flowing submontane rivers. The presence of the demon shrimp was not confirmed from these watercourses (Altermatt et al. 2019). Thus, the most probable invasion vector, in this case, is an unintentional introduction associated with recreational activities such as sailing or diving, as it was reported for D. villosus by Bącela-Spychalska et al. (2013). High potential for such dispersal was experimentally evidenced also for the demon shrimp (Bącela-Spychalska 2016).

The small sample of D. haemobaphes from the Moskva River in Russia, also situated outside the native range of studied species (ca. 3000 km up from the Volga delta), seems to be a sister group to all the other GDUs within lineage B found in the Black Sea as well as in Central and Western Europe (Fig. 3). A single available COI sequence mined from GenBank (GenBank accession number: AY529049) and labelled as coming from the Caspian Sea (Cristescu and Hebert 2005) represented the same haplotype as eight out of ten sequences from the site in the Moskva River. This suggests that the Caspian Sea served as a donor for this population, which most likely arrived there via the Volga River belonging to the Northern Corridor of invasion. This scenario may be supported by the fact that D. haemobaphes has been recorded all along the Volga River as well as in its tributary Oka (Zadin 1964) and finally from the Moskva river (Lvova et al. 1996). Intentional introductions and extensive shipping accelerated its invasion process to a large number of water reservoirs built on the Volga River (Karpevich 1975; Slynko et al. 2002). The lack of other molecular data from the Caspian and Azov seas does not allow to track the invasion process in the Northern Corridor in detail. However, it has to be underlined that the shores of the Azov Sea were covered by our sampling and no demon shrimp were found there.

The divergence of the Black and the Caspian Sea populations was observed in the cases of several Ponto-Caspian crustaceans (Cristescu et al. 2003, 2004) and similar situations may be expected for the presently studied species. Although we have found some private haplotypes occurring only in the invaded range of D. haemobaphes, they are all closely linked with the most common haplotypes found in the Dnieper River and its liman. No haplotypes of unknown origin were recorded in Western Europe, so the Caspian Sea may be excluded as a likely source of invading population along the Central and Southern invasion corridors.

Although a weak signal of IBD was recorded, all other measures expressing changes in the level of genetic diversity along the Central Corridor of the invasion were not significant, showing no bottleneck effect associated with the invasion of D. haemobaphes along the invasion corridor. The detected isolation by distance may come from the fact that a relatively high number (four) of private haplotypes existed at some sites in the invaded area (Fig. 6). It could lead falsely to suggest isolation of the invading population from the source one. Generally, it is expected that the molecular diversity in the invaded range might be reduced due to a relatively small number of newcomers. However, recent studies show that such a phenomenon does not have to be as common as previously predicted (Roman and Darling 2007). No genetic diversity loss was observed either during the invasion of D. villosus invasion in Europe (Wattier et al. 2007; Rewicz et al. 2015) or in case of dreissenid mussels and fish in the North American Great Lakes area (Stepien et al. 2005). The phenomena that can explain the lack of genetic diversity loss may be associated with rapid colonization by a large population of invaders or subsequent multiple invasion waves (Stepien et al. 2005; Roman and Darling 2007). In the present case, historical data suggest that the first alternative is more likely.

In conclusion, our results further confirm that the geological history of the Ponto-Caspian region stimulated divergence and speciation in the local freshwater and oligohaline organisms. Some of the resulting closely related lineages cannot be separated morphologically and, at the present stage, have to be treated as cryptic ones. Some of them, such as the lineage A of D. haemobaphes, remained in the area of origin, while others, such as the lineage B, successfully conquered new areas. The colonization of Western Europe by the demon shrimp followed mainly the Central Corridor, while another successful invader, the killer shrimp, used the Southern Corridor predominantly. These differences may be explained by interspecific competition with the killer shrimp (Kobak et al. 2016) that could force D. haemobaphes to leave the native range to reduce the competitive pressure. Another reason may be competitive niche partitioning, as observed by Borza et al. (2017) and perhaps different habitat characteristics of the rivers making up the two invasion corridors. Finally, we conclude the demon shrimp extended its range rapidly and that the expanding population was large as evidenced by the lack of loss of genetic diversity in the populations studied along the Central Corridor.