In the foothill zone—Sabanejewia balcanica (Karaman 1922), in the lowland zone—Sabanejewia bulgarica (Drensky, 1928): Myth or reality?

Abstract The status of golden loaches (genus Sabanejewia) in the region of Central Europe and Balkans is still ambiguous. The greatest controversy is caused by species Sabanejewia balcanica and S. bulgarica. Both species are characterized by a wide spectrum of morphological variability and overlapping of distinguishing features, which then lead to difficulties in their determination. Previous phylogenetic studies aimed on the resolving of their taxonomic status did not include samples from their type localities and so led to a lack of their true distribution in this region. Therefore, the main aim of this study was to identify taxonomic status of golden loaches populations in the region of the middle Danube basin and adjacent areas on the model territory of Slovakia. For this purpose, we used novelty approach (morphological, molecular, and microhabitat) and we also included the missing samples from the type localities of both species. Based on mtDNA all the Slovakian samples reflected haplotype richness revealed on the type locality of S. bulgarica, although the genetic distances from other representatives of the genus Sabanejewia occurring are not significant. Within the morphology, we have revealed a great measure of variability in studied populations, which is largely caused by different habitat conditions and thus representing a phenotypic plasticity of these fish.


| INTRODUC TI ON
Systematics of loaches of the genus Sabanejewia actually include 10 fish species (Kottelat, 2012), of which eight occur in Europe (Kottelat & Freyhof, 2007;Marešová et al., 2011) and two remaining are widespread in southwestern Asia (Sayyadzadeh, Abbasi, & Esmaeili, 2018). However, in the beginning representatives of this genus were assigned to the related genus Cobitis. Until Vladykov (1929) performed a detailed morphological analysis and said: "It's Sabanejewia!" But, from of its establishment the genus by itself was questioned. As generally accepted among scientists, the validity of common name Sabanejewia has met with recognition until paper published by Nalbant (1963), who acknowledged the Vladykov's claims of significant morphological difference of this genus as justified. For a long time, taxonomy of individual representatives of Sabanejewia genus was also complicated. Almost all populations of golden loaches in Europe were perceived as polytypic species Sabanejewia aurata (Filippi 1863) (Bănărescu, Nalbant, & Chelmu, 1972). Subsequently, several of its subspecies were described by other classic morphologists (Drensky, 1928;Economidis & Nalbant, 1996;Jászfalusi, 1951;Karaman, 1963;Nalbant, 1957;Vasiľeva & Vasiľev, 1988;Witkowski, 1994).
For these reasons, we decided to examine these irregularities on the model territory of Slovakia, where the status of Sabanejewia fish has also been ambiguous. The occurrence of both morphological forms together with their intergradation forms in some localities was recorded in this area. Novomeská and Kováč (2016) state that there is more than one species of Sabanejewia occurring in this country.
However, these claims are not supported by any further information.
Based on the variability of cytochrome b gene in samples taken from six Slovakian rivers, Bartoňová et al. (2008) has included them into the sublineages III and IV of the DB complex (Perdices et al., 2003) and concluded that only species S. balcanica (Karaman 1922) occur in this territory. However, individuals resembling species Sabanejewia bulgarica (Drensky, 1928) by their pigmentation and physical proportions have been recorded in the catchment area of lowland streams in Eastern Slovakia (unpublished data). Some of literature sources (Csipkés & Stündl, 2015;Movchan, 2011;Szepesi & Harka, 2013) report the occurrence of this species near this territory. Kottelat and Freyhof (2007) even mention its occurrence in the Tisza basin and in the Danube itself up to Bratislava (capital city) in Slovakia. This investigation aimed at helping to clarify taxonomic issues, but surely it did not enable complete clarification. Consequently, the main objective of this study was to identify status of fish of the genus Sabanejewia in region of the middle Danube basin and adjacent areas on the model of Slovakia simultaneously based on morphological, microhabitat, and molecular approach, which has not been carried out up to present.

