Ethics statement

This investigation was conducted entirely in accordance with ethical standards and Croatian legislation. The Ethical Committee of the Faculty of Science, University of Zagreb, approved the work. Permission for sampling was issued by competent authorities: for sampling on localities within Croatia (Jadova, Suvaja, Ričica, Krbava, Čikola, Cetina, Vinelić, Studenac, Sušik, Jasenak field, Šmit lake, Gradaščica, Bregana and Konavosko field) permit was issued by the Ministry of Nature Protection; for Greek localities (Pavllo, Acheloos, Louros, Pinios, Evinos, Kalamas, Alfios, Kifissos) by the Greek Ministry of Environment, Energy and Climate Change; for sampling on localities in Montenegro (Cijevna, Morača, Zeta) permission was issued by Ministry of Agriculture and Rural Development; Italian samples (localities Bevera, Staffora, Arno, Tiber, Volturno, Po) were obtained based on permissions from following authorities: ‘Amministrazione Provinciale di Cuneo–Servizio Caccia e Pesca’, ‘Citta Metopolitana di Torino–Servizio Tutela Fauna e Flora’, ‘Parco del Monviso’, ‘Ente di gestione delle Aree Protette del Po e della collina Torinese’, ‘Ente Tutela Pesca del Friuli Venezia Giulia’, ‘Administrazione Provinciale di Arezzo–Servizio Caccia e Pesca’; on French localities (Argens, Herault, Lergue and Var) permission from ‘Conseil Superieur de la Peche’ was obtained; and for localities within Bosnia and Herzegovina permit was issued by ‘Republički zavod za zaštitu kulturno-istorijskog i prirodnog nasljeđa, Banja Luka’.

We have collected samples from 40 localities covering the whole geographic range of the genus Telestes (Fig 1, Table 1) and representing all 14 recognized species. For species presented with more than one population, namely those distributed on more than one locality that are not connected with permanent ground connections, samples were obtained from more than one locality. Sampling was conducted by electrofishing. Small piece of fin tissue was taken from each individual and deposited in ethanol or frozen for further analyses. Two molecular markers were employed in this investigation: the mitochondrial gene for cytochrome b (cyt b) and the nuclear RAG1 gene with altogether 2490 base pairs (bp) analyzed. The cyt b gene is considered as the most useful marker in recovering phylogenetic relationships among closely related taxa [10], [11] and its suitability for taxonomic and population genetic studies was proven in several previous investigations on various vertebrates (e.g. [10], [12], [13]). RAG1 gene, on the other hand, due to its lower mutation rate, is adequate for phylogenetic investigation on higher level, namely discovering relationships among more distantly related species. Inclusion of both, maternally inherited and biparentally inherited markers shall enable investigation of phylogenetic structure and evolutionary history of this genus, but also point out to eventual hybridization.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Sampling localities and number of specimens belonging to certain species collected at each locality.

https://doi.org/10.1371/journal.pone.0187366.t001

Total genomic DNA was extracted from fin tissue samples using a standard extraction product (DNeasy tissue kit, Qiagen). Polymerase chain reaction (PCR) amplifications were performed using HotStarTaq Master Mix (Qiagen) and primers GluF and ThrR [14] for cyt b gene, and RAG1F1 and RAG1R1 [15] for RAG1. Amplification of cyt b gene was conducted through 35 PCR cycles with 45 s of denaturation at 92°C, 90 s of annealing at 48°C and 1 min and 45 s of extension at 72°C. PCR protocol for RAG1 gene comprised 5 cycles with 40 s at 94°C, 1 min at 60°C, 2 min at 72°C, followed by 35 cycles with 30 s at 95°C, 1 min at 56°C and 2 min at 72°C. Sequencing was carried out by Macrogen Service Centre (Amsterdam, Nederland) with internal primer pairs CB4-GLU (5’CCTGAAAYATYGGYGTRGT3’) and PHOX-THR (5’AGGAGGAARTGRAATGCGAA3’) (Doadrio and Perea, personal communication) for cyt b, and RAG3F (5’GGGTGTGTCAGYGAGAAGCA3’) (Quenouille et al. 2004 in [2]) and RAG6R (5’ATGGCTTTCCGCTCTGCTAC3’) (Doadrio and Perea, personal communication) for RAG1. We have obtained cyt b sequences from 319 samples and RAG1 sequences from 152 samples.

