Cryptic invasion of Baltic lowlands by freshwater amphipod of Pontic origin

Gammarus varsoviensis is morphologically close to G. lacustris, with which it is often misidentified. Geographic range of G. varsoviensis includes Germany, Poland, Lithuania, Latvia, Belarus and Ukraine. Such a distribution pattern led us to the assumption that the species might have originated in the Black Sea drainage area. From there, as early as the nineteenth century, it could possibly have migrated to the Baltic basin through the Pripyat-Bug canal. Thus, the goals of this study are: (1) to indicate the level of genetic divergence of G. varsoviensis from the morphologically closest species – G. lacustris and (2) to investigate the possibility of the Pontic origin of G. varsoviensis and its range expansion across the Black Sea/Baltic Sea watershed to Central Europe through the artificial canal network. Altogether 128 partial 16S rDNA sequences of Gammarus varsoviensis from 19 localities were gained. They were analysed in conjunction with a sequence of G. lacustris obtained in this study from the Dnieper system and sequences of G. lacustris, G. pulex and G. fossarum available in GenBank in order to estimate the relationships among the species. GenBank accession numbers for all the haplotypes defined within this study are: G. varsoviensis – from JN641868 to JN641875; G. lacustris – JN641876. The genetic distance within and between the species was calculated, as well as phylogenetic relationships among haplotypes, which were inferred with Neighbor-Joining method. The haplotype relationships were analysed with the Minimum Spanning Network. Also mismatch distribution of the haplotypes were tested under sudden expansion model and sequence deviations from selective neutrality. Neighbor-Joining analysis revealed that G. varsoviensis and G. lacustris haplotypes formed separate well defined clades. Mean genetic diversity between the two species was ca. 15× higher than the intraspecific distance and similar to the value obtained for G. fossarum/G. pulex species pair. Therefore, G. varsoviensis can be definitely considered as a distinct species from G. lacustris. Among 128 sequences obtained from G. varsoviensis, 8 haplotypes were identified and grouped into two clades: one found only in the lower Dnieper (two haplotypes) and the second one encompassing the rest of haplotypes observed in the upper Dnieper and Baltic Sea basin. Only one haplotype was found in the Baltic Sea drainage area. A mismatch distribution curve as well as selective neutrality tests demonstrated sudden expansion model. Our findings suggest that G. varsoviensis is an alien gammarid that originated in the Pontic area. Its expansion in Central Europe apparently started soon after the opening of the artificial waterways joining the Black and the Baltic Sea drainage basins.


