Diversity of spined loaches from Asia Minor in a phylogenetic context (Teleostei: Cobitidae)

Accurate determination of species diversity in areas of high endemicity, particularly those lacking comprehensive systematic knowledge, represents a challenge for both taxonomists and conservationists. This need is particularly evident in areas greatly affected by anthropogenic disturbances such as the Eastern Mediterranean and its freshwater environments. To improve our knowledge of Eastern Mediterranean freshwater fishes, we phylogenetically studied Western Palearctic Cobitis species, focusing on those found in Turkey. Overall, our results provide a robust framework to assess the number of species of Cobitis. Phylogenetic reconstructions based on mitochondrial (cyt b) and nuclear (RAG1) sequences show seven major clades (Clades 1–7) grouping all Western Palearctic Cobitis species, except C. melanoleuca. In general, each major clade comprises Cobitis species that inhabit geographically close areas and have similar secondary sexual characters. Multiple divergent lineages were identified in our analyses, some of which were highly divergent such as the ones inhabiting Turkish freshwaters. Moreover, in some analyses, several of the identified lineages were incongruent with a priori defined species. Furthermore, our analyses identified eight potentially new candidate species, six that had been suggested in previous studies and two that are reported here for the first time. Our results reveal Turkey as the area with the greatest diversity of spined loaches in the Mediterranean.


Introduction
Accurate assessment of the number of species and correct species identification are fundamental to systematic biology. Correct species identification and knowledge of endemism patterns are particularly important for biodiversity studies, endangered species lists and the prioritization of regions for conservation [1,2]. Precise species recognition is especially challenging in groups exhibiting very low morphological variation but is crucial for adequate conservation strategies, particularly of geographical regions rich in endemic species under strong anthropogenic pressure, such as those in the Eastern Mediterranean [3].
Renewed interest in assessing biodiversity levels has been driven by the implementation of molecular COI DNA barcoding [4] and the emergence of new analytical methods for species ethanol in the field. Permission for sampling in Turkish waters was issued by the Ministry of Food of the Republic of Turkey.

Taxon sampling and data analyses
We identified species according to their morphology and geographical distribution following previous studies [13], [15], [17], [21]. We used external secondary sexual characters for Cobitis identification, specifically relying on a character known as the lamina circularis or scale of Canestrini, a thickening of the pectoral rays on Cobitis males. The absence, presence or duplication of this character provided the basis for identification. In the time since [14] published a global checklist of loaches, four new species have been described for the area. Therefore, in this study, we recognized a total of 46 valid species for the Western Palearctic [19], [21], [23] and molecularly analysed 41 of these species (Fig 1, S1 and S2 Tables). For phylogenetic analyses, we used C. choii as the primary outgroup [29] and Sabanejewia aurata, Oxynoemacheilus hanae, O. gyndes, Pangio pangia and Kottelatlimia pristes as more distant outgroups. If possible, specimens were sampled from multiple localities across the known range of each species. In order to place Eastern Mediterranean Cobitis species into a broader phylogenetic context, we first analysed the 41 selected species with other spined loaches of the Cobitis sensu lato group that are included in the Northern Clade defined by [29] (S1 Fig). The full-length mtDNA cyt b (1140 bp) and nuclear RAG1 (877 bp) genes were analysed as described by [30]. These genes were selected for phylogenetic reconstructions as previous analyses revealed both genes to be informative at different phylogenetic depths [28]. We combined data obtained from previous studies [27][28][29] with new data obtained in this study from Western Asia, particularly Turkey. When possible, we isolated both genes from the same specimen; otherwise, any available sequences were used for the analyses. A total of 211 new sequences were obtained and deposited into GenBank (for accession numbers, see S1 and S2 Tables).
Total DNA was extracted from fins using the ChargeSwitch gDNA Micro Tissue Kit (Invitrogen, Carlsbad, CA, USA). The conditions and primers used for PCR amplifications of cyt b and RAG1 were as previously described by [30,31]. Purified PCR products were sequenced directly using the primers used for amplification.