Homologous regions of both genes were aligned manually against previously published sequences. Chromatograms and alignments were checked visually and were found to contain no gaps or stop codons. In order to determine nuclear haplotypes in heterozygous individuals, Bayesian statistical method implemented in PHASE 2.1 software [16], [17] was conducted. The analysis was run five times with different values of the seed of the pseudo-random number generator. Each run consisted of a burn-in-period of 100 followed by 1000 iterations. In order to test whether all mutations were selectively neutral, statistical tests D* and F* [18] and Tajima’s test [19] were conducted using DnaSP v5 [20]. The same software was employed to estimate the recombination parameter, R [21], and the minimum number of recombination events, RM [22], for nuclear data set.

Phylogenetic reconstruction was conducted with the aim of revealing the phylogenetic structure inside the genus Telestes, i.e. confirming the position of already investigated Telestes populations and revealing phylogenetic relationships for populations analysed for the first time. It was based on two methods of phylogenetic inference: maximum parsimony (MP), implemented in PAUP (version 4.0b10; [23]), and Bayesian inference (BAY), implemented in MrBayes (version 3.1.2; [24]). For MP analysis, a heuristic search mode was used, with randomized input orders of taxa, and TBR branch swapping with all codon sites and nucleotide substitutions types weighted equally. Nonparametric bootstrapping (1000 pseudo-replicates, 10 additional sequence replicates) was used to assess branch support (BS). Each BAY analysis consisted of two simultaneous runs. For each, Markov Chain Monte Carlo was run four times for three million generations with trees sampled every 100 generations. The first 20% of the sampled trees were discarded and Bayesian posterior probabilities (BPP) were estimated from the 50% majority-rule consensus tree of the retained trees. For rooting of cyt b phylogenetic tree we used a sequence of Squalius cephalus, more distant representative of the same subfamiliy [2], and of Cyprinus carpio, that belongs to the same family, but to a different subfamily. A sequence of Squalius squalus was used for rooting the RAG1 phylogeny. As the most suitable model of sequence evolution, TN93 model was identified using MEGA6 software and used for BAY, as well as *BEAST analyses (see later).

Additional method, median-joining approach (MJ) using Network 4.5.1.6. software (Fluxus Technology Ltd.) was employed on RAG1 data set. In a resulting phylogenetic network it is possible to notice horizontal gene transfers, so this approach is especially adequate for phylogenetic reconstruction based on nuclear genes.

Divergence times between investigated species, as well as between phylogenetic lineages, were estimated by a Bayesian MCMC coalescent method, using BEAST 2.4.7software [25]. The analysis was conducted on cyt b data set and on the concatenated data set (cyt b and RAG1). Concatenation was achieved using Mesquite 2 software [26]. The rate homogeneity across phylogenetic lineages was assessed by the log-likelihood ratio test (LRT), comparing the likelihood of phylogenetic trees (reconstructed using maximum likelihood approach) with and without molecular clock enforcement in PAUP*. Since the likelihood scores were the same in both cases for cyt b data set, we applied a strict molecular clock. Since we aimed to reconstruct species tree and gene tree at the same time, while preparing an input file for analysis we used StarBeast template in BEAUti software. This analysis is based on a multispecies coalescent model, which allows inclusion of an algorithm of coalescence in the speciation process. The number of MCMC steps (the length of chain) was twenty millions. The molecular clock calibration was conducted based on the divergence rate of cyt b gene of 0.4% per lineage per million years [2]. For concatenated data set we employed relaxed molecular clock and divergence rate of 0.2% per lineage per MYA, based on [2].

