Genetic structure of three invasive gobiid species along the Danube-Rhine invasion corridor: similar distributions, different histories

Ponto-Caspian gobiids have expanded their ranges throughout Europe since the 1990s. While genetic studies have been widely used to assess the invasion history of gobiids in North America, complex genetic studies involving multiple sites and species have been less common in Europe, severely limiting our understanding of invasion processes along navigable rivers and their tributaries. In this study, we used both nuclear and mitochondrial markers to assess genetic diversity and structure in native and non-native Western tubenose goby Proterorhinus semilunaris, round goby Neogobius melanostomus and bighead goby Ponticola kessleri sampled from the main areas of their joint distribution, i.e. the lower Danube, middle Danube and lower Rhine. Additionally, we describe expansion into Danubian tributaries and provide early genetic data for N. melanostomus from the River Elbe. Our data revealed i) a founder effect in non-native P. semilunaris, ii) an increase in genetic diversity in non-native N. melanostomus samples from the Rhine and Elbe, and iii) no genetic structuring in P. kessleri. This suggests greater initial propagule pressure in P. kessleri; strong propagule pressure with introductions from multiple sources followed by admixture for N. melanostomus in the Rhine and Elbe; and one or very few introduction events for P. semilunaris. We provide further support for the Danubian origin of all three goby species in the Rhine and document lower genetic diversity in fish colonising non-navigable tributaries. Our results illustrate how the ranges of invasive species can become sympatric, despite clear differences in their invasion histories.


Introduction
Five different gobiid species have greatly increased their ranges in recent decades.From their native area around the Black and Caspian Seas and the lower stretches of large rivers such as the Dnieper and Danube, they have colonised numerous fresh-and brackish water habitats in Europe and North America (Kornis et al. 2012;Roche et al. 2013).The speciesspecific spread patterns of these five species appear to differ, however, for reasons that are currently not well understood.Of the five species, the monkey goby, Neogobius fluviatilis (Pallas, 1814), and the racer goby, Babka gymnotrachelus (Kessler, 1857), have spread mainly through the Rivers Dnieper and Vistula, although patchily-distributed populations extend to the middle Danube (Grabowska et al. 2010;Szalóky et al. 2015), and N. fluviatilis has also colonised the lower Rhine (Borcherding et al. 2011).The bighead goby Ponticola kessleri (Gunther, 1861) has colonised parts of the Dnieper, the whole of the Danube and almost all of the Rhine (Roche et al. 2013).The Western tubenose goby, Proterorhinus semilunaris (Heckel, 1937), and the round goby, Neogobius melanostomus (Pallas, 1814), have similar nonnative ranges to P. kessleri but have also spread into other European waters (e.g.P. semilunaris in the Vistula and N. melanostomus in the Baltic Sea; Grabowska et al. 2008;Skóra and Stolarski 1993).Both P. semilunaris and N. melanostomus have been introduced into the Laurentian Great Lakes Basin of North America (Jude et al. 1992).
Inference of invasion pathways based on direct observations of presence and absence data can be problematic.However, indirect methods, based on the distribution of genetic variability, have enabled precise reconstructions of the invasion process (Handley et al. 2011), especially when combined with relevant geographical and ecological data (Croucher et al. 2013).The use of molecular markers for reconstructing a species' invasion history is now well established (Estoup and Guillemaud 2010), having been successfully applied in identifying the source(s) of introduction, invasion pathways, type(s) of vector involved, number of independent introductions, and pathways of secondary spread for a wide range of invasive species (reviewed in Cristescu 2015).
Using such methods, the invasion histories of North American gobiids have been well documented, with a number of studies describing the origin and spread of gobiids in the Great Lakes Basin (Dougherty et al. 1996;Dillon and Stepien 2001;Brown and Stepien 2008;Brown and Stepien 2009;Neilson and Stepien 2009a, b).For example, a recent study examining genetic changes in the Great Lakes since invasion found that while populations have stayed genetically consistent, they show considerable between-population divergence, thus retaining the signatures of their different sources over time (Snyder and Stepien 2017).
The invasion histories of introduced European populations are less clear, as genetic studies have often been restricted to a comparison of relatively few sites or populations (for exceptions, see Ohayon andStepien 2007 andGrabowski et al. 2016).This is particularly so for the expansion of P. semilunaris, N. melanostomus and P. kessleri along European rivers.All three species appear to have spread simultaneously along the Danube-Rhine corridor (an important route for Ponto-Caspian species migrating west; Tittizer 1997) over the past 20 years or so; yet our knowledge of the pathways through which these species have achieved such a significant overlap in distribution remains relatively scarce.Ondračková et al. (2012) suggested the lower Danube as the source for P. kessleri and N. melanostomus colonisations in the middle Danube.Cerwenka et al. (2014), who have provided the sole comparison between Rhine and Danubian gobies to date, found no significant difference between P. kessleri and N. melanostomus from the Rhine and upper Danube.Finally, Adrian-Kalchhauser et al. (2016) suggested the lower Rhine as the source (or one of the sources) for the P. kessleri population in the upper Rhine.
In this multi-species, multi-site study, we aim to reveal potential inter-specific differences in source area and invasion history for three gobiid species now found all along the Danube-Rhine corridor: P. semilunaris, N. melanostomus and P. kessleri.To this end, we compare genetic diversity and structure in samples obtained from sites (both indigenous and non-native) where all three species co-occur, using a combination of nuclear (10 microsatellites) and mitochondrial (cytochrome b sequences) markers.Additionally, by sampling P. semilunaris and N. melanostomus from tributaries of the Danube (the Rivers Morava and Dyje), we aim to clarify the processes governing dispersal into tributaries from main stem rivers.Finally, we analyse samples from a newly established N. melanostomus population in the upper Elbe (Roche et al. 2015), thereby providing early insight into their genetic structure and the putative source(s).

