Population Genetic Patterns of Threatened European Mudminnow (Umbra krameri Walbaum, 1792) in a Fragmented Landscape: Implications for Conservation Management

The European mudminnow (Umbra krameri) is a Middle Danubian endemic fish species, which is characterised by isolated populations living mainly in artificial habitats in the centre of its range, in the Carpathian Basin. For their long term preservation, reliable information is needed about the structure of stocks and the level of isolation. The recent distribution pattern, and the population genetic structure within and among regions were investigated to designate the Evolutionary Significant, Conservation and Management Units (ESUs, CUs, MUs) and to explore the conservation biological value of the shrinking populations. In total, eight microsatellite loci were studied in 404 specimens originating from eight regions. The results revealed a pronounced population structure, where strictly limited gene flow was detected among regions, as well as various strengths of connections within regions. Following the results of hierarchical structure analyses, two ESUs were supposed in the Carpathian Basin, corresponding to the Danube and Tisza catchments. Our results recommend designating the borders of CUs in an 80–90km range and 16 clusters should be set up as MUs for the 33 investigated populations. How these genetic findings can be used to better allocate conservation resources for the long term maintenance of the metapopulation structure of this threathened endemic fish is discussed.


Introduction
Although large floodplain rivers and their associated habitats (e.g., dead arms, backwaters, wetlands) represent a small fraction of the total land area, they nonetheless harbor an outstanding the mudminnow is a strictly protected species in Hungary, only few studies have carried out a comprehensive analysis of the distribution of the mudminnow, and are in need of a more recent update [32]. In order to promote the conservation of the remaining stock, captive breeding and rearing techniques have been developed [40,41], along with some pilot habitat revitalization and population translocation (transferring residual stocks from threatened habitats to newly established nature-like habitats) experiments [42,43]. Nevertheless, we still do not have a comprehensive overview of the genetic structure of the remaining populations distributed across shrinking and strongly fragmented habitat patches in the Carpathian Basin. According to the findings of field surveys, many of the known populations now consist only of very few individuals, which indicates a high probability of inbreeding and the extreme risk of their extirpation [32,42,44]. We do not know, however, how the total genetic diversity of the mudminnow is distributed among populations. Regarding the rapid decrease of mudminnow populations and the limited amount of resources available for their conservation, it is especially Fig 1. (a) Recent river network of the Carpathian Basin, with 33 sampling sites across eight sampling regions. Sites indicated by different colours belong to different regions. Open circles indicate further known (but not analysed) stocks of mudminnow. Periodically and permanently flooded areas before the beginning of river regulation works (in the mid. 19th century) are indicated by light and dark grey patches respectively. For detailed information, refer to the text and Table 1. (b) Distribution area of European muddminow (red coloured area). The location of the Carpathian Basin in Europe is indicated by a dotted rectangle.
important to explore population genetic patterns across the distribution area (i.e., toward the identification of ESUs), and to identify those populations containing most of the genetic diversity in order to set priorities for an effective action plan. For instance, imperative species conservation actions strongly demand information on which populations may be used for supplementing (i.e., serving as a pool of parent fish in artificial breeding programs) endangered populations and stocking reconstructed habitats in different geographic areas and on which populations have the genetic uniquity requiring them to be preserved without mixing them with other populations (i.e., the determination of conservation and management units, MUs and CUs) [17,45].
Therefore, the aims of the current study were: (i) to explore the recent distribution of mudminnow populations in the Hungarian part of the Carpathian Basin, the area inhabited by the majority of the known stocks; (ii) to reveal spatial patterns in the population genetic structure and among sites' migratory rates; (iii) to assess the conservation biological value of extremely small populations; and (iv) to designate ESUs, CUs and MUs.

