Strong versus weak population genetic differentiation after a recent invasion of gobiid fishes ( Neogobius melanostomus and Ponticola kessleri ) in the upper Danube

Approximately ten to 15 generations after first inoculation, two invasive goby species Neogobius melanostomus and Ponticola kessleri have dispersed and established rapidly the upper Danube River. Population genomic amplified length polymorphism (AFLP) data show that the genome of the more recent newcomer, i.e. the globally invasive N. melanostomus, is significantly differentiated to a comparatively large degree (~ 5%) and exhibits pronounced small-scale population structure along a recently invaded 200 km river section. MtDNA haplotype identity over N. melanostomus samples suggests that an admixture of phylogenetically strongly differentiated source populations is unlikely. Fine-scaled local genetic population structure of N. melanostomus as deduced from Bayesian assignment tests suggest a trisection of the upper Danube instead of a clinal pattern: one downstream sample is assigned together with distant upstream samples to one population cluster. A second cluster comprises central samples, whereas two samples from the margins of this central region appear to have mixed ancestry. AFLP genome scan results indicate this population structure is strongly correlated with extrinsic (geographic) parameters, i.e. migration barriers of anthropogenic origin. However, divergence of at least one AFLP locus correlates positively with a proxy for trophic differentiation, i.e. variation of white muscle δN stable isotope signature. In contrast to N. melanostomus, no significant population differentiation was detectable in P. kessleri along the analyzed invasion pathway. In genome scans of P. kessleri, variation of a single locus is strongly positively correlated with an extrinsic parameter combination but not with any ecological parameter.