Phylogenetic analyses
We used the complete cyt b (N = 199) and complete RAG1 (N = 158) datasets for gene tree reconstructions. Phylogenetic trees were inferred for the complete datasets using Maximum Likelihood (ML) and Bayesian inference (BI) approaches. Robustness of inferred trees was assessed by bootstrapping (1000 replicates) in ML analyses and posterior probability (pp) values in BI analyses. The ML method was implemented in RAxML [32] through its graphical interface RAxMLGUI 1.3 [33], using the GTR+I+G model of evolution with 1000 bootstraps (BT). All trees were used to construct a 50% consensus tree in PAUP � v. 4.0a146 [34] for each dataset. The BI method was implemented in MrBayes v.3.1.2, partitioning the dataset by codon position for cyt b and under a GTR+I+G substitution model (nst = 6 and pinvGamma). We ran four simultaneous Monte Carlo Markov Chains (MCMC) for 2 million generations, sampling every 1000 generations, chain temperature 0.2. Log-likelihood stability was attained after 10,000 generations, and we excluded the first 100 trees as burn-in. The remaining trees were used to compute a 50% majority rule consensus tree in PAUP � . We performed a ML analysis in RAxMLGUI using a concatenated dataset of the two markers (N = 154), and the GTR+I+G model of evolution with 1000 bootstraps (BT). A BI analysis was also performed on the combined dataset, independently by gene, using the parameters described above. For the broader phylogenetic analysis of Cobitis species, we performed a ML analysis in RAxMLGUI (GTR+I+G,1000 BT) using an extended (N = 215) concatenated dataset (cyt b and RAG1) that contained species included in the Cobitis sensu lato group [29].
All taxa have the same number of codons for cyt b and RAG1 with no stop codons when translated to amino acid sequences. No gaps or ambiguous alignments were found, and all positions were used in the analyses. Sequences were checked and aligned with Sequencher ver. 4.8 (Gene Codes Corp., Ann Arbor, MI, USA).

Species delimitation and K2P sequence divergence
To identify potential species boundaries within our dataset, we used the Poisson tree processes (PTP) model proposed by [35], [36], and the maximum likelihood (ML) phylogeny generated with the complete cyt b dataset. We used two versions of the PTP model: bPTP, which adds Bayesian support (pp) values to branches that delimited species in the input tree, and the refined multi-rate mPTP. Both bPTP and mPTP are single locus species delimitation methods that use non-ultrametric phylogenies in which the number of substitutions are directly reflected by branch lengths. A key assumption is that the number of substitutions between species is significantly higher than the number within species. The RAxML gene tree was then used as the input for the PTP analyses as implemented on the PTP server, with 100,000 MCMC generations, a thinning of 100 and a burn-in of 0.1.
We used the complete cyt b dataset comprised of 199 spined loach sequences to calculate Kimura 2 parameter (K2P) distances [37] for direct comparison with the barcode (COI) reference database of Mediterranean freshwater fishes that includes multiple Cobitis species [10] [we excluded one C. bilineata specimen that exhibited a highly discordant COI sequence (Gen-Bank ID = KJ553116)]. Mean K2P distances (range and standard deviation) using the "pairwise deletion option" were calculated in Sequencer 6.1 (written by B. Kessing). We estimated K2P sequence divergence among the major lineages recovered in the RAxML gene tree and between species delimited with PTP methods. These analyses will allow us to compile values for species identified on the basis of morphology and for molecular lineages delimited by analyses of cyt b sequences that will then be used in subsequent studies to confirm potentially new Cobitis species based on COI sequences [10].