Sample collection
Fishes were collected during the gobiid expansion from 2005 to 2015, with all fish from each site being sampled in the same one or two year period (Table 1).Original samples from 2005-2006 were obtained to compare gobiid populations from the Austrian and Bulgarian Danube.Samples from other sites were added following the N. melanostomus invasion of the Rhine or Danubian tributaries less than five years after the first report of the species (Table 1).A total of 118 P. semilunaris, 146 N. melanostomus and 94 P. kessleri were collected, mostly from three sites along the Danube-Rhine corridor where all three species are now present: 1) the Bulgarian Danube near the town of Vidin (43º57′35″N; 22º53′16″E; all three species indigenous); 2) the Austrian Danube near the village of Orth an der Donau (48º07′23″N; 16º42′46″E; P. semilunaris indigenous [present for 100+ years], other two species non-native), close to the first recorded site for P. kessleri outside of its native range (Zweimüller et al. 1996); and 3) the lower Rhine near the town of Rees (51º45′42″N; 06º20′21″E; all three species non-native) in the navigable German stretch (Figure 1, Table 1).Samples of non-native P. semilunaris and N. melanostomus were also collected from the River Morava (Czech Republic; 48º39′11″N; 16º58′03″E), a non-navigable Danubian tributary.Additional non-native P. semilunaris were sampled from the River Dyje, a tributary of the Morava (48º38′35″N; 16º55′27″E), close to the site where P. semilunaris was first detected (Lusk and Halačka 1995) in the Czech Republic (Figure 1).Samples of N. melanostomus were obtained from the upper Elbe (Czech Republic; 50º39′58″N; 14º06′04″E) in 2015, the year it was first recorded (Figure 1, Table 1).
Fish were sampled from a natural stretch of the Bulgarian Danube (gravel beach with sparse rocks), while fish from all other sites were sampled from the  rip-rap banks of the main channelised stem of the rivers.The only exceptions were P. semilunaris from the Rhine and middle Danube, which were taken from backwaters as species presence in the main river channel was too sparse (Borcherding et al. 2011).Fish sampling was conducted with respect to the habitat sampled, i.e. electrofishing was used along rip-rap, beach seining along the natural river stretch and a combination of both methods in backwaters.