Introduction
Invasive species, by definition, arrive, establish and spread in novel environments within very short time frames (Keller et al. 2011).The study of the evolutionary dynamics of invasions may yield inferences about the causal factors leading to invasion success, a key topic in invasion biology.In addition, invasions may reveal otherwise intractable insights into evolutionary processes, e.g. the first and difficult-to-observe steps towards population differentiation and speciation (Hendry et al. 2000;Reznick and Ghalambor 2001).Whereas the effects and consequences of invasive species on natural communities have been studied in great detail (e.g.Gozlan et al. 2010; for a summary see Sanders et al. 2003), the evolutionary biology of alien species themselves is increasingly receiving attention (e.g.Sakai et al. 2001;Lambrinos 2004;Hastings et al. 2005;Dlugosch and Parker 2007;Hänfling 2007).
Contemporary evolution and invasion success is hypothesized to be shaped to a large degree by intrinsic characteristics of the source populations (e.g.number and genetic constitution of introduced specimens) (Lambrinos 2004;Björklund and Almqvist 2010), but also effects of propagule pressure (Allendorf and Lundquist 2003;Lockwood et al. 2005;Colauti et al. 2006), inbreeding (Nei et al. 1975;Young and Seykora 1996), phenotypic plasticity (Parker et al. 2003), life history traits (Tsutsui et al. 2000) and migration and dispersal abilities (Sakai et al. 2001;Phillips et al. 2006).In particular, genetic bottlenecks and founder effects can promote but also restrict the speed of adaptive evolution (Tsutsui et al. 2000;Colautti et al. 2004;Stepien and Tumeo 2006;Prentis et al. 2008).These effects are not apparent in every invasion (Stepien et al. 2005) and the genetic variability of invasive populations can even exceed the one of the source population (Lockwood et al. 2005).Invasive populations evolve under novel and diverse extrinsic conditions, which may differ not only in their ecology but also in the degree of natural or anthropogenic habitat fragmentation and connectivity (Ricciardi and MacIsaac 2000;Lambrinos 2004;Bronnenhuber et al. 2011).The interplay between population-intrinsic and environmental factors may fuel or delay the rate of spatial and/or adaptive diversification by changing locally divergent standing genetic variation (Kolbe et al. 2008;Stepien et al. 2005;Mitchell-Olds and Schmitt 2006;Novak 2007).Admixture of different native stocks in a single novel inoculation site or after secondary contact of previously allopatric invasive populations may lead to increased standing genetic variation upon which natural selection might act and may lead to local adaptation (Verhoeven et al. 2011).Thus, range expansion, admixture and/or population genetic diversification are often concurrent (Kolbe et al. 2008;Olivieri 2009) and common in invasions.Especially in human-mediated introductions (Kolbe et al. 2004;Therriault et al. 2005;Roman and Darling 2007), and at frontend expanding sites (Price and Sol 2008) differentiation is widespread and can occur within short time frames.
Over the last two decades, the upper Danube River in Germany (Figure 1) has been invaded by numerous invasive species mostly originating from the Ponto-Caspian region (Gollasch and Nehring 2006).Among those, invasive goby species (Teleostei: Gobiidae) have reached the region probably in ballast water of freight vessels commuting between the lower Danube (Black Sea) region and the lower Rhine (North Sea) along the Rhine Main Danube canal (RMDcanal) (Wiesner 2005).Bighead goby, Ponticola kessleri (Günther, 1861) was first recorded in the central part of the upper Danube close to the city of Straubing in 1999 (Seifert and Hartmann 2000).Five years later in 2004, round goby, Neogobius melanostomus (Pallas, 1814), was detected simultaneously next to the port town of Passau in the lower reach of the Danube in Germany and again in Straubing (Paintner and Seifert 2006).Both, N. melanostomus and P. kessleri have a similar ecology (Eros et al. 2005) and their expansion has been fast and successful in terms of fish densities (Brandner et al. 2013a, b) although they have low natural migration rates and small home ranges (Ray and Corkum 2001;Brownscombe andFox 2012, 2013).
Neogobius melanostomus is a globally invasive species, which has expanded its range rapidly with and without anthropogenic support (Stepien and Tumeo 2006;Bronnenhuber et al. 2011;Kornis et al. 2012), whereas P. kessleri is restricted to central and eastern Europe (Ahnelt et al. 1998;Borcherding et al. 2011;Brandner et al. 2013a;Kalchhauser et al. 2013).On a local scale in the upper Danube River, both species have colonized the whole 200 km stretch between Passau and the most recently (2010) invaded uppermost area at the junction of the Danube with the RMD-canal, hereby providing a link between invasive populations of the Rhine and the Danube River (Brandner et al. 2013c;Cerwenka et al. 2014).
Population genetic analyses of N. melanostomus from its native Ponto-Caspian as well as from its invasive Eurasian and North American ranges identified two major native mtDNA-lineages from the Caspian and the Black Sea drainages, each with a high intralineage genetic diversity according to microsatellite results (Stepien et al. 2005).Invasive populations from the middle Danube River in Serbia and Slovakia appeared most closely related to a population sample from Odessa (Black Sea drainage).On average, invasive populations exhibited comparatively low levels of genetic diversity, except for the upper Volga population, which contained haplotypes from both divergent lineages (Brown and Stepien 2008).Populations from the upper Danube River or the River Rhine have not yet been investigated, and population genetic data for P. kessleri are not available yet.
The initial phase of invasions is short and few population genomic studies have assessed the correlation of invasive population differentiation with spatial and environmental factors and genetic admixture (Sakai et al. 2001;Lee 2002;Kolbe et al. 2004).A comparative approach assessing patterns of population genomic differentiation of two or more invasive species under identical eco-geographical settings could elucidate the relative contribution of intrinsic versus extrinsic ecological and/or geographical factors to invasion dynamics.The present study was designed to compare the dynamics of two simultaneous goby invasions along a small-scale two-dimensional river continuum.The two species are sympatric throughout the investigated river stretch with both recent ("leading edge") and comparatively old ("established") inoculations.The Danube River is intersected by migration-barriers potentially facilitating a rapid built up of locally adapted populations.Using population genomic (Amplified Fragment Length Polymorphism (AFLP)) and to a smaller extent mtDNA data, we investigate baseline population genetics, ancestry and admixture of the upper Danube goby populations and assess the general hypothesis that intraspecific differentiation of two sympatric invasive goby species has developed on a small geographical scale in about ten generations after first introduction.We further hypothesize that the globally less successful species, P. kessleri is characterized by a less pronounced local population structure and is therefore less potent for the evolution of local genetic adaptation.In contrast, we expect the highly invasive N. melanostomus to exhibit increased local population differentiation correlated with both, barriers to gene flow and ecological parameters.

Materials and methods
Field sampling, environmental and specimen data 1053 goby specimens (471 P. kessleri and 582 N. melanostomus) were sampled from October 2009 to October 2011 at ten river stretches (stretch 1 (downstream) to stretch 10 (upstream) along the upper Danube River, Figure 1) and a single site at the River Rhine (near Rees) (Table 1).All specimens were collected and processed using an electro-fishing gear following the procedure described in Brandner et al. (2013a, b).Pectoral fin clips were preserved in 96% ethanol p.A. for genetic analyses.All specimens and tissue vouchers were stored in the ichthyology collection of Bavarian State Collection Munich (ZSM).On each sampling site and of each species, two males and two females were collected with a selected total length (L T ) of 8-12 cm.In the field, L T was measured to the nearest 0.1 cm and sex was determined externally and later verified in the laboratory (for a detailed description, see Brandner et al. 2013a).For each specimen (i) environmental and (ii) specimen-specific parameters were recorded.Environmental parameters were: (i) stretch 1-10 (i.e.sampled localities which despite of missing information on origin and gene flow are subsequently referred to as "populations"; Figure 1), (ii) interdam population (ip) 1-6 (i.e.combination of stretches separated by artificial dams limiting upstream migration; ip 1 = stretch 1, ip 2 = stretches 2 to 6, ip 3 = stretch 7, ip 4 = stretch 8, ip 5 = stretch 9, ip 6 = stretch 10; Figure 1), (iii) the distance measured in km from the lowermost stretch (river-km), (iv) left or right bank of the Danube River (bank side), (v) gravel or rip-rap substrate (habitat type), and (vi) densities of con-and heterospecific gobies at each sampling site (competitors).Specimen parameters were (i) the number of acanthocephalan parasites (parasite), (ii) the stable isotopic signatures of muscle tissue (δ 15 N and δ 13 C, determined as in Brandner et al. 2013a) and (iii) body shape measurements (i.e.principal components (PC) of geometric morphometric analyses calculated as in Cerwenka et al. 2014; for N. melanostomus PC 1-9 and for P. kessleri PC 1-10).Table 1 summarizes the number of specimens sampled from each stretch differentiated according to specimen and ecological parameters.The input data for analyses of N. melanostomus and P. kessleri specimens are provided in the Supplementary Material (Appendix 1).