Ethics Statement
This study was carried out following relevant national and international guidelines pertaining to the care and welfare of fish. Our surveys (e.g., tissue sampling during field surveys) are not qualified as animal tests by the operative Hungarian law (N o : 40/2013.(II.14.). However, since the studied species is strictly protected in Hungary, any procedure to be applied in connection with it is controlled by the order of government N o : 348/2006.(XII.23.), and is subject to authorisation from the National Inspectorate for Environment, Nature and Water. The collection and storage of mudminnow tissue samples was authorized by this bureau (permission numbers: 14/881/5/2011, 14/678-9/2012). During the sampling procedure, all efforts were made to minimize suffering. Field studies did not involve other fish species that were endangered according to the IUCN Red List of Threatened Species v. 2014.3 (www.iucnredlist.org).

Sample collection, DNA extraction and microsatellite amplification
Studied areas and sites were selected based on the last published distribution data of mudminnow [32]. Specimens were collected by electrofishing between 2011 and 2013, from each region of the Carpathian Basin where the presence of this species had previously been noted (Table 1). To minimize suffering, individuals collected for this study (N = 5-20 per site, depending on the stock size) were narcotized using clove oil. The active agent of clove oil is eugenol (4-allyl-2-ethoxyphenol), which is considered non-carcinogenic, non-mutagenic, and a "Generally Recognized As Safe" (GRAS) substance by the U.S. Food and Drug Administration [46,47]. The aqueous emulsion of clove oil was deposited into the dangerous refuse containers of the Balaton Limnological Institute after use. For the genetic survey, an approximately 2mm 2 anal fin clip was sampled from each specimen and stored in 96% ethanol at -20°C until DNA extraction. After fin tissue sampling, fish were released to the habitat they were captured from. DNA was isolated with DNeasy Blood and Tissue kit (Qiagen, Germany), using 10-20 mg of tissue, as per the manufacturer's instructions. Quality and quantity of the extracted DNA were verified using a NanoDrop 2000c Spectrophotometer (Thermo Scientific, USA).
Nine microsatellite markers previously developed by Winkler and Weiss [48] were investigated by multiplex PCR. 50 ng of DNA was amplified in a multiplex PCR reaction with three primer pairs in a PCR reaction using Type-it Microsatellite PCR Kit (Qiagen, USA), and fluorescently labelled forward primers according to the protocol of the kit (details and multiplex combinations are provided in Table 2). The final concentration of each primer was 0.2 μM, and the annealing temperature was 60°C in each PCR. Amplified products were detected using an ABI 3130 sequencer (Applied Biosystems, USA). Peak Scanner Software v1.0 was used to analyse the data. The functionality of all primer pairs and the reliability of multiplex PCR were tested in singleplex reactions using the same PCR conditions and reagents.

Statistical analyses
Genetic variation. Microsatellite markers ( Table 2) and primers published by Winkler and Weiss [48] were checked in GenBank before PCR and microsatellite analysis. Based on the  [49] for each locus using a Markov chain of 10,000 dememorization steps, 20 batches, and with 5,000 iterations per batch. MICRO-CHECKER 2.2.3. [50] was used to estimate the frequency of null alleles for each locus. Number of alleles, and locus specific F st values, were calculated using GenAlEx 6.5 [51]. Similarly, this software was used to estimate the mean of allelic richness (Ar), Shannon's Information Index (I), observed and expected heterozygosity (H o , H e ), fixation index (F) and list the private alleles for each population.
Population structure. The Lynch & Ritland [52] estimator was used to calculate betweenspecimen pairwise relatedness. From this semimatrix, mean within-population pairwise values were calculated with 999 permutations, and 1,000 bootstraps in GenAlEx 6.5 [51]. Because in some cases null alleles were detected, their effect on pairwise F st values were corrected using the ENA procedure of Chapuis and Estoup [53,54]. Genetic distances among populations were calculated using the Cavalli-Sforza and Edwards [55] estimator after INA correction [53], and presented in PCoA ordination. GenAlEx 6.5 [51] was used to perform the analysis of molecular variance (AMOVA), with 999 permutations for population differentiation and hierarchical partitioning of genetic variation among and within regions and populations (F-statistics: F rt -H 0 : individuals are shuffled among regions, F sr -H 0 : individuals are shuffled within regions, F st -H 0 : populations are shuffled among regions, F is -H 0 : individuals are shuffled within populations, and F it -H 0 : individuals are shuffled in the whole sample).
Genetic population structure was inferred using the hierarchical approach [56] of the STRUCTURE analysis [57] to estimate the most probable number of genetic groups (clusters, K) for all analysed individuals. Namely, the STRUCTURE analyses was first run including all samples, then samples were separated by river drainage basin, and finally by region. (A, B, C, Table 2. Primer sequences, fluorescent dyes and combinations of primers in three (1., 2., 3.) multiplex PCR, number of alleles, locus-specific F st values, number of populations out of the 33 in which deviations from Hardy-Weinberg Equilibrium (HWE) were detected. (*UkrTet2 UkrTet3 sequences are identical, therefore locus UkrTet2 was omitted from the further analyses, for more details refer to text), and Mean ± SD estimated frequency of null alleles. etc.) Values of K were investigated from 1 to 20, with a burn-in period of 100,000 followed by 100,000 MCMC iterations and 10 runs for each K using an admixture model with correlated allele frequencies. Results of these Bayesian statistics were evaluated by STRUCTURE HAR-VESTER [58], implementing the (deltaK) Evanno method [59]. Results of the 10 repetitions were combined using the software CLUMPP 1.1.2. [60]. For the genetic assignment of the studied individuals, Bayesian cross validation tests [61] were carried out on drainage basin, regional, and population levels using GeneClass2 [62] software. Cross validation tests were carried out similarly on the clusters defined at various levels by the hierarchical STRUCTURE analyses.
Migration and spatial analyses. Migration rate estimation was performed by MIGRA-TE-N 3.2.15 [63,64], using the maximum likelihood method under the following parameterization: 20 short chains (500 trees used out of the sampled 10,000) and 5 long chains (5,000 trees used out of the sampled 100,000). Missing data were not included. Theta values were generated from the F st -calculation. Mutation rates among loci were estimated from the data. Migration was estimated for the original regional groups designated by catchment basin, and based on the results of hierarchical STRUCTURE analyses at the three levels.
Spatial genetic structure was assessed at the population level using Mantel tests [65], through the comparison of pairwise F st data and pairwise straight line geographic distances (GGD) with 999 randomisations for the whole dataset and at drainage basin levels. On the individual level, two kinds of spatial autocorrelation [66,67] computation were made. Spatial autocorrelation computation is suitable to accurately identify the scale at which population genetic structure is detectable. In this case, "r" coefficients were calculated using multiple distance class analyses. This method plots "r" as a function of increasing distance class sizes. The last distance class for which "r" is significant is considered to be the limit of detectable isolation by distance (IBD). In the second case, autocorrelation coefficients, "R", were plotted as a function of discrete distance classes, partitioned so as to achieve a similar number of pairwise comparisons for each class [68]. Positive and significant "R" values indicate IBD and the "x" intercept provides an estimate of the extent of IBD. For more details see Peakall et al. [69]. All of the spatial autocorrelation computations, including the GGD semimatrix calculation from geocoordinates, were made in GenAlEx v6.5 [51].

Results and Discussion Results
Distribution patterns. As a result of our broad field investigations, European mudminnow was observed at more than 40 sampling sites from the Middle Danubian catchment (Fig  1). A total of 404 mudminnow specimens were sampled for population genetic research, from 33 sites, originating mainly from artificial habitats (i.e., ditches, canals, and ponds), across eight regions (Fig 1, Table 1). In several locations, no more than 5 individuals were able to be collected at one time, despite the intensive sampling effort.
Data quality and Genetic variation. Results of Fisher's exact tests did not reveal evidence for large allele drop-out for any locus, or for linkage disequilibrium between any pairs of loci. All microsatellite loci showed significant deviation from the Hardy-Weinberg equilibrium state if all 404 specimens were analysed as one entire population (404 samples), and we calculated high F is values (>0.08) at all loci among the 404 samples. These results suggest strong genetic substructure caused by geographical isolation and/or non-random mating in subgroups. Characteristics of each microsatellite locus used in our work are provided in Table 2. The total number of alleles per locus varied between 6 and 25 (mean: 17.6). Locus specific F st values ranged between 0.148 and 0.482. Chi-Square Tests for Hardy-Weinberg Equilibrium for each loci at each population showed significant values in 17 cases (see Table 2 and S1 Table).
A total of 136 alleles were observed for the eight loci used in the analysis. The whole raw dataset used for further analyses is indicated in S2 Table in GENEPOP format. Average allelic richness, such as observed heterozygosity, showed high-level differences among populations. These ranged from 2.75 to 10 and from 0.44 to 0.85 respectively, within the studied populations. Shannon's Information Index showed a high level of variability as well, ranging between 0.8 and 1.97 (Table 1). Deviation from HWE was significant (p<0.05) in two studied populations (E8, G2).
As MICROCHECKER analyses proved the presence of false homozygotes (null alleles) in some cases (i.e., for UkrTet1 at sites B3, D2, E3, F1, for UkrTet3 at site E3, for UkrTet4 at site B2), therefore null allele correction was made prior to further analyses (see S3 Table).
13 of the detected 136 alleles were unique to a single location ( Table 3). The frequency of the unique alleles ranged between 0.031 and 0.071. Eight of the unique alleles were observed in sites of the Tisza River drainage basin, and five from the Danubian part of the study area. No unique alleles were found for small stocks (where N = 5), in populations A3, C1, E5 and F2.
Mean within-population pairwise correlation relatedness values, "r", ranged between 0.063 and 0.651, and were significant in all cases (Fig 2). These values correlated negatively with the observed heterozygosity and with Shannon's Information Index (Spearman's rho: -0.7681, p<0.01, and -0.956, p<0.001 respectively). Therefore, higher "r" values refer to higher separation and lower genetic diversity (these features may indicate inbreeding). In some cases (e.g., populations G2 = 0.651, C2 = 0.498, D2 = 0.467), high values were found, but pronounced regional patterns were not. In almost all sampling regions (except regions A and B), one or more populations formed by strongly related individuals were observed.
The PCoA plot based on pairwise population genetic distances (Fig 3) showed clear separation along the "x" axis of eastern (Tisza drainage basin) and western (Danube drainage basin) sites. Moreover, the populations in the Tisza drainage basin showed higher similarity. In the  case of these populations, the separation along the second axis is also important. Pronounced separation of neighbouring populations (i.e., populations from the same region) within a region can only be seen in the case of population E1. This population belongs to the Danube drainage basin Middle Hungarian region, but shows stronger similarity to the Mura region or to the C1, C2, and D4 populations of the Tisza drainage basin. AMOVA analysis showed that among-population and among-regional differences represented 12.8% and 7.8% of the total variation respectively, while within-population differences accounted for 78.9% of the total genetic variance. Results of the regional F-statistics, except F is , showed significant within-and among-regional differentiation and population level differences.
Bayesian hierarchical clustering separated Danubian and Tiszaian mudminnow stocks in the first step (Fig 4). The only exception is the E1 population, which was assigned clearly to the Tiszanian clade. The Tiszanian group was split into 4 clusters; the E1 population, Tápió region (C1-2), and sites from NE Bihar region (D1-3) were separated from the others. Stocks form the Upper Tisza region (A1-4), Borsodi-mezőség region (B1-3), and SW of Bihar plain (D4-5) were separated only in the third clustering step, and the extent of separation is not as unambiguous as in the previous steps. The Danubian branch split to two larger groups; the Middle Hungarian region together with the Hanság-Szigetköz region (E+F) formed a group separated from the Balaton and Mura (G+H) regions. The E+F cluster split further into four minor groups (E2 +E8; E3+E4+E5; E6+E7; F1-F4). Populations in the G+H cluster formed five minor groups. Distinct populations were found in this cluster, except G3-G4-G5 stocks, which all originated from the Kis-Balaton marshland area and formed one population.
Bayesian cross-validation showed that 98.3%, 89.6%, and 71.8% of the individuals were grouped correctly on drainage basin, region and population levels respectively (Fig 5). Most of the reclassifications were made within the same region, mainly between neighbouring sites (e.g., C1-C2).
Result of the cross validation procedure showed negligible linkage between the Danubian and Tisza drainage basins. Only three specimens (i.e., 0.7% of the total sample) were misclassified between the two drainage basins. Within a drainage basin, the percentage of misclassified cases is higher; five specimens of the total 156 (i.e., 3.2%), and 13 specimens of the 248 (i.e., 5.2%) were classified to other regions on the Tisza and Danube drainage basins respectively. The highest number of misclassified cases was found between the Middle-Hungarian and Hanság-Szigetköz regions (E and F regions) where 4.9% and 10% of the stocks were misclassified.
Population assignment results from GeneClass2 strongly support the structure revealed on the three hierarchical levels identified by Bayesian STRUCTURE analyses. 97.8, 96.0 and 94.6% of the individuals were classified correctly on the three hierarchical levels respectively (see: Fig  4). The third level of the hierarchical STRUCTURE analyses was therefore found to have classified the individuals much more reliably than the original population based classification.
The between-region (A-H) migration ranged between 0.28 and 2.31 individuals per generation. Mutation-scaled immigration rate is less than one individual per generation in most cases, therefore negligible gene flow is assumed among the eight studied regions in the Carpathian Basin, especially in the case of regions C and H.
Migration computations between the two clusters (corresponding to the Danube and Tisza catchments, designated on the first level of hierarchical STRUCTURE analysis) showed that   On the second level of STRUCTURE analysis, the immigration rate (M) between the four clusters (Figs 4 and 6) found in the Tisza catchment ranged between 0.75 and 1.52. The highest level of mutation-scaled immigration rate was found between the "a" and "c" clusters. In the case of the two Danubian clusters (see e and f clusters on Fig 4) the migration rates were less than one in both directions, resulting in a 1.36 and 1.48 individual per generation migration rate. These groups are therefore practically separated. Pairwise between-region and betweenclusters immigration rates are presented in Tables 5-9. A Mantel test conducted on the entire dataset at the population level showed no significant correlation between pairwise F st and GGD semimatrices (Rxy = 0.061, p>0.05). The inverse was found for the two drainage basins, where the genetic differences showed significant correlations both for the Tisza (Rxy = 0.4196, p<0.01) and the Danubian drainage basins (Rxy = 0.2489, p<0.01) as well. Therefore, population level IBD can be detected only within the Danubian and Tisza drainage basins (see Fig 7) Spatial autocorrelation computation carried out on the entire dataset at the individual level showed that IBD was not appreciable beyond an approximate 340-360 km distance. There was a significant positive correlation for the first two distance classes, which intercepted the "x" Estimated population structure as inferred by three rounds of hierarchical structure analyses. Each individual is represented by a thin horizontal line, which is partitioned into K-coloured segments representing individual's estimated membership fractions in K clusters. Black lines separate individuals from different sampling sites. The most probable K for a sample, given by the arrows, is based on the results of StructureHarvester using the Evanno method. Codes of sampling sites correspond to those of Table 1. Red codes represented on the three hierarchical levels correspond to the codes used in Fig 6. doi:10.1371/journal.pone.0138640.g004 axes at 124 km (Fig 8A and 8B) and was therefore significant within this range. Autocorrelation calculation at the drainage basin level showed similar results in both cases. The IBD is still detectable to approximately 150-160 km, but it is significant only within an 80.1 and 86.7 km range for the Tisza and Danube drainage basins respectively (see Fig 9).