Ancestral distribution ranges of lineages were reconstructed using Statistical Dispersal-Vicariance Analysis implemented in S-DIVA software [27]. This method reconstructs the ancestral distribution in a phylogeny by optimizing a three-dimensional cost matrix, in which extinctions and dispersals ‘cost’ more than vicariance, and it determines statistical support for ancestral range reconstruction [27]. Reconstruction of ancestral geographic ranges was conducted based on the concatenated data set in order to investigate bi-parental origins. An input set of trees was obtained by BAY analysis. Altogether six recent geographic ranges for Telestes were denoted (marked in the Fig 1).

Pairwise comparisons of uncorrected sequence divergence (p-distances) of both genes were analysed using MEGA version 3.1 [28].

The level of intraspecific genetic diversity was estimated by calculating several measures of DNA polymorphism for each genetic marker, employing DnaSP v.5: number of polymorphic sites (S), haplotype diversity (Hd), average number of nucleotide differences (K) and nucleotide diversity (π). Furthermore, the frequency of each haplotype was calculated as percentage of a certain haplotype in a population. In the intraspecific genetic diversity estimation we did not include sequences of introgressed individuals, in order not to overestimate measures; all measures are based on sequences belonging to a single species, as recognized by current taxonomy.

Phylogenetic structure of the genus Telestes

Even though species of the genus Telestes have been included in numerous previous phylogenetic investigations (see Introduction), this is the first attempt to achieve a comprehensive depiction of its present phylogenetic structure, as well as evolutionary mechanisms that shaped it. Similar as in previous investigations, the basal lineage inside Telestes cannot be determined without doubt. However, the majority of the phylogenetic analyses denoted two lineages as the first ones to diverge; one comprising T. souffia (revealed by analyses based on RAG1 sequences) and the other containing T. fontinalis and T. croaticus (proposed by the majority of reconstructions of the cyt b data). We speculate that the reason for the discrepancy is fast and old separation event of both mentioned lineages. A long period of isolation of these ancient lineages led to their high variation compared to the remaining lineages. Further, a majority of the analyses corroborated an old origin and the distinctiveness of southern Balkan lineage containing T. pleurobipunctatus and T. beoticus. The long period of independent evolution is also corroborated by morphological features of the species belonging to the oldest lineages: T. croaticus and T. fontinalis can be distinguished from the remaining Telestes species by gleaming coloration (personal observation), while T. pleurobipunctatus can be differentiated from the remaining species by having pharyngeal teeth in one row only [29], [3].

The remaining phylogenetic lineages form a common cluster, in which the phylogenetic relationships do not always follow a geographic pattern, e.g. T. croaticus and T. fontinalis are very distant from the remaining species distributed in northern Dinarids; T. montenigrinus seems to be related to T. muticellus, even though their recent distribution ranges are apart. This investigation, furthermore, rejected the presumption of a closer relatedness of T. souffia, T. muticellus and T. pleurobipunctatus, which were sometimes considered to form a single complex [3].

Recent relationships and intraspecific structuring was better resolved with mitochondrial marker, due to its faster mutation rate, whereas nuclear phylogeny seems to be adequate for describing speciation. Significant intraspecific structuring and the existence of more than one lineage are present in T. souffia, T. muticellus, T. croaticus and T. pleurobipunctatus. Divergence of two T. souffia lineages is of Pliocene origin and is concordant with findings of previous authors [30], [9] with one lineage distributed in France (rivers Argens, Hérault, Lergue and Var as the type locality) and a second lineage inhabiting localities in the Danube basin, as well as the Soča River. Two T. souffia lineages were denoted also by [5], whereas [4] recognized four distinct lineages inside T. souffia, one of which is concordant with our second lineage. Previous reports regarding T. muticellus were contradictory (e.g. [8] described low genetic diversity of that species, whereas [6] presented five distinct haplotype groups), and sometimes not reliable due to the inclusion of T. souffia samples (e.g. in the investigation of [6]). Based on our results, T. muticellus comprises at least three highly distinct lineages: I. Volturno population, whose distinction was already observed by [31], II. populations from the northern part of T. muticellus distribution range, III. Arno and Tiber populations. Two lineages of T. croaticus diverged in the late Pliocene or early Pleistocene. However, haplotypes belonging to both lineages are present in the Jadova R., implying their secondary contact. Exceptionally pronounced structuring is also present inside T. pleurobipunctatus.