DNA-extraction, AFLP-genotyping and mtDNA sequencing
DNA from 0.4-0.6 cm 2 pectoral fin tissues of samples from 2009 and 2010 was extracted using the Genomic DNA from Tissue Kit (Macherey-Nagel) and of samples collected in 2011 using the DNeasy® Blood and Tissue Kit (QIAGEN).Subsequently, AFLP were detected following Vos et al. (1995) modified by Herder et al. (2008).
In order to control for systematic errors due to sample position on 96 well microtiter plates (seven plates for N. melanostomus and six plates for P. kessleri specimens) in downstream genetic analyses, samples were processed plate by plate.Each plate contained specimen DNA of one goby species from all eleven different localities with the full range of environmental variation (bank side, habitat type).The season and the year of sampling were present on each plate as far as possible.Samples were AFLP-genotyped with six restrictive amplifications using an ABI 3130 capillary sequencer (PE Applied Biosystem, Foster City, CA, USA) and ROX 500 XL as internal size standard.The primer combinations were EcoAGG/MseCTG (Albertson et al. 1999), EcoACA/MseCAA (Barluenga et al. 2006;Albertson et al. 1999), EcoACA/MseCTG (Barluenga et al. 2006), EcoACT/MseCAA (Geiger et al. 2010), EcoAGG/MseCTC (Geiger et al. 2010) and EcoACC/MseCTA (Geiger et al. 2010), fluorescently labeled with HEX and FAM.Bin sets were created with Peak Scanner TM Software Version 1.0 (Applied Biosystems) and peaks were automatically selected and scored using tinyFLP (Arthofer 2010) with modified adjustments following Geiger et al. (2010).Four N. melanostomus and six P. kessleri specimens were replicated on each well of the corresponding species after the DNA-extraction.
The complete cytochrome b gene (cytb: 1138 base pairs (bp)) and partial sequence of the threonine tRNA gene (66 bp) was amplified and partially sequenced for a representative subset of N. melanostomus samples (n = 28) from all stretches (four samples from stretch 6; three samples from stretches 10, 8, 7, 5, and 4; and two samples from stretches 3-1).Primers L14724 (Meyer et al. 1990) and H5 (Iwata et al. 2000) were used to PCR amplify a single fragment in 10 µl volume with 5 µl Multiplex Mix (QIAGEN), 1 µl genomic DNA, 0.8 µl of each primer (2.5 nmol), Q-Solution (QIAGEN) and HPLC water.The PCR temperature profile was: 94°C initial denaturation (120 s); 35 cycles with denaturation at 94°C (45 s), annealing at 52°C (30 s) and extension at 72°C (60 s); final extension at 72°C (180 s).PCR products were purified using ExoSAP-IT (USB) and were diluted in 10-20 µl HPLC water.Cycle sequencing was performed using Big Dye 3.1 (Applied Biosystems) with the internal sequencing primer L15066 (Brown and Stepien 2008), and products were electrophoresed and read using an ABI 3130xl DNA sequencer (Applied Biosystems).Sequences were edited using BioEdit v.7.05.3 (Hall 1999) after a preliminary alignment using default settings of ClustalW algorithm (Larkin et al. 2007).Resulting cytb haplotypes were compared for sequence identity with 81 cytb haplotypes from the native and introduced range of N. melanostomus (Brown and Stepien 2008).Analyses for cytb haplotypes were not conducted for P. kessleri.