Introduction
The Baltic Sea watershed, including the catchments of the Elbe, Oder, Vistula, and Nemunas rivers, is a relatively young system shaped mostly by the last episode of Pleistocene glaciations.As a result, freshwater fauna of the Baltic lowland rivers, streams and lakes is a product of remaining periglacial aquatic species, Holocene re-colonization process and recent human-mediated introductions of exotic species (Hewitt 2004;Bhagwat et al. 2008;Karatayev et al. 2008;Provan and Bennett 2008).The rivers are interconnected, through a net of man-made canals, resulting in merging the North, Baltic and Black Sea catchment areas (Jazdzewski 1980;Panov et al. 2009).In consequence this system has served as migration route for a number of alien species and has been defined by Bij de Vaate et al. (2002) as the Central Invasion Corridor -one of the most important invasion routes for aquatic organisms in Europe.Thus, among 19 amphipod species (excluding subterranean Niphargidae) inhabiting fresh waters of the Baltic basin, 10 are recognised as aliens (Jazdzewski and Konopacka 1995;Eggers and Martens 2001;Grabowski et al. 2007).As many as eight of them are of Ponto-Caspian origin.The prominent example is Dikerogammarus villosus (Sowinsky, 1894), a successful coloniser of all Central and Western European large rivers -nowadays listed among 100 of the worst invasive species in Europe (DAISIE 2009).
Among native amphipods, Gammarus lacustris G.O. Sars, 1863 is one of the most widely distributed aquatic species, occurring in Europe, Asia and North America (Vainio et al. 2008).Quite recent findings of Vainio and Vainola (2003) and Ianilli et al. (2004) on the molecular diversity of G. lacustris in Europe suggest that in a large part of its range, the species is a post-glacial relic.In southern Europe the species is present only in alpine lakes, while in Central and Northern Europe it thrives in lowland post-glacial lakes and in some typically lowland, slow flowing rivers.A thorough study on G. lacustris from Poland, Germany and Belarus done by Jazdzewski (1975a) revealed not very conspicuous yet stable morphological differences on intra-and interpopulational level, which in consequence led to the description of a new species G. varsoviensis Jazdzewski, 1975.Identity of the latter taxon was then verified by Vainio et al. (1995), who concluded that genetic differentiation of G. varsoviensis vs. G.lacustris detected with allozyme electrophoresis is of the level found between other previously defined Gammarus species.As no conclusions regarding phylogenetic relationships between the species could be drawn by employment of the allozymes, the authors suggested a need for further studies upon the species with other, more accurate, genetic markers.Gammarus varsoviensis is a typical lowland species -the highest sampling site was situated 213 m asl.It inhabits mainly slow running rivers and streams.It is also present in lakes close to the inflow or outflow of rivers, as well as in flooded meadows, old river beds and canals (Jazdzewski 1975a;b;Konopacka 1988;Meijering et al. 1995).Analysing the geographical distribution (Figure 1, Appendix 1) of the newly described species, Jazdzewski (1975a) and Arbaciauskas (2008) reported that within the Baltic basin its range spreads from the Elbe system on the west through the Notec River, lower Vistula drainage, the Narew and Bug rivers to the Nemunas and Daugava rivers on the east.In addition, Jazdzewski (1975a) mentioned one locality of G. varsoviensis in the lower Pripyat River belonging to the Black Sea basin.Based on these data and taking into account the fact that all the above rivers are, or were historically, interconnected by navigable canals, Jazdzewski (1980) as well as Vainio et al. (1995) noticed that the G. varsoviensis distribution range does not conform to any general zoogeographical pattern among native European inland water animals (e.g.Illies 1978).The above authors hypothesized that it might be a recent immigrant from sources outside central Europe, expanding its range through the canal network.Unfortunately, due to the general lack of data on the species distribution outside the already studied area, no further conclusions could be derived.In recent years, during the amphipod collection in the inland waters of the Black Sea basin, we surprisingly found several populations of G. varsoviensis inhabiting tributaries of the middle and the lower Dnieper River (Figure 1, Appendix 1), doubling the already known geographic range of this species, and shifting its limits some 830 km south-eastwards (Grabowski et al. 2012).The Dnieper basin was already formed during the Miocene (some 11.6 million years ago) and survived through Pleistocene glaciations (Matoshko et al. 2002).Taking that into account and following the above cited opinions of Jazdzewski (1980) and Vainio et al. (1995), we came to the supposition that the Dnieper basin may be the missing origin area/glacial refugium of the species, from which in recent times it migrated north-and westwards.Having access to a large collection of G. varsoviensis from various parts of its range, we attempted to fill the above mentioned question marks.
The analysis of mitochondrial DNA (mtDNA) polymorphism is used frequently to investigate the intraspecific genetic structure and to derive phylogenetic relationships among closely related species (Avise 1986;Harrison 1989;Galtier et al. 2009).The advantage of mtDNA is its fast mutation rate (Brown et al. 1979) and its nonrecombining character due to maternal inheritance.Thus, the entire mtDNA molecule can be considered as a single genetic unit with numerous alleles, which can be used to trace maternal phylogenies (Taberlet 1996).16S rDNA is one of the mitochondrial markers suitable for testing phylogeographic hypotheses or for use as a DNA barcoding tool (eg.Müller 2000;Müller et al. 2002;Vences et al. 2005;Lefebure et al. 2007).In the present study, a short fragment of 16S rDNA was analysed to investigate: (1) the level of genetic divergence of G. varsoviensis from the morphologically closest species -G.lacustris; (2) the possibility of the Pontic origin of G. varsoviensis and its range expansion across the Black Sea/Baltic Sea watershed to Central Europe through the artificial canal network.* -site from which a reference sample of G. lacustris was taken.Sampling sites coded as in Table 1.