| Study area and samples collecting
Fish from nine sites in the Slovak territory and one site near the town of Vidin, Bulgaria (type locality of Sabanejewia bulgarica), were sampled for this study (Table 1). In addition, 14 voucher specimens (catalogue numbers NPM P6V 85,299,85,310,and 85,317)   .
Identification of specimens from Slovakia was based on external morphological characters and coloration patterns as reported by Kottelat and Freyhof (2007).
During the sampling, selected microhabitat parameters were recorded using point sample method (Copp & Peňáz, 1988) modified according to Pekárik, Koščo, and Švátora (2012). At each sample point, where Sabanejewia specimen was present, four microhabitat variables were recorded: water depth to the nearest centimeter; wetted width; average velocity taken in 5-s interval measured 5 cm above the bottom using of flow probe (Valeport Flow Meter, Valeport Ltd.) and substratum type classified to categories as follows: silt, mud, clay, sand, gravel, pebbles, cobbles, boulders, and bedrock according to Pekárik et al. (2012). Due to the low abundance of Sabanejewia specimens at some sampling sites (Bodrog, Kysuca and Latorica Rivers), fish from previous samplings without evaluating the microhabitat parameters were also included to replenish the material for morphological studies.
Immediately after capture fish were anaesthetized, individually labeled and fin clip was taken and stored in 96% ethanol for later molecular analyses. The specimens were placed in labeled plastic bottles and preserved in 6% of formaldehyde solution. Voucher specimens are stored at the Department of Ecology of the University of Prešov (Slovakia).

| Morphological analyses
Since a preservation can cause deformations on the body shape and hence to affect final morphological analysis (Sotola et al., 2019), all measurements were taken at least after 3 months of their preservation. To minimize any ontogenetic differences and conservation bias, only well preserved sexually identified adult specimens (SL > 55 mm) ; own findings) were used for our study. Before each measurement, fish were placed into the cold water for at least 24 hr. Then, a total of 26 morphometric characters (including SL and TL) were measured on the left side of body ( Figure 2) to the nearest 0.01 mm using a digital caliper. To avoid any bias, all measurements were made point to point by one author. In order to minimize the resulting measurement error, each measurement was repeated three times and subsequently averaged (Morinaga & Bergmann, 2017).
Despite of fact that several significant differences occur between males and females of the genus Sabanejewia (Bohlen, 2008;Nalbant, 1963;Vasiľeva & Vasiľev, 1988), there was considerable overlap between both sex groups character ranges.
Moreover, we have assumed a significant impact of local habitat conditions on body shape independent of sex. Therefore, sexual dimorphism was not expected to affect the results.
In addition to morphometric measurements, 12 meristic parameters were counted (Table 6). Fin rays were counted under the light microscope with sufficient zoom. The last two unbranched rays in dorsal and anal fin, which articulate on single pterygiophore, were counted as "1 1 / 2. "

| Molecular analysis
Total genomic DNA was extracted from a small piece of the pectoral fin by a commercial kit (GT300, Geneaid). The entire sequences of cytochrome b (1,140 bp) were amplified by polymerase chain reaction (PCR) with primer pair GluDG.L (Palumbi, 1996) and H16460 (Perdices & Doadrio, 2001 Marešová et al., 2011;Perdices et al., 2003) were added for comparison with our samples. Detail list of all studied taxa, their haplotype classification to sampling sites, haplotype frequencies, and GenBank accession numbers are shown in Table S1.

| Data analysis
Multivariate normality was tested by visualization of morphometric variables (MVs) through the histograms and Mahalanobis multivariate QQ-plot. Before analyzing, morphometric dataset was standardized by arcsine square root transformation in terms of its percentage character. All statistical analyses were performed in R statistical software ver. 3.5.2 (R Core Team, 2019) using functions of packages morphoTools (Koutecký, 2014), MASS (Venables & Ripley, 2002), vegan (Oksanen et al., 2013), and pairwiseAdonis (Arbizu, 2017).
Visualization of PCA scatterplot was conducted by functions of package ggplot2 (Wickham, 2016).
For phylogenetic reconstructions and delimitation of boundaries within the DB complex, all forward and reverse sequences were assembled, edited, and aligned using the Seqman module (Lasergene v15) and also were checked by eye. Furthermore, as a final quality control, cyt b sequences were translated to verify that they were free of stop codons, frame-shifts, and gaps. The genetic dataset was analyzed by Bayesian inference (BI) using MrBayes 3.1.2 (Ronquist et al. 2012), the maximum-likelihood (ML) method using PhyML (Guindon et al., 2010), and neighbor-joining algorithm (NJ) using PAUP* 4.0B.10 (Swofford, 2002). The best-fit model of molecular evolution was determined for mitochondrial dataset using the Akaike Criterion (AIC) in Modeltest ver. 2.1.4 (Posada, 2008).
MrBayes was run with six substitution types (nst = 6) and considered gamma-distributed rate variation and the proportion of invariable positions (GTR + G + I). For BI, we ran four simultaneous Monte Carlo Markov Chain (MCMC) for two million generations and sample frequency every 100 generations. The first 5,000 trees were excluded as burn-in. The remaining trees were used to compute a 50% majority rule consensus tree. For ML analysis, we conducted heuristic searches under a GTR + I + G. For NJ analysis, DNA distance was calculated using MEGA 7 (Kumar, Stecher, & Tamura, 2016).