Population genomic analyses
Plate-specific effects were reduced by the following pairwise comparisons of peak frequencies after binning.Histograms were computed to visualize differences in frequencies between fragments with the same number of bp using PAST 2.15 (Hammer et al. 2001).All fragments showing higher values of differences than indicated by the first minimum of the according Kernel density were excluded from the following comparisons and the subsequent analyses.Furthermore, following Collin and Fumagalli (2011), all fragments not present on replicated individuals were excluded.This deletion decreases plate-specific effects and increases the likelihood of detecting potentially masked divergence (see Geiger et al. 2010).This procedure results in a comparatively large number of weakly amplified low-frequency AFLP loci being excluded, but is conservative with regard to controlling for type II error based on platespecific systematic error.Mariette et al. (2002) and Singh et al. (2006) proposed ~200 bands to be sufficient for measuring population genetic variation and differentiation (for a summary see Bonin et al. 2007).Population based pairwise genetic differentiation was measured using analysis of molecular variance (AMOVA) in GenAlEx 6.41 (Peakall and Smouse 2006) and partitioned in among and within population differentiation.Pairwise ΦPT values (after sequential Bonferroni correction, 9999 permutations) are analogous to F ST -values but applicable to haploid markers, and indicate levels of hierarchical genetic differentiation among populations (Huff et al. 1993).Single pairwise comparisons showing negative values for ΦPT were converted to zero.Differentiation according to river stretch and ip was assessed in additional AMOVAs.
To test for isolation by distance (IBD) versus a potential influence of anthropogenic barriers to migration (i.e.dam), three independent approaches were used.Population differentiation (stretch and ip) was compared using (i) AMOVAs, (ii) pairwise ΦPT values and (iii) matrix comparisons of pairwise ΦPT values and associated geographic distances between sampling stretches (river-km) using Mantel tests in PAST 2.15.The longest continuous free-flowing section of the upper Danube River includes five river stretches (i.e.stretches 2-6) not disrupted by dams.Pearson's correlation coefficients (R) between ΦPT (population differentiation) and river-km were calculated within this stretch as well as within our down-and upstream ("sliding window") stretch sections of similar length but disrupted by one or more dams (i.e.stretch 1-5, 3-7, 4-8, 5-9, 6-10).Levels of significance for correlation coefficients were computed using 10000 random permutations.These values should be comparable among all six stretch-sections in case of only IBD determining population differentiation, whereas values of stretch-sections interrupted by dams should be higher if dams contributed to population differentiation in addition to geographic distance.
Population genetic structure and AFLP loci linked to genomic regions potentially under divergent selection were further examined using a combination of logistic regression and F SToutlier based methods.First, F ST -outlier loci were identified using the DFDIST algorithm (Beaumont and Balding 2004) as implemented in the workbench MCHEZA (Antao and Beaumont 2011), as well as with BAYESCAN (Foll and Gaggiotti 2008), both preferentially used for AFLP data (Pérez-Figueroa et al. 2010).In a second step, Bayesian factors were calculated for every marker using BAYESCAN which is based on logistic regressions.Following Mattersdorfer et al. (2012) the threshold to reject the null hypothesis of log 10 (BF) was set to a value smaller than 0.5, all other adjustments were applied using default settings.Candidate loci identified by DFDIST and BAYESCAN were compared with each other as well as with the selection of candidate loci identified by logistic regressions associated to environmental and specimen-specific parameters using MatSAM Version 2Beta (Joost et al. 2008).Here, potential genetic differentiation was correlated with population-genetic independent environmental or specimen-specific parameters shaping divergent selection; it was tested locus by locus, based on likelihood ratios using "G" and "Wald" tests and the "Cumulated test" (Hosmer and Lemeshow 2000).All p-values were adapted using the conservative Bonferroni correction.The more robust Cumulated test was only performed if both, G and Wald were significant (Joost et al. 2007).
An assumption free, individual based Bayesian algorithm (Falush et al. 2003) implemented in STRUCTURE 2.3.4 (Pritchard et al. 2000) was used to identify genetically distinct population clusters.Individuals were assigned to K populations without prior information on their origin.K = 1 to K = 11 were assessed, each with nine independent replicates with 400000 MCMC-iterations and a burn-in-value of 200000.STRUCTURE uses a model-based multivariate analysis and a Bayesian approach under the assumption of Hardy-Weinberg or linkage disequilibrium within each population.However, STRUCTURE appears robust with regard to violations of this assumption (Falush et al. 2003).Calculations were performed using the Bioportal computer service of the University of Oslo (http://www.bioportal.uio.no;Kumar et al. 2009).The most likely number of populations (K) was estimated following Evanno et al. (2005) using STRUCTURE HARVESTER (Earl and vonHoldt 2012) and the method proposed by Pritchard et al. (2000).To display results graphically CLUMPP version 1.1.2(Jakobsson and Rosenberg 2007) and distruct version 1.1 (Rosenberg 2004) were used.To reveal the potential influence of selection on spatial genetic structuring, all calculations and graphical illustrations were performed once using all loci and once without the candidate loci detected by MatSAM, DFDIST and BAYESCAN.
To test for "surfing allele" candidates (sensu Manel et al. 2009) allele frequency distributions were screened under two criteria: (i) they should be outliers detected by DFDIST but not by BAYESCAN and MatSAM (their frequency deviation should not correlate with extrinsic parameters), and (ii) allele frequencies at "leading edge" populations from stretch 8, 9 or 10 according to Brandner et al. (2013a) should be significantly increased as pairwise compared to "established" populations.Band presence frequencies between ten Danubian "populations" were tested for differences by multiple pairwise nonparametric Kruskal-Wallis tests (Bonferroni corrected).

Neogobius melanostomus
All 28 partially sequenced cytochrome b haplotypes (862 bp) as well as the three nearly completely sequenced (1,115 bp) haplotypes were identical to each other on the 862 bp stretch as well as to the most common "Black Sea Table 2. Number of clusters (K) of N. melanostomus (upper part of every cell) and of P. kessleri (lower part of every cell) including all loci (all) and excluding the loci potentially under selection (without) detected by DFDIST, BAYESCAN and MatSAM.Clusters are inferred by Ln P(X|Y) with its SD over 9 runs following the method of Pritchard et al. (2000) and ∆K, the second order rate of change of Ln P(X|Y) proposed by Evanno et al. (2005)  basin" haplotype 1, i.e. there is no indication of mtDNA admixture of phylogeographically different groups in the upper Danube River.
After correction for potential plate-specific effects, the number of detected polymorphic AFLP-bands was 189.Individual band frequency ranged between 7 and 100% (mean = 19%).No fragment occurred with frequencies lower than 5% and 29 fragments were present in more than 95% of all individuals.
STRUCTURE analyses for K = 1 to K = 11 revealed two or seven genetically distinct clusters being most likely.Applying the method of Evanno et al. (2005) it was K = 2 for the complete dataset, whereas after excluding the three candidate loci potentially under divergent selection detected by MatSAM, DFDIST and BAYESCAN, the number increased to K = 7 (Table 2, Figure 2).One of the two major clusters combined individuals from disjunct regions, i.e. the uppermost three populations (stretch 8-10), the River Rhine (stretch 11) and the lowermost population (stretch 1); in contrast, Table 3. Pairwise ΦPT of goby populations at river stretches according to species.N. melanostomus are indicated in the upper right part of the table, whereas P. kessleri in the lower left part.Stretches 1 to 10 are part of the Danube River and stretch 11 of the River Rhine.Nonsignificant differences are indicated by n.s., significant differences according to p < 0.05 by *, p < 0.01 by ** and p < 0.001 by ***.Comparisons with negative ΦPT were converted to zero.most individuals of the central populations (stretch 2-7) where assigned to the second cluster.Removal of the three candidate loci resulted in a shift of individual assignments to populations of the central stretch.Individuals of "populations" at stretch 1, 8 and 9 were almost entirely assigned to the central cluster, and about 20% of the individual genetic constitution of the uppermost Danubian "population" at stretch 10 was assigned to the cluster from River Rhine (before removal they shared almost 100% of the loci, Figure 2).The majority (95%) of the genetic variance was explained by within Danubian "population" structure (i.e. by stretch), whereas 5% was explained by among stretch.Overall ΦPT was 0.05 (Danubian specimens only: ΦPT = 0.045) which is significantly different compared to the variance calculated from randomly generated data (AMOVA all: p < 0.05, Danubian "populations" only: p < 0.001).Individuals from River Rhine (mean ΦPT = 0.176) showed highest levels of differentiation (mean ΦPT of several pairwise stretch-based comparisons ranged between 0.036 at stretch 2 and 0.085 at stretch 9).Without the "population" from the River Rhine mean ΦPT was smaller and ranged between 0.025 at stretch 2 and 0.076 at stretch 9).All pairwise comparisons were significant (with and without specimens from River Rhine, using Bonferroni corrections and 9999 permutations: all p < 0.05), after exclusion of comparisons showing negative ΦPT values and with exception of comparisons between stretch 6 and stretch 3 to 5 and of comparison between stretch 4 and stretch 6 (all p > 0.05).Detailed results are given in Table 3.
Candidate loci potentially under selection were identified using three population genetic analysis tools.BAYESCAN identified a total of 15 candidate loci under the criterion of log10(BF) greater than 0.5.Following Jeffreys' interpretation of the Bayes factors, five of those loci were under "strong" selection, one under "very strong" and nine were indicated to be under "decisive selection".
DFDIST detected outliers potentially being under positive selection and outliers possibly having a balancing effect by comparing the calculated F ST -values with simulated ones under neutral conditions.For N. melanostomus the overall calculated F ST was 0.045, the simulated one was 0.044.DFDIST suggested 28 loci as candidates under balancing and 9 under divergent selection.
MatSAM logistic regressions of extrinsic parameters provided an independent evaluation of BAYESCAN and DFDIST outliers.Seven loci were detected to be potentially under selection.Five loci were associated with parameters related to spatial heterogeneity and barriers to migration: one locus correlated with the parameters stretch, river-km and interdam-population, and four loci to interdam population.In addition one outlier was assigned to isotopic signature of δ 15 N and one to δ 13 C.
The total number of outliers, thus being candidate loci under divergent selection recognized by all three methods (MatSAM, DFDIST and BAYESCAN) was three.One locus was found to be correlated with all three tested spatial parameters stretch, river-km and interdam population, one single with the spatial parameter interdam population, and one single with δ 15 N isotopic signature.
A single AFLP-locus potentially surfing at the leading edge of N. melanostomus of the upper Danube River (i.e.stretch 10) was detected being present in 32% of all analyzed individuals at stretch 10.This locus was not present in individuals at stretch 3-5, 7 and 8 and a significantly lower number of individuals had this fragment at stretch 6 (p < 0.001) and stretch 9 (p < 0.05).Higher frequencies were observed at stretch 2 (12.5%) and at stretch 1 (15%).