Discussion
Tockner et al. [4] noted that fragmentation, water stress, land use change, and the introduction of non-native species are the most important threats to freshwater biodiversity. All these factors currently threaten mudminnow stocks, and their impacts are reflected in present day distribution and population genetic patterns. In this study, population genetic structure and dynamics of genetic mixing among populations in a strongly altered and fragmented landscape were analysed, and revealed some basic principles of optimal species conservation strategies.

Distribution patterns
The field surveys of the current work have provided new information about the present day distribution of the European mudminnow at the centre of its range. In the "C" and "F" regions, stable and large assemblages were found, contrary to earlier communication which noted that the mudminnow was extremely rare and had been extirpated from most habitats in these regions [32]. At the same time, due to the expansion of Amur sleeper (Perccottus glenii) [70], European mudminnow has disappeared from many locations in the Hungarian Upper Tisza region (i.e., region "A" of this study), where it had previously been abundant [44]. At present, these shrinking populations seem to be the most threatened ones within the Carpathian Basin.
In almost all regions, the largest assemblages were found in artificial habitats. Of the 33 sampling sites, only seven were in a near-pristine state (Table 1). Consequently, similar to other threatened stagnophilous fish, it seems that some of the human-altered or man-made habitats (e.g., irrigation and drainage canals) with rich macrophyte coverage and permanent water supply might be the last chance to provide refuge for the shrinking European mudminnow populations [27,71]. Genetic and spatial structure This study presents the first data on the population genetic structure of threatened and endemic European mudminnow for the centre of its distribution range. In general, remarkable genetic differentiation was found among the 33 populations originating from eight regions. The spatial structure revealed seems to be quite fixed with negligible migration among the studied regions ( Table 5). The genetic structure corresponded clearly to past and recent geographical migratory routes (i.e., it reflected natural distribution patterns). The only identified anthropogenic alteration within this pattern is a result of a documented stock transfer from Middle Hungary to the Hanság-Szigetköz region (i.e., from region E to region F) in April 2005 [42]. Results of Bayesian cross validation tests (Fig 5), confirmed by migration tests (Fig 6), suggest that 10% and 5% of the specimens occurring here are closely related to Middle Hungarian (E) and Hanság-Szigetköz (F) populations respectively (Figs 5 and 6; Tables 5-9.). An alternative explanation for the mixed genetic composition of the Hanság-Szigetköz (F) region population could be the dispersion along the Danube valley route, as has been proven for many other species [72,73]. In any case, SRUCTURE analysis showed (Fig 4) that the aforementioned stocking activity had only a moderate effect on the genetic pattern of mudminnow assemblages inhabiting this region. Some of the studied populations showed complete isolation and a high level of relatedness (Figs 2 and 5). Populations which are characterised by special environmental requirements and restricted expansion ability may hold unique pieces of genetic diversity, and their loss can therefore cause a significant reduction in genetic diversity in general [74]. Therefore, to conserve the genetic diversity of this species, special attention must be paid to the preservation of these separated stocks (e.g., B2, D1, and G2 on Fig 2).
Although European mudminnow has a restricted dispersal ability [28], and historical processes therefore mainly underlie the pronounced spatial pattern in the genetic structure, our analyses also captured some human-induced effects related to the anthropogenic modification of the Carpathian Basin floodplains. Historical hydrological data [75,76] suggest that prior to Table 9. Within "f" cluster theta values (Θ) and immigration rates (M). Clusters defined in the third round of hierarchical STRUCTURE analysis. + receiving population (For cluster codes see Fig 4).  Population Genetic Patterns of Hungarian Umbra Stocks river regulation there were no notable barriers among some sub-drainage basins in the floodplain of the Tisza River (Fig 1). Periodic floods inundated almost the entire floodplain, sometimes for more than 3 months at a time (see: [22]), and migration could therefore take place freely among A-B-D regions in the Tisza drainage basin. As a result of river regulation, which began in 1846 and ended in the 1930s [77], the connections between regions became much more restricted. The strong distinctiveness of populations inhabiting the Tápió region (i.e., the "C" region here), outside the main floodplain, also supports this hypothesis (Fig 3). According to data from the literature (cited in: [28,39]), mudminnow may spawn at the age of one year for the first time. The separation of some Tiszanian regions could therefore have begun more than 100 generations ago. The separation of the Danubian and Tiszanian systems, as well as the mudminnow stocks, is presumably much older. The within-region separations are not complete in many cases, however, because the continuous flood-protection and regulation activity (i.e., dredging, flow direction shifts within the drain-canal systems) could have caused the connections between stocks to be re-established from time to time [78]. E1 is the only population among the 33 studied which seems to be misclassified by both the STRUC-TURE and the pairwise population genetic distance analyses (Figs 3 and 4). Although this site currently belongs to the Danubian drainage system, the E1 population seems to be more related to those of the Tisza catchment. This relationship can be explained by paleohydrological changes caused by the emergence of Gödöllő Hilly Region in the Quaternary period [79], when the river system of this area was partly disjointed and the flow direction of the upper sections of the streams were diverted to the NW (i.e., to the Danube). At the same time, the two catchment areas were separated only by marshy areas referred to as "intra-valley drainage basin divides" [80], which these assemblages can traverse. This assumption is supported by the fact that the closest populations (C1-C2) to the Tápió system showed the highest similarities and bidirectional gene flow was found between the two areas. However, it should be noted that both assemblages showed high levels of isolation and inbreeding, and the connection between them is therefore presumably strictly limited. Similarly, contact between the G and H regions is plausible, because the isolation of the Middle and South Danubian catchments started only at the end of Pleistocene (by the formulation of the Lake Balaton basin) and this separation is still not complete [81]. Our results correspond with the findings of Brauer et al. [11], who investigated a South Australian endemic fish species with very similar environmental requirements to those of European mudminnow, the pygmy perch (Nannoperca obscura Klunzinger, 1872), and found similarly high levels of isolation, and designed CUs in very similar ranges. These results therefore show that preferred habitat type and dispersion abilities may have major roles in the formulation of these features rather than the geographic range (i.e., Australia vs. Europe) or the taxonomic position (i.e., Perciformes vs. Esociformes) of the studied species.

