1. Introduction
The freshwater systems of the Dinaric karst region maintain exceptionally rich fish communities with a high portion of endemic species [
1]. Several events and phenomena have been proposed as possible explanation for such diversity [
2,
3,
4,
5]: geographic isolation of rivers and their underground connections, geological history, ice-age refugia and specific mutation rates. Allopatric speciation, as a consequence of long-term isolation of water basins, has been considered the most important divergence mode in southern Europe (e.g., [
6]), including the Dinaric karst region [
5]. Evidence of parapatric speciation is very rare (e.g., [
6]) and disputed, while hybridization events, although locally recorded (e.g., [
7]) were not considered as important for the speciation of freshwater fishes in this area. Nevertheless, for two cyprinid species from the Iberian Peninsula,
Iberocypris alburnoides (Steindachner, 1866) and
I. palaciosi Doadrio, 1980, hybridization profoundly affected their evolutionary histories [
8]. Hybridizations and mixed genetic material in a single population, originating from two species, was also found in the
Squalius svallize Heckel & Kner, 1858 population from southern Croatia [
9]. Furthermore, several cases of mitochondrial introgressions have been reported among cyprinid fishes in southern Europe (e.g., [
9,
10,
11]). However, data on the evolutionary history and recent genetic composition of species distributed in the Dinaric freshwaters remain scarce.
The genus
Delminichthys Freyhof, Lieckfeldt, Bogutskaya, Pitra & Ludwig, 2006 is endemic to the freshwater systems of the Dinaric karst in Croatia and Bosnia and Herzegovina. Its closest relative, the genus
Pelasgus Kottelat & Freyhof, 2007 [
3], is distributed in the southern Balkans, primarily in Greece with one species in Lake Ohrid [
8]. Unlike most genera within the Leuciscidae family,
Delminichthys inhabits a very restricted geographic area [
12]. It consists of only four endemic species:
D. krbavensis (Zupančič & Bogutskaya, 2002), restricted to the Krbavsko karstic field in Croatia;
D. jadovensis (Zupančič & Bogutskaya, 2002), inhabiting the small Jadova River in Croatia;
D. adspersus (Heckel, 1843), found in the Imotski Field, Matica River (in the Jezero Field), Norin River, Kuti and Baćinska Lakes in Croatia, as well as in the Trebižat River and Krenica Lake in Bosnia and Herzegovina; and
D. ghetaldii (Steindachner, 1882), found in the Popovo, Ljubomirsko, Dabarsko and Fatničko karstic fields in Bosnia and Herzegovina, and in the Ombla River in Croatia. All four
Delminichthys species are endemic species with very restricted distribution ranges and all four species are considered endangered.
Delminichthys jadovensis and
D. krbavensis are considered critically endangered (CR), based on the IUCN Red List, whereas
D. adspersus and
D. ghetaldii are listed as vulnerable species (VU).
Though
Delminichthys represents an ideal model system for investigating the evolution in the karstic river systems of the Dinarides, it has not yet been comprehensively studied. Ref. [
3] suggested that the
Delminichthys/
Pelasgus lineage is one of the basal lineages inside Leuciscidae, and its origin was estimated to around 14 MYA (million years ago). Ref. [
12] also recognized
Delminichthys as the remnant of the first great Leuciscidae radiation. Ref. [
4] investigated the population genetic structure of
D. adspersus and found evidence of its underground migrations. The lack of the stygobiontic characters was explained by the limited time this species annually spends in the underground environment. On the other hand, ref. [
3] reported the presence of reductive characters in both
Delminichthys and
Pelasgus.
With an aim of contributing to the understanding of evolutionary peculiarities in Dinaric water systems, this study examines the population genetic structure and demographic history of each
Delminichthys species. Considering that recent genetic structure and polymorphism bear witness to the past geology and evolutionary history of a certain species, we wanted to look for evolutionary responses and adaptations to specific conditions of karstic freshwaters. Another goal was to further investigate underground migrations in the genus
Delminichthys. Namely, if adaptations to underground life are characteristic to genus, this should be visible in the intraspecific structure of
D. ghetaldii (another species occurring at more than one locality), and perhaps also in the genetic polymorphism of
D. jadovensis and
D. krbavensis. The two latter species are located in an area that was affected by glaciations led to bottleneck effects in previously investigated freshwater fishes [
13].
2. Materials and Methods
Samples of all four
Delminichthys species were collected (, ) and used for phylogenetic and population genetic analyses. A total of five molecular markers were analysed: mitochondrial cytochrome
b (cyt
b) and control region (CR), and nuclear recombination activating gene 1 (RAG1), gene for interphotoreceptor retinoid binding protein (IRBP) and gene for beta-actin (BA). All but CR are protein coding genes and CR is a non-coding DNA with function in RNA and DNA synthesis.
Total genomic DNA was extracted from fresh and deep-frozen fin tissue using a standard extraction product (DNeasy tissue kit, Qiagen, Hilden, Germany). Polymerase chain reaction (PCR) amplifications were performed in a 50 µL reaction volume containing 25 µL AmpliTaq Gold Master Mix (Applied Biosystems, Waltham, MA, USA), 2 µL of each primer and 4 µL of template DNA. Amplification protocols and primers, as well as sequencing primers are presented in the . Sequencing was carried out by Macrogen Europe Service Centre (Amsterdam, The Netherlands).
. Map of the study area indicating the distribution ranges of <i>Delminichthys</i> species and sampling localities. J—Jadova River, K—Krbavsko Field, D—Vrljika River, I—Prološko blato Lake (in the Imotsko Field), C—Sija, A—Rastočko Field, M—Matica River (in the Jezero Field), B—Baćinska Lakes, E—Popovo Field, F—Bilećko Lake, O—Ombla R. Abbreviations: CRO—Croatia, BIH—Bosnia and Herzeovina, ME—Montenegro; RS—Republic of Serbia, ITA—Italy.
. Number of samples from each locality used for analyses, as well as haplotype codes and GenBank accession numbers. Country codes: CRO—Croatia, BIH—Bosnia and Herzegovina.
. PCR protocols, amplification and sequencing primers used to amplify and sequence for investigated genetic markers.
Homologous regions of the investigated genes were aligned manually. Chromatograms and alignments were checked visually and were found to contain no gaps or stop codons. Haplotype variants of nuclear genes in heterozygous individuals were reconstructed with the Bayesian statistical method implemented in PHASE 2.1 software [
16,
17]. To test whether all mutations were selectively neutral, statistical tests D
∗ and F
∗ (proposed by [
18]) and that of Tajima [
19] were conducted using DnaSP v5 [
20].
Phylogenetic relationships of the genus
Delminichthys were investigated using two methods of phylogenetic reconstruction: maximum parsimony (MP) and maximum likelihood (ML) implemented in PAUP (version 4.0b10; [
21]). For MP analyses, a heuristic search mode with 100 replicates was used, with randomized input orders of taxa, and tree bisection-reconnection (TBR) branch-swapping with all codon sites and nucleotide substitutions types weighted equally. ML analyses were performed under the heuristic search option using the TBR branch-swapping algorithm. Branch support (BS) was assessed by nonparametric bootstrapping (1000 pseudoreplicates, ten additional sequence replicates). Each gene was analysed separately, due to differences in their phylogenetic performance (). Sequences of species belonging to other Leuciscidae lineages (as defined by [
3]) were also included in the phylogenetic reconstruction with the aim of describing the position of the genus
Delminichthys within the subfamily. List of included sequences with their GenBank accession numbers and references is presented in Supporting information (Table S1). Sequences of
Cyprinus carpio were used as outgroups. Namely, since in the analyses we have also included representatives of other Leuciscidae species, the rooting could be accomplished by a species that is outside Leuciscidae family, but inside Cypriniformes and for which sequences of all investigated genetic markers are available, so
Cyprinus carpio was chosen as an outgroup.
. Length and phylogenetic performance of genetic markers.
Estimation of divergence times was conducted with the Bayesian MCMC coalescent method, using BEAST 1.7.0 software [
22]. The rate homogeneity across phylogenetic lineages was assessed by the log-likelihood ratio test. The likelihood of phylogenetic trees (reconstructed by maximum likelihood approach using PAUP) with and without molecular clock enforcement were compared. Since they were the same in both cases, a strict molecular clock was applied. Branch rates were drawn from an uncorrelated log-normal distribution and a Yule speciation prior with random starting tree. The applied substitution model was HKY with Gamma site heterogeneity model. We used default prior distributions for kappa, frequencies and alpha, whereas substitution rate parameters were unlinked across codon positions. The number of MCMC steps (chain length) was three million. Molecular clock calibration was conducted based on the divergence rate of cyt
b gene in Leuciscidae, of 0.4% per lineage per million years [
3]. Since nuclear genes did not enable resolution of
Delminichthys species, they were not used for divergence time estimation or for ancestral ranges reconstruction. Likewise, due to the lack of adequate data set for investigation of evolutionary history based on CR, this genetic marker was not included in the divergence times estimations.
In order to reconstruct ancestral geographic distributions, we performed Statistical dispersal-vicariance analysis implemented in the S-DIVA program [
23]. 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 determines statistical support for ancestral range reconstruction [
23]. A total of 11 geographic localities inhabited by
Delminichthys were denoted (concordant with the sampling sites shown in ). A set of trees used for ancestral distributions reconstruction was obtained by BI analysis of cyt
b haplotypes. This, thereafter, resembles the maternal ancestral ranges, while further study is needed on more nuclear markers to describe ancestral ranges of paternal populations. Similarly as with divergence times estimations, this analysis was based on cyt
b data set.
Intrapopulational and intraspecific genetic diversity was described using several measures of genetic polymorphism calculated using the DNAsp software [
21]: number of haplotypes (h), number of polymorphic sites (S), haplotype diversity (Hd), average number of nucleotide differences (k), and nucleotide diversity (π).
In order to detect possible gene flow among populations of the same species (
D. adspersus and
D. ghetaldii, because those two species comprise more than one population) and to test the hypothesis of underground migrations, we examined the interactions among populations using maximum likelihood approach [
24,
25] implemented in MIGRATE 3.2.1 [
26]. The computer programme MIGRATE calculates the profile likelihoods of the population parameters using genetic data. We determined the immigration rates between populations, as well as the number of migrants in one generation. Immigration rates (more specifically—mutation scaled immigration rates) are defined as a measure of how much more important immigrations are over mutations in bringing new variants into populations. In MIGRATE, this parameter is considered as a long-time average over the genealogy of the individuals in the sample [
26]. Immigration rates were calculated both as mutation-scaled effective immigration rates and as the number of immigrants per generation. Due to low resolution of nuclear genes and their inability to clearly distinguish
Delminichthys species, migration rates were calculated based only on the mitochondrial genetic markers (cyt
b and CR).
Changes in past effective population sizes were analysed using Bayesian skyline plots (BSP), implemented in Beast 1.7.0. The BSP model was employed on the cyt
b data set and it estimates the history of changes in effective population size from the variability among sampled haplotypes, assuming a mutation rate of 0.4%/MY. The BSP settings were the same as in the divergence time estimations with the tree prior set to Coalescent: Bayesian Skyline.
3. Results
A total of 84 cyt
b sequences of a length of 1140 base pairs (bp), 71 CR sequences of a length of 320 bp, 120 RAG1 sequences (from 60 individuals) of a length of 1300 bp, 82 IRBP sequences (from 41 individuals) of a length of 849 bp and 118 BA sequences (from 59 individuals) of a length of a 460 bp were obtained. No stop codons, insertions or deletions were recorded; sequences of the same gene were of the same length in all
Delminichthys species. Number of sequences was determined by the availability of samples in the field sampling campaigns.
Of the 60 samples sequenced for RAG1, 20 (33%) were heterozygous. The number of heterozygous individuals for IRBP gene was 23 (56%), whereas in the BA data set there were 30 heterozygous individuals (51%). The samples provided 43 unique cyt
b haplotypes, 18 CR haplotypes, 25 RAG1 haplotypes, 23 IRBP haplotypes and 21 BA haplotypes. Neutrality tests conducted on all data sets suggested that generally all species are in mutation-drift equilibrium. Only exception was
D. adspersus for which neutrality tests implied some deviation from mutation-drift equilibrium since calculated statistics (Fu and Li’s D
∗ and F
∗ statistics, as well as Tajima’s D) for its cyt
b and BA data sets turned out to be statistically significant (
p < 0.05), whereas statistics for CR and IRBP genetic markers implied mutation-drift equilibrium. The number and proportion of variable and parsimony informative characters was much higher in mitochondrial genetic markers than in nuclear genes, with the exception of BA (). Even though investigated BA data set provides high amount of variable and parsimony informative characters (comparable with mitochondrial markers), those characters are mostly present in other species included in the data set and not in
Delminichthys, yielding phylogenies in which positions of
Delminichthys haplotypes are not well resolved.
The phylogenetic reconstruction of cyt
b sequences revealed clear resolution of the four investigated species () and corroborated the monophyletic origin of all
Delminichthys species, which are most closely related to
Pelasgus prespensis (Karaman, 1924). Based on MP phylogenetic reconstruction, two sublineages can be observed within the
Delminichthys lineage, each comprised of two geographically more closely located species (where
D. adspersus is a sister species to
D. ghetaldii, while
D. jadovensis and
D. krbavensis form a second species pair). Strong support for the separation of two species pairs is suggested by the 100% bootstrap value. ML analysis, however, implied an earlier separation of
D. ghetaldii and its distinct position. There is, however, no obvious intraspecific structuring that could be related to different geographic localities or pertinence to separate populations. Similar situation was observed in CR phylogenies (), but with the exception of
D. ghetaldii position. Three
Delminichthys species,
D. adspersus,
D. jadovensis and
D. krbavensis form separate lineages and are monophyletic, with
D. jadovensis and
D. krbavensis forming sister pair. Position of
D. ghetaldii was not resolved well by CR.
Delminichthys ghetaldii sequences form soft polytomy in the CR phylogenies, event though all of them clearly belong to
Delminichthys genus ().
The RAG1 gene tree, on the other hand, implied a polyphyletic origin of
Delminichthys (). Samples from all localities clustered together, with the exception of the Ombla River, as a sister lineage to
P. prespensis. Haplotypes from the Ombla River (
D. ghetaldii) were included into different lineage and appear to be more closely related to other Leuciscidae, in particular to
Telestes pleurobipunctatus (Stephanidis, 1939), than to the remaining
Delminichthys, including other
D. ghetaldii populations. It is important to mention that 75% of the Ombla samples are heterozygous for RAG1, but both alleles of all samples clustered with
T. pleurobipunctatus and not with the remaining
Delminichthys samples. Sequences of the second investigated
D. ghetaldii population (from Bilećko Lake) clustered with the remaining
Delminichthys species and no structuring or species differentiation can be noted within that lineage. Moreover, the same haplotypes were found in samples belonging to different species.
The remaining nuclear genes (phylogenetic trees presented in the Supporting information, Figure S1) did not enable clear resolution of the position and relationships of four
Delminichthys species inside the Leuciscidae phylogenetic tree. Yet they resulted in soft polytomies of
Delminichthys haplotypes. Moreover, some of the haplotypes are present in more than a single species (). Nevertheless, these markers clearly pointed out monophyletic position of all
Delminichthys species and no traces of hybridogenetic origin or introduction of DNA belonging to other species could be noticed from phylogenies based on IRBP and BA. Interestingly, in IRBP and BA data sets sharing of haplotypes among various species was noticed. In both data sets, a haplotype with the highest frequency (IrbpJAD1 and BAJAD1) was found in the Jadova River and in both
D. ghetaldii populations (Bilećko Lake and Ombla River). BA most widespread haplotype (BAJAD1) was found besides in
D. jadovensis and
D. ghetaldii, also in
D. adspersus.
Due to the observed situation, intraspecific structures, divergence time and ancestral ranges estimations were calculated based only on the mitochondrial data sets.
. MP phylogenetic tree of <i>Delminichthys</i> cyt <i>b</i> sequences. Numbers at nodes represent MP and ML branch support. Asterisk denotes full branch support (100).
. MP phylogenetic tree of <i>Delminichthys</i> CR sequences. Numbers at nodes represent MP and ML branch support. Astersik denotes full branch support (100).
. MP phylogenetic tree of <i>Delminichthys</i> RAG1 sequences. Numbers at nodes represent MP and ML branch support.
The origin of the genus
Delminichthys, i.e., its separation from
Pelasgus, probably already occurred in the Upper Miocene (). The intrageneric divergences, however, are much younger, with two sublineages (
D. jadovensis/
D. krbavensis and
D. adspersus/
D. ghetaldii) separating in the Lower Miocene, and the species are of Pliocene origin. Even though the ancestral distribution analysis () did not yield clear results for earlier nodes, it seems likely that the ancient
Delminichthys population, as the ancestor of all the present day species, was distributed in the western and middle Dinarides region, concordant with today’s Gorski Kotar, Dalmatian hinterland and southern Bosnia and Herzegovina. The Matica River in the Jezero Field most likely harboured the ancestral population of
D. adspersus, while the origin of
D. ghetaldii might be connected with the karstic fields of eastern Herzegovina.
. Divergence time estimations within the genus <i>Delminichthys</i> based on cyt <i>b</i> sequences. The timing of splitting events is presented by mean values and the 95% credibility range (in brackets), in million years ago. Reconstruction of ancestral distribution ranges is marked by the nodes.
BSPs were estimated to analyse the changes in past population sizes of
Delminichthys species ().
Delminichthys adspersus seems to presently be in a reduction period, which was preceded by significant growth that lasted for about 220,000 years, while the BSPs of
D. ghetaldii and
D. krbavensis imply stability and minimal growth in the last period of their evolution.
Delminichthys jadovensis also appears to be stable.
. Bayesian skyline plots based on cyt <i>b</i> sequences. Changes in the effective population size (millions of individuals on a log scale; Y-axes) are depicted over time (in million years; x-axes). The black central lines represent the median values of Ne, while the purple lines represent the 95% highest posterior density of the Ne estimates.
Analysis of the intraspecific structure based on the cyt
b gene revealed small migrations among
D. adspersus populations (from Baćinska Lakes to Matica River, M = 5809.32, Nm = 4.5), and migrations among
D. ghetaldii populations (from Ombla River to Bilećko Lake, M = 341.44, Nm = 4.5; from Ombla River to Popovo Field intensive migration with M = 24019). Estimates for migrations among the remaining population pairs were not significantly different from 0. Even though control region data set was employed for gene flow estimations, no reliable results were obtained.
The measures of genetic polymorphism based on cyt
b, as the most relevant genetic marker, revealed that moderate to high levels of genetic diversity are characteristic for all species and populations ().
Delminichthys ghetaldii and
D. krbavensis contain higher genetic diversities than the other two species, whereas genetic polymorphism of
D. jadovensis is the lowest within this endemic genus. CR pointed to similar conclusions, but overall levels of genetic diversity are lower and differences among populations and species even more pronounced. Based on CR,
C. jadovensis express low level of genetic diversity and the remaining species contain moderate levels of genetic polymorphisms. As expected, the RAG1 gene yielded smaller values of genetic polymorphism, with populations and species containing small to moderate levels of genetic diversity. Contrary to the situation with mtDNA,
D. krbavensis contains a very small polymorphism of RAG1 (only 3 haplotypes found in the sample of 24 sequences, Hd = 0.163). Despite the low genetic diversity of
D. adspersus based on RAG1, it is moderate in the Vrljika population. Higher values of nuclear polymorphism of
D. ghetaldii are visible in the number of polymorphic sites, haplotype diversity and average number of nucleotide differences. The remaining nuclear markers, IRBP and BA expressed highly variable results for different populations and species. The haplotype diversity in the IRBP data set ranged between 0 (recorded in
C. jadovensis despite significant sample number) and 1 (recorded in the Bilećko Lake population of
D. ghetaldii). Generally, species with distribution ranges located more to the south (
D. adspersus and
D. ghetaldii) have higher genetic polymorphisms based on nuclear markers, whereas
D. jadovensis express lower or no genetic diversity in nuclear markers. RAG1 is an exception because, besides enabling better phylogenetic resolution, implied different results regarding genetic polymorphism measures (high in
D. ghetaldii, low to moderate in the remaining species).
. Measures of genetic polymorphism of the investigated populations and species. As it is not possible to calculate genetic diversity measures for populations where a single individual was sampled, those localities are not included in this table. However, when calculating the genetic diversity of a species, all individuals belonging to that species were included.
4. Discussion
4.1. Evolutionary History of Delminichthys
Although mito-nuclear discordance was observed, mitochondrial DNA corroborated a clear separation and monophyly of the four
Delminichthys species. As already proposed by [
3], the closest relative of the genus
Delminichthys is
Pelasgus. Assuming a lineage divergence rate of 0.4%, the origination of the genus
Delminichthys can be dated back to the Oligocene/Miocene boundary, which is earlier than proposed by [
3], but similar to the results obtained by [
12]. This was a period of increased continentalization in Europe and significant tectonic activity in the Mediterranean region [
27], which, together with a warm climate, might have induced ancient divergences within freshwater systems. However, intrageneric divergences are of a much younger origin. More than 15 “million years of silence” presents a period with events that did not leave a trace in the recent genetic structure of this genus. A period from the lower Miocene or early Pliocene until the middle Pleistocene has been proposed as important for divergences of freshwater fishes in the Adriatic basin, and this period also brought divergences within the
Delminichthys genus. Connection of the ancestor of recent
Delminichthys species with the central Dinarides region was corroborated by the ancestral ranges’ reconstruction. Separation of two sister species pairs is concordant with the separation of
Cobitis bilineata Canestrini, 1865 (from the Zrmanja River) and the ancestor of
C. dalmatina Karaman, 1928 and
C. narentana Karaman, 1928 (Cetina and Neretva River basins) [
13], which was likely induced by changes in the Dinaric Lake system (DLS). This system of tectonic freshwater lakes originated in the early Miocene within the western Dinaric belt [
28,
29,
30]. In accordance with reports that evolution of DLS during favourable Miocene climatic conditions enabled diversification of several snail genera [
29], as well as the freshwater fishes of the genus
Cobitis [
13] and
Telestes [
5], it might have shaped early divergences within the genus
Delminichthys. Furthermore, ancestral geographic ranges fall within the geographic range of the DLS. It seems that the spectacular Miocene radiation [
29] comprised even more animal taxa than previously considered and might be one of the most important evolutionary events for freshwater fishes in the Adriatic basin. Divergence of
D. adspersus and
D. ghetaldii occurred more than a million years earlier than the separation of
T. miloradi Bogutskaya, Zupančič, Bogut & Naseka, 2012 and
T. metohiensis (Steindachner, 1901) [
5], species that are distributed in similar areas. The youngest divergence event, the separation of
D. jadovensis and
D. krbavensis, occurred in the Middle Pliocene. Interestingly, their separation is somewhat younger than the separation of
Telestes croaticus (Steindachner, 1866) and
T. fontinalis (Karaman, 1972), distributed in the same area [
5]. Thereafter, it can be concluded that although several geological events generally influenced the evolutionary history of cypriniform fishes recently distributed in the Dinaric region, it is also evident that the evolutionary history of the endemic
Delminichthys genus is not completely concordant with other fishes.
4.2. Peculiar Case of Nuclear Introgression
Hybridization phenomenon, leading to the pertinence of genetic material of two species in a single individual or population, has been repeatedly observed in fishes (e.g., [
31,
32,
33,
34]). There is even evidence that in some cases, significant proportions of genomes are derived by hybridization [
34] or that hybridization led to speciation [
35]. Even though hybrid biotypes (considered as any taxonomic or ecologic units originated by past or current hybridization events) blur the species borders and obscure the definition of species as taxonomic units, they certainly present an amazing part of biodiversity that is still vastly unknown and without adequate protection. With only a few reported exceptions, hybrid biotypes that succeeded in reproducing and persisting for a longer period originated as a result of the interbreeding of closely related species (e.g., [
32,
36]). Furthermore, although widely observed in the central and northern European freshwater basins, hybrid biotypes are the exception in the Adriatic watershed, which has been explained by the long-term isolation of these river basins (see [
13]).
Based on the genetic structure of
Delminichthys species revealed here, the existence of hybrid lineages in the Adriatic watershed can be hypothesized and the possibility that, in addition to allopatric speciation, hybridization might also played yet unnoticed, but nevertheless important role in the evolution of biodiversity in this ara. The mito-nuclear discordance observed within
Delminichthys, with mitochondrial markers implying monophyletic origin of the four
Delminichthys species and RAG1 discovering polyphyletic origin of
D. ghetaldii in the Ombla River, whereas nuclear markers in general did not enable species distinction, is the probable result of two phenomena. Incomplete lineage sorting is the likeliest explanation that the nuclear markers, having lower mutation rates, did not enable resolution of intrageneric structuring of
Delminichthys, yet all species (with the exception of
D. ghetaldii from the Ombla River, based on the RAG1) share the same or similar nuclear haplotypes. However, extensive nuclear introgression of RAG1 observed in the Ombla population invokes the possibility of its hybridogenetic origin. Mitochondrial haplotypes of the Ombla population are closely related to those from Bilećko Lake and, together with morphological determination, corroborate the pertinence of the Ombla population to
D. ghetaldii. On the other hand, both RAG1 alleles of all specimens from the Ombla River are very distinct from the remaining
Delminichthys sequences. They are the most similar, though not entirely, to
T. pleurobipunctatus RAG1 haplotypes. Even in heterozygous individuals, both alleles clustered with
T. pleurobipunctatus, which is very unusual since
T. pleurobipunctatus is currently distributed in Greece. The time, locality and pathway of their connection are not known. The possiblet explanation for this peculiar case of nuclear DNA capture is ancient hybridization that caused mosaicism in the nuclear genome, although further, preferably genomic investigations are required to distinguish between possible ancient hybridization and incidental cross-generic gene flow. Noteworthy, pattern of mitochondrial isolation with nuclear introgression is highly unusual. The reversed situation of mitochondrial introgression has been reported among different taxa of freshwater fishes and occurs much more frequently than the observed nuclear introgression (e.g., [
10,
11]). Namely, cytoplasmically inherited genes are more likely to cross species boundaries than nuclear ones [
37,
38]. Accordingly, there are only several reports on nuclear introgression causing cytonuclear and/or nuclear/nuclear discordance [
39,
40]. It is thought that nuclear genes may be prone to introgression in cases with different incompatibilities between species, depending on which species is the female parent [
39], or that some loci are selectively advantageous in the “wrong” species and not closely linked to genes causing intraspecific incompatibility so they might more easily cross species boundaries [
40]. Since there are no evidences for nuclear introgressions of IRBP and BA genes, we can suspect that hypothesis on selective advantage of certain nuclear genes is more likely. The hypothesys that the RAG1 pattern is a retention of an ancestral polymorphism is not supported by our phylogenetic analysis that revealed distant relations among the two parental taxa. On the other hand, presence of the same haplotype in various species, as observed in nuclear markers, is a likely example of ancestral polymorphism retention. This is also corroborated by higher frequencies of those haplotypes, which is usually characteristic for ancestral haplotypes.
4.3. Subterranean Migrations
One of the most astonishing adaptations of freshwater fishes for life in karstic rivers might be their ability to survive in and migrate through underground waters. Namely, due to the carbonate terrain that makes up the karstic surfaces in the Dinarides Mountains area, surface water flows are scarce. The hydrographic network of the Adriatic watershed is fragmented, comprised of short, isolated rivers. The underground water network is, on the other hand, well developed, complex and widespread [
41], providing possible shelter grounds, migration and/or colonization routes for fishes. Nevertheless, evidence of underground migrations have only been provided for
D. adspersus [
4]. No significant underground migrations among
D. adspersus populations investigated here were detected (with migration observed only between the Matica River and Baćinska Lakes that are connected by an artificial tunnel), while migrations were found among
D. ghetaldii populations that have no surface connections. Adaptation to underground life, thereafter, seems to be characteristic to this genus, developed during the specific evolutionary history in the Dinaric karst freshwater systems. The revelation of morphological characters responsible for this adaptation should assist in explaining the evolutionary history of this species. Results of the migration between the Ombla River and Popovo Field populations are not likely precise, due to low sample size from the Popovo Field. Nevertheless, the indication of the Ombla population as a source of immigrants is very important for two reasons. First, since individuals migrate upstream, it provides evidence that underground migrations within this genus are active. Secondly, this knowledge is important as it indicated the need for immediate, effective conservation of the Ombla population.
4.4. Genetic Diversity
No significant reduction of genetic polymorphism was observed in species distributed in areas affected by glaciations (
D. jadovensis and
D. krbavensis), which might also be a consequence of its ability to survive in underground shelters. High genetic diversity contained within all
Delminichthys species is most likely a consequence of their old origin. Recent genetic diversity presents a reservoir for coping with future environmental changes and a guarantee of the viability of a species. It does seem possible that hybridization led to an increased mutation rate [
42,
43] and higher genetic variability. Namely, genetic polymorphism of the RAG1 gene within the Ombla River population is unexpectedly high (more than threefold higher than in the conspecific Bilećko Lake population, higher than the cyt
b polymorphism of this population), which is contrary to what was expressed by the cyt
b gene. Lower genetic diversities observed in IRBP and BA are probable consequences of their lower mutation rates and ancestral polymorphism retention phenomenon. High genetic diversity of the majority populations and species, and the stability observed by BSPs are favourable when discussing the extinction risk of the
Delminichthys species. Their restricted distribution ranges and high level of endemicity, nevertheless, require maximum care and effective protection.
Supplementary Materials
The following supporting information can be found at: https://www.sciepublish.com/article/pii/285, Table S1: Accession numbers and references of the sequences retrieved from the GeneBank; Figure S1: MP phylogenetic trees of Delminichthys IRBP and BA sequences. Numbers at nodes represent MP and ML branch support.
Acknowledgments
We are very grateful for suggestions and corrections made by two reviewers, that helped us improve this paper!
Author Contributions
Conceptualization: I.B.; Data curation: I.B., Z.M., M.Ć. and R.Š.; Formal analysis: I.B., S.R.; Investigation: I.B., Z.M., M.Ć., R.Š. and S.R.; Methodology: I.B.; Resources: I.B. and R.Š.; Software: I.B. and S.R.; Writing—original draft: I.B., Z.M. and M.Ć.; Writing—review & editing: I.B., Z.M., M.Ć., R.Š. and S.R.
Ethics Statement
This study was conducted entirely in accordance with the ethical standards and Croatian legislation, and was approved by the Ethics Committee of the Faculty of Science, University of Zagreb. Sampling permits were issued by the competent authorities: for localities within Croatia, the permit was issued by the Ministry of Nature Protection; and for localities within Bosnia and Herzegovina the permit was issued by the Republic Department for the Protection of Cultural, Historical and Natural Heritage in Banja Luka.
Informed Consent Statement
Not applicable.
Funding
This research received no external funding.
Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
References
1.
Ćaleta M, Buj I, Mrakovčić M, Mustafić P, Zanella D, Marčić Z, et al. Endemic Fishes of Croatia; Croatian Environment Agency: Zagreb, Croatia, 2015.
2.
Guinand B, Durand JD, Laroche J. Identifying main evolutionary mechanisms shaping genetic variation of
Leuciscus cephalus L. 1758 (Cyprinidae) in Western Greece: Discordance between methods.
Comptes Rendus de l'Académie des Sciences-Series III-Sciences de la Vie 2001,
324, 1045–1060. doi:10.1016/s0764-4469(01)01361-0.
[Google Scholar]
3.
Perea S, Bohme M, Zupancic P, Freyhof J, Sanda R, Ozulug M, et al. Phylogenetic relationships and biogeographical patterns in Circum-Mediterranean subfamily Leuciscinae (Teleostei, Cyprinidae) inferred from both mitochondrial and nuclear data.
BMC Evol. Biol. 2010,
10, 265. doi:10.1186/1471-2148-10-265.
[Google Scholar]
4.
Palandačić A, Matschiner M, Zupančić P, Snoj A. Fish migrate underground: The example of
Delminichthys adspersus (Cyprinidae).
Mol. Ecol. 2012,
21, 1658–1671. doi:10.1111/j.1365-294X.2012.05507.x.
[Google Scholar]
5.
Buj I, Marčić Z, Ćaleta M, Šanda R, Geiger MF, Freyhof J, et al. Ancient connections among the European rivers and watersheds revealed from the evolutionary history of the genus
Telestes (Actinopterygii; Cypriniformes).
PLoS ONE 2017,
12, e0187366. doi:10.1371/journal.pone.0187366.
[Google Scholar]
6.
Dubut V, Fouquet A, Voisina A, Costedoat C, Chappaz R, Gilles A. From Late Miocene to Holocene: Processes of Differentiation within the Teleost Genus (Actinopterygii: Cyprinidae).
PLoS ONE 2012,
7, e34423. doi:10.1371/journal.pone.0034423.
[Google Scholar]
7.
Gilles A, Costedoat C, Barascud B, Voisin A, Banarescu P, Bianco PG, et al. Speciation pattern of
Telestes souffia complex (Teleostei, Cyprinidae) in Europe using morphological and molecular markers.
Zool. Scr. 2010,
39, 225–242. doi:10.1111/j.1463-6409.2010.00417.x.
[Google Scholar]
8.
Kottelat M, Freyhof J. Handbook of European Freshwater Fishes; Kottelat & Freyhof: Berlin, Germany; Cornol, Switzerland, 2007.
9.
Buj I, Marčić Z, Čavlović K, Ćaleta M, Tutman P, Zanella D, et al. Multilocus phylogenetic analysis helps to untangle the taxonomic puzzle of chubs (genus
Squalius: Cypriniformes: Actinopteri) in the Adriatic basin of Croatia and Bosnia and Herzegovina.
Zool. J. Linn. Soc. 2020,
189, 953–974. doi:10.1093/zoolinnean/zlz133.
[Google Scholar]
10.
Freyhof J, Lieckfeldt D, Pitra C, Ludwig A. Molecules and morphology: Evidence for introgression of mitochondrial DNA in Dalmatian cyprinids.
Mol. Phylogenet. Evol. 2005,
37, 347–354. doi:10.1016/j.ympev.2005.07.018.
[Google Scholar]
11.
Perea S, Vukić J, Šanda R, Doadrio I. Ancient Mitochondrial Capture as Factor Promoting Mitonuclear Discordance in Freshwater Fishes: A Case Study in the Greece
Squalius (Actinopterygii, Cyprinidae) in Greece.
PLoS ONE 2016,
11, e0166292. doi:10.1371/journal.pone.0166292.
[Google Scholar]
12.
Freyhof J, Lieckfeldt D, Bogutskaya NG, Pitra C, Ludwig A. Phylogenetic position of the Dalmatian genus
Phoxinellus and description of the newly proposed genus
Delminichthys (Teleostei: Cyprinidae).
Mol. Phylogenet. Evol. 2006,
38, 416–425. doi:10.1016/j.ympev.2005.07.024.
[Google Scholar]
13.
Buj I, Ćaleta M, Marčić Z, Šanda R, Vukić R, Mrakovčić M. Different Histories, Different Destinies-Impact of Evolutionary History and Population Genetic Structure on Extinction Risk of the Adriatic Spined Loaches (Genus
Cobitis; Cypriniformes, Actinopterygii).
PLoS ONE 2015,
10, e0131580. doi:10.1371/journal.pone.0131580.
[Google Scholar]
14.
Brito RM, Briolay J, Galtier N, Bouvet Y, Coelho MM. Phylogenetic Relationships within Genus
Leuciscus (Pisces, Cyprinidae) in Portuguese Fresh Waters, Based on Mitochondrial DNA Cytochrome
b Sequences.
Mol. Phylogenet. Evol. 1997,
8, 435–442. doi:10.1006/mpev.1997.0429.
[Google Scholar]
15.
Waap S, Amaral AR, Gomes B, Coelho MM. Multi-locus species tree of the chub genus Squalius (Leuciscinae: Cyprinidae) from western Iberia: New insights into its evolutionary history.
Genetica 2011,
139, 1009–1018. doi:10.1007/s10709-011-9602-0.
[Google Scholar]
16.
Stephens M, Smith NJ, Donnelly P. A New Statistical Method for Haplotype Reconstruction from Population Data.
Am. J. Hum. Genet. 2001,
68, 978–989. doi:10.1086/319501.
[Google Scholar]
17.
Stephens M, Scheet P. Accounting for Decay of Linkage Disequilibrium in Haplotype Inference and Missing-Data Imputation.
Am. J. Hum. Genet. 2005,
76, 449–462. doi:10.1086/428594.
[Google Scholar]
18.
Fu YX, Li WH. Statistical Tests of Neutrality of Mutations.
Genetics 1993,
133, 693–709.
[Google Scholar]
19.
Tajima F. Statistical Method for Testing the Neutral Mutation Hypothesis by DNA Polymorphism.
Genetics 1989,
123, 585–595.
[Google Scholar]
20.
Librado P, Rozas J. DnaSP v5: A software for comprehensive analysis of DNA polymorphism data.
Bioinformatics 2009,
25, 1451–1452. doi:10.1093/bioinformatics/btp187.
[Google Scholar]
21.
Swofford DL. PAUP*: Phylogenetic Analysis Using Parsimony (*and Other Methods), Version 4 [Computer Software and Manual]; Sinauer Associates: Sunderland, MA, USA, 2002.
22.
Drummond AJ, Suchard MA, Xie D, Rambaut A. Bayesian Phylogenetics with BEAUti and the BEAST 1.7.
Mol. Biol. Evol. 2012,
29, 1969–1973. doi:10.1093/molbev/mss075.
[Google Scholar]
23.
Yu Y, Harris AJ, He X. S-DIVA (Statistical Dispersal-Vicariance Analysis): A tool for inferring biogeographic histories.
Mol. Phylogenet. Evol. 2010,
56, 848–850. doi:10.1016/j.ympev.2010.04.011.
[Google Scholar]
24.
Beerli P. Estimation of migration rates and population sizes in geographically structured populations. In Advances in Molecular Ecology; NATO-ASI Workshop Series; Carvalho G, Ed.; IOS Press: Amsterdam, The Netherlands, 1998; pp. 39–53.
25.
Beerli P, Felsenstein J. Maximum likelihood estimation of a migration matrix and effective population sizes in n subpopulations by using a coalescent approach.
Proc. Natl. Acad. Sci. USA 2001,
98, 4563–4568. doi:10.1073ypnas.081068098.
[Google Scholar]
26.
Beerli, P. How to use Migrate or why are Markov chain Monte Carlo programs difficult to use? In Population Genetics for Animal Conservation, Volume 17 of Conservation Biology; Bertorelle G, Bruford MW, Hauffe HC, Rizolli A, Vernese C, Eds.; Cambridge University Press: Cambridge, UK, 2009; pp. 42–47.
27.
Rögl F. Mediterranean and Paratethys. Facts and hypotheses of an Oligocene to Miocene paleogeography (short overview).
Geol. Carpath. 1999,
50, 339–349.
[Google Scholar]
28.
Jiménez-Moreno G, Mandic O, Harzhauser M, Pavelić D, Vranjković A. Vegetation and climate dynamics during the early Middle Miocene from Lake Sinj (Dinaride Lake System, SE Croatia).
Rev. Paleobot. Palyno. 2008,
152, 237–245. doi:10.1016/j.revpalbo.2008.05.005.
[Google Scholar]
29.
Harzhauser M, Mandic O, Latal C, Kern A. Stabile isotope composition of the Miocene Dinaride Lake System deduced from its endemic mollusc fauna.
Hydrobiologia 2012,
682, 27–46. doi:10.1007/s10750-011-0618-3.
[Google Scholar]
30.
Harzhauser M, Mandic O. Neogene lake systems of Central and South-Eastern Europe: Faunal diversity, gradients and interrelations.
Palaeogeogr. Palaeoclimatol. Palaeoecol. 2008,
260, 417–434. doi:10.1016/j.palaeo.2007.12.013.
[Google Scholar]
31.
Bohlen J, Ráb P. Species and hybrid richness in spined loaches of genus Cobitis (Teleostei: Cobitidae) with a checklist of European forms and suggestions for conservation.
J. Fish Biol. 2001,
59 (Suppl. A), 75–89. doi:10.1006/jfbi.2001.1751.
[Google Scholar]
32.
Janko K, Culling MA, Ráb P, Kotlík P. Ice age cloning—Comparison of the Quaternary evolutionary histories of sexual and clonal forms of spiny loaches (
Cobitis; Teleostei) using the analysis of mitochondrial DNA variation.
Mol. Ecol. 2005,
14, 2991–3004. doi:10.1111/j.1365-294X.2005.02583.x.
[Google Scholar]
33.
Saitoh K, Chen W-J, Mayden RL. Extensive hybridization and tetraploidy in spined loach fish.
Mol. Phylogenet. Evol. 2010,
56, 1001–1010. doi:10.1016/j.ympev.2010.04.021.
[Google Scholar]
34.
Cui R, Schumer M, Kruesi K, Walter R, Andolfatto P, Rosenthal GG. Phylogenomics reveals extensive reticulate evolution in
Xiphophorus fishes.
Evolution 2015,
67, 2166–2179. doi:10.1111/evo.12099.
[Google Scholar]
35.
Mallet J. Hybrid speciation.
Nature 2007,
446, 279–283. doi:10.1038/nature05706.
[Google Scholar]
36.
Seehausen O. Hybridization and adaptive radiation.
Trends Ecol. Evol. 2004,
19, 198–207. doi:10.1016/j.tree.2004.01.003.
[Google Scholar]
37.
Whittemore AT, Schaal BA. Interspecific gene flow in sympatric oaks.
Proc. Natl. Acad. Sci. USA 1991,
88, 2540–2544. doi:10.1073/pnas.88.6.2540.
[Google Scholar]
38.
Matos JA, Schaal BA. Chloroplast evolution in the
Pinus montezumae complex: A coalescent approach to hybridization.
Evolution 2000,
54, 1218–1233. doi:10.1111/j.0014-3820.2000.tb00556.x.
[Google Scholar]
39.
Abe TA, Spence JR, Sperling FAH. Mitochondrial introgression is restricted relative to nuclear markers in a water strider (Hemiptera: Gerridae) hybrid zone.
Can. J. Zool. 2005,
83, 432–444. doi:10.1139/Z05-030.
[Google Scholar]
40.
Di Candia MR, Routman EJ. Cytonuclear discordance across a leopard frog contact zone.
Mol. Phylogenet. Evol. 2007,
45, 564–575. doi:10.1016/j.ympev.2007.06.014.
[Google Scholar]
41.
Bolić J, (Ed.). Croatian Waters—Monograph on Waters and Water Management in Croatia; Ministry of Agriculture, Forestry and Water Management: Zagreb, Croatia, 1992; p. 215. (In Croatian)
42.
Soltis DE, Soltis PS. The dynamic nature of polyploid genomes.
Proc. Natl. Acad. Sci. USA 1995,
92, 8089–8091. doi:10.1073/pnas.92.18.8089.
[Google Scholar]
43.
Salzburger W, Baric S, Strumbauer C. Speciation via introgressive hybridization in east African cichlids?
Mol. Ecol. 2002,
11, 619–625. doi:10.1046/j.0962-1083.2001.01438.x.
[Google Scholar]