Ponticola kessleri
After correction for potential plate effects, the number of detected polymorphic AFLP-fragments was 372 and individual band frequency ranged between 1 and 100% with a mean of 15%.280 fragments occurred with frequencies lower than 5% and 30 were present in more than 95% of all individuals.
In P. kessleri no clusters were recognizable for K = 1 -11 in STRUCTURE analyses.As the method of Evanno et al. (2005) cannot be used to estimate the number of clusters for the K extremes (1 and 11) the log probability of Ln P(X|K) was used following Pritchard et al. (2000).However it indicated K = 10 to be most probable (Table 2).Excluding the single candidate locus potentially under selection revealed by the three methods used (see below) the highest likelihood for the number of clusters was K = 2, applying the method of Evanno et al. (2005) and K = 3 using the log probability of Ln P(X|K) but K = 2, 6, 8, 9 and 11 showed also comparatively high values of ∆K (Table 2).In both cases with and without consideration of the single locus under selection no apparent spatial population structure could be detected (data not shown).
The main part of genetic variance (99%) was explained by "populations" structure (i.e.within stretch), whereas 1% was explained by among stretch.An exclusion of low sample size populations (n = 8 individuals from River Rhine and n = 5 from stretch 10) did not change the molecular variance of the remaining populations.ΦPT of all specimens was 0.008 and did not indicate significant population differentiation (AMOVA: p > 0.05).However, differentiation was significant when regarding Danubian specimens only (ΦPT = 0.007, AMOVA: p < 0.001).Mean ΦPT values of pairwise comparisons were highest for the River Rhine "population" (mean ΦPT = 0.034) and lowest for the one from stretch 4 (mean ΦPT = 0.006).In total, 32 ΦPT values of pairwise stretch comparisons were not significant (all: Bonferroni corrected, 9999 permutations, p > 0.05).Detailed results are given in Table 3.
BAYESCAN identified four candidate loci potentially under selection using the criterion of log10(BF) greater than 0.5.Following Jeffreys' interpretation of the Bayes factors, three of those loci were under "substantial" and one under "decisive" selection.DFDIST identified 79 loci potentially being under selection, 18 under positive and 61 under balancing selection.Logistic regressions using MatSAM detected a single locus indicating spatially controlled divergence, i.e. for the factors stretch, river-km and ip.All three methods identified this particular locus.
No locus could be identified to potentially surf at the leading edge of P. kessleri.
Cytochrome b haplotype analyses of native P. kessleri populations are still lacking and thus invasive specimens from the upper Danube River were not analyzed.