Implications for conservation and management
Although pairwise F st computation showed significant differentiation in most cases, none of the small stocks differed considerably from their neighbouring populations (Table 5; Fig 5). Moreover, these small stocks did not carry any private alleles (Table 3). It is therefore likely that these small stocks belonged to a larger metapopulation system, and could have been isolated only recently. Therefore, special consideration for their preservation is necessary only if there is no larger population nearby presenting the same genetic features.
Based on the current data, it is recommended that any reconstructed or newly established habitats should be stocked by individuals originating from the same CU, namely from populations within an 80-90 km range. The 16 clusters separated in the third step of the hierarchical STRUCTURE analysis can be accepted as MUs for conservation, instead of populations. In many cases, the MUs agree with some strongly isolated populations (e.g., populations C1, E1, G1, and G2). These separate entities need special attention to prevent substantial genetic diversity loss. The ex situ preservation of these populations could also be considered, preferably by artificially recruiting and translocating them to fish-free revitalized or newly established nature-like habitats (see e.g., [43]).

Conclusions
The investigations reported here reveal unexpectedly high genetic diversity of this endemic fish species, despite its declining and fragmented populations. The studied stocks show high levels of genetic variability, and this pattern seems to have been only slightly influenced by anthropogenic impacts (i.e., resettlements) so far. Continued habitat loss, along with the invasion of non-native competitors, strongly threaten this habitat specialist, endemic fish. Climate change increases the probability and severity of dry periods in the area, and represents a high risk for European mudminnow with its shrunken and fragmented population structure [82]. Consequently, there is a strong need for implementing comprehensive conservation management programs. Nevertheless, the extremely small and/or isolated populations (e.g., in region "C") would require special attention, and thus should be prioritized for conservation.
Supporting Information S1