Methods
In order to cover a broad geographic range of Gammarus varsoviensis, specimens (including a reference sample of G. lacustris) used in the study were drawn from a large collection of the species housed in the Department of Invertebrate Zoology and Hydrobiology, University of Lodz collected during several surveys in Poland, Ukraine, and Belarus (Table 1, Figure 1).All material was preserved in 96% ethanol.

Molecular methods
Altogether 128 individuals of Gammarus varsoviensis from 19 sites, both in the Black (7 sites) and the Baltic Sea basins (12 sites), including the locus typicus of the species, were used for DNA extraction.Additionally, four individuals of G. lacustris from the locality where it co-occurred with G. varsoviensis, were analysed as a reference sample.Around 3 mm 3 of the muscle tissue was taken with a sharp edged tweezers from each individual and incubated over night at 55°C in a 1.5ml tube containing 200µl of Queen's lysis buffer with 5µl of proteinase K (20 mg/ml) (Seutin et al. 1991).
DNA was then extracted based on a standard phenol-chloroform method as described by Hillis et al. (1996).Air-dried DNA pellets were eluted in 100µl of TE buffer, pH 8.00, stored at 4°C until amplification, and subsequently stored long-term at -20°C.As no detailed phylogenetic inference requiring longer DNA sequence was a study requirement, only short (309 bp) easily amplifiable fragment of 16S rDNA was obtained from all individuals with primers designed by Müller (2002) (forward: LR-J-GAM: 5'- ATTTTAATTCAACATCGAGGTTGC-3', reverse: LR-N-GAM: 5'-TTTAACGGCTGCG GTATTTTGAC -3').Each polymerase chain reaction had a total volume of 20µl and contained 10µl of DreamTaq Master Mix (2x) Polymerase (Fermentas), 1.6µl of each primer (concentration 5µM), and 2µl of DNA template.The PCR was performed on MaxyGen Gradient Thermocycler in thermal regime consisted of initial denaturation at 94°C for 3min, and then 34 cycles of denaturation at 94°C for 20s, annealing at 50°C for 45s, and elongation at 65°C for 60s, followed by a final extension for 2min at 65°C.Then 3µl of PCR product was visualized in MidoriGreen-stained (Nippon Genetics) 1.0% agarose gels to verify the quality of the amplification process.Sequences were obtained using BigDye technology by Macrogen Inc., Korea.Sequences were adjusted and aligned using BIOEDIT v 7.0 (Hall 1999), and deposited in GenBank under accession numbers (G.varsoviensis: JN641868, JN641869, JN641870, JN641871, JN641872, JN641873, JN641874, JN641875; G. lacustris: JN641876).

Statistical analysis
Part of the analysis, including calculation of within and between species genetic distance (model T92+G), as well as a rough estimation of phylogenetic relationships among haplotypes was performed with Mega 5 (Tamura et al. 2011).The following reference sequences from GenBank were used in the analysis: G. pulex from the Netherlands -GenBank Accession number EF582877 (Hou et al. 2007) Analysis of the haplotype affinities (Minimum Spanning Network), mismatch distribution as well as deviations from selective neutrality were tested with Tajima's D and Fu's F approaches in the package ARLEQUIN v3.11 (Excoffier et al. 2005).The neutrality tests were used as an indication of recent population expansion when the null hypothesis of neutrality was rejected due to significant negative values.The distribution of pairwise nucleotide differences (mismatch distribution) was calculated as an additional test for demographic expansion with the sum of square deviations (SSD) and Harpending's test computed to estimate whether the observed distributions deviated significantly from those expected under the population expansion.