DNA extraction and genotyping
Fin clips were taken from each fish immediately upon capture and stored in 96% ethanol for subsequent genetic analysis, the material being deposited in the collection of the Institute of Vertebrate Biology of the Czech Academy of Sciences in Brno (Czech Republic).DNA was extracted from the fin clips using the Genomic DNA mini kit (Geneaid, Taiwan), following the manufacturer's protocols.We genotyped all samples at 10 microsatellite loci previously used successfully in studies of gobiid invasion, using primers and annealing temperatures designed by Vyskočilová et al. (2007), Feldheim et al. (2009) and Dufour et al. (2007).Amplifications were performed on 10 µl samples using the Qiagen Multiplex kit (Qiagen, the Netherlands), following the manufacturer's protocol.The design of the three multiplexes used for fragment analysis was the same for all three species, with Mix 1 containing loci NG71, NG111, NG215; Mix 2 containing NG70, NG115, NG135; and Mix 3 containing Ame01, Nme5, Nme7 and Nme8.A multiplex polymerase chain reaction (PCR) was performed on both Mix 1 and Mix 2. For Mix 3, however, PCR was performed separately for each locus, the PCR products being mixed later for fragment analysis.The amplified products were subjected to standard fragment analysis on an ABI PRISM 3130 Genetic Analyser (Life Technologies, Canada) using a 500 LIZ Size Standard.Electrophoretograms were scored in GENEMAPPER v. 3.7 (Life Technologies, Canada).Cytochrome b (CYTB) was chosen as a mitochondrial marker as it is widely used in phylogeographical studies, it allows detailed analysis of matrilineal genetic structure, and numerous sequences are freely available in GenBank.CYTB was amplified using the AJG15 and H5 primers (Akihito et al. 2000; targeting 1138 bp in P. kessleri and 1140 bp in P. semilunaris and N. melanostomus).Amplifications were performed on 10 µl samples using PPP Master Mix (Top-Bio, Czech Republic), following the manufacturer's protocol.Annealing temperature (T a ) was 53 °C during the first five PCR starting cycles, followed by 30 cycles with T a = 58 °C.PCR products were Sangersequenced bi-directionally in a commercial laboratory.New CYTB sequences were submitted to GenBank under accession numbers KJ605186-KJ605189 (P.kessleri) and KJ605190-KJ605212 (P.semilunaris).