Phylogenetic analyses
In the extended phylogenetic analysis, Cobitis species described for the Western Palearctic were included in Subgroup II of the Cobitis sensu lato group recovered by [29]. This subgroup contains species that are widespread across Europe Asia Minor, the Black Sea and the Caucasus, with some lineages related to species restricted to East Asia (S1 Fig). All analyses resolved a Western Palearctic clade that included 40 previously described species of Cobitis from this region; the clade is further organized into seven major phylogenetic clades (Clades 1-7) with unresolved basal relationships (Figs 2-4, S1 Fig). Cobitis melanoleuca does not belong to the Western Palearctic clade but rather is more closely related to East Asian species. All seven clades were well supported in phylogenies derived from the complete cyt b and combined datasets (>87BT, >0.94pp). Only partial groups were supported in the more conserved nDNA phylogeny based on the complete RAG1 dataset (Table 1 and Fig 3).
Each clade was comprised of geographically close Cobitis species that display similar secondary sexual characters, namely in the number of lamina circularis (Fig 2), congruent with previous studies that had less geographic coverage of species [27], [29], [30], [38].
Nuclear gene analyses. The consensus nuclear phylogeny showed moderate to high support (97BT, 0.67pp) for a Western Palearctic clade but was generally less informative for major groups (Fig 3). Two of the seven major clades identified in the phylogenies based on mitochondrial and combined datasets (Clades 3 and 4) were recovered as monophyletic in ML and BI analyses of RAG1 with moderate to high support (>70BT, 1 pp). Clade 2 was not only recovered as monophyletic in any of these analyses. Clades 1, 5, 6 and 7 also were not  7), were recovered as independent lineages. Unlike in the mtDNA gene tree (see Fig 2), many of the Cobitis species resolved as non-monophyletic groups in the RAG1 gene tree. However, some of the highly divergent lineages recovered in the mitochondrial analysis, namely Cobitis sp. Aksu-Eğirdir (76BT, 1pp) and Cobitis sp. Menderes (99BT, 1pp), were also highly divergent in the nuclear phylogenies, supporting the idea that these lineages are potentially undescribed species.

Species delimitation and K2P genetic divergence
Although we obtained a wide range of K2P values with cyt b among the seven major clades, mean values between these clades were always greater than 8% (mean 11.5 ± 1.24%, range 8.92-16.20%), much higher than mean values between species, which were always less than 3% (mean 0.59 ± 0.52%, range 0-2.61%). In four of the major clades, we observed higher cyt b distances between certain lineages than between recognized species of Cobitis (K2P >2-3%).
None of these lineages clustered with any of the a priori defined species (Fig 2). Species delimitation analysis using the bPTP approach, and the phylogeny derived from the complete cyt b dataset, identified at least 49 weakly to strongly supported lineages (0.02-1pp), whereas the more conservative mPTP analysis identified 33 lineages. In many cases, individual lineages detected in the bPTP analysis corresponded to one of the morphologically identified species. In contrast, in the mPTP analysis, many of the a priori morphologically defined species grouped together within individual lineages.
Most of the highly differentiated lineages obtained in our study (i.e. Cobitis sp. Aksu-Eğirdir, Cobitis sp. Ceyhan, Cobitis sp. Gölbasi, Cobitis sp. Menderes, Cobitis sp. Sapanca, and Cobitis sp. Soysalli,) had been previously identified by [10] as potential species candidates. Our bPTP analysis strongly supported these lineages with pp values between 0.81 and 1. Other highly differentiated lineages found in our study, namely Cobitis sp. Marmara-Black Sea and Cobitis sp. Manyas, were not analysed by [10]. These two lineages were moderately to strongly supported with pp values of 0.23 and 0.85, respectively (Fig 2).
We also found some discrepancies between our bPTP results and some of the previously described species. Specifically, some species considered as Central European (C. elongatoides, C. satunini and C. tanaitica), Adriatic (C. ohridana) or Greek (C. arachthosensis, C. hellenica, C. stephanidisi and C. vardarensis) were not well supported as monophyletic lineages (<0.36 pp) or grouped together in lineages according to geographic location (Fig 2 and Table 1).

Discussion
Our mitochondrial and combined results based on the analysis of~89% of the recognized Western Palearctic Cobitis species identified seven highly distinctive clades, Clades 1-7, which, in some cases, were corroborated by the nuclear data. These major clades, and the species comprising them, are largely allopatric and correspond to groups recovered in previous studies with restricted sampling [16], [27], [29]. Our attempt to use molecular data to assess the number of Cobitis species in the Western Palearctic provides evidence for some incongruence  Table 1. List of Cobitis species analysed. Support values (BT/pp) are also provided for the different analyses: bPTP, ML and BI using the complete cyt b, complete RAG1 or combined datasets. NO refers to species not supported in the bPTP analysis or COI barcode study [10]. COI refers to values found by [10]; NO/NO refers to species not supported in ML/BI analyses, respectively.-indicate species not included in the analyses; 1 indiv = 1 specimen analysed. Lamina circularis (scale of Canestrini), a thickening of the pectoral rays on Cobitis males, refers to its absence (0), presence (1), or duplication (2). between a priori morphologically defined species and genetically identified lineages. The Western Palearctic loaches exhibit a complex phylogenetic pattern with multiple divergent lineages that not only supports the existence of the six candidate species suggested by [10] on the basis of COI sequences but also two additional candidate species, which were detected due to the extended sampling in Turkey conducted in this study.