Genetic divergence of Gammarus varsoviensis from G. lacustris
The partial 16S rDNA sequences (309 bp long) were obtained from 128 individuals of Gammarus varsoviensis collected from a total of 19 localities (seven within the Black Sea basin, twelve within the Baltic Sea basin) covering most of its currently known range.Among them  2. Neighbor-Joining (NJ) analysis revealed that G. varsoviensis and G. lacustris haplotypes grouped according to the species and formed well defined clades with, respectively, 100% and 98% bootstrap values (Figure 2).Within species, genetic diversity of the identified G. varsoviensis haplotypes was comparatively similar to that of G. lacustris, 0.012±0.004and 0.009±0.004,respectively (Table 2).Mean genetic diversity between the two species was 0.146±0.027,around 15× higher than the intraspecific distance.That is relatively similar to the value obtained for another pair of morphologically close species -G.fossarum/G.pulex (0.230 ± 0.033) and lower than between G. varsoviensis or G. lacustris and any of the other two species (Table 2).The above genetic distance data and NJ tree topology support well the distinctness of G. varsoviensis from G. lacustris, yet they show their closer relatedness than with the other Gammarus species from European lowlands.Haplotypes are coded as in Table 2.

Pontic origin and range expansion of Gammarus varsoviensis
As seen from the Minimum Spanning Network (MSN) analyses, the G. varsoviensis haplotypes could be divided into two groups (Figure 3a).The first one groups 6 haplotypes (Gvar1 to Gvar6) differing pairwise by up to 3 transitions, with Gvar1 or Gvar 2 being probably the ancestral ones.These variants can be found in the north-western part of the species distribution range (Figure 3b).All of them occur in the upper part of the Dnieper catchment belonging to the Black Sea basin.Only Gvar1 may be found in the Baltic Sea drainage area, including the Nemunas River as well as in all the localities within the Vistula and Oder river systems.The second group comprises only two haplotypes (Gvar6, Gvar7), differing from each other with just one transition, yet separated from Gvar1 with 4-5 substitutions (1 transversion, and 3-4 transitions).That group was found only in the Konka River, being a left tributary of the lower Dnieper system, in the south-eastern part of its distribution range (Figure 3b).As this locality was distant from all the others in geographical terms, and as the local haplotypes were sitespecific and quite divergent from the rest, it was excluded from the subsequent mismatch distribution analysis.
The population inhabiting the upper part of the Dnieper catchment area and the Baltic basin clearly showed an unimodal pattern of mismatch distribution curve, as expected in cases of population expansion (Figure 4).The sum of squared deviations (SSD) of mismatch distribution is very low, indicating that the observed curve fits well the sudden expansion model tested.Although the Harpending's raggedness index value is higher, it still does not differ significantly from the expected one, further indicating that the curve fits well for the sudden expansion model tested.Additionally, values of both selective neutrality tests, Tajima's D and Fu's Fs, are negative and statistically significant suggesting a departure from neutrality that may be related to population demographic expansion.The haplotype frequency pattern is congruent with these results (Figure 3b).If, based upon haplotype distribution and results from MSN analysis, the upper part of the Dnieper catchment area and the Baltic basin are treated as one demographic unit, then the overall frequency of Gvar1 in this area is 0.89 (Table 3).Haplotype Gvar1 was found there in all analysed localities, solely represented in all 12 samples from the Baltic basin and in three samples from the upper Dnieper area.In the remaining localities from the latter area, Gvar1 had a frequency from 0.36 to 0.83 (Figure 3b).

Genetic divergence of G. varsoviensis from G. lacustris
The performed NJ analysis clearly indicated that G. varsoviensis is a well defined species.All haplotypes were grouped according to the species, regardless their geographic origin -eg.G. lacustris from the Pripyat River was closer to its conspecifics from Asia and North America than to G. varsoviensis from the same sampling area.This supports well the taxonomic features differentiating G. varsoviensis from otherwise morphologically similar G. lacustris, as defined by Jazdzewski (1975a).It also agrees with results obtained by Vainio et al. (1995), separating the two species based on allozyme polymorphism.The above authors concluded that genetic differentiation of G. varsoviensis vs. G.lacustris detected with allozyme electrophoresis is at the level found in other previously defined Gammarus species.The genetic distance between G. varsoviensis and G. lacustris (0.146±0.027) was similar to the divergence level between two other Central European species -G.fossarum and G. pulex (0.230 ± 0.033).Such a level of divergence corresponds well to the median value of 0.037 for 16S rDNA marker found among congeneric, morphologically well-defined, species within Crustacea (Lefebure et al. 2006).Also the branching order supported by high bootstrap values and branch lengths on the NJ consensus tree suggest a much closer phylogenetic relationship of G. varsoviensis with G. lacustris than with any other of the above species.The same phylogenetic relationship among these species may be seen in the paper by Hou et al. (2011) showing divergence among 115 Holarctic Gammarus species, based on four mitochondrial and nuclear genes (altogether ca. 5 kbp long) -however, no intraspecific diversity was taken into account by these authors.