In order to investigate more detail intra- and interpopulational, as well as within lineages structure and gene flow, further investigations are advisable that should focus on a certain lineage, but include more genetic markers.

Evolutionary history scenario

The proposed evolutionary scenario (Fig 7) relates to estimates based on both cyt b gene and nuclear data set, that resulted in similar estimation of timing of divergence events and known data on the historical biogeography of Europe ([32], [33], [34], [35]). Reconstruction of ancestral ranges, even though not reliable for older events, is also concordant with the proposed scenario, as explained below.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Colonization routes of the Telestes lineages, based on the proposed evolutionary scenario.

*—Telestes ancestor in 7a and 7b, ancestor of T.croaticus/T. fontinalis in 7c and 7d; #—ancestor of T. souffia; x–ancestor of the species cluster (T. metohiensis/T. dabar/T. miloradi/T. muticellus/T. montenigrinus/T. ukliva/T. turskyi/T. karsticus/T. polylepis); o–ancestor of T. pleurobipunctatus/T. beoticus.

https://doi.org/10.1371/journal.pone.0187366.g007

In the first stage of the Miocene, southern Europe was an archipelago, with an ancestor of Telestes probably persisting on one of the islands (Fig 7A), which it might have reached from the European mainland or Anatolia during earlier landmass connections. However, by the end of the Burdigalian (16 MYA), some of the islands became joint into a single landmass, forming a continental bridge between western Europe and Anatolia [33] and providing an opportunity for colonization of eastern Balkan region (concordant with the divergence of T. pleurobipunctatus/T. beoticus lineage) (Fig 7B). In the same period the northern part of the Balkan was connected with central Europe and the Dinarids were not formed yet, explaining the connection of T. souffia with the Mediterranean species. Phylogenetic reconstructions and divergence time estimations imply the hypothesis that already in the Burdigalian, when connections between central Europe, Balkan and Anatolia were possible [33], the ancestor of T. croaticus and T. fontinalis was isolated. In the Langhian stage (15.9–13.8 MYA), following the Middle Miocene transgression, the majority of central and southern Balkan became separated from the European landmass, with several land fragments also persisting in the Parathetys [33]. The ancestor of T. croaticus and T. fontinalis did not come into a secondary contact with the remaining Telestes lineages, leading to the clear distinctiveness between those two species and the remaining ones, even though some of them (e.g. T. polylepis and T. karsticus) are presently located in their immediate geographic proximity. By the Langhian stage (Fig 7C) T. souffia was also isolated from the remaining lineages, and its further evolutionary course was determined by the biogeography of the Black Sea watershed and central Europe. Nevertheless, on the already mentioned southern island the diversification of Telestes continued and consequently led to the high species richness. During the Langhian phase the ancestor of Herzegovinian and southern Dalmatian T. metohiensis, T. dabar and T. miloradi became isolated from the remaining species belonging to already mentioned species cluster. The divergence of T. muticellus, distributed on the Apennine peninsula can be dated to the Serravallian stage (12.6–7.8 MYA, Fig 7D). The colonization of the Apennine peninsula might be provoked by the closing of the Slovenian corridor. Furthermore, communications between the Black Sea and the Adriatic watershed also existed in that period, at least until the early Tortonian (11.6 MYA) when lineages of T. karsticus/T. polylepis and T. turskyi/T. ukliva diverged. Separation of species inside the mentioned lineages occurred in the Messinian stage (7.2–5.3 MYA), and is probably connected with Dinaric Lake System [34], [35] and/or intense tectonic activity. The Pliocene epoch had a final impact on Telestes diversity, because divergence of several species (T. metohiensis, T. miloradi) or lineages within T. muticellus, T. souffia and T. pleurobipunctatus can be dated back to that time.