Genetic diversity and population structure from microsatellite data
No microsatellite pair was in linkage disequilibrium and all loci but one were in Hardy-Weinberg (HW) equilibrium.The only exception was locus Ame01 in four of the five N. melanostomus samples, disequilibrium probably being caused by null alleles.As the increased frequency of null alleles has a negligible effect on measures of genetic variability and structure, we present the results based on the complete dataset (for more details see Appendix 1).
Intra-population genetic diversity was compared using the number of different alleles (Na), number of private alleles (Np) and unbiased expected heterozygosity (uHe), all of which were calculated in GenAlEx v. 6.5 (Peakall and Smouse 2006).As these measures can be influenced by sampling effort, we also calculated allelic richness (AR) corrected for sample size (using the rarefaction method in FSTAT v. 2.9.3.2;Goudet 2002) based on a minimum sample size of 15, 29, and 22 diploid individuals of P. semilunaris, P. kessleri and N. melanostomus, respectively.A relatively low degree of polymorphism and low number of loci made results of common statistical comparisons of genetic diversity (e.g.Friedman test) extremely conservative (see Appendix 2 for test results and more details).Hence, we base our interpretations on a simple visual comparison of the calculated values (all trends in the interpretation are supported by detailed loci-specific comparisons presented in Appendix 2).Detailed information on allele occurrence, fish sex and fish size is provided in Appendix 3.
Genetic differentiation between sampling sites was quantified by computing pairwise estimators of F ST according to Weir and Cockerham (1984), significance being tested with 1000 permutations in GENETIX v. 4.05 (Belkhir et al. 1996(Belkhir et al. -2004)).The Benjamini-Hochberg procedure (Benjamini and Hochberg 1995) was applied to control for false discovery rate in multiple comparisons.The genetic relationship between all genotyped individuals for each species were obtained by factorial correspondence analysis (FCA) implemented in GENETIX v. 4.05.FCA detects the best linear combination of variables (allele frequencies at different loci) and describes the variation between observations (individuals in this case).
The Bayesian clustering procedure implemented in STRUCTURE v. 2.3.4 (Falush et al. 2003) was used to infer the number of distinct genetic populations represented by the samples, and assignment of individuals to these genetic clusters.The Bayesian model assumes K (unknown) populations with different allele frequencies at a set of independent loci.The program was run with five independent simulations for each value of K from 1 to the maximum number of populations sampled (P.semilunaris = 5, P. kessleri = 3, N. melanostomus = 5), each with 10 6 iterations following a burn-in period of 10 5 iterations.In all simulations, we used an admixture ancestry model and an independent allele frequency model (with default parameters).The likelihood of K, i.e.Ln Pr(X|K), was used to infer the most likely number of real populations in the dataset, using the delta K method (Evanno et al. 2005).The results of five replicate runs for each value of K were combined and evaluated in CLUMPAK (http://clumpak.tau.ac.il).Summary outputs for the most supported K value were displayed graphically using DISTRUCT v. 1.1 (Rosenberg 2004).

Genetic variability
Patterns in genetic variability clearly differed between species, particularly the non-native Rhine samples.Samples from both Danubian sites showed similar genetic diversity in all three species (similar AR in all three species; similar uHe in P. semilunaris and P. kessleri; Table 2; Appendix 2).

Population structure -microsatellites
All pairwise F ST values but one were significantly higher than zero, though values showed high variability (Table 3).The most distinct P. semilunaris sample was that from the Dyje, whilst that for N. melanostomus was from the Elbe.These results were supported by those from FCA and STRUCTURE.For P. semilunaris, the best STRUCTURE model indicated three genetic clusters (K = 3).According to both STRUCTURE and FCA, the two Danubian P. semilunaris samples were almost identical, suggesting they belonged to a single population.The Dyje sample appeared most distinct, with FCA indicating it as a subset of the Morava sample.The Morava sample itself showed possible admixture between the Dyje and Danube.Both STRUCTURE and FCA showed that the Rhine sample appears to be distinct from Danubian samples (Figures 2A and  3A), though there was a partial overlap in FCA with individuals from the Austrian Danube (Figure 3A).
The best model in STRUCTURE also indicated three genetic clusters (K = 3) for N. melanostomus.Despite a slight shift in genetic composition along the Danube from Bulgaria to Austria and on up the Morava (note the gradual increase in the proportion of the "blue" cluster in the STRUCTURE barplot; Figure 2B), all three sites (DAN-BG, DAN-AU, MOR) showed high overlap on the first two FCA axes (Figure 3B).On the other hand, the Rhine and Elbe samples both appear to differ genetically (with either "red" or "green" clusters prevailing at each site; Figure 2B) and to form separate groups in FCA (Figure 3B), only partially overlapping with the Danubian and Moravian samples.
The best STRUCTURE model for P. kessleri indicated no structure between the three samples (i.e.model K = 1 had the highest likelihood; not shown).This was further indicated by a non-significant pairwise F ST for the Austrian Danube and lower Rhine, and a relatively uniform cloud of genotyped specimens in FCA, with especially high overlap between the Austrian Danube and the Rhine (Figure 3C).FCA uniformity was not complete, however, and some differences were visible between the groups.Note, however, that the major differences occurred along the second factorial axis, which explained just 20% of the data variability, while the first axis, for which there was major overlap, explained 80% of the variability.

Haplotype networks
Based on the CYTB marker, P. semilunaris appeared more diverse than the other gobiid species examined.We recorded 23 haplotypes clustered in two distinct haplogroups (Figure 4).Haplogroup A was dominated by central haplotype Pro9, which was present in all sampling regions.Pro9 matched a haplotype previously reported from the Danube in Serbia (Neilson and Stepien 2009b; Pro10 was also found at the same Serbian locality; Appendix 3).Some of the peripheral haplotypes around Pro9 were not present in indigenous Danubian samples.The Rhine and Dyje samples contained haplotypes from Haplogroup A only.The second grouping, Haplogroup B, showed high diversity for the Bulgarian Danube, though two haplotypes were also common in the Morava and one (Pro68) in the Austrian section of the Danube (Figure 4).No Haplogroup B sequences were found in available GenBank sequences (for more details see Appendix 3).
All N. melanostomus possessed haplotype ame1 (Appendix 3), the dominant central haplotype for the Black Sea basin.This haplotype is widespread and abundant in all non-native samples in the Danube and in the Gulf of Gdansk in the Baltic Sea (Brown and Stepien 2008).
Four CYTB haplotypes were recorded for P. kessleri, in addition to the two downloaded from GenBank (Figure 5

Discussion
Our data revealed (i) a founder effect in non-native P. semilunaris samples, (ii) relatively high genetic diversity and private alleles in non-native (Rhine and Elbe) N. melanostomus samples, and (iii) lack of genetic change during P. kessleri expansion (though there was slight differentiation based on FCA).In combination with data on population dynamics and first records of the species, these results allow us to develop and compare putative invasion history theories for these three sympatric gobiid species.

Colonisation of the middle Danube
For all three species, there was relatively high similarity between the two Danubian samples (Bulgarian and Austrian).This was expected for P. semilunaris, the only invasive gobiid that has historically been found upstream of the Iron Gates Dam (the barrier separating the lower and middle Danube; see Figure 1).As such, it is now considered indigenous in the Danube up to the Austrian stretch where we obtained our samples (Roche et al. 2013).The high genetic similarity between the Austrian and Bulgarian samples confirms a continuous pan-Danubian population.
The indigenous distribution of both N. melanostomus and P. kessleri is restricted to the lower Danube (up to the Iron Gates for P. kessleri and around 70 river km downstream of the Iron Gates for N. melanostomus [our Bulgarian site in Vidin]).Both species colonised the middle Danube at approximately the same time, first appearing near our sampling site in 1994 (P. kessleri) and 2000 (N.melanostomus).The lack of any records between the lower Danube and Austria at that time makes transport in ship ballast water the most plausible vector for these new populations (Polačik et al. 2008;Roche et al. 2013).Austrian samples of both species showed high AR, suggesting high propagule pressure, and were indistinguishable from Bulgarian samples using FCA (also STRUCTURE analysis for P. kessleri), suggesting the lower Danube as a putative source for the middle Danube occurrences.The same conclusion was also reached by Ondračková et al. (2012) when comparing genetic diversity and parasite assemblages.
While there is strong support for the lower-tomiddle Danube invasion pathway for P. kessleri and N. melanostomus, evidence for direct transition from the Bulgarian to Austrian Danube is equivocal.Austrian samples of both species possess private microsatellite alleles that are absent in the Bulgarian samples (Table 2).Moreover, mitochondrial haplotype Hap_1, common in both middle Danubian and Rhine N. kessleri, was not found in our Bulgarian sample, though was observed in another sample from the lower Danube from Prahovo, Serbia (Neilson and Stepien 2009b).Nonetheless, based on our analyses the putative source for middle Danubian P. kessleri is likely to be another stretch of the lower Danube.A putative non-Danubian Black Sea basin source for P. kessleri is unlikely, given that the two previously reported non-Danubian Black Sea haplotypes (Neilson and Stepien 2009b;Medvedev et al. 2013) were not found in any of the samples covered by this study.
As for N. melanostomus, Brown and Stepien (2009) found high similarity between two non-native samples (middle Danube [Slovakia] and lower Danube [Prahovo]) and a sample from Odessa (northwest Black Sea).In doing so, they detected haplotype ame7 in the Prahovo sample, which has since only been detected in samples from the port of Kherson on the Dnieper (a putative source for North American N. melanostomus).These observations may suggest that both the Prahovo and middle Danubian samples originated through an introduction from Ukraine.On the other hand, Brown and Stepien (2009) did not include an analysis of indigenous Danubian samples.It is possible that haplotype ame7 is also present at low frequencies within the indigenous range of N. melanostomus in the lower Danube.As Prahovo lies just downstream of the Iron Gates, with no transversal barrier between Prahovo and the species' indigenous range, the population could also have originated through natural upstream migration.The lower Danube thus remains the most likely source for middle Danubian N. melanostomus.

Colonisation of the Rhine
All three species were introduced into the lower Rhine between 1999 and 2007 (Table 1), with introduction via cargo vessel ballast strongly implicated (Roche et al. 2013).As for the middle Danube, genetic diversity for both N. melanostomus and P. kessleri in the Rhine remains relatively high.In contrast, P. semilunaris from the Rhine showed relatively low diversity compared to Danubian samples, suggesting just one or few introduction events linked with a loss of genetic diversity (i.e. a founder effect).
Our data suggest that Danubian fish played an important role in colonisation of the Rhine by all three species.This is particularly visible with P. kessleri, where the Rhine sample was indistinguishable from Danubian samples.The same conclusion was reached by Cerwenka et al. (2014) when comparing Rhine P. kessleri with those from the upper Danube.This maintenance of genetic diversity during the invasion process suggests multiple introduction events, or at least strong initial propagule pressure, as observed following the recent confirmed introduction of P. kessleri via shipping in the upper Rhine (Adrian-Kalchhauser et al. 2016).
Haplotype analysis also suggests a Danubian origin for Rhine P. semilunaris (Figure 4).Although FCA results indicate a degree of dissimilarity between Rhine and Danube samples, this could be attributable to genetic drift.Mombaerts et al. (2014) also suggested the Danube as the original source for P. semilunaris in Belgian rivers and canals (connected with the lower Rhine).They identified the Serbian Danube as the putative source based on the sample of Neilson and Stepien (2009b), which included the Serbian haplotype Pro9.Our results, however, demonstrate that Pro9 is common throughout the native stretch of the Danube and, as such, the Rhine P. semilunaris could have originated from any part of the Danube up to Vienna.
The true source(s) of Rhine N. melanostomus remains unresolved.All lower Rhine samples analysed to date (see Cerwenka et al. 2014;Tserkova et al. 2015; this study) possess haplotype ame1 only.While this excludes introduction from the Caspian drainage, and implicates the Black Sea, haplotype uniformity prevents any finer-scale mtDNA-based location of origin.Neogobius melanostomus could have entered the Rhine from the Black Sea basin by several putative pathways.A first option is introduction from the Baltic Sea.The species was first introduced into the Baltic in early 1990 (Skóra and Stolarski 1993), probably originating from one of the Black Sea ports as ame1 is the only haplotype that has ever been detected in the Baltic Sea population (Brown and Stepien 2009).Since the 1990s, the species has spread along the Baltic coastline and is now found in several major ports, from where transport to the lower Rhine is possible.Rhine N. melanostomus could also have originated from Ukrainian coastal waters.Mombaerts et al. (2014) recently detected haplotype ame18, typical for N. melanostomus from the Sea of Azov, the River Southern Bug and the northern Black Sea (Brown and Stepien 2009), in gobies from the River Scheldt, the Albert Canal and the Gent-Terneuzen Canal, all of which have a connection to the Rhine.Despite this, the chances of "Ukrainian" lower Rhine gobies being the source for the Belgian samples is relatively low.The relatively high frequency of the ame18 haplotype in Belgian samples (five of 34 samples), and its absence in the Rhine (Cerwenka et al. 2014;Tserkova et al. 2015;this study;58 samples in total) or River Waal (which connects the Rhine with the Belgian sites; Mombaerts et al. 2014), suggests that Belgian gobies originated from an independent introduction via one or more of the large Belgian ports.Finally, a Danubian origin for Rhine N. melanostomus cannot be excluded.According to our STRUCTURE analysis, gobies from the lower Rhine appear to share a cluster distribution with those from the lower Danube, whereas Cerwenka et al. (2014) noted similar clustering between N. melanostomus from the lower Rhine and the upper Danube.Large differentiation between the Danubian and Rhine samples (as determined by FCA in our study) and the presence of private alleles in Rhine gobies, however, argues against a purely Danubian origin.
Probably the most likely scenario, supported by the presence of private alleles in the Rhine samples, would be the admixture of two or more source populations.Admixture arises through interbreeding between two or more previously isolated populations.This process usually requires multiple or parallel introduction events, which are common and ongoing in N. melanostomus and P. kessleri (Brown and Stepien 2009;Cerwenka et al. 2014;Adrian-Kalchhauser et al. 2016).Admixed populations may benefit from increased genetic variation, the creation of novel genotypes and the masking of deleterious mutations (Verhoeven et al. 2011), all especially beneficial during range expansion (Keller et al. 2014).The Rhine ports of Rotterdam (Netherlands) and Duisburg (Germany) are among the largest and busiest inland ports in the world and, as such, are good candidates for sites where admixture could have taken place.Large ports are not only natural hot spots for admixture of invasive aquatic species (Adrian-Kalchhauser et al. 2016); they also play an important role in the onward spread of non-native admixed populations to novel locations.

Colonisation of Danubian tributaries and the Elbe
The upper Elbe N. melanostomus population became established in 2015 (Roche et al. 2015), almost certainly originating via ship transport from the river mouth population established just a few years earlier near the port of Hamburg (Hempel and Thiel 2013).Both STRUCTURE and FCA show differences between samples from the upper Elbe and the Danube, with Elbe samples also displaying private alleles.While our data do not allow us to identify the definitive source of the Elbe sample, likely candidates are the Baltic or an admixture from more than one source.Given that the Elbe sample was taken from a very young and isolated population, its relatively high genetic variability is unusual, either mirroring high genetic variability in the source (supporting the admixture hypothesis) and/or high propagule pressure (less probable, given the relatively infrequent shipping traffic to the Czech stretch of the Elbe).
Colonisation of the non-navigable Morava by both P. semilunaris and N. melanostomus most probably took place through natural upstream dispersion (i.e.swimming) from the middle (Austrian) Danube, as supported by the similarity between Moravian and Danubian samples in FCA (Figures 3A and 3B).An apparent decrease in genetic diversity suggests that the number of Austrian colonisers entering the Morava is relatively small.This would also explain the slight genetic differences between the Morava and middle Danube discerned by STRUCTURE (mostly in P. semilunaris) and F ST , which can thus be attributed to genetic drift.
The Dyje sample showed the greatest genetic difference in the STRUCTURE analysis, formed uniform clusters in both STRUCTURE and FCA, had the lowest genetic diversity and showed no private alleles.These results strongly suggest a founder effect (i.e.genetic distinctiveness caused by random genetic drift), which is supported by the species' known invasion history.The first occurrence of P. semilunaris in the Dyje basin (1994) was restricted to a small isolated population in the Nové Mlýny reservoir, 42 km from the Dyje's confluence with the Morava and more than 70 km from the "natural" invasion front that year (Lusk and Halačka 1995;Prášek and Jurajda 2005).In fact, the population is believed to have originated from an angler's bait-bucket introduction, comprising just a few individuals from the nearby Danube or Morava (our FCA suggests the Dyje to be a subset of the Morava sample).

Summing up and looking forward
While the distribution patterns of the three species examined are relatively similar, they each originated through different invasion histories.Both Rhine and Elbe N. melanostomus and Rhine P. kessleri appear to have experienced multiple introductions and strong propagule pressure, which was clearly not the case for Rhine P. semilunaris or with gobiids that have colonised tributaries.Adult P. semilunaris are much smaller than adult N. melanostomus and P. kessleri so are better adapted to streams, backwaters and other habitats away from the main river stem.Hence, there is less chance for the species to be taken up with the ballast water of large cargo ships, and thus a reduced probability of multiple introductions via this vector.Such transfers can and do occur, however, as documented by their introduction into the Rhine and/or North American Great Lakes.On the other hand, given the species' habitat preferences and inferior competitive ability (Valová et al. 2015), arrivals from subsequent introductions would be unlikely to communicate with those already established.
In contrast to P. semilunaris, N. melanostomus has established stable and abundant populations in harbours and main river stems, and frequently dominates local fish assemblages.The probability of onward transport via cargo ship is therefore high.This is clearly reflected in the species' invasion history, which strongly suggests multiple and/or parallel introductions and frequent transmission between ports (see also Brown and Stepien 2009;Cerwenka et al. 2014), along with genetic admixture, especially in the Rhine and Elbe.While N. melanostomus has successfully colonised non-navigable tributaries like Morava and Dyje, also reaching high abundances and commonly dominating local fish assemblages (see e.g.Janáč et al. 2016), our results show that natural colonisation of these rivers may have resulted in genetically less diverse (i.e. more vulnerable) populations than those introduced via navigation to the large rivers.
Maintenance of high genetic diversity (compared to the source population) in P. kessleri from the middle Danube and Rhine reflects strong propagule pressure.In contrast to N. melanostomus, P. kessleri shows no evidence of population structure or admixture.There are two possible reasons for the absence of admixture.First, all non-native populations studied to date have apparently originated from the same lower Danubian source.Second, newly introduced P. kessleri tend to exhibit high abundance in the first years of expansion only, numbers dropping rapidly over the following years (Cerwenka et al. 2014;Adrian-Kalchhauser et al. 2016;Borcherding et al. 2016).As such, the probability that ships regularly transport P. kessleri between non-native "populations" is low, as is the chance for gene flow.
The invasion histories for P. semilunaris, N. melanostomus and P. kessleri put forward in this paper are based on our best interpretation of the observed genetic patterns, dates of first record and known population characteristics and habitat preferences.Similar patterns could potentially arise via other means.For example, differences in genetic diversity could be attributable to differences in species plasticity, or strong selection after arrival causing genetic uniformity instead of bottleneck events during introduction.In the absence of any hard evidence for such alternatives, however, we believe our interpretations to be the most likely under the circumstances.
Range expansion of invasive Ponto-Caspian gobies in European freshwaters is an ongoing process, with colonisation of new locations being reported annually.In addition to identifying their origins, future studies should focus on the possible effects of admixture on invasion success.All introduced N. melanostomus sampled so far in Europe only contain haplotype ame1, typical for all populations from the Black Sea basin (Brown and Stepien 2008).This tends to hamper any effort to reveal the invasion history of introduced populations using mtDNA markers.To further elucidate the source(s) of Rhine N. melanostomus, we suggest analysis of samples from North Sea and Baltic Sea ports alongside samples from the Danube, and comparison of specific scenarios for pure and admixed colonisation, e.g. by using approximate Bayesian computation approaches (see Konečný et al. 2013).In general, the use of additional marker types, such as singlenucleotide polymorphisms (e.g. by RAD-sequencing) would help clarify the invasion histories of all invasive gobiids.

Figure 1 .
Figure 1.Map of localities sampled.DAN-BG = Bulgarian Danube, DAN-AU = Austrian Danube, RHI = Rhine, MOR = Morava, DYJ = Dyje, ELB = Elbe, PS = Proterorhinus semilunaris, NM = Neogobius melanostomus, PK = Ponticola kessleri.Blue species labels indicate the species is indigenous to the site, while red labels indicate non-native species.The red triangle marks the position of the Iron Gates dam.Dotted lines indicate channel connections between rivers (i.e. between the Danube and Rhine and the Rhine and Elbe).

Figure 2 .
Figure 2. Genetic population structure for K = 3 estimated in STRUCTURE from 118 Proterorhinus semilunaris (A) and 146 Neogobius melanostomus (B) from five localities.Each individual is represented by a vertical line partitioned into three coloured segments, the length of each colour being proportional to the estimated membership coefficient (Q).White lines separate different sample locations.DAN-BG = Bulgarian Danube, DAN-AU = Austrian Danube, RHI = Rhine, MOR = Morava, DYJ = Dyje, ELB = Elbe.Numbers in parentheses designate the number of years after first detection of the species (N = native population).Note that we do not provide a plot for P. kessleri as no population structure was observed in this species (estimated K = 1).

Figure 5 .
Figure 5. CYTB haplotype network for Ponticola kessleri.The diameter of the circles corresponds with sample size and different colours represent different sampling locations.DAN-BG = Bulgarian Danube, DAN-AU = Austrian Danube, RHI = Rhine.For more details (including GenBank accession numbers), see Appendix 3. Note that the best model in STRUCTURE indicated no structure at all (i.e.K = 1).
; Appendix 3).Two haplotypes from this study were recorded in the Austrian Danube and Rhine samples, with Hap_1 dominant and Hap_5 less frequent in both areas.Interestingly, Hap_1 was not found in the Bulgarian sample (native area), though it contained the other three haplotypes.Hap_1 has previously been recorded in the Serbian Danube (GenBank no.FJ526770; Neilson and Stepien 2009b).Haplotypes Hap_2 (previously reported from the Dniester; Neilson and Stepien 2009b) and Hap_6 (previously reported from the Odessa region; Medvedev et al. 2013) were not recorded in our study.

Table 1 .
Overview of samples used for genetic analysis.Sample size = number of fish genotyped for microsatellites and cytochrome b (in parentheses); Years after detection = number of years between first report of the species and the year sampled.
Compared to samples from the Danube, non-native P. semilunaris from the Rhine, Morava and Dyje and N. melanostomus from the Morava showed lower diversity overall, with lower AR, lower uHe and no private alleles (Table2, Appendix 2).In contrast, non-native N. melanostomus from the Rhine and Elbe displayed higher genetic diversity than those from the Danube, especially as regards AR (Table2, Appendix 2).Both the Rhine and Elbe samples also possessed private alleles.Rhine P. kessleri showed similar genetic diversity to those from the Danube, with similar values for both AR and uHe, and private alleles in all three samples (Table2).