Pontic origin and range expansion of Gammarus varsoviensis
The diversity hotspots of freshwater amphipods in Europe are observed in the middle latitudes, that were free of ice during the Pleistocene glaciations.The formerly glaciated areas have been subsequently colonised from such areas (Vainola et al. 2008).The distribution of Gammarus varsoviensis (Figure 1) revealed by our latest findings (Grabowski et al. 2012) encompasses both the Black Sea basin (Dnieper system) as well as the Central European plains belonging to the Baltic basin (Daugava, Nemunas, Vistula, Oder systems).Its most western localities reach the right tributaries of the Elbe River, belonging to the North Sea basin (Zettler et al. 1999).The plains inhabited by G. varsoviensis between Elbe and Daugava rivers underwent several glaciation events in Pleistocene.In contrary, most of the Black Sea basin was free of ice during that time (Ehlers and Gibbard 2004).There is evidence that the Dnieper system had already originated in the Miocene and survived relatively intact through the Pleistocene (Matoshko et al. 2002).Thus the most possible scenario is that it served as the origin or a glacial refugium area for the species from which it has colonized Central Europe.
Such history is well evidenced by our results showing relatively high haplotype diversity within the Dnieper system if compared to no diversity detected over the wide area of Baltic plains.Minimum Spanning Network (Figure 3a) indicates the presence of two haplotype groups.The one from the southern part of the Dnieper system is represented in our results by two closely related haplotypes (Gvar7, Gvar8) restricted only to the Konka River.The other group consists of six closely related haplotypes (Gvar1 to Gvar6) inhabiting the upper part of the Dnieper system and the Baltic plains.An obvious bottleneck effect may be observed between the population from the tributaries of the middle or lower Dnieper and the populations inhabiting the tributaries of Pripyat or the Baltic rivers, where only one 16S haplotype (Gvar1) was found.Our analysis of mismatch distribution curves as well as neutrality indices (Tajima's D, Fu's F) strongly support such an expansion event.A question arises whether the colonisation progressed naturally through trans-watershed periglacial lakes existing at the beginning of Holocene (Rolik and Rembiszewski, 1987) or in historical time through artificial navigable canals joining the two watersheds (Jazdzewski 1980).The answer cannot be concluded from the molecular data, yet we may try to derive it from some indirect evidence.Firstly, the present distribution of G. varsoviensis outside the Black Sea basin (Figure 1) strictly follows the artificial canal network from the Dnieper system to the Elbe.Daugava was joined to the upper Dnieper in 1805, Nemunas and Vistula were connected to the Pripyat River in 1784, Oder to Vistula in 1774, and Elbe to Oder in 1746 (Jazdzewski 1980;Bij de Vaate et al. 2002;Kubijovye 1984).As far back as the records are known, G. varsoviensis occurs predominantly in the rivers connected by canals, or in their direct affluents (Jazdzewski and Konopacka 1995).That clearly indicates colonisation through the canals.The first record of the species near Berlin (the western border of its range) was in 1898, suggesting that it was one of the first successful alien colonizers of Central European waters from the Pontic region, along with the zebra mussel (Dreissena polymorpha) and Chelicorophium curvispinum (see Mordukhaj-Boltovskoj 1964;Jazdzewski 1980;Bij de Vaate et al. 2002).What could promote such a colonisation?Since 1970 our team has collected specimens of G. varsoviensis from 158 sampling locations in Poland and in Ukraine (Jazdzewski 1975a;own unpublished data).In most cases it was the only amphipod recorded in the site.Only on 26 sites did it co-occur with other gammarids (in 18 sites with natives: G. fossarum, G. lacustris; in 8 with invasives: G. roeselii, Dikerogammarus haemobaphes, D. villosus, Chaetogammarus ischnus).This trend may be tentatively explained in several ways.One is that this species colonised water bodies avoided or abandoned by other amphipods, e.g.due to ecosystem degradation.Another possibility is that G. varsoviensis outnumbered formerly occurring species.If that is true, then G. varsoviensis should be considered not only as alien but also as an invasive species in Central Europe, according to the definition from IUCN (2000).There are several factors listed as promoting invasion of a nonindigenous gammaridean species.Among them are: relatively short generation time, early sexual maturity, high fecundity, size of an animal, euryoeciousness and effective predation (Bij de Vaate et al. 2002;Grabowski et al. 2007).There is only one study done by Konopacka (1988) upon the life history of G. varsoviensis.Grabowski et al. (2007) made multiple factor comparison of life history traits of native vs. invasive gammaridean species.Interestingly, G. varsoviensis grouped together with G. lacustris G.O. Sars, 1863 and G. roeselii Gervais, 1835.Gammarus lacustris is a typical limnic cosmopolitan species widely distributed in the whole Holarctic (Holsinger 1976;Karaman and Pinkster 1977), while G. roeselii is an alien species of the Balkan origin, spreading to Northern and Western Europe most probably in the early twentieth century (Karaman and Pinkster 1977;Jazdzewski 1980;Jazdzewski and Roux 1988).All these species are relatively large gammarids and are characterized by quite high fecundity compared to other freshwater Gammarus spp. in Europe.According to Grabowski et al. (2007) G. varsoviensis is the most fecund species with the brood size on average 1.5 times higher than observed in the native species that could have occurred previously in the ecosystems: G. fossarum, G. pulex and G. lacustris.On the other hand G. varsoviensis, has a univoltine life cycle and the reproductive period is very short -it lasts only 5 months (Konopacka 1988).There is no more data available on its ecology or physiological response to ecological stresses such as high temperature or salinity.From the field data we may conclude that this species can occur in freshwaters with relatively high ionic content up to 0.6 PSU (Grabowski et al. 2009).Such feature could promote the species expansion.
Concluding from the above, G. varsoviensis recognised until now as native to Central Europe is apparently an alien gammarid originated in the Pontic area.It has colonised Central European lowland rivers through the artificial navigable canals that had already joined the Dnieper, Vistula and Elbe catchment areas in the eighteenth and nineteenth centuries.This colonisation apparently happened soon after opening of these waterways.Therefore, the species should become an object of thorough studies on both its ecology and impact upon local aquatic biodiversity in the future.

