Toward the conservation of the endemic monotypic fish genus Aulopyge from the Balkan Dinaric karst: Integrative assessment of introduced and natural population

Abstract The complex biogeographical history of the Balkan Peninsula caused remarkable freshwater fish diversity and endemism, among which Cyprinidae fish dominate. The Dinaric karst was a Pleistocene refugium and it harbors ancient and endemic cyprinids, including Aulopyge huegelii, a sole representative of its genus. Being highly distributionally restricted, it faces various threats that promote a critical decline in population abundance and even population extinction. Phenotypic and molecular diversity of the introduced (Šator Lake, Šator Mountain) and natural (Studena River, Duvanjsko Polje) populations of Dalmatian barbelgudgeon from Bosnia and Herzegovina was studied by using two mitochondrial genes and morphometric traits (linear and geometric morphometrics). Nonparametric ANOVA showed that two analyzed populations significantly differed in six linear measurements, except snout length and postorbital head length. Contrary to centroid size, two populations were found to be significantly different in body shape. Deformation grids indicated that individuals from Studena River are characterized by wider and slightly shorter body comparing to individuals from Šator Lake. Incongruence in cytochrome c oxidase subunit I (COI) and cytochrome b (cyt b) mitochondrial DNA (mtDNA) variation was observed since a common COI haplotype was observed, while four and three cyt b haplotypes were registered in Šator Lake and Studena River, respectively. Since it was demonstrated that cyt b mtDNA was a faster evolving gene, we encourage its use in intraspecies studies, especially for evaluating the connectivity of fragmented populations and for studying the evolutionary footprint of the processes incorporated into the distinctive evolution of Aulopyge. Finally, findings herewith provide a firm basis for designing a long‐term sustainable conservation strategy for endemic species in Dinaric karst.