Discussion
Invasive organism structure evolves dynamically and may leave different signatures resulting from intrinsic and extrinsic factors.Population genomic AFLP data of invasive gobies in the upper Danube River show that the more recent newcomer, i.e. the globally invasive N. melanostomus, is significantly differentiated to a comparatively large degree (~ 5%) and exhibits pronounced small-scale population structure along a 200 km river section.Local genetic population structure of N. melanostomus suggests a trisection: one downstream sample is assigned together with distant upstream samples to a first population cluster, the central samples to a second one, and two samples from the margins of the central region appear to have mixed ancestry.Divergence of at least one locus correlates with a proxy for trophic differentiation, i.e. variation of white muscle δ 15 N stable isotope signature in this species.No significant population differentiation of P. kessleri is detectable, and in genome scans, variation of only one single locus was strongly correlated with an extrinsic, geographic parameter combination.
The comparison between P. kessleri and N. melanostomus in the upper Danube River highlights that rapid population differentiation in invasive organisms can be different under identical extrinsic settings.Apparently, the interplay of in-and extrinsic factors e.g. the number of inoculation events, propagule pressure, origin of invaders, and/ or potential genetic admixture acts differentially resulting in interspecifically different evolutionary responses.In addition, intrinsic factors as phenotypic plasticity or different levels of standing genetic variation may act and change the population genetic constitution and therefore the different population genomic basis of non-native species.