| RE SULTS
The mean values and standard deviations of morphometric characters expressed in relations to SL, H, and c for studied populations are listed in Table 2. Broken-stick model detected first component axis to be suitable for PCA interpretation, since its percentage of explained variation was higher than broken-stick percentage. The first principal component (PC1) accounted for 49.5%, while the second (PC2) for 10.5% of the total variance explained. Morphometric characters with the highest absolute correlation were lengths of lb1, lb2, and lb3 to the first and hc, io, and aA to the second axis, respectively (Table 3). The

| Meristic and coloration
The number of fin rays did not show any significant differences between studied populations. Their number was almost constant with only minimal differences (Table 6). Based on coloration, two main groups of fish were formed. Populations from larger lowland streams (Bodrog, Danube, Latorica) were set aside, where the number of lateral and dorsal spots was significantly lower than in others.
However, great differences in number of spots were also found in individuals from the same populations (Table 6).

| Delimitation of golden loaches clades
The global cyt b dataset was analyzed using the STRUCTURE, bPTP, and mPTP to ascertain the DB complex structure. The uppermost hierarchical level of structure was two clusters at K = 12 and K = 14 suggested STRUCTURE Harvester analysis (Figures 7 and 8). At both K, this analysis indicated nine distinct groups (Figure 9; Figure S1 and Table S4) in agreement with mitochondrial network ( Figure 5).
The species delimitation methods bPTP and mPTP recognized the same number of candidate species in agreement with sublineages designation (Figures 5 and 6 (Table S3). Graphic representation of the mutual relations within DB complex (Figures 5 and 6) is a majority consensus based on the results of the network reconstruction, phylogenetic, and delimitation methods, and therefore, group 4 of sublineage III is not supported more (the hatched network design).

| D ISCUSS I ON
In general, the variation in Slovakian populations observed by ordination analysis ( Figure 3) and compared to both samples from type localities for S. balcanica (Treska estuary in Skopje, MK) and S. bulgarica (Danube River in Vidin, BG) reflects great morphological variation within the genus Sabanejewia distribution. Similar largescale variation of populations referred to as species S. balcanica has been observed in Romania (Bănărescu, 1966;Bănărescu et al., 1972;Iftime, 2002) and Croatia ). Most of the morphometric and meristic and coloration traits exhibit wide range of variability.

TA B L E 7 Importance of environmental variables used in RDA analysis
Relatively distant position of our populations from Laborec and Vlára Rivers in PCA scatterplot (Figure 3) is mainly due to very short length of barbels of these specimens. These sites were the only ones, where boulders substrate type was dominating. Similarly, short barbels in relation to faster water velocity and stony substrate type were reported by Vasiľev (1988, 2019) for population of Sabanejewia kubanica in Kura River (Russian Federation).
The remaining morphometric characters used in our study have not been shown to be of significant use in distinguishing individual populations. However, the character loadings of PCA (Table 3) revealed several similar identifying features for bulgarica and balcanica-like populations as reported in several previous studies (Bănărescu et al., 1972;Iftime, 2002;Oliva et al., 1952;Sivkov, 1991;Vasiľeva & Vasiľev, 1988). Toward lowland populations, head length (c), preanal postdorsal (pD) distance are increasing toward foothill morphotype populations. However, in our study we did not confirm the significant difference in body depth (H) reported by several authors (Bănărescu et al., 1972;Iftime, 2002;Kottelat & Freyhof, 2007;Sivkov, 1991) as one the main discriminatory morphometric features. Surprisingly, the highest value of this character was observed in populations from Kysuca and Treska Rivers, that is, typically balcanica-like morphotype (Table 2). In our study, the character of body depth was constantly measured at the origin of dorsal fin.
The typical bulgarica-like "hump-backed" appearance is most pronounced on the body at the level of behind the head. Iftime (2002) however reported that this "hump-backed" appearance is also considerably variable and is related to breeding conditions. By author, ovigerous females also present distend abdomen, which adds to the overall appearance of body depth. In our case, most of the specimens from Kysuca River were sampled at the beginning of summer, which marks the spawning period for Sabanejewia sp. (Juchno & Boroń, 2012), while the other populations were mostly sampled in postspawning period. Therefore, the idea of spawning period impact on the body depth can be explained. Track changes in this and other characters between pre-and postspawning period should be a subject of further observations.
In terms of fin rays, our results correspond to previous published data of their numbers within the Central European (Mišík, 1958;Oliva et al., 1952) or Balkan populations (Bajrić et al., 2018;Buj et al. 2008;Sivkov, 1991;Šumer & Povž, 2000). Their number is almost constant in all observed populations, and small deviations between results of individual studies may be due to different counting methods and techniques. The only one more significant difference was observed in a few specimens (Ipeľ, Laborec, Torysa Rivers), in which up to four spines in anal fin were found. So far this number has been reported only by Witkowski (1994) in S. baltica. In this case, it is necessary to emphasize the need to use a microscope with a sufficient zoom as well as the need of skin disruption at the location of the fin origin. Some of the spines are of a very short length and also hidden in the skin, making them difficult to observe.