Figure 1 .
Figure 1.Distribution of Gammarus varsoviensis in Europe (Appendix 1).Numbered circles indicate sites sampled for DNA analysis.*-site from which a reference sample of G. lacustris was taken.Sampling sites coded as in Table1.

Figure 3 .
Figure 3. (a) Minimum spanning network of 16S rDNA haplotypes identified in Gammarus varsoviensis.Number of mutational steps among haplotypes are indicated.Ts -transition, Tv -transversion.(b) Geographical distribution of G. varsoviensis haplotypes in the Black and the Baltic Sea basins.Numbered circles represent sampling localities coded as inTable 1. Dashed line indicates watershed boundaries.Haplotypes are coded as in Table 2.

Figure 4 .
Figure 4. Observed and expected mismatch distribution under population expansion model, sums of squared deviations (SSD), Harpending's Fu's Fs and Tajima's D neutrality test value for Gammarus varsoviensis 16S rDNA haplotypes.

Table 1 .
Locations sampled and number of Gammarus varsoviensis individuals used for DNA analyses.Geographic coordinates given in decimal degrees.m asl -meters above sea level, n -sample size, * -location from which reference sample of G. lacustris was taken.

Table 2
Mean intra-and interspecific genetic distance ± standard error (model T92+G) calculated for selected Gammarus species, sample size in brackets.NA -not available

Table 3
Frequency of Gammarus varsoviensis haplotypes in the studied regions.GenBank accession numbers and sample size in brackets.