Origin of invasive genomic diversity
Invasion success and rapid population differentiation appear to act synergistically, and might enhance the speed of invasion (Grosholz 2002;Lee 2002;Björklund and Almqvist 2010).The genomic constitution of native population(s) potentially contributes to the success of invasive species (Mitchell-Olds et al. 2008;Geneva and Garrigan 2010).Theoretically, population differentiation may be enhanced, if allopatrically differentiated strains amalgamate into a new (invasive) population, characterized by an instantaneously elevated standing genetic variation (Lucek et al. 2010).Despite being significantly differentiated, mtDNA variation of N. melanostomus in the upper Danube is zero, as all analyzed individuals carried the same Black Sea basin haplotype.Therefore, rapid differentiation of this species is most likely not caused by an admixture of phylogenetically strongly distinct source populations from the Caspian and Black Sea basin (i.e. by a Wahlund effect; Björklund and Almqvist 2010).This has been suggested for an invasive N. melanostomus population of the Volga region (Brown and Stepien 2008).Investigations of native P. kessleri populations are still lacking.
However, since this result is based on matrilinearly inherited and comparatively slowly evolving mtDNA only, an admixture of related populations, even of a male Caspian contribution, cannot be excluded completely.Therefore it remains open, whether the observed rapid population differentiation in N. melanostomus is the result of multiple introductions of closely related but nevertheless pre-differentiated populations, as it has been shown by several studies in invasion biology summarized by Prentis et al. (2008) and Vellend et al. (2007).Support for this scenario comes from the interspecific comparison, because N. melanostomus showed considerably higher overall genetic variability than P. kessleri.An alternative explanation, i.e. decreased genetic variability due to genetic bottlenecks in P. kessleri potentially restricting differentiation (Stepien and Tumeo 2006), cannot be ruled out without comparative data for both species from different invasion regions as well as from source populations.Under this scenario, low effective population size and low levels of immigration (Fitzpatrick et al. 2011) would have contributed to lower values of population differentiation in P. kessleri.A higher vulnerability to inbreeding (Frankham 2005) might then explain the observed numerical decrease of P. kessleri over the years (own data), independent from effects of interspecific competition.

Factors correlated with invasive population differentiation
Population differentiation may be favoured by barriers to gene flow and by IBD (see Meldgaard et al. 2003).In invasive N. melanostomus and P. kessleri, genetic differentiation was shown across very short geographic distances, and factors correlating with population structure were mostly of geographic nature.Hence, local population structure evolved over very short time spans and few generations, and appears to be supported by extrinsic factors in the upper Danube River.(Anthropogenic) barriers to gene flow seem to be decisive for locally different success of both gobiid species, even if different source populations would be a major reason for population variation and differentiation.In conjunction with established subpopulations, population genetic effects may promote the rapid evolution of population structure, e.g."allele surfing".A rapid increase of previously low allele frequencies in expanding frontend populations (Klopfstein et al. 2006;Excoffier and Ray 2008) is hypothesized to be typical for invasive populations.Gobies, having an extended spawning period with males guarding nests aggressively and females spawning several times (Charlebois et al. 1997;Corkum et al. 2004;Groen et al. 2012), should be prone to allele surfing through the "Hedgecock effect" (see Hedrick 2005), where a low number of parental individuals have a high number of offspring.Allele surfing probably is likely to have contributed to local population structure of N. melanostomus and was detectable in a strongly differentiated, but yet comparatively young "leading edge" sampling site, i.e. at the uppermost stretch of the Danube River, which was seeded after the year 2009 or even 2010 (Brandner et al. 2013a, c).In contrast, no surfing allele could be detected in P. kessleri, underlining the lower genetic and phenotypic (Cerwenka et al. 2014) variability and the lower population differentiation in this species.

Proxies of trophic and morphometric differentiation as indicators of genomic adaptation
Natural selection may trigger local adaptation even at small geographic scales, but it is often difficult to identify single factors driving local adaptation (Collin and Fumagalli 2011).In addition, local adaptation can be confounded with genetic signatures of introgression if only F ST -based genome scans are used to identify candidate loci, because allochthonous introgressed alleles in a subsample may mimic selectively favored high allele frequencies (Gagnaire et al. 2011;Mattersdorfer et al. 2012;Gosset and Bierne 2013).However, in N. melanostomus genomic differentiation at a candidate locus has been identified not only on the basis of F ST genome scans but also by logistic regression against a proxy for niche segregation (δ 15 N isotopic signature), indicating genomic adaptation to alternative trophic niches.Nitrogen stable isotopes provide a temporally integrated, quantitative perspective on individual diet and are indicative of the relative trophic position of an individual.Although both goby species are generalistic omnivores (Borcherding et al. 2013;Brandner et al. 2013a), the differential variation of N. melanostomus in the upper Danube matches with results showing that this species exhibits a greater feeding niche width, and a lower degree of specialization than P. kessleri; this possibly reflects a higher degree of individual adaptation to available prey in comparison to N. melanostomus (Brandner et al. 2013a).
Nevertheless, the individual trophic niche positions of N. melanostomus specimens could be the result of a phenotypically plastic response.Neogobius melanostomus is indeed known for its broad diet and high feeding versatility, which could indicate phenotypic plasticity, as deduced from observations in different ecosystems (Gaygusuz et al. 2007).Phenotypic plasticity, where different phenotypes are expressed under different environments (Fitzpatrick et al. 2012), can be important in successfully invaded environments (Agrawal 2001).In N. melanostomus of the upper Danube River phenotypic plasticity is highly probable (Cerwenka et al. 2014).In addition, population differentiation based on geometric morphometric data revealed a similar geographic trisection into an upper, central and lower part of the River as in the genetic analysis.Habitat parameters and body shape variation were not identified as significant in logistic regression analysis, i.e. no single allele frequency was detected that significantly corresponded to proxies of morphometric PCs or habitat type.This renders phenotypic plasticity more probable, which is considered as an important "jump-starter" directly after inoculation (Collyer et al. 2007).It may drive subsequent directional evolution and it may facilitate rapid adaptive evolution leading to rapid success in invasive species.