F I G U R E 5
The unrooted TCS haplotype network for the sublineages I-VI of the Danubian-Balkanian complex based on sequences of the cyt b. The haplotype numbers refer to Table S1. The node sizes are proportional to haplotype frequencies. Haplotype numbers from type locality in Vidin, Bulgaria, are highlighted in red The results of our study prove that the variability within morphology does not reflect groups created from molecular analyses.
On the contrary, one of the most important factors affecting the body shape of these small bottom-dwelling fish is likely represented by local habitat conditions, which are a result of long-term hydrological conditions at a given site. Therefore, the wide spectrum of morphological variability within the Sabanejewia populations in Danube basin could also be understand as a phenotypic heterogeneity among populations caused by diverse environmental characteristics. After analyzing several populations of Sabanejewia in Croatia, Buj et al. (2008) came to a conclusion that similar ecological factors are most likely a reason for a morphometrical similarities between populations. The specimens from rivers forming parts of different watersheds but having similar habitat conditions were more uniform than the others.
Our results indicate that morphotype of S. bulgarica is bound by its occurrence to larger and deeper lowland rivers with slow velocity and fine substrate bottom. Comparable results have also been reported from the Romania, Bulgaria, or Hungary (Bănărescu et al., 1972;Iftime, 2002;Sivkov, 1991;Stefanov, 2007) as well as from lower courses of the rivers in Central Asia, where specimens of species S. aurata have also some lowland morphotype features (reduced body pigmentation and deeper body) similar to S. bulgarica description (Bănărescu et al., 1972). The position of specimens F I G U R E 6 Bayesian consensus tree resulting from the analysis of the cyt b data in studied golden loaches taxa with Bayesian posterior probabilities/ML bootstrap/NJ bootstrap values listed near the nodes. Only values > 75% are shown. Haplotype numbers from type locality in Vidin, Bulgaria, are highlighted in red  Sabanejewia baltica from Treska River close to lowland type populations in PCA analysis ( Figure 3) may be due to the nature of microhabitats on this site. The sampling locality on this river was situated near the estuary to Vardar River (Marešová et al., 2011), which is relatively large river in this area. Hence, the local ecological conditions on this site can be similar to the lowland streams, where typically bulgarica-like morphotype occurs. Therefore, we consider it appropriate to confirm this theory also through a comprehensive study of Sabanejewia populations in the Vardar basin. These conclusions also lead us to claim that body shape of several Sabanejewia populations reflects only phenotypic adaptation to diverse habitats. Generally, fish morphology as a manifestation of phenotypic plasticity is a well-known phenomenon due to diversity of environmental factors (Keeley, Parkinson, & Taylor, 2006;Laporte, Claude, Berrebi, Perret, & Magnan, 2016;Ramler et al. 2016;Senay, Boisclair, & Peres-Neto, 2014). Phenotypic variability among populations may arise without major genetic differentiation when they occupy heterogeneous habitats across their distribution range (Cheng et al., 2017;Colihueque, Corrales, & Yáñez, 2017).
When comparing two main species concerned of this study (S. balcanica vs. S. bulgarica) based on molecular analyses, it is necessary to point out the fact that most of previous studies focused on resolving the taxonomic status did not include samples from their terra typica (Bartoňová et al., 2008;Perdices et al., 2003;Buj et al. 2008). Our results comprising samples from both of these species have shown that the haplotypes of S. bulgarica population from the type locality are spread across almost all haplotype groups in Slovakia and they are also clustered with most of the samples from Danube basin previously considered as species S. balcanica Halačka, Muška, Mendel, & Vetešník, 2017;Perdices et al., 2003). All phylogenetic and delimitation methods used reliably differentiated the two species and, at the same time, drew our