Species delimitation
As described above, we used PTP-based methods to identify putative Cobitis species. However, to fully discuss whether identified lineages qualify as strong species candidates, other criteria must also be considered. To do this, we applied criteria successfully used in other freshwater fishes, such as cyprinids or siluriforms [39][40], including i) prior established taxonomy, ii) K2P distances, iii) phylogenetic congruence between mitochondrial and nuclear genes and iv) geographic congruence and range disjunction. Three lineages with high support in the bPTP analysis and moderate to high K2P distances for cyt b (> 2-3%) were highlighted as potentially new taxa: Cobitis sp. Gölbasi, Cobitis sp. Sapanca, and Cobitis sp. Soysallı. However, each of these lineages consist of a single individual. As they have already been identified as candidate species by [10], we will not discuss them further here. Although we observed a broad range of molecular divergence values, all analyses showed higher mean K2P cyt b divergences among major Cobitis clades than within species. Based on the 49 Cobitis lineages identified by bPTP, we found a mean K2P value of 0.59% for cyt b intraspecific divergence, similar to the value found for COI divergence (0.82%) [10] (Table 1). In our conservative interpretation, which does not consider lineages comprised of single individuals as new species, Clades 2, 3, 6 and 7 each contain at least one potentially undescribed species: Cobitis sp. Manyas, Cobitis sp. Marmara-Black Sea, Cobitis sp. Menderes, Cobitis sp. Ceyhan and Cobitis sp. Aksu-Eğirdir. Of these, Cobitis sp. Aksu-Eğirdir, Cobitis sp. Ceyhan and Cobitis sp. Menderes were moderately to well supported in the bPTP analysis, fairly well supported in the mitochondrial and nuclear phylogenies and showed K2P cyt-b distances between 3 and 5%. Cobitis sp. Aksu-Eğirdir is restricted to the Aksu river drainage, a current outflow of Lake Eğirdir in Southern Anatolia. Cobitis sp. Ceyhan and Cobitis sp. Menderes also inhabit particular drainages in restricted areas or lakes and fulfil all criteria to be recognized as new species. Cobitis sp. Manyas occurs in sympatry with C. puncticulata. Although K2P distances revealed the distinctiveness of Cobitis sp. Manyas, the nuclear phylogeny did not support it as a distinct lineage. Unlike the first three potentially new Cobitis species mentioned above, Cobitis sp. Manyas is not as geographically restricted and, therefore, does not satisfy many of the criteria to be considered a new species. Indeed, all available evidence suggests that Cobitis sp. Manyas represents a divergent population of C. puncticulata. Similarly, Cobitis sp. Marmara-Black Sea is widely distributed across drainages that terminate in the Marmara and Black seas. In spite of its mitochondrial differentiation, it was only weakly supported in the bPTP analyses and not supported at all in the nuclear phylogenies. Therefore, the Cobitis sp. Marmara-Black Sea lineage also does not satisfied all criteria to be considered a candidate new species. Our phylogenetic analyses suggest that this lineage is most closely related to C. taenia, a Central European species well known to consist of populations that easily hybridize with other Cobitis species [41,42]. A study of the phylogeography of C. taenia populations revealed an old lineage inhabiting drainages of the northern Black Sea that did not contribute to the recolonization of Europe [43]. However, our results highlight a new southern Ponto-Caspian lineage that may provide insights on European Cobitis recolonization. The high levels of genetic variation and the complex phylogenetic relationships found in the Turkish Cobitis species suggest the need for further phylogeographic evaluation.
In addition to the potentially new species, we also found some discrepancies between our criteria for species delimitation and a priori defined species (Table 1). In most cases, morphologically defined species satisfied some of the adopted criteria, with nuclear phylogenetic support being the most commonly violated criterion. Of the previously defined species, C. elongatoides and C. tanaitica lacked nuclear support and had low bPTP values (<50%). These species are well known for their tendency to hybridize and have not accumulated large genetic diversity [44][45]. Similarly, C. ohridana is also considered a potential species of the hybrid complexes formed by C.elongatoides and C. tanaitica [30,46]. In the case of C. vardarensis and C. stephanidisi, which recovered as sister species in the mitochondrial phylogeny, none of the other criteria support their delimitation as distinct species. However, in this case, a hybrid complex has not been described for this species. Indeed, C. stephanidisi has a highly restricted distribution: it is exclusively found in a single karstic spring in Greece (Kefakovriso, in Velestino village) [15]. Given its distribution, this species likely represents an offshoot of the more widely distributed C. vardarensis. The time elapsed since their split is likely not long enough to show monophyly with the nuclear genes. However, in some cases, such as with the Greek species C. hellenica and C. arachthosensis, none of the criteria were achieved (Figs 2-4 and Table 1). Therefore, the independent specific status of C. hellenica and C. arachthosensis is unsupported, in agreement with findings reported in previous studies [27,29]. The pigmentation pattern, especially in young specimens, represents a major difference between the two species [15]. The adjacent distribution of the species and their genetic similarity suggest a recent split and/or recent gene flow leading to unsorted lineages. Mediterranean peninsulas (Italian and Iberian, Clades 1 and 5) and Central European drainages (Clade 2) hold lower numbers of spined loaches when compared to the Balkan and Eastern Mediterranean drainages. The early isolation of the Iberian and Italian peninsulas and the lack of subsequent connection with the geographically closest ichthyofauna from Central Europe or North Africa most likely promote their high endemicity and low number of species observed in the peninsulas. Our data confirm that the Adriatic region and the Italian and Iberian peninsulas represent long-term persistence areas of spined loaches; however but contrary to what has been observed for other vertebrates [47], these areas do not represent refugia or sources of evolutionary assemblages that colonized Northern Europe at any period. The long-term isolation of the ichthyofauna inhabiting these areas favoured their genetic differentiation and promoted their high level of endemicity. In general, the major Cobitis clades, and the species comprising them, are allopatric with the exception of some of the Eastern Mediterranean species inhabiting Turkey and Greece. The observed distribution of species found in Turkish drainages indicate that these habitats are not continuous or linear as the distribution pattern of some current Cobitis species appear to correspond to ancient geological units. Drainages occupying geologically active areas, particularly boundary zones between plates such as the river Ceyhan [48], contain different taxa i.e. C. evreni and Cobitis sp. Ceyhan. Furthermore, the observed high endemicity of spined loaches of Clade 7 (7 species and 2 newly proposed species) indicate that the lacustrine system of Central Anatolia as a geologically complex area. Within this are, Lake Beysheir is one the places with the highest number of Cobitis species. The phylogenetic pattern found in our study supports the importance of geological history for the speciation of Cobitis in the Western Palearctic. Congruent spatial genetic subdivisions have been reported for other vertebrates, including water frogs and lizards, in the Eastern Mediterranean [49][50][51]. Our results, therefore, provide further evidence of a common pattern of co-distributed taxa within this region.

Major Cobitis clades in the Western Palearctic
Overall, our biodiversity assessment indicates that Turkey has the highest level of genetic diversity and endemicity of spined loaches in the Western Palearctic followed by the Balkan region. The remarkable richness observed for this group, which has also been observed for other Eastern Mediterranean groups of freshwater fishes, including cyprinids and aphaniids [13], [52], [53], and has highlighted Turkey as a present-day biodiversity hotspot.  Table. List of new specimens analysed in this study with corresponding collection information. Group refers to major Clades 1-7 obtained in this study. (XLS) S2 Table. List of specimens obtained from GenBank with corresponding collection information. Group refers to major Clades 1-7 obtained in this study. (XLS)