Inter-population variability within Sabanejewia populations from the western Balkans, and their phylogenetic position in comparison to other European populations were investigated. Of 79 samples analysed, 51 unique haplotypes were identified. Network analysis divided thirteen populations from five river basins into two clusters: cluster I was composed of populations from the Vardar drainage and tributaries of the neighbouring River Morava (River Danube basin), while cluster II contained the River Timok (eastern Serbia) and all the River Sava populations. The only locality that housed haplotypes of both clusters was the River Kutinska reka in the upper Morava basin. When compared with the haplotypes reported in former studies, both clusters fell within the ‘Danubian-Balkanian complex’. Cluster II was included in the most heterogeneous sub-lineage S. montana — S. bulgarica — S. balcanica (III), while cluster I was related to the sub-lineages S. doiranica — S. balcanica (II) and S. balcanica (VI). Recently published haplotypes from the Croatian Sava (rivers Petrinjčica and Rijeka) and Drava (rivers Drava and Voćinska) basins, as well as Czech and Slovak individuals from the Danube and Tisza river basins were included. The River Drava showed the same population subdivision as the River Kutinska reka.
Loaches of the genus Sabanejewia, Vladykov, 1929 are small, economically unimportant species, widespread in the western Palaearctic. They occur in the Ponto-Caspian basins, some Baltic Sea rivers and in several river basins of the Aegean and North Adriatic Sea watersheds (Kottelat & Freyhof 2007). Altogether, eight species are now recognised in Europe, of which only two species are widespread (Kottelat & Freyhof 2007). Sabanejewia baltica Witkowski, 1994 occurs in northern Black Sea river drainages and in southern Baltic Sea rivers. Sabanejewia balcanica (Karaman, 1922) is known from the Danube basin and from several Aegean Sea river drainages, namely Marica and rivers from Gallikos to Pinios, including the River Vardar. Three other species occur in the River Danube basin: Sabanejewia bulgarica (Drenski, 1928) in the lower Danube, Sabanejewia romanica (Bacescu, 1943) in several Danube inflows in Romania and Sabanejewia vallachica (Nalbant, 1957) in the River Ialomita and the lower Siret basin (Kottelat & Freyhof 2007).
Sabanejewia balcanica mainly inhabits middle and upper parts of medium-sized rivers and the populations are often isolated (Ludwig et al. 2000, Deli´ et al. 2003, Pekárik et al. 2008). However, its presence in rather large rivers is also known (Mrakovčić et al. 2000). This species has attracted attention during the last decade as a model for phylogenetic and phylogeographic studies (Ludwig et al. 2000, Perdices et al. 2003). Genetic diversity and phylogenetic relationships of S. balcanica and related species have been studied in populations mainly from Romania and Greece (Ludwig et al. 2000, Perdices et al. 2003), Czech Republic and Slovakia (Bartoňová et al. 2008) and Croatia (Buj et al. 2008).
Based on mtDNA analysis Perdices et al. (2003) defined inside the genus Sabanejewia six evolution lineages. One of these, the so-called Danubian-Balkanian complex, contains all samples of S. balcanica, S. vallachica and S. bulgarica from most parts of the Danube basin and Balkan Peninsula drainages. This lineage is further subdivided into six sub-lineages: (I) S. vallachica — S. balcanica from the eastern Romania, (II) S. doiranica — S. balcanica from western Greece, (III) S. montana — S. bulgarica — S. balcanica from most tributaries of the middle Danube basin, (IV) S. radnensis — S. balcanica from the Mures system, (V) S. thracica — S. balcanica from Greece, and (VI) S. balcanica from the Mur system. Perdices et al. (2003) suggested Pleistocene origin of the Danubian-Balkanian complex with ongoing vicariant separation and a recent range expansion of the single lineage S. balcanica due to its presence in all sub-lineages. Nevertheless, genetic diversity and relationships of S. balcanica populations from some parts of its distribution area remain unknown. Sabanejewia balcanica is widespread in the western Balkans part of the Danube drainage (unpublished data), and is also present in the neighbouring Vardar drainage (Grupče & Dimovski 1973, Economidis & Nalbant 1996). Perdices et al. (2003) included into their analysis samples from several drainages in Greece (Evros, Pinios and Lake Doirani, which belongs to Vardar basin), while Buj et al. (2008) analysed several populations from western Croatia. Other sites from the western Balkans have not been investigated. Gaining an understanding of the relationships among Sabanejewia populations from the western Balkans and the proximate Vardar drainage is important for identification of the level of isolation, the taxonomic status and proper conservation management of the species in this area.
The aims of the present study were to define this genetic diversity, and identify and evaluate relationships among S. balcanica populations from the western Balkans and among previously identified sub-lineages in the Danubian-Balkanian complex sensu Perdices et al. (2003).
Material and Methods
A total of 79 individuals were newly sequenced from 14 localities in two river drainages: Danube drainage (Sava, Morava and Timok river basins) and Vardar drainage (Fig. 1, Table 1). Locality numbers, corresponding to those in Fig. 1 and Table 1, are given in parentheses in the text after the locality name, wherever necessary. Fish were collected by portable battery electroshocker. Small pieces of fin tissue were preserved in 96 % ethanol for DNA analysis, while voucher specimens were labelled, preserved in 5 % formaldehyde and later removed to 70 % ethanol. Voucher specimens are deposited in the ichthyological collection of the National Museum, Prague, Czech Republic, under catalogue numbers NMP P6V 83745, 83747, 83749-51, 83832-5, 83837-9, 84006-7, 84710-5, 84846-7, 84936-7, 84940-2, 84944-5, 85018-20, 85023-7, 85132, 85134-5, 85137-8, 85297-8, 85303, 85305-8, 85874, 85876-7, 85879, 85882-5, 85888-9, 86293-4, 86296-9, 86301-3, 86306, 86456-60, 86693-4 and 86697-9.
DNA extraction, PCR amplification and sequencing
Genomic DNA was extracted from fin tissue by use of the phenol-chloroform-isoamyl alcohol method of Sambrook et al. (1989) or JETQUICK Tissue DNA Spin Kit (GENOMED) following manufacturer's instructions. The mitochondrial protein-coding gene cytochrome b was amplified in its entire length of 1140 bp. Two pairs of primers were used: L14735 and H15149 under the conditions described in Ludwig et al. (2000); GluF and ThrR according to the conditions described in Machordom & Doadrio (2001). The PCR products obtained were purified via precipitation with PEG/Mg/NaAc. Sequencing was carried out on an ABI 3730XL DNA sequencer by Macrogen (Seoul, Korea) and ABI PRISM 310 by IVB (Brno, CZ) with the same primers or with internal primers L-Cyp and H-Cob (Buj et al. 2008). All PCR products were multiply sequenced in both directions to ensure high quality reads. Unique haplotypes were submitted to GeneBank under accession numbers HQ291787– HQ291838.
Overview of samples analysed, collection localities (Croatia (CRO), Serbia (SRB), Bosnia and Herzegovina (BIH), Macedonia (MK)) and haplotypes found (numbers in brackets represent occurrence of haplotypes in case of more than one individual).
Sequence and phylogenetic analyses
All sequences were aligned and edited manually in programme MEGA 4 (Tamura et al. 2007) and translated to amino acid sequences for verification of correct alignment. No stop codons were found within the sequences. To detect the presence of saturation, a plot of uncorrected against corrected distances was constructed. Intra-specific relationships between investigated populations were estimated, with inclusion of nine previously published sequences from Buj et al. (2008), by network analysis using the Median Joining algorithm in Network 22.214.171.124 (Bandelt et al. 1999). According to the structure of the network 27 haplotypes were selected and employed in phylogenetic analysis, together with sequences determining the Danubian-Balkanian complex (16 sequences from Perdices et al. 2003), with other previously published Balkan sequences from mainland Croatia (9 sequences from Buj et al. 2008) and sequences from the Czech Republic and Slovakia (20 sequences from Bartoňová et al. 2008). Furthermore, fourteen sequences of six other Sabanejewia species from Perdices et al. (2003) were used as close outgroups to root the Danubian-Balkanian complex, while Misgurnus fossilis (AF263097; Perdices & Doadrio 2001) was used as a distant outgroup to root Sabanejewia. Phylogenetic trees were constructed in PAUP* ver. 4.0b10 (Swofford 2002) using maximum parsimony (MP) and neighbour joining (NJ) algorithms. For NJ the best model of sequence evolution TrN + I + G (Tamura & Nei 1993) was selected according to Akaike information criterion (AIC) implemented in Modeltest 3.7 software (Posada & Crandall 1998). All analyses were undertaken using heuristic search with the TBR branch swapping (tree bisection reconnection) algorithm with all codon sites and with transitions/transversions weighted equally. Nonparametric bootstrap analysis was performed to evaluate robustness of topology (1000 replicates for NJ; 500 for MP) with a cut off level of 50 %. For Bayesian analysis four Monte Carlo Chains were run simultaneously in parallel MrBayes version (Altekar et al. 2004) for 5000000 generations with a sampling frequency of 100. The first 25 % trees were discarded, and the resulting trees taken for 50 % majority consensus tree. All analyses were conducted on a computation cluster at IVB AS CR (Brno, CZ).
Some statistical parameters (number of haplotypes, haplotype diversity, nucleotide diversity, Strobeck's statistic, number of polymorphic sites) were estimated using DnaSP v5 (Librado & Rozas 2009). Genetic distances between populations in the main river basins (Sava, Morava, Vardar and Timok in this study, and Drava, Tisza and Danube with data from Buj et al. 2008 and Bartoňová et al. 2008) were computed with MEGA 4 using the evolution model reported above.
Molecular diversity of examined populations (number of specimens (n), haplotypes (h) and polymorphic sites (S); nucleotide diversity (Pi), haplotype diversity (Hd), Strobeck's statistic (St); -, impossible to estimate the value).
Genetic distances between drainages and river basins computed from TrN + I + G. Values below the diagonal are genetic distances; those above the diagonal are standard errors.
In 79 newly sequenced individuals, 7.1 % (81) sites were variable and 4.7 % (54) parsimony informative from a total of 1140 characters. Variability in the protein level was 4.2 % (16 variable amino acids of 380), and 1.6 % (6) parsimony informative. Sequence variation was mainly due to third position substitutions and transitions. Nucleotide composition (T = 30.2 %, C = 27.3 %, A = 26.0 %, G = 16.5 %) showed anti-G bias typical for the mitochondrial genome. A plot of uncorrected (p-distance) against corrected distances did not reveal any trend towards saturation (plot not shown).
Out of 79 sequences 51 unique haplotypes were identified, giving an overall Hd of 0.9864 ± 0.005. Thirty-three haplotypes were each found in a single individual, but only six haplotypes were present in more than three individuals (Table 1). Nucleotide diversity Pi (π) varied within populations between 0.0001 ± 0.0031 to 0.0130 ± 0.0028. High values of haplotype diversity were obtained in all examined populations (Table 2).
Intra-population genetic variability did not exceed 1 % except in the River Kutinska reka (12). Interpopulation genetic distances were computed both at population (not shown) and river basin levels (Table 3). The largest distance was revealed between Vardar drainage and Timok river basin (1.85 %) and the smallest between Sava and Timok river basins (1 %), with that between Sava and Danube basins resolving in between (1.03 %). The closest populations within the Sava basin were, as expected, the River Ilova (2) and its tributary, the River Rijeka (1) (0.43 %). In contrast, populations from rivers Petrinjčica (3) and Usora (7) proved to be the most distant (1.53 %).
Firstly, we employed in the analysis only 79 newlysequenced individuals, among which two distinct clusters, differing in nine substitutions, were revealed (Fig. 2). The first cluster included all samples from Vardar drainage and Zapadna Morava river basin and half of the individuals from the Južna Morava river basin (three out of seven samples from the River Kutinska reka (12)). Cluster II included all populations from Sava river basin (1–9) and from the River Timok (10). Some populations from the River Sava tributaries formed separate branches (Jadar — Drina (8), Usora — Bosna (7), Vrbanja — Vrbas (6)), while the remaining populations were mixed, creating three other subgroups in which the samples were equally distributed. The most divergent position in the River Kutinska reka (12) was identified, being the only population whose haplotypes occurred in both clusters in the network.
Subsequently, previously published Croatian haplotypes from Sava (the River Petrinjčica, KUP; the River Rijeka, RI1, RI2, RI3) and Drava basins (the River Voćinska, VO1, VO2; the River Drava, DRA1, DRA2, DRA3) (Buj et al. 2008) were entered into network analysis. All individuals from the Sava basin were related to specimens of cluster II, and with two haplotypes from the River Drava (DRA2, DRA3). Interestingly, one haplotype DRA1 and both haplotypes from the River Voćinska were included in cluster I (Fig. 2).
Based upon the structure obtained from the network, 27 individuals were chosen for phylogenetic purposes. To identify the phylogenetic position of the studied populations within the context of other European localities sequences from these populations were compared with those determining the Danubian-Balkanian complex sensu Perdices et al. (2003). In addition, populations from the Czech Republic and Slovakia, forming groups 1–5 from Bartoňová et al. (2008) (Danube basin: rivers Vlára (CZ), Malý Dunaj and Ipel' (both SK); Tisza basin: rivers Olšava/Hornád, Olšava/Topl'a, Ublianka and Ondava (all SK)), and sequences mentioned above (Buj et al. 2008) were included. All methods for constructing phylogenetic trees were congruent. MP analysis generated a single MPT with L = 725, retention index (RI) = 0.816 and consistency index (CI) = 0.618. A Bayesian phylogenetic tree, with values of Bayesian posterior probabilities/NJ bootstrap/MP bootstrap is avalible on http://www.ivb.cz/download/strom_1.pdf.
All populations of cluster II were included within the heterogeneous sub-lineage S. montana — S. bulgarica — S. balcanica (III), that related specimens mainly from middle Danube tributaries. Additionally, most Czech and Slovak specimens (groups 1–4) were included in sub-lineage III. Individuals from cluster I were included in two other sub-lineages. Vardar drainage and Zapadna/Južna Morava river basin samples clustered within sub-lineage S. doiranica — S. balcanica (II) from western Greece, and the River Drava samples were included within sub-lineage “S. balcanica” (VI) from the River Mur flowing from Austria. The remaining group (5) from Tisza river basin was most related to the S. radnensis — S. balcanica (IV) sub-lineage described from the Mures system in Romania.
In the present study low to high levels of nucleotide diversity (0.0001–0.0130) and high values of haplotype diversity (0.778–1.000) were observed.
The Hd parameter was most diverse in the River Ilova (1.000 ± 0.500), probably a result of the small number of individuals included in analysis (also a likely reason for the extremely low value of nucleotide diversity). In the case of two populations (rivers Crni Timok (10) and Lepenac (13)), it was not possible to estimate the values due to absence of polymorphic sites (by Pi) or presence of only one haplotype (Hd). The most diverse population was found to be the River Kutinska reka (12) with the highest value of nucleotide diversity and most polymorphic sites. Its division into subpopulations is reflected in a low value for the Strobeck (1987) statistic (S = 0.325), which expresses the probability of obtaining a sample with equal or fewer haplotypes than the number observed. This finding was confirmed by network analysis: three samples were present in cluster I and four in cluster II. The same pattern was observed in the River Drava population, where samples ranked in both cluster I and cluster II. The inclusion of specimens from the same populations into different sub-lineages was also observed in the River Olšava/Topl'a from the Tisza basin (Bartoňová et al. 2008). Presence of haplotypes belonging to different sub-lineages, often geographically distant, at one locality, indicates recent genetic flow or would be caused by incomplete lineage sorting, i.e. the persistence of ancient polymorphism which leads to incomplete isolation of species. Most probably the population discontinuities of this rather cold water adapted species in the lowlands during the time with not optimal conditions would have disrupt the gene flow for a certain period, causing diversification of the isolated populations (Englbrecht et al. 2000, Šedivá et al. 2008). These populations went later again into contact. According to the position of the shared haplotypes in the phylogenetic tree, the gene flow seems to be more probable (Joly et al. 2009). If the shared haplotypes would result from incomplete lineage sorting, their position should be basal to the rest of haplotypes in the lineage (Joly et al. 2009). On the other hand, the presence of variability between populations from different drainages or river basins is evident from the network analysis, with separate branches (Vardar drainage, Morava and Sava river basins). Within the Sava river basin, subgroups corresponding to partial catchment areas (Drina, Bosna and Vrbas river basins) were formed. Although gene flow between clusters I and II appears to have occurred, it has not been strong enough to remove variability between river basins. The same pattern was found by Perdices et al. (2003). A single lineage S. balcanica is included in all mitochondrial lineages of the Danubian-Balkanian complex, which suggests ongoing vicariant separation but the elapsed time since the establishment of the different sublineages (I–VI) was insufficient to consolidate genetic differences (Perdices et al. 2003). This is important for conservation management of S. balcanica, as any stocking should use populations originating from nearby localities.
Our results support adouble origin of Aegean populations. We ascertained the western Greek populations from Pinios and Vardar have close affinities to the River Morava population while Perdices et al. (2003) found out the Evros population is closely related to the River Vit population (Danube inflow) from Bulgaria. In this regard, there is geological evidence of an ancient direct connection between Danubian tributaries and the River Vardar (Economidis & Bănărescu 1991). On the other hand, there is no evidence for connection between the River Vit (or any other of Danubian tributaries) and the River Evros (Angelova 2003). However, river capture of the headwater stream can explain the observed close relationships of Sabanejewia from both river basins.
Inter-population variability in the S. balcanica complex varied between 0.4 % and 2.3 %. Within the Sava catchment area we obtained similar values (0.4 %-1.5 %) to those obtained by Bartoňová et al. (2008) for the middle-Danubian populations (0.6 %-1.5 %). Phylogenetic analysis integrated specimens from cluster II, together with groups 1–4 of the Czech and Slovak samples, into sub-lineage III of Perdices et al. (2003), indicating a close relatedness (low genetic distance of 1.03 % between rivers Danube and Sava) of Danube and Balkan ichthyofauna, previously observed in cobitids (Perdices & Doadrio 2001) and cyprinids (Zardoya et al. 1999).
Inclusion of populations from the western Balkan area and central Europe into four of six sub-lineages of Perdices et al. (2003) (Sava river basin, sub-lineage III; Morava river basin, II, III; Vardar drainage, II; Drava river basin, III, VI; Danube drainage, III, IV) would imply recent or extant gene flow between these sub-lineages, resulting into genetic homogenisation of the species across Europe. This is in accordance with Perdices et al. (2003) who suggested recent genetic homogenisation during Pleistocene glaciations because of low sequence divergence among all mtDNA sublineages in the Danubian-Balkanian complex. Further detailed phylogeographic research throughout the whole distribution area of S. balcanica is necessary to clarify the taxonomic status of the numerous nominal taxa and to identify evolutionary significant units for conservation purposes.
We are grateful to Ivan Bognt from University of Mostar for his help with obtaining permissions to work in Bosnia and Herzegovina and to the anonymous referees for their remarks, which contributed to improvement of this manuscript. We thank Iain Wilson for language correction o f the article. This research was supported by project no. DE06P04OMG008 of the Czech Ministry of Culture and grant numbers VaV-SM/6/3/05 from the Ministry of Environment of the Czech Republic and 1/2360/05 VEGA. This work was further supported by Ministry: of Science and Technological Development of the Republic of Serbia (Grant No. 173045). Sequencing was supported within the framework of research project M200930901 supported by the programme o f internal support for international collaborative projects of the Academy of Sciences of the Czech Republic.