F I G U R E 8
The ΔK plot describing the rate of change in the log probability of the data between successive K values from 1 to 17. The modal value of this distribution is the true K, or the uppermost level of hierarchical structure  III3  V  I I  V I  I II1  III4  III2  I  I Table S1 attention to new areas of their occurrence. The new description of distribution of the haplotypes of both the species is in contradiction to the general hypothesis of the dominant position of S. balcanica in region of middle Europe and Balkans (Marešová et al., 2011;Perdices et al., 2003). However, populations containing haplotypes of Lineage II (sensu S. balcanica) typical for Aegean Sea basin can also be found in the peripheral part of the Danube basin (Marešová et al., 2011) ( Figure 10). More precise determination of the border line of occurrence or confirmation of hybrid individuals of both the species will require further investigation especially that performed using the nuclear marker analysis.
We believe that recent dispersion of variety of mtDNA haplotypes from the type locality of S. bulgarica throughout the Danubian corridor has took place probably during cyclical cold and warm periods in Pleistocene glaciations as reported by Perdices et al. (2003) for the whole DB complex clade. However, much more detailed phylogeographical analysis must be performed for determination of various parameters of distribution, for example, in how many waves, in which numbers of individuals, etc., but this goes beyond the extent of this study. These glaciations played an important role in secondary recolonization from the Danube refuge (Seifertová, Bryja, Vyskočilová, Martínková, & Šimková, 2012;Sommerwerk et al., 2009) leading to low genetic homogenization of freshwater species in this region (Perea et al., 2010). This fact is also most probably the cause of low genetic distances (Table 8) and simultaneous presence of haplotypes of different sublineages of the DB complex at some localities within the Danube basin (Bartoňová et al., 2008;Buj et al. 2008). At present, the relatively short elapsed time from forming the current state of the Danube basin (approximately 700,000 years ago) (Hsü, 1978) and since the establishment of DB complex within Sabanejewia genus (Pleistocene period) (Perdices et al., 2003) was not enough to make the genetic distances between lineages more pronounced. However, the blending of haplotypes from type locality of S. bulgarica occurring only within the most diverse sublineage III of DB complex is the basis of claim that populations of golden loaches previously referred as species Sabanejewia balcanica (Karaman 1922) within Central Europe and Balkan region are closer to naming Sabanejewia bulgarica (Drensky, 1928). This is also underlined by fact that morphotype of these fish is very diverse, strongly dependent on local habitat conditions and thus does not allow unambiguous determination based on external morphological features.

| CON CLUS ION
Our results demonstrated a high degree of morphological variability among the studied populations of the genus Sabanejewia, which is mainly caused by the adaptation of these fish to the ecological conditions on a given habitat. The body shape and coloration pattern in diverse environments reflects local microhabitat conditions and is thus a manifestation of significant phenotypic plasticity. From a phylogenetic point of view, this issue can be characterized as a previously mentioned complex (Perdices et al., 2003) that is currently still in the process of evolution and clear allocation of its species is difficult.
We confirmed that none of the Vardar haplotypes (representing species S. balcanica) have been found among Slovakian or other samples included in the sublineage III. Oppositely, haplotypes from Vidin (type locality for S. bulgarica) occurred within the sublineage III of Danubian-Balkanian complex (Perdices et al., 2003) as well as Slovak samples. All these findings form the basis of the claim that populations of golden loaches within the middle part of Danube basin and adjacent regions are closer to name S. bulgarica. However, taxonomically there is also Vladykov's description of Sabanejewia montana from the mentioned area (Šanda, Vukić, & Švátora, 2010), whose validity could also be reassessed on the basis of further analyses.
In further studies, we suggest a comparison of the biological indicators such as growth differences, fecundity, or more complex molecular studies (nuclear or microsatellite markers) of the DB complex.
These could lead to further important knowledge and clarification of this complex issue.

ACK N OWLED G M ENT
The study was implemented within the projects VEGA 0916/17 and GaPU 32/2018. The authors would like to thank (in alphabetical order) Juraj Hajdú, Karel Halačka, Yuliya Kutsokon, Sebastián Pavlinský, Tihomir Stefanov, and Ľubomír Šmiga for their assistance in the field or during laboratory work. Especially, we would like to thank Radek Šanda for lending us the voucher specimens of Sabanejewia balcanica from the collection of National Museum in Prague.

CO N FLI C T O F I NTE R E S T
The authors declare no competing interests.