Population expansion and spreading
Population structure of non-native species may differ according to expansion mechanisms (e.g.Currat et al. 2008;Bronnenhuber et al. 2011).Results of population structure analyses revealed assignment of P. kessleri specimens along the upper Danube River to only a single population which corresponded to comparatively low levels of overall genomic variability.This could either suggest a bottleneck situation at or after inoculation, or an already depauperate native P. kessleri stock.The N. melanostomus population trisection along the upper Danube River suggests disjunct inoculations from multiple founder populations.The most likely mode of inoculation is by transportation of eggs or larvae in ballast water vessels rather than by active migration since this species has small home ranges and limited migration rates in adults (Bronnenhuber et al. 2011;Gutowsky and Fox 2011;Brownscombe et al. 2012;Kornis et al. 2012).Population expansion proceeded in upstream direction, with dams acting as barriers to gene flow, in both species.Nevertheless, passive downstream drift of juvenile gobies, may explain the existence of a genetically intermediate population (stretch 7), which is located between the central (stretch 3 to 6) and the upper part of the upper Danube (stretch 8 to 10).Drift has been shown to be significant for invasive gobies (Hensler and Jude 2007;Hayden and Miner 2008;Björklund and Almqvist 2010;Janáč et al. 2013), but its importance is most likely underestimated as compared to active dispersal.Nevertheless, single N. melanostomus are known to move long distances at least in upstream direction as described by Kornis et al. (2012), Bronnenhuber et al. (2011) and Brandner et al. (2013c), and hereby could have contributed to an admixture of genetic clusters at least at stretch 2 and stretch 7. Multiple inoculations in combination with subsequent downstream drift and active dispersal may thus have contributed to population admixture and possibly to invasive success.
In conclusion, population differentiation and expansion, as well as factors correlating with it are clearly species-specific in our case.Despite highly similar invasion histories of N. melanostomus and P. kessleri in the upper Danube River their invasive populations respond differentially to spatial and ecological parameters.The species having a higher variability in life-history traits, phenotype and nutrition (N.melanostomus), responds to its novel non-native area by rapid population genomic differentiation on a local level.Neogobius melanostomus is by far the most successful of invasive goby species in terms of fastly establishing high density populations on a global scale.The correlation between rapid responses to locally different environments suggests a significant contribution of genomic adaptability to invasion success.Barriers to gene flow conducting to a subdivision of non-native populations increase rather than decrease the potential for local adaptation in "plastic invaders" as N. melanostomus.Apparently the less plastic and more inertially invader, i.e.P. kessleri is responding less flexible at the genomic level to these extrinsic factors.

Figure 1 .
Figure 1.River stretches 1 to 10 in the study area of the upper Danube River (within the Danube drainage, upper right part).Triangles crossing the river indicate hydroelectrical dams.Populations combined in some analyses to "interdam populations" (ip 1 to ip 6) are given in grey shaded circles.The beginning of the RMD-canal at the city of Kelheim is marked by an arrow.

Figure 2 .
Figure 2. Population structure of N. melanostomus identified using STRUCTURE Bayesian assignment analysis, were K is the number of clusters predefined.The uppermost graphs show the population and the lowermost the individual matrix, respectively.A: K = 2 including all loci (n = 189).B: K = 2 and C: K = 7, both without loci potentially under selection (n = 186) identified by MatSAM, BAYESCAN and DFDIST.Numbers of river stretches are indicated at the uppermost and interdam populations (ip) at the lowermost part of the figure.

Table 1 .
Stretches (i.e.names of localities with the corresponding River system, number of sampling locality and GPS-coordinates of downstream sampling site boundaries) and numbers of analyzed goby specimens separated according to specimen-specific parameters and ecological data.Data on N. melanostomus (n = 582) are given at the upper part of every cell and data on P. kessleri (n = 471) at the lower part, respectively.