in the Qinghai-Tibetan Plateau (QTP) about 19 MYR (Benovics et al., 2017). The orogenesis of the QTP was tightly linked with the radiation of Cyprininae (29-18 MYR) and their first migration route, which was followed by the divergence of European Barbus sensu stricto (Palaearctic tetraploids) and the tetraploid Aulopyge lineage (around 16.6-15.5 MYR) (Machordom & Doadrio, 2001;Tsigenopoulos & Berrebi, 2000;Wang et al., 2013). As a result of ancient origin, Aulopyge is a basal lineage of Barbus sensu stricto and the sister group of Barbus and Luciobarbus subgenera (Machordom & Doadrio, 2001). The next migration wave of cyprinids from Asia to the peri-Mediterranean region occurred during Plio-Pleistocene, and it was associated with complex geomorphological and hydrogeological events that lead to further diversification and biogeographical structuring of barbels (Wang et al., 2013).
Intensive biogeographical history of the Balkan Peninsula caused remarkable freshwater fish diversity and endemism, among which Cyprinidae dominate (Oikonomou et al., 2014). The Dinaric karst, as the area of karstic springs, caves, and subterranean hydrological network, was a Pleistocene refugium and it harbors ancient and endemic cyprinids, including A. huegelii (Oikonomou et al., 2014). Aulopyge huegelii possess a unique morphological assemblage that includes autapomorphic and synapomorphic traits with other cyprinids (Bless & Riehl, 2002;Howes, 1987). Regarding its phylogenetic relictness, evolutionary uniqueness, and endemism, Dalmatian barbelgudgeon is considered as a species of high priority in conservation of the notable, but sensitive and thus vulnerable, Dinaric karst region. However, population structure of the species is still poorly understood, which implies further research related to spatial and temporal distribution of molecular and phenotypic variation. This is of fundamental importance since Dalmatian barbelgudgeon is seriously threatened due to habitat degradation and loss, water pollution, unsustainable water extraction, and introduction of alien invasive species (Crivelli, 2006). Hence, A. huegelii faces various threats, leading its populations to critically small abundance or even extinction (http://www.iucnr edlist.org/detai ls/61350/ 0). Thus, the species is assessed as Endangered B1ab(iii,v) according to IUCN (http://www.iucnr edlist.org/detai ls/61350/ 0). As a consequence, urgent habitat protection and population monitoring of the species were proposed (Ćaleta et al., 2009). Moreover, A. huegelii fulfilled both criteria for assessing biodiversity conservation priorities -endemism and phylogenetic diversity as the measures of the amounts of evolutionary history (Isaac et al., 2007). More importantly, deeper insights into genetics of the species might uncover ancient history of and Croatia (Krka. Cetina, Čikola, Ruda and Rumin Rivers) (Ćaleta et al., 2015Mrakovčić & Mišetić, 1990;Vuković & Ivanović, 1971). Livanjsko Polje is the largest karst field in the world (Ritter-Studnićka & Grgić, 1971), characterized by extraordinary karstic features circled by mountain ranges, including Šator Mt. One of the geomorphological and hydrological phenomena of this karstic region is a glacial Šator Lake situated on the Šator Mt. Šator Lake has been noted as one of two localities (another one being Blidinje Lake) where Dalmatian barbelgudgeon was introduced (Delić et al., 2005).
Unlike natural habitats of the species that belong to the Adriatic basin, Šator Lake belongs to the Black Sea basin through the Una and Danube rivers. Delić et al. (2005) regarded Šator Lake as a locality at the highest altitude of Aulopyge populations. Consequently, the translocated population of the unknown origin faces harsh environment which favors specific adaptations (Delić et al., 2005).
In addition, unlike other populations, Aulopyge from Šator Lake is supposed to spend its whole life cycle in aboveground water (Delić F I G U R E 1 Linear morphometric measurements: BD, body depth; CPD, caudal-peduncle depth; HL, head length; OL, orbit length; PHL, postorbital head length; SL, standard length; SNL, snout length; TL, total length. The white dots represent 13 homologous landmarks used for geometric morphometry (landmarks 1-13). The dot marked with number 14 is used only for the "unbend" procedure et al., 2005). Another large Dinaric karst field, Duvanjsko Polje located in Southeastern part of Bosnia Herzegovina is a 20 km long and 12 km wide. Duvanjsko Polje represented by paleogenic limestones and dolomites from Jurassic and Cretaceous, while newer sediments origin from Miocene (Radoš et al., 2012). One of the geological unique features is the subterranean Studena River that flows through Duvanjsko Polje, while in southwestern corner of the field the river goes underground within the main estaville and further continuing through karstic subterranean network. Then, the Studena River partially re-emerges at the source of Ricina and finally empties into the reservoir of Buško Blato (Radoš et al., 2012). Hence, Studena River is a natural home of the Dalmatian barbelgudgeon that faces both underground and belowground environments. Therefore, the primary objective of our study was to evaluate intra-and interpopulation phenotypic and molecular diversity of the Dalmatian barbelgudgeon in order to provide invaluable data for conservation management of this unique evolutionary cyprinids clade. In this study, we complementary used phenotypic traits (linear and geometric morphometrics) and two mitochondrial genes (cytochrome c oxidase I -COI and cytochrome b -cyt b) of the two populations, one of them was recently introduced in Šator Lake (Šator Mt), while the second was autochtonous population from Studena River (Duvanjsko Polje). The data on variation of the two mitochondrial genes was contrasting to evaluate their usefulness in conservation and forensics of this endangered fish species. Indeed, choice of appropriate molecular marker is a prerequisite for population discrimination and identification of divergent units. Given the power of DNA taxonomy in the biodiversity assessment, our additional goal was to fill the gap of DNA barcode database of the vulnerable ichthyofauna in the Mediterranean Hotspot region. Furthermore, regarding that phenotypic variation of A. huegelli has been studied based on linear morphometrics (Dekić et al., 2016;Mihinjač, 2018;Mušović et al., 2018) we choose to contrast traditional approach (linear measurements) with geometric morphometrics (body size and shape). Hence, study of the two populations, one faces to different selection regimes of under-and aboveground ecosystems (a sinking Studena River) and another population adapted to aboveground environment (Šator Lake), further provides comprehensive information on the phenotypic variation. It is important from the ecological perspective since pattern of subtle intraspecific phenotypic (e.g., body shape) is tightly linked to environmental variables (Collin & Fumagalli, 2011). Finally, findings herewith further contribute to understanding cryptic diversity on the Balkan Peninsula and provide a firm basis for designing a long-term sustainable conservation strategy for endemic species in Dinaric karst.   (Kottelat & Freyhof, 2007;Vuković & Ivanović, 1971). During the sampling of fish in Šator Lake, water parameters (water temperature, pH, dissolved oxygen, oxygen saturation, conductivity, and turbidity) at five points were also collected (  individuals were used-cyt b mtDNA: 24 from Šator Lake and 12 from Studena River; COI mtDNA: 23 from Šator Lake and 10 from Studena River. In addition, to study cyt b mtDNA and COI mtDNA variability, three cyt b mtDNA sequences (two from Buško Lake in Bosnia and Herzegovina and one from Krka River in Croatia) and four COI mtDNA sequences (one from Livno drainage in Bosnia and Herzegovina and three from Cetina River in Croatia) were downloaded from GenBank ( Figure 2; Table S2). Therefore, the total sample size included 39 cyt b mtDNA and 37 COI mtDNA sequences. Due to the lack of sexual dimorphism revealed by preliminary analysis, female and male specimens were not considered separately.

| Morphometric analyses
A total of 38 specimens from Šator Lake (25 individuals) and Studena River (13 individuals) were analyzed. Digital images of the right lateral side were photographed with a Nikon D7100 digital camera with 50-mm f/1.4D objective appointed on a tripod stand and positioned vertical to the surface of the object. Additionally, since the linear measurements and precise position for landmarks could not have been precisely determined on right-side images, for six specimens left lateral side was captured and use in morphometric analyses. The scale factor is determined for each photo according to the millimeter scale before the digitization process.

| Linear measurements
Linear morphometric characteristics were measured (in mm) on digital images using measurement mode in tpsDig 2.30 (Rohlf, 2017).
According to the standard for cyprinoid fishes (Armbruster, 2012), the total of eight linear measurements used for analysis were as follows: total length (TL), standard length (SL), orbit length (OL), snout length (SNL), postorbital head length (PHL), head length (HL), body depth (BD), and caudal-peduncle depth (CPD) (Figure 1). In order to remove allometry effect, linear measurements were transformed into (a) the percentage ratios of measured distance and SL (for OL, SNL, HL, BD, and CPD) and HL (for OL, SNL, and PHL), and (b) sizeindependent adjusted measurements. The second one was calculated in a way the size component was eliminated from a multivariate data set of measured distances following Elliot et al. (1995) procedure using formula: where M is the original morphometric measurement, M adj the size-independent measurement, L s the overall mean of SL for all fish in both population, L o the SL of the fish, b the slope (allometric coefficient) of regression of logM on logL o calculated for both population. Correlation coefficients between adjusted size-independent variables and SL were calculated in order to confirm that the size effect was eliminated and significant correlation was not found.
Prior to statistical analysis, a normal distribution test was performed (for each measure separately) on percentage ratios and size-adjusted measurements data using both Shapiro-Wilk test (for univariate normality). Also, for matrix of size-adjusted measurements, Doornik & Hansen omnibus test for multivariate normality was done.
Nonparametric ANOVA (with 10,000 permutation runs) was used to determine differences in variances among two populations.
Furthermore, principal component analysis (PCA) based on size-adjusted measurements was performed. All statistical analysis on linear measures was executed in PAST 4.03 (Hammer et al., 2001) software.

| Geometric morphometrics
In total, 13 homologous landmarks ( Unbending procedure was implemented in tpsUtil 1.74 (Rohlf, 2017) and four landmarks were aligned (1, 6, 10, 14). Landmark (14) was removed from following analysis. Procrustes fit procedure was used for superimposing landmarks' coordinates (Dryden & Mardia, 1998;Klingenberg & McIntyre, 1998) and information about shape variable (Procrustes coordinates) was obtained. Before unbending procedure was applied, centroid size (size variable) was also calculated. To test for the presence of allometry (the relationship between size and shape), a multivariate regression of Procrustes coordinates against centroid size on pooled within-group (pooled by population) variation was conducted. Permutation test with 10,000 iterations was used for checking significance of the allometry.
To examine body size differences among defined groups, nonparametric ANOVA with 10,000 permutaions was used. In order to determine and visualize differences in body shape between two a priori defined groups (populations), the matrix of shape variables (Procrustes coordinates) was subjected to discriminant function analysis (DFA). The dependability of the discrimination among groups is assessed using leave-one-out (jackknifing) cross-validation procedure. Finally, to estimate the extent of morphological variation of body shape, Procrustes distances between pairs of individuals and morphological disparity were calculated for each population.

| Length-weight relationship (LWR)
To analyze the length-weight relationship (LWR) of fish body we measured the weight of 11 and 13 individuals from Šator Lake and Studena River, respectively, using digital weighing scale KERN 440-33 to the nearest 0.01 g. The relationship between standard length (in cm) and weight was estimated according to the formula (Froese, 2006): where W -weight of fish, L -standard length (SL) of fish, a -scaling coefficient, b -length exponent (shape parameter for the body form).
Model of LWR was transformed into linear type of data using logarith- To assess whether the obtained b values statistically differed from the isometric value (b = 3), t test was used (Froese, 2006). Statistical analysis was performed in Microsoft Excel 2010 and PAST 4.03 (Hammer et al., 2001) software.  (Folmer et al., 1994). PCR reactions were performed using an illustra PuReTaq Ready-To-Go PCR Beads kit (GE Healthcare Life Sciences). PCR conditions for cyt b mtDNA and COI mtDNA amplification were described in Palandačić et al. (2012) and Milankov et al. (2009), respectively. To check the success of reactions, amplification products were separated on a 2% agarose gel. PCR products were then purified using ExoSAP-IT™ PCR Product Cleanup Reagent (Thermo Fisher Scientific), and bidirectionally sequenced on ABI3730XL by Macrogen.

| Data analyses
Chromatograms obtained by mtDNA sequencing were edited in Chromas 2.6 (Tehnelysium Pty Ltd) for erroneously called bases, while sequence alignment was performed in BioEdit 7.2.5 (Hall, 1999).
Haplotype networks were constructed in Network 10.1.0.0. (Fluxus Technology Ltd.) using a median joining approach. Haplotype divergences (p-distances) were obtained using MEGA X 10.0.5 (Kumar et al., 2018) by dividing the number of nucleotide differences by the total number of nucleotides compared. In addition, we determined total and private haplotype numbers per geographic sample.
Estimates of θ (where θ = 2Nu, N is the effective population size and u is the average mutation rate per locus per generation), were based on the infinite-allele model implemented in Arlequin version 3.11 (Excoffier et al., 2005). Using both mtDNA markers, four θ estimates were calculated: θ K obtained from the distinct number of haplotypes (K), θ H obtained from the observed homozygosity (H), θ S obtained from the observed number of segregating sites (S = number of polymorphic sites) and θ π obtained from the mean number of pairwise differences (π). θ K was estimated from the infinite-allele mutation model equilibrium relationship between the expected number of haplotypes (K), the sample size (n) and θ using Ewens' sampling formula. θ S Watterson's estimator is based on the infinite-site mutation model relationship between the number of segregating sites (S), the sample size (n) and θ. Tajima's estimator of θ (θ π ) is also based on the infinite-site mutation model, but on the relationship between the mean number of pairwise differences (p) and θ.
To assess the correlation of genetic distance (p-distances obtained from cytb mtDNA sequences) and phenotypic distances (Euclidian distance calculated from Procrustes-fitted landmark coordinates) Mantel's test with 10,000 permutations was applied.  (Table 1).

| Linear morphometric analysis
Nonparametric ANOVA performed on all types of linear measurements showed that two analyzed populations significantly differed in all but SNL(%HL) and PHL(%HL) ( Table 1). PCA on size-adjusted measures revealed that majority of the variation was explained with the first principal component (PC1 97.4%). Scatterplot for the first two PCs axes showed that those two populations clearly separated along PC1 ( Figure 3). The measure with highest loading on PC1 was BD adj (0.78).

| Geometric morphometrics
Normality test showed that the size variable (centroid size) deviated from the normal distribution (Shapiro-Wilk test, p < .05).
Therefore, nonparametric ANOVA on centroid size was applied and showed no statistically significant differences between two populations (F = 0.026, p = .88). Multivariate regression of shape variables on centroid size showed no significant allometry (p = .34) and accounted for only 2.98% of the overall shape variation.
Discriminant function analysis conducted on Procrustes coordinates showed significant body shape difference between populations (Procrustes distance = 0.0449, p < .0001) that were separated with no overlap along discriminant axis (Figure 4). Percentage of correct classification was 100% (92% after cross-validation).
Deformation grids indicated that displacement of landmarks 1, 7, 12, 13 had main influence on shape changes; individuals from Studena River are characterized by wider and slightly shorter body comparing to individuals from Šator Lake (Figure 4). Disparity analysis showed that Procrustes distance between the most divergent specimens within population was higher in Šator Lake (0.077) than in Studena River (0.066). Procrustes variances of two populations (Šator Lake: 0.00098; Studena River: 0.00106) did not differ significantly (p = .62).

| Length-weight relationship
The values of r 2 for both population were higher than 0.93. The obtained b value for Šator Lake is 2.490 (95% confidence limits: 1.987-2.994) and for Studena River is 2.922 (95% confidence limits: 2.443-3.401). The t test revealed statistically significant deviation from 3 for b value for population from Šator Lake (p < .05).

| Molecular diversity
A total of seven cyt b mtDNA haplotypes (around 1,110-bp long) were found among the 39 analyzed individuals (sequences obtained in this study will be submitted to GenBank after the manuscript acceptance) (Figure 1; Table S2). There were seven variable positions defining cyt b mtDNA haplotypes (Table S3), diverging up to three bases from each other (Table S4). Four haplotypes (HI, HIII-HV) were found in Šator Lake, three in Studena River (HI, HVI, HVII), whereas two and one were registered in Buško Lake (HII, HIII) and Krka River (HV), respectively. Unique haplotypes were found in Šator Lake (HIV), Studena River (HVI, HVII) and Buško Lake (HII) (Table S2; Figures 2 and 5). Intraspecific sequence divergences (uncorrected p-distance as a percentage) ranged from 0.09% to 0.38% (Table S4).
All the 23 individuals from Šator Lake, ten from Studena River, and two from Cetina River shared the same haplotype (HA), while unique haplotypes B and C were registered in Cetina River and Livno drainage, respectively (Table S2; Figures 1 and 5).
The genetic heterogeneity and discordance in variation of the two mtDNA genes in Šator Lake and Studena River is in agreement with the θ values. Contrary to cyt b mtDNA that possess four (Šator Lake) and three (Studena River) distinct haplotypes, COI mtDNA showed a lack of variation in both populations. Hence, the highest θ values were observed for cyt b mtDNA in the Šator Lake sample (Table S6).
Using the Mantel's test no significant correlation between morphology (body shape) and genetics (p-distance) was found (R = .04, p = .32).

| D ISCUSS I ON
In this study, we tested the usefulness of molecular markers (COI and cyt b mtDNA), and linear and geometric morphometrics in F I G U R E 3 Plots of the first two principal component (PC) scores from PCA on size-independent adjusted measurements of the Aulopyge huegelii populations. The percentage of explained variance of each PC is in parentheses F I G U R E 4 DFA histogram (above) and deformation grids of body shape differences between populations of Aulopyge huegelii (below). Numbers in the deformation grids refer to landmarks shown in Figure 1 quantification of variation in populations with different evolutionary histories (natural -Studena River vs. introduced -Šator Lake) and selection regimes (under-and aboveground ecosystems of sinking Studena River and strictly aboveground environment of Šator Lake) of A. huegelii, the relict cyprinid evolutionary lineage.
Contrary to cyt b mtDNA which possess four and three haplotypes in Šator Lake and Studena River, respectively, a lack of variation was observed at the barcoding fragment of COI mtDNA. Indeed, we revealed cyt b mtDNA as a faster evolving gene in this species, which is in concordance with the variation observed for A. huegelii from Buško Blato, Bosnia and Herzegovina (Mušović, 2016), and other endemic cyprinid species of Telestes (Buj et al., 2017;Gilles et al., 2010;Ketmaier et al., 2004;Perea et al., 2010), Delminichthys (Palandačić et al., 2012;Perea et al., 2010), and Phoxinus genera (Palandačić et al., 2015Perea et al., 2010;Vucić et al., 2018) from the Dinaric karst. However, contrary to this study, COI mtDNA expressed variation in autochthonous populations of Telestes metohiensis and T. dabar, Delminichthys and Phoxinus genera (Francuski et al., 2019;Perea et al., 2010). Significant discordance in variation of two mtDNA genes once again advocates implementation of specific approaches and molecular markers for studying any evolutionary entity, which has already been proposed in the DNA barcoding debate (e.g., Krishnamurthy & Francis, 2012;Moritz & Cicero, 2004;Rubinoff, 2006;Sheth & Thaker, 2017).
Interpopulation differentiation of the two genes was observed in our study in spite of the limited DNA data that was available.
On the contrary, COI mtDNA as a slowly evolving gene expressed lower variation within A. huegelii. For instance, individuals from Šator Lake and Studena River shared the same haplotype (HA), which was also registered in Cetina River (Geiger et al., 2014), while HC found in Livno drainage (Geiger et al., 2014) slightly differed from HA (p-distance: 0.15%) and HB (p-distance: 0.31%).
Contrasting variation of the two molecular markers was shown by θ values as well. The genetic heterogeneity at cyt b mtDNA expressed by θ K , θ H and θ S was slightly higher in Šator Lake compared to Studena River. However, introduced population possessed more cyt b haplotypes (I, III, IV, V) than natural populations (I, VI, VII), which is in contrast to the number of private haplotypes (Šator lake: IV; Studena River: VI, VII). Since the Šator lake population shares the same haplotypes with Studena River (I), Buško Blato (III), and Krka River (V) we assumed that more than once Dalmatian barbelgudgeon has been introduced to the novel environment. Indeed, data on origin and number of introduced fish is still unknown (Delić et al., 2005), further step of spatial analyses of molecular and phenotypic variation that will include more individuals and samples would provide better understanding of this subject.
Furthermore, by testing linear and landmark-based geometric morphometric approaches in quantifying phenotypic variation between the two populations of A. huegelli we found contrasting results. Phenotypic divergence was found to be statistically significant considering both raw and transformed (allometry-free) linear measurements. So far, linear measurements were used in study of phenotypic variation of the Dalmatian barbelgudgeon from Buško Blato (Guzina, 2000;Mušović et al., 2018) and Ždralovac canal in Bosnia and Herzegovina, and Čikola River, Visovac and Torak Lakes in Croatia (Mihinjač, 2018), but spatial pattern across the entire species area remains uncover. Comparing measurements used in these studies (see Mihinjač, 2018;Mušović et al., 2018) and sample from F I G U R E 5 Spanning network of cyt b and COI mitochondrial DNA sequence haplotypes in Aulopyge huegelii. Each circle represents one haplotype, the size of the circle is proportional to the overall number of individuals with that haplotype. Cyt b mtDNA: HI-Šator Lake, Studena River; HII-Buško Lake; HIII-Šator Lake, Buško Lake; HIV-Šator Lake; HV: Šator Lake, Krka River; HVI and HVII: Studena River COI mtDNA: HA-Šator Lake, Studena River, Cetina river; HB-Cetina river; HC-Livno drainage. Thick marks on lines connecting haplotypes represent single-nucleotide substitutions Studena River (Duvanjsko Polje), we found similarity in mean values as well as that ranges (min-max values) were greatly overlapped for all measurements we reported. Contrary to already analyzed natural (autochthonous) populations (Mihinjač, 2018;Mušović et al., 2018) and Studena River (data herewith), some individuals of the introduced sample attain smaller body size influencing smaller mean and minimal values of the studied traits of the Šator Lake sample. The larger body dimensions of fish in the river habitat compared to those in the lake we found in this study contradict the assumption that in lake habitats, due to uniform environmental conditions throughout the year and available food, larger individuals develop (Mihinjač, 2018). Indeed, Mihinjač (2018) reported higher body length for A. huegelli from lake systems than from Čikola River. However, standardized linear measurements independent to SL obtained through allometric transformation revealed that fish from Čikola River are characterized by longer head length in relation to lake individuals (Mihinjač, 2018) which is consistent with our findings. In addition, but not for Šator Lake where upper limit of 95% confidence interval was close to 3 (1.987-2.994). Lower b values we found and deviation from isometric growth for population from Šator Lake are possibly due to small sample size and standard length range covered by sample (Froese, 2006). Also, variability of the parameter b in some species is in relation with different sampling seasons and spawning period (Treer et al., 2008). It should take into account that this lack of consistency with published data regarding the difference in body Contrary to almost all studied linear measurements, no significant difference in centroid size obtained by geometric morphometric approach was observed. However, phenotypic dissimilarity was found regarding body shape, suggesting that changes were mainly associated with displacement of landmarks which influenced the body depth and head length in individuals. As far as we know, this is first study which implemented geometric morphometrics in understanding pattern of phenotypic variation of the Dalmatian barbelgudgeon.
However, geometric morphometrics has been approved as a valuable tool for uncovering subtle body shape divergence among conspecific populations of European freshwater fish underlying the influence of heterogeneous ecological environments (e.g., Bajić et al., 2018;Collin & Fumagalli, 2011Francuski et al., 2019;Marić et al., 2015;Ramler et al., 2017;Zaccara et al., 2019;Závorka et al., 2020). Namely, it was found that differences in body shape associated with environmental factors result in a morphologically optimized phenotype for a given habitat such as a more streamline (slender) body shape and larger head of lake fish compare to stream ones  as we also found. The observed shape differences are likely closely linked to swimming and feeding performances, and predatory pressure in different habitats . In addition, based on a very few available ecological data (physico-chemical parameters) for the Šator Lake and Studena River (Table S1) difference between these two sites, primarily in water temperature and dissolved oxygen concentration, is noticed. It has already known that oxygen availability, temperature as well as elevation could influence morphological changes of cyprinids (Collin & Fumagalli, 2011). Therefore, it is (again) important to highlight that recognizing morphological differentiation of the adaptive traits, such as shape divergence obtained in this study, mirrors effect of divergent local natural selection on populations from heterogeneous habitats. Considering recent introduction of Dalmatian barbelgudgeon to Šator Lake, which flows into the Unac River, and thus, belongs to the Una River and the Black Sea basin, this population has specific evolutionary history, biogeography, and habitats. Indeed, during over 40 years since its introduction, A. huegelii has adapted to the extreme environment (Šator Lake is situated at the highest altitude within the distribution area of the species, 1,488 m a.s.l.), unlike the introduced trout species that did not manage to survive (Delić et al., 2005).
Finally, given that the uniqueness of genetic, morphological, ecological, and life-history traits of A. huegelli provide valuable information on its distinct evolutionary history, broader study of the species across the distribution area is essential for setting conservation priorities and understanding the evolutionary history and phylogenetic diversity of the family Cyprinidae. This is of high relevance since assessment of taxonomic diversity is a prerequisite for indentifying priority areas of the conservation interest such as the Balkan karst region, which according to Reed et al. (2004) represents one of the most famous karstic regions of the world. As such, it is highlighted as the hotspot of freshwater fish biodiversity.

ACK N OWLED G M ENTS
We thank the two anonymous reviewers for their constructive comments which improved the manuscript quality. We would like to