Our estimates for early divergences inside this genus are older than proposed by previous authors [7], [1], [8], [4], which is probably a result of molecular clock calibration and smaller sample sizes included in earlier investigations. In disagreement to previous investigators [1], [4], our results do not point out the Messinian Salinity Crisis (MSC) as critical point for the divergence of Mediterranean Telestes species, supporting findings of [2] for the Leuciscinae, but also for Adriatic spined loaches (genus Cobitis; [36]). It is also clear that the separation of T. souffia and T. muticellus did not happen simultaneously, as previously proposed [5], but evolutionary history of Telestes can rather be described as containing sequences of gradually occurring events intermittent with longer periods that did not leave a trace in recent genetic structure.

Even though widely observed inside Leuciscinae, hybridization did not play an important role in the evolutionary history of Telestes, although it was sporadically noticed. Previous investigations [37], [5] revealed recent gene flow between T. souffia and T. muticellus, which was explained as hybridization by [4]. Since we have investigated both nuclear and mitochondrial DNA, we are able to conclude that interspecific connections include both nuclear hybridization and introgression of mtDNA, but are very restricted based on our sampling.

Contrary to the opinion of [1] that relationships of Telestes species mostly follow the geographic pattern, it seems that the complex geological history of southern Europe led to complicated biogeographical pattern with genetically distinct species distributed close to each other, most likely a consequence of multiple colonization and crossing events between rivers and water basins. For example, in the Adriatic watershed in Croatia and Bosnia and Herzegovina species belonging to three lineages are present, implying three colonization events of this area by different lineages and in different time periods, explaining the high diversity observed there.

Molecular diversity of Telestes species

The distribution range of Telestes comprises the Mediterranean and central European area. Inside its range there are smaller regions with exceptionally high Telestes diversity both, at the species and genetic levels. The Eastern Adriatic coast is the species’ ‘hot spot’, but also contains high levels of genetic diversity, especially in the region of central Dalmatia (Krka and Cetina River drainages). Even though the number of recognized species is smaller in Greece, it contains the highest genetic diversity of Telestes at all.

Elevated levels of genetic polymorphisms based on the mitochondrial marker (Hd>0.95) were found in T. pleurobipunctatus, T. muticellus and T. croaticus, three species with pronounced intraspecific structuring; but also inside T. turskyi and T. ukliva (Table 4). The latter two species have smaller distribution ranges and their high genetic diversity is a probable consequence of their old origin and unconstrained evolutionary history, due to the absence of glaciations in their distribution ranges that might have provoked bottleneck events. The main difference in the genetic polymorphism measures between these two groups of species is in the values of nucleotide diversity and the average number of nucleotide differences that are significantly higher in well-structured species (34.9 in T. pleurobipunctatus, 21.9 in T. muticellus comparing to maximum of 5.8 in species with homogenous structure). Likewise, T. souffia, another species comprising more than one genetically distinct unit, also contains high average number of nucleotide differences and high nucleotide diversity, although its haplotype diversity is smaller. The nuclear RAG1 gene also pointed out T. croaticus, T. muticellus, T. souffia, T. ukliva, T. fontinalis and T. pleurobipunctatus as genetically highly diverse. Very high values of nucleotide diversities and average number of nucleotide differences contained inside T. pleurobipunctatus, T. muticellus and T. souffia, are not characteristic for a single species and imply the possibility that they represent cases of species complexes, which is concordant also with other results and some previous reports. Similarly as obtained in this investigation for cyt b gene, inadvertent inclusion of cryptic species in samples presumably belonging to a single species has been considered as one of the most powerful causes of bias in haplotype and nucleotide diversity estimations based on mitochondrial cytochrome c oxidase sequence data [38].