Genetic Heterogeneity Reveals On-Going Speciation and Cryptic Taxonomic Diversity of Stream-Dwelling Gudgeons (Teleostei, Cyprinidae) in the Middle Danubian Hydrosystem (Hungary)

Although stream-dwelling gudgeons (Cyprinidae, genus: Gobio) are widespread in Central Europe, the taxonomy of this group and the distribution of its species are still unexplored in detail. The aims of our study are to ascertain taxonomic composition and distribution of the former Gobio gobio superspecies in the inner area of the Carpathian Basin. Since the presence of cryptic species is suspected in this area, we examined the taxonomic and phylogenetic relationships of Central European Gobio taxa by sequencing the mitochondrial DNA control region (mtCR). Additionally, we characterized the genetic structure of 27 stream-dwelling gudgeon populations of this area by Amplified Fragment Length Polymorphism (AFLP). Results of mtCR analysis proved the presence of three species already known as G. obtusirostris (dominant in NW-Hungary), G. gobio (sporadic) and G. carpathicus (sporadic). Additionally, the analysis revealed the existence of one doubtful taxon, G. sp1 (dominant in NE-Hungary), and a new isolated haplogroup (dominant in SW-Hungary). Although Network analysis showed significant detachment among haplogroups, their genetic distances were quite small. Therefore Bayesian phylogenetic analysis showed weak nodal support for the branching pattern both for newly described haplotypes, and for the already accepted species. AFLP data showed distinct population structure and a clear pattern of isolation was revealed by distance of stocks. At the same time, level of separation was not affected by the altitudinal position of sites. Moreover we found three major clusters of populations which were separated according to hydrographic regions, and corresponded to the findings of mtCR analysis. Our results suggest the on-going speciation of gudgeons in the Carpathian Basin, however the separation of haplogroups seems to only be an intermediate phase. The discovered natural pattern seems to be only slightly influenced by anthropogenic impacts. Additionally our results put into question the suitability of the recently accepted within Gobio genus taxonomy.


Introduction
Freshwaters are exceedingly diverse ecosystems, but at the same time they are extremely sensitive to habitat degradation and pollution [1]. Accurately quantifying their taxonomic and functional diversity is a fundamental requirement of conservation biological, ecological, biogeographical, and macroevolutionary research [2], [3], [4]. However, when species are not clearly distinguishable by the conventional methods using ecological and morphological traits, have highly similar environmental needs and reveal a high level of phenotypic plasticity, then estimates of biodiversity and ecosystem functioning will be biased [5], [6], [7], [8].
With the increasing use of molecular techniques, it has become evident that many species formerly believed to have widespread geographical distribution can in fact be divided into numerous more or less discrete entities -so called cryptic or sibling species by Mayr [9]-or represent genetic gradients between separating species (i.e. on-going speciation) [10], [11], [12], [13], [14]. Such phenomena are more likely to occur in organisms with limited dispersion ability and/or in organisms living in separated or narrowly connected habitats [15], [16], [17]. Many streamdwelling fishes have specific environmental needs and therefore form discrete populations, not only between geographical areas with separated catchment systems, but also between closely related sites of the same catchment [18], [19], [20], [21]. These isolated fish populations may then genetically differentiate with time, although they may still maintain their similar morphological appearance and ecological function [22]. However, this process is not yet fully understood in seemingly well connected catchment systems.
The type species of the Gobioninae subfamily (Fam: Cyprinidae), the common gudgeon Gobio gobio Linnaeus (1758), was known as the most widely distributed lentic gudgeon species in West Eurasia [23]. However, the high between and within population morphological variability of this taxon [24] resulted a long-standing debate regarding its taxonomical status [25], [26]. As a result, numerous forms/varietas/subspecies were described, nearly all from larger catchment areas within its range (e.g. [23]). Recent genetic studies have raised some of the former subspecies to species level and several new species have been described as well [27], [28] [29], [30], [31], [32]. However, the taxonomic and genetic status of gudgeon populations in the Carpathian Basin is still unclear. Gobio gobio was considered a common species in the waters of the Carpathian Basin for a long time [33], [34]. On the contrary, recent studies [31], [35] excluded the Carpathian Basin from the potential range of G. gobio and suggested the occurrence of the Danube gudgeon, Gobio obtusirostris Valenciennes (1842) in the western region of the basin, and the Carpathian gudgeon, Gobio carpathicus Vladykov (1925) in the drainage system of River Tisza. Furthermore, Mendel et al. [31] indicated the presence of a ''species'' (sic!: Gobio sp1) from the catchment of River Tisza, but the taxonomic position of this newly described taxon has not yet been firmly established. It is also important to note that the findings of Kottelat and Freyhof [35] are based exclusively on data from the literature. Moreover, although the study of Mendel et al. [31] is the most comprehensive genetic study on the Middle European Gobio species to date, it included only a very limited number of samples from the Carpathian Basin and furthermore all of those samples originated from the edges of this region. The above mentioned Gobio taxa (i.e. G. gobio, G. obtusirostris, G. carpathicus and G. sp1) are morphologically very similar [35], [36] and thus their distribution and ecology cannot be explored without molecular justification [24]. Furthermore, recent studies on G. gobio and the related species [37], [38], [39], [40] found remarkably high genetic and morphologic variability between and within catchment areas, supporting the likelihood of the presence of cryptic species.
This study aims to ascertain taxonomic composition and distribution of the former G. gobio superspecies in the inner Carpathian Basin. Specifically, we (i) examine the taxonomic and phylogenetic relationships of Gobio taxa and the presence of cryptic species, and (ii) characterize the genetic structure of gudgeon populations with special attention paid to the effects of hydrological distance and elevation as possible forces facilitating genetic separation.
To unravel the taxonomic relationships of stream-dwelling gudgeons inhabiting the central area of the Carpathian Basin we use the same methodology and molecular marker (sequencing the mitochondrial Control region) as was used by Mendel et al. [31], thus our results are comparable with their findings. Moreover, we also screened for Amplified Fragment Length Polymorphisms (AFLP's) as a supplementary method of analysis for genome-wide genetic variation [41].

Ethics Statement
This study was carried out following relevant national and international guidelines pertaining to the care and welfare of fish. Collections were made by electrofishing, partially from sampling sites which are located within protected areas. Moreover, each species within the Gobioninae subfamily is protected in Hungary. Electrofishing in protected areas and any procedure (collection and storage of tissue samples) to be applied to protected species are subject to authorisation in Hungary. Fin tissue collection and storage were approved by the National Inspectorate for Environment, Nature and Water, Hungary (permission numbers: 14/ 3714-2/2009, 14/1237/2/2010, 14/881/5/2011, 14/678-9/ 2012). Fish collected for this study were narcotized using clove oil. After fin tissue sampling, when they regained consciousness, they were returned to the wild. Field studies did not involve fish that were endangered (The IUCN Red List of Threatened Species v. 2013.1; www.iucnredlist.org).

Study Area
All sampling sites were situated in the inner area of the Carpathian Basin, which belongs to the drainage system of the Danube River. Based on its topographic characteristics, the Hungarian part of Middle Danubian hydrosystem can be divided into two larger catchments and ten smaller sub-catchments ( Table 1, Figures 1 and 2). The hydrography of this area is biaxial. From the western region, all watercourses empty into the Danube River. The eastern part of the country belongs to the drainage system of the Tisza River, which is the largest tributary of the Danube (157 000 km 2 catchment area). The structure of this drainage system is dendritic, with the Tisza River forming the central axis. All the studied watercourses connected to the middle section of the River Tisza [42], therefore this region was not differentiated further. The hydrography of the Danubian system is more complicated. This system consists of three comparatively isolated subsystems: North, Middle and South Danubian regions (Table 1). North Danubian region is formed by the drainages of River Rába (Raab), River Ipoly (Ipel) and by the drainages of some direct inflowing streams. Middle Danubian catchment originally joined to the River Danube through a marshy area. Until the construction of the Sió canal at the end of the 19th century, there was only intermittent connection between the Danube and the Lake Balaton drainage system. Therefore this subdrainage had been hydrologically isolated to some degree from the others until the last century. Waters from South Danubian region flow into the River Dráva (Drau), which empties to the Danube River at Osijek (Croatia).

Fish Sampling
Fish samples were collected between 2009 and 2012 by electrofishing from 27 sites across five sub-catchments of the Danube River and five sub-tributaries of the Tisza River (Table 1, Fig. 1). Since individuals of the Gobio genus show notable phenotypic plasticity, we investigated only adult (.60 mm standard length) specimens characterised by Gobio gobiolike morphological features, such as dispersedly spotted dorsal and caudal fins, and with no epithelial crests on the scales situated on the predorsal region of the body [35].

Molecular Methods
DNA extraction and purification. Fin clips of 241 specimens were sampled and stored in 96% ethanol at 220uC until DNA extraction. DNA was isolated with a DNeasy Blood and Tissue kit (Qiagen, Germany), using 10-20 mg of fin tissue as per the manufacturer's instructions. Quality and quantity of the extracted DNA were verified using a NanoDrop 2000c Spectrophotometer (Thermo Scientific, USA).
Mitochondrial sequence data. DNA of 168 out of 241 individuals (111 from the Danube River and 57 from the Tisza River catchments) were used for the amplification of the mitochondrial control region (mtCR). The sequences of mtCR were amplified by polymerase chain reaction (PCR) using the primers CR159 (CCCAAAGCAAGTACTAACGTC) and CR851 (TGCGATGGCTAACTCATAC) ( [33]). PCRs were carried out using 0.2 ml of 5 U/ml Taq DNA polymerase (Fermentas), 2.5 ml of 10X Taq buffer, 1.7 ml MgCl 2 (25 mM), 0.2 ml dNTPs (10 mM), 0.3 ml of each primer (20 mM), 2.0 ml template DNA, and 17.8 ml purified and distilled water in a final volume of 25 ml. Reactions were performed in a MJ Research PTC-200 Peltier Thermal Cycler under the following conditions: 95uC for 1 min, followed by 37 cycles of 94uC for 45 s, annealing at 52uC for 30 s, and an extension temperature of 72uC for 45 s, followed by a final extension at 72uC for 8 min. PCR products were purified from 1% agarose gel using the Millipore Ultrafree-DA DNA extraction kit. PCR products were sequenced on an ABI 3730XL sequencing machine by MWG-Biotech AG (http://www. mwg-biotech.com). Sequences were edited manually and aligned using the program Geneious 5.4 [43] and ClustalX 2.0.11 [44]. Newly described haplotypes have been deposited in the Gene-Bank. Calculation of sequence polymorphism and haplotype detachment was performed using DnaSP 5.10 software [45]. Sequence divergence was calculated with net nucleotide divergence (Da) in MEGA5 [46].
AFLP. To verify the results of the mitochondrial CR sequencing, a complementary method, Amplified Fragment Length Polymorphism (AFLP), was carried out; which is a reproducible, PCR-based molecular genetic method [41]. Altogether, 241 specimens were surveyed according to the following protocol. 200 ng DNA extracted from caudal fin tissue was digested at 37uC for 2 hours in a final volume of 10 mL with 2.5 U MseI, 5 U EcoRI enzymes and 2 mL NEBuffer4 (New England BioLabs, USA). Enzymes were then inactivated at 65uC for 20 min. Adaptor ligation was carried out at 24uC for 16 hours in 20 mL final volume containing the total digestion mixture, 0.25 mM EcoRI, 2.5 mM MseI adaptors, 200 cohesive end unit T4 Ligase and 1 6 T4 Ligase Reaction Buffer (New England BioLabs). After heat inactivation at 65uC for 10 min, 10 mL of digested, ligated mixture was diluted 10 fold with nuclease free water. Pre-selective PCR was carried out with AmpliTaqGold 360 Master Mix (Applied Biosystems, USA) in a final volume of 20 mL containing 0.5 mM Eco-A (59 GACTGCGTACCAATTCA 39), 0.5 mM Mse-C primer (59 GATGAGTCCTGAGTAAC 39) and 5 mL diluted, digested, ligated DNA. The PCR was started at 94uC for 2 min followed by 20 cycles of 94uC for 30 sec, 56uC for 1 min, 72uC for 1 min and a final elongation at 72uC for 7 min. Selective PCR was performed in a final volume of 20 mL containing AmpliTaqGold 360 Master Mix, 0.1 mM fluorescently labelled Eco-ACT primer (59 6FAM GACTGCGTACCAATT-CACT 39), 0.25 mM Mse-CTT primer (59 GATGAGTCCT-GAGTAACTT 39), 2 mL PCR product from the pre-selective PCR. Cycling conditions of the touchdown PCR were as follows: enzyme activation at 94uC for 2 min followed by 13 cycles for 30 sec at 94uC, for 30 sec at 65uC and decreased by 0.7uC in each cycle, and for 1 min at 72uC, then 23 cycles for 30 sec at 94uC, for 30 sec at 56uC and for 1 min at 72uC, followed by 5 min at 72uC. Digestion, ligation and PCRs were carried out in a Gene Amp PCR System 9700 (Applied Biosystems). Fragment analysis was performed on an ABI 3130 sequencer (Applied Biosystems, USA) and data were analysed with Peakscanner v1.0 (Applied Biosystems). Electropherograms were automatically analysed with tinyFLP [47], by the following scoring parameters: min. height: 90, max. width: 1, min. size: 50, max. size: 500, range (+/2): 0.5,  Table 1. doi:10.1371/journal.pone.0097278.g001 Table 1. Name, location, code, geographic position, and altitude (alt.) of sampling sites, the number of individuals (N) analysed by the two methods (mtCR, AFLP).  Cryptic Speciation of Stream-Dwelling Gobio PLOS ONE | www.plosone.org min peak-peak dist.: 1, peak height difference: 0, min. freq.: 0.1, max. freq.: 90. From the 553 peaks detected in total, 154 selected bands were retained. After further evaluation (e.g. specimens with a small number of peaks were excluded from the analysis) a dataset of 229 specimens was used for further statistical analyses.

Data Analysis
Mitochondrial sequence data. To shed light on the taxonomic relationships, alignment of all haplotypes found in this study and the previously published Gobio haplotypes (source: [31]) described from the neighbouring regions (e.g. Central Europe, Balkan Peninsula, and Anatolia) was performed ( Table 2). Originally the sequences revealed in this study were 651 bp long, but for the Network analyses we had to align them to the GenBank sequences of these closely relative Gobio species. Therefore for the Network analyses a 652 bp long dataset were used. Network was constructed using the median-joining algorithm in Network v. 4.6. [48]. Similar haplotypes were classified arbitrarily into haplogroups (see ''boxes'' in Fig. 2). Differentiation within and among haplogroups was tested by analysis of molecular variance (AMOVA; [49]) with 9999 permutations.
To construct the phylogenetic tree, the (652 bp long) sequence set analyzed in the Network analysis was aligned against further haplotypes used as outgroup taxa of varying putative phylogenetic depths (sources: [50], [51], [31], and Mendel et al.  [52]. The best-fitting models of DNA substitution were selected for analysis using Akaike information criterion (AIC) implemented in the jModelTest 0.1.1 [53], [54]. jModelTest indicated that Hasegawa, Kishino and Yano substitution model [55] with gamma-distributed rate heterogeneity (a = 0.5710) (HKY+G) was the best fitting. We conducted Bayesian tree construction with 6 chains, 2 independent runs and 7 million generations. Trees were sampled every 1000th generation. The first 10000 generations were discarded as burn-in. We plotted the log-likelihood scores of sample points against generation time using Tracer 1.5 [56] to ensure that stationariness was achieved after the first 10000 generations by checking whether the loglikelihood values of the sample points reached a stable equilibrium plateau. We used the remaining trees with average branch lengths to create a 50% majority-rule consensus tree with the sumt option of MrBayes. Posterior probabilities were obtained for each clade.
AFLP analysis. Higher level differentiation of Gobio assemblages was assessed using STRUCTURE 2.3.3 [57] to estimate the most probable number of genetic groups (clusters, K) for all analysed individuals. Values of K were investigated from 1 to 10, with a burn-in period of 10000 followed by 100,000 iterations and 10 runs for each K using an admixture model with correlated allele frequencies. Results of the Bayesian statistics were evaluated by Structure Harvester [58] implementing the Evanno method [59]. To characterise the standard measures of population genetic diversity, the percentage of polymorphic loci (%), mean unbiased heterozygosity, and unbiased Nei's gene diversity [60] were calculated. Within population genetic distance (GD) was calculated using the following equation:  GD~n 1{ 2n xy 2n where 2n xy equals the number of shared character states and n is the total number of binary characters. Population genetic structure was characterised by hierarchical AMOVA [49] with 9,999 randomisations. Isolation by distance was estimated by a Mantel Figure 3. Bayesian consensus tree derived from the analysis of the mtCR sequence data. Haplotypes revealed in this study are marked with their CR codes (see Table 3). Bayesian posterior probabilities are assigned on nodes. For GenBank accession numbers see the text and  Table 2. GenBank haplotypes used for the network computation. Numbers (N o -) of haplotypes correspond with the numbers displayed in Figure 2.
GenBank accession numbers code taxon name by GenBank

Mitochondrial DNA Sequence Analysis
Aligned sequences of 651 39-end CR mtDNA were obtained from 168 individuals grouped into 17 haplotypes. Sequence data of 15 previously undescribed haplotypes are highlighted in bold in Table 3. These have been deposited in the GenBank database  under Accession Nos. KC757328-42. The sequences of the CR03 and CR14 haplotypes had already been identified from the Danube catchment in the Czech Republic and in Slovakia [33], and were demonstrated by megaBLAST [66] to be Gobio gobio (100% similarity with the specimen: EU131542) and G. carpathicus (100% similarity with the specimen: EU131559) respectively. From the North Danubian region only three and from the Tisza River catchment a total of five haplotypes were displayed. The Middle and South Danubian regions were the richest in haplotypes. Moreover nine out of the 10 haplotypes found were unique to these regions ( Table 3). Results of the median-joining Network analysis showed that haplotypes described from Hungary were classified into five haplogroups (A, B, C, D, F 'boxes' in Fig. 2). Through the AMOVA analysis, 85% of the total genetic variance was explained by among haplogroup differences, and in addition significant (p,0.001) differentiation was found in each  (Síkhegyi, 2002)  pairwise comparison of haplogroups, confirming the arbitrary classification pattern.
With the exception of haplogroup 'B', each one corresponds to a previously described ''species'' (Table 4). In the Network analysis, 46 investigated specimens were identified as G. obtusirostris (haplogroup 'A'), 55 as G. sp1 described by Mendel et al. [31] (haplogroup 'C'), two as G. gobio (haplogroup 'D') and two as G. carpathicus (haplogroup 'F') (Fig. 2, Table 4). Altogether 63 specimens originating from the Middle and South Danubian region form a distinct, currently unidentified group (haplogroup 'B'), which was in a transitional position between G. sp1 and G. obtusirostris (Fig. 2). The results of Bayesian phylogenetic analysis (Fig. 3) were similar to those obtained through the Network analysis. However, the posterior probabilities of nodes showed a high level of uncertainly (weakly supported branching) at both ''lower'' (i.e. species) and at ''higher'' (i.e. genus) levels as well.
AFLP. The mean of the estimated Ln probability values from STRUCTURE analysis of the final matrix strongly increased between K = 1 and K = 3 and then consolidated at higher K values (Fig. 4a). The comparative statistics [67], [68] supported three major clusters (Fig. 4b, c). A triangle plot of the results (Fig. 4d) showed that individuals classified into the Cluster 1 originated only from the Danubian region while the individuals in Cluster 2 originated completely from the Tiszanian catchment area. The Cluster 3, consisting of mainly Middle and South Danubian specimens, shows a transitional position (continuous transition) between the two aforementioned clusters. This pattern was similar to the hydrological distances between the sampling sites (see Fig. 1c).
For population genetic analyses, data were used of those 21 populations (196 individuals) where N$8. The average number of bands per specimen was 40.569.4 (ranging between 16.0 and 62.0). Base population genetic features as: P. loci %, UHe, F is , F st are given for each population where N$8 in Table 3. Kruskall-Wallis tests revealed that the North, Middle and South Danubian and Tiszanian groups of populations did not show any significant differentiation in terms of their population genetic features (Table 3).
Within population genetic distance (GD), pairwise Wpt and pairwise Unbiased Nei Genetic Distances data are displayed in Table 5. GD ranged between 25.4 and 38.8 (av. 6SD = 34.863.1) and neither differed significantly by subregion, nor correlated with the altitudinal position of the collection site. Mean Unbiased Nei Genetic Distances ranged between 0.007 and 0.100 (av. 6SD = 0.03560.020). AMOVA analysis showed that among group differences accounted for 12% of the total genetic variance, and 193 out of the 210 pairwise comparisons (93%) showed significant (p,0.05) differentiation. Pairwise Wpt data ranged between zero and 0.321 (av. 6SD = 0.11660.06). This markedly strong population separation was verified by the results of DAPC. Assessing the AFLP dataset, 86% and 99% of the individuals were grouped correctly on population and region levels respectively in multidimensional space based on the cross-validation procedure within DAPC (Fig. 5). Results of Spearman rank correlations supported that the population genetic variables (P. loci %, UHe, F is , and F st ) were not significantly correlated (p,0.05) with the altitudinal position of the sites.
According to the results of Mantel tests, hydrological distances correlated significantly and positively both with Unbiased Nei genetic distances (Rxy = 0.499, p,0.01) and with the pairwise Wpt values (Rxy = 0.551, p,0.01) of the selected (N$8) populations.

Discussion
Our results only partially support the earlier assumptions of the taxonomic composition and distribution of Gobio species in the inner Carpathian Basin [31], [35], and show that Gobio reveal high diversity, both taxonomically and from the population genetic standpoint, within this relatively small and well-connected catchment area. At the same time, the presence of more cryptic species is indicated and the fine-scale separation of the identified genetic lineages between sub-catchments supports the existence of ecological barriers and on-going speciation in Hungarian drainage systems.

Taxonomic and Phylogenetic Features
All of the Gobio haplotypes found in the inner Carpathian basin can be classified into the north European clade and all except one haplotype (G. gobio) belong to the north-eastern European subclade described by Mendel et al. [31]. At the same time, our results indicate that the taxonomic status of the stream-dwelling gudgeons inhabiting the inner area of the Carpathian Basin is more complex than was previously presumed. Although no remarkable differences were detected in the morphological and meristic traits of the specimens analyzed, haplotypes of three previously described species (G. obtusirostris, G. gobio, and G. carpathicus), a doubtful taxon (G. sp1) and an additional, distinct haplogroup were distinguished from the study area.
Although the haplotype of G. carpathicus occurred in the Tisza River catchment, this area was dominated by the haplotypes of G. sp1 (Fig. 2, Table 4). Gudgeon stocks inhabiting the Danubian part of the country showed greater taxonomic complexity. Contrary to the earlier hypotheses [35], we found the haplotype of G. gobio in the Carpathian Basin (in Cuhai-Bakony-ér). Furthermore, haplotypes of G. obtusirostris proved to be dominant only in the North Danubian region. A distinct, but highly diverse haplogroup (haplogroup 'B') was dominant in the Middle and South Danubian regions (Fig. 2, Table 4).
Differentiation of haplogroups 'A' and 'B' (Table 4) may be attributed to a population split caused by paleohydrographic changes that took place in the geologic recent past. Namely, approximately 140,000 years ago a new watershed with an eastwest direction formed, separating the North and Middle Danubian region [69]. This changed the flow direction of the Zala River, originally flowing northward following the current channel of the Marcal River, southward to the Dráva River. Separation of the Middle and South Danubian regions began only at the end of the Pleistocene, 14-16,000 years ago by the formation of Lake Balaton [70], [71]. Namely, the Zala River, and all the smaller streams flowing southward until then, turned toward this newly formed depression. The phylogenetic effect of the relatively old watershed separating the North and Middle Danubian region and the incomplete splits -''intra-valley drainage divides''- (Fig. 6) between the Middle and the Southern regions [72] was proven by the fact that (1) the Network analysis showed significant differentiation between haplogroup 'A' and 'B', but (2) haplotypes did not differ between the drainages of Lake Balaton, River Kapos-Sárvízsystem, and River Dráva (Fig. 2, Table 4).
Although some specimens characterised by CR01 (G. obtusirostris) haplotypes also occurred in the Middle Danubian region, their presence was restricted to only those sites that were in the vicinity of fish ponds and/or where trout (Oncorhynchus and/or Salmo) stocking occurred [73]. Therefore, we assume unintentional G. obtusirostris introductions with gamefish (trout) stocking and thus a secondary, anthropogenic contact between the phylogenetically separated haplogroups. Table 3. Details of sample locations, population code and the number of individuals (N) analysed by the two methods (mtCR, AFLP).  Haplogroup 'B' showed similar genetic distances from the G. obtusirostris (1.31%60.31%) and from G. sp1 (1.61%60.17%), as an interspecific genetic divergence among some already accepted gudgeon species such as G. skadarensis and G. carpathicus (1.29%60.18%) or G. skadarensis and G. ohridanus (1.38%60.17%) ( Table 4). These differences make species level detachment of the haplogroup 'B' or G. sp1 taxon as well. At the same time, the results of Bayesian phylogenetic tree analysis showed weakly supported differentiation among haplogroup 'B' and some already accepted Gobio species (e.g. G. obtusirostris), and did not clearly support the recently accepted within-genus taxonomy. Our results indicate that the reproductive isolation of these entities may have only began in the geohistorical recent past, presumably in the middle Pleistocene. Contrary to the findings of the Network analysis, Bayesian phylogenetic computations in most cases query the species level differentiation within the Gobio genus. This premise is supported by the fact that in the case of other fish species, e.g. topmouth gudgeon -Pseudorasbora parva (Temminck & Schlegel, 1842), stone loach -Barbatula barbatula (Linnaeus, 1758), and grayling -Thymallus thymallus (Linnaeus, 1758), a similar or higher degree of differentiation among haplogroups is considered to be not more than subspecies level detachment [20], [74], [75]. In addition, some authors [76], [77] have suggested that the genetic distance of haplogroups must be greater than or equal to ten times the level of within-haplogroup differences to distinguish separate species. In our study, the G. insuyanus is the only taxon which fulfils the above mentioned criteria (Table 4).

Population Genetic Variables
The values of basic population genetic parameters (P. loci%, UHe, F is F st , and GD) did not show significant differences among subdrainage basins. Moreover, none of these variables showed significant correlation with the altitudinal position (as a possible marker for population isolation due to differences in the habitat use of fish) of the sampling sites. In the case of basic population genetic parameters, the local environmental conditions and the degree of hydrographic isolation are likely to be more important than either the altitudinal position or the taxonomic arrangement of the inhabiting specimens.
The population genetic features of gudgeon assemblages inhabiting the northwest region of Hungary differed notably from the other studied Hungarian assemblages. It is the only area where statistical analyses suggested considerable gene flow (Fig. 5, Table 5). This may be attributable to the species level differences. At the same time, landmarks of the river systems characterising this area assure convenient migration routes among sites. Similarly, the occurrence of Gobio gobio (haplotype) may be attributed to the role of Danube River. At the other Danubian Regions the population structure was much more differentiated. For the Middle Danubian Region, notable differences were found among closely situated sites. Results of the mtCR analyses revealed the existence of a different haplotype (CR01) in addition to the assumed ''native'' haplotype group from this area ( Table 3). Therefore these differences may be caused by accidental G. obtusirostris introductions to this area.
Results of STRUCTURE analysis are in accordance with the results of mtCR sequence analyses. Both inferred the existence of three larger clusters/haplogroups within the Carpathian Basin. Furthermore, both analyses indicated the transitional position of Cluster 2/haplogroup 'B' between the Cluster 1/haplogroup 'A' and the Cluster 3/haplogroup 'C'. The two specimens identified as G. gobio and two specimens identified as G. carpathicus by mtCR analyses did not show notable separation from the others by AFLP analysis (Fig. 4D) which may suggest interspecific hybridisation in Table 4. Estimates of evolutionary divergence over sequence pairs within and between haplogroups (see Fig. 4) based on 652b long mtCR sequence data.
Mantel test results revealed a clear pattern of isolation by hydrographical distance. Taxonomic and population genetic differences of the studied Gobio stocks were simultaneously changed by the growing hydrographical distances. This natural pattern is just slightly diminished by anthropogenic impacts.
Consequently, the genetic analyses confirmed the results of former analyses, which were based on mainly morphologic/ morphometric variables [23], since they revealed that the Middle European Gobio ''species'' form an extremely diverse and variable group. At the same time explanation of the phylogenetic relationships and within-genus taxonomic features are still partly unresolved. Our results showed that because of the casual immigration and/or the accidental introductions, and the sympatric occurrences, the location of the collection site is not a convenient feature to discriminate these ''species'' occurring in Hungarian waters.
Our results indicate that these cryptic Gobio entities form a relatively young phylogenetic group and that the genetic differences among them are not strong enough in most cases for species level differentiation. Moreover, considering the possibilities of interspecific and intergeneric hybridisation, the recent taxonomic partitioning of the Gobio genus needs re-evaluation. However, the discovered genetic diversity is probably very vulnerable. Since the separation of haplogroups seems to be only an intermediate phase of an on-going speciation and streamdwelling gudgeons have specific environmental needs and a restricted habitat area at present, habitat alteration and accidental stocking may easily damage the integrity of haplogroups. Consequently, conservation actions should be implemented to preserve the exceptional diversity of this fish group.