Unique , Ancient Stygobiont Clade of Hydrobiidae ( Truncatelloidea ) in Bulgaria : the Origin of Cave Fauna

Knowledge about subterranean aquatic gastropods is still limited and species inhabiting caves are typically only known from empty shells found at the surface. Phylogenetic analyses of such fauna are extremely rare. In the present work, we analysed seven species of five genera from 13 cave localities in Bulgaria, using molecular (mitochondrial COI and three nuclear markers: 18S rRNA, 28S rRNA and H3) techniques. All analysed DNA sequences from cave specimens formed one, highly supported, monophyletic lineage, divided into six subclades. According to molecular clock analyses, this subterranean clade originated 7-6.75 Mya. Our data support the “climatic-relict” hypothesis for the speciation of these taxa, assuming that surface ancestors colonized the caves and, after extinction of the surface population, cave populations evolved in allopatry. Extinction of the surface population was most probably caused by the climatic shift associated with the Messinian Salinity Crisis (5.96–5.33 Mya). A new model of climate-driven separation promoting speciation is proposed for the origin of the cave fauna.

The caves, intuitively known as rare and isolated habitats, harbouring specific fauna, whose origin, ecology, taxonomy and adaptation are still poorly known, are interesting objects to study.Their putative spatial isolation, simple community structure, and habitat specialization make them excellent models for studying evolutionary, biogeographic, and ecological issues.Worldwide there are over 350 described species of stygobiont (obligate subterranean-dwelling aquatic species) gastropods, 97% of them, and over 50 genera represent high flow, especially during spring (HAASE 1995;RICHLING et al. 2016).At a few localities, more living specimens can sometimes also be found in springs and streams at the surface.Collection of stygobiont Truncatelloidea is therefore difficult, especially in adequate numbers for any study of population genetic structure.There are also subterranean snails whose shells can only be found deep in caves, to say nothing of the extreme scarcity of living specimens (GEORGIEV 2011(GEORGIEV , 2012(GEORGIEV , 2013)).In fact, nothing but shell morphology is known for the vast majority of the cave gastropods.Animals adapted to cave environments are also thought to be highly geographically isolated because of their limited dispersal ability, resulting from limited physiological tolerances and, especially in the case of snails, physical limitations of their locomotion (PURCHON 1977;TRUEMAN 1983).However, this isolation is not obvious.In fact, there are many subterranean habitats which are not caves, but also harbour eyeless, depigmented animals, especially unconsolidated sediments bordering and underlying streams and rivers.They are parts of the interstitial habitat, neither rare nor discontinuous, making thus possible migration between caves (LAMOREAUX 2004;CULVER et al. 2009).
The limited available material has been, perhaps, the reason why molecular data are scarce.Thus, neither species boundaries nor phylogenetic relationships with the snails inhabiting surface environments are well understood.KANO andKASE (2004, 2008) studied neritopsid Neritilia cavernicola Kano et Kase, 2004, from anchialine (connected with sea) caves on two islands in the Philippines.Cytochrome oxidase mtDNA sequences reflected no genetic isolation, most probably due to freeswimming veliger larval dispersal.In Georissa (Neritimorpha) -a land snail from a limestone cave in Borneo, 16S mtDNA sequences analysis (SCHILTHUIZEN et al. 2005) showed that its ancestor was the local epigean population of Georissa saulae (van Benthem-Jutting, 1966), inhabiting the rainforest directly at the cave entrances.The shell morphometry revealed that the troglobite has diverged morphologically, but there are populations of intermediate morphology in the twilight zone of the cave.Thus, ongoing gene flow between the epigean and subterranean (sub)populations, via the intermediate ones, and the troglobite divergence without prior genetic isolation was hypothesized.
As concerns the stygobiont Truncatelloidea, BICHAIN et al. (2007) studied Bythinella in France, also in two caves (Padirac and Folatière).COI (mtDNA) sequences indicated several subterranean lineages, still belonging to the epigean species, thus confirming that colonizations of hypogean habitats by epigean individuals are not rare events in the genus, as already earlier suggested (BOETERS 1979;GIUSTI & PEZZOLI 1982;BOLE & VELKOVRH 1986;BERNASCONI 2000).BICHAIN et al. (2007) found in sympatry or parapatry with those epigean lineages also distinct, exclusively hypogean phylogenetic lineages, indicating thus multiple invasions of the caves and/or radiation already in caves.
COI sequences were used to estimate the time of divergence between Heleobia dobrogica (Grossu et Negrea 1989) from the Movile Cave in Dobrogea, Romania and epigean H. dalmatica (Radoman, 1974) (FALNIOWSKI et al. 2008). SZAROWSKA et al. (2014) demonstrated that all the European brackishwater Heleobia Stimpson, 1865 belong to the same species, slightly polymorphic, thus H. dalmatica may represent the (hypothetical) epigean ancestor.The time 2.172±0.171Mya coincided with the period marking the beginning of the fall in temperature and precipitation that predated the Pleistocene, when H. dobrogica found a safe shelter within a warm cave, and the epigean populations became extinct in Dobrogea, but not in much less cold Dalmacija. FALNIOWSKI and SARBU (2015) described new species of Iglica and Daphniola from the Melissotripa Cave in Greece, inferring their molecular relationships with epigean species.OSIKOWSKI et al. (2015) studied the genetic and morphological differences between Bulgarian surface-and cave-dwelling populations of the genus Bythinella, also belonging to the Truncatelloidea.RYSIEWSKA et al. (2016) have studied 16 populations of the Bulgarian Pontobelgrandiella.
Current phylogeographic studies confirm theoretical assumptions that animal cave taxa are often cryptic -morphologically undistinguishable -and have highly restricted geographical distributions, despite potential gene flow from surface populations (JUAN et al. 2010).Caves are relatively stable, long-lasting environments and individual ones often have an island character with no subterranean connections to any others.However, as already mentioned above, such connections may exist for stygobiont fauna.In some cases, particular caves can be characterized by endemic taxa with long, independent, evolutionary histories (e.g., FALNIOWSKI et al. 2008) that differ strongly from their sister taxa occurring outside caves.
The competing "climatic-relict" and "adaptiveshift" hypotheses have been proposed to explain the origins of cave organisms (JUAN et al. 2010).According to the "climatic-relict" model, preadapted ancestors took refuge (or simply colonized) in caves, and gradually adapted to the cave environment.Later, the surface climate was altered by glaciation or aridification, which caused local extinctions of the surface populations, leaving each subterranean relict population to evolve allopatrically in an isolated cave system.The "adaptiveshift" model assumes that preadapted ancestral species actively colonized caves to exploit novel resources and diverged under a gene flow scenario.Divergent selection between surface and subterranean habitats gradually overcame the homogenizing process of gene flow and eventually led to parapatric speciation.
Molecular phylogenies are helpful in testing those hypotheses, since they allow reconstruction of sister relationships between cave and surface taxa (RIVERA et al. 2002;JUAN et al. 2010).Phylogenetic relationships -especially sister-clade ones -as well as spatial distribution of genetically divergent cave and surface populations: either surface and cave populations in allopatry, or existence of only subterranean relict species, or parapatric distribution of both surface and cave populations, may be criteria of distinction between the two hypotheses (RIVERA et al. 2002;LEYS et al. 2003;JUAN et al. 2010).Allopatric distribution of surface ancestor vs subterranean descendant, without any zone of contact, or existence of the relict subterranean species with absence of surface one, may indicate 'climatic-relict' model.On the other hand, a troglobiont/stygobiont which highly diverged from parapatric epigean ancestor may indicate "adaptive-shift" hypothesis.Distinction between the two models may be more difficult when some of the ancestral, surface population have become extinct.The relative time of branching of the troglob-ites, and possible correlations with geological events, like changes of climate, are also helpful.
The aim of the present study is: (i) to test the monophyly of the Bulgarian cave gastropods and to infer their possible origin and evolutionary history; and (ii) to infer the model of origin of Bulgarian cave snails ("climatic-relict" or "adaptive-shift").To achieve this, mitochondrial COI and three nuclear genes: H3, 18S, 28S were analysed.

Sample collection and fixation
Snails were collected from 13 cave localities in Bulgaria (Fig. 1, Table 1) either by hand from stones and rocks or sieving sediments and bacterial mats with a 0.5 mm mesh.The snails (Fig. 2) were washed in 80% ethanol and left to stand in it for about 12 hours.The ethanol was then changed twice during 24 hours.After a few days, the samples were transferred to 96% ethanol and stored at -20°C prior to DNA extraction.

DNA extraction and sequencing
DNA was extracted from foot tissue using a Sherlock extracting kit (A&A Biotechnology) and dissolved in 20 µl TE buffer.From the mtDNA we amplified a fragment of the cytochrome oxidase subunit I (COI), using the following primers: LCOI490 Fig. 1.Localities of the studied cave-inhabiting gastropods in Bulgaria; for details see Table 1.Colours correspond to COI clades (see Fig. 3).
The PCR conditions were as follows: COI: initial denaturation step of 4 min at 94°C, followed by 35 cycles of 1 min at 94°C, 1 min at 55°C, 2 min at 72°C; H3: initial denaturation step of 2 min at 94°C, followed by 35 cycles of 30 s at 94°C, 30 s at 50°C, 1 min at 72°C; 18S: initial denaturation step of 4 min at 94°C, followed by 40 cycles of 45 s at 94°C, 45 s at 51°C, 2 min at 72°C; 28S: initial denaturation step of 2 min at 94°C, followed by 35 cycles of 30 s at 94°C, 30 s at 50°C, 90 s at 72°C; after all cycles were completed, there was an additional elongation step of 4 min at 72°C.
Ten µl of the PCR product were run on 1% agarose gel to check the quality.The PCR product was purified using Clean-Up columns (A&A Biotechnology).The purified PCR product was sequenced in both directions using BigDye Terminator v3.1 (Applied Biosystems), following the manufacturer's protocol and using the primers described above.The sequencing reaction products were purified using ExTerminator columns (A&A Biotechnology), and the sequences were read using an ABI Prism sequencer.

Molecular Data Analysis
Sequences were initially aligned in the MUSCLE (EDGAR 2004) program in MEGA 6 (TAMURA et al. 2013) and then checked in BIOEDIT 7. 1.3.0 (HALL 1999).Basic sequence statistics, like nucleotide divergence, were calculated in DnaSP 5.10 (LIBRADO & ROZAS 2009) and MEGA 6.The presence of substitution saturation was determined with the program DAMBE 5 (XIA 2013), using the saturation test of XIA et al. (2003).
In the phylogenetic analysis, we used 19 reference sequences of the COI and 9 reference sequences of 18S genes from GenBank (Table 2).The sequences represented the main lineages of Sadlerianinae and Hydrobiinae, to infer the relationships of the Bulgarian stygobiont clade.To demonstrate that the Bulgarian clade does not represent any Bythiospeum (in the literature reported from Bulgaria so far), in analysis we included also three species of the Moitessieriidae, and Truncatella scalaris as an outgroup (Table 2).T. scalaris was used as an outgroup for three nuclear genes as well.The data were analysed using an approach based on Bayesian inference (BI) and maximum likelihood (ML).Prior to phylogenetic analyses, we calculated the best fit nucleotide substitution models under the Bayesian information criterion in PartitionFinder 1.1.1(LANFEAR et al. 2012) testing for all possible partitioning schemes.The concatenated dataset (mtDNA and three nuclear genes) was partitioned by genes, with specific model for each partition (GTR+Ã+I for COI and HKY+I for all nuclear genes).
The Bayesian analyses were run with MrBayes ver.3.2.3 (RONQUIST et al. 2012).The default value of four Markov chains was used (one cold and three heated) and each chain was started from a random tree.The Monte Carlo Markov chain length was 50,000,000 generations and trees were sampled every 10,000 generations.Trees generated before the steady phase were discarded as burn-in (first 25% of trees sampled).The analyses were summarized on a 50% majority-rule tree.Convergence of the Bayesian runs was checked with Tracer v1.6 (RAMBAUT et al. 2014) and AWTY (http://ceb.csit.fsu.edu/awty)(WILGENBUSCH et al. 2004).
A maximum likelihood (ML) approach was performed with RAxML v8.0.24 (STAMATAKIS 2014), with GTR+Ã+I model.The sufficient number of bootstrap replicates was automatically detected by the analysis.RAxML analyses were performed using the free computational resource CIPRES Science Gateway (MILLER et al. 2010).
To estimate the time of divergence, molecular clock was used for COI, applying the calibration point of WILKE (2003), the separation of two hydrobiids, Peringia ulvae (Pennant, 1777) and Salenthydrobia ferreri Wilke, 2003during Messinian Salinity Crisis (WILKE 2003).The divergence  (TAJIMA 1993) was performed in MEGA.As Tajima's RRTs and the LRT test rejected the equal evolutionary rate throughout the tree for snails, time estimates were calculated using a penalized likelihood method (SANDERSON 2002) in r8s v.1.7 for Linux (SANDERSON 2003).For the comparison with penalized likelihood analysis, we conducted Bayesian evolutionary analysis for COI sequences using the BEAST package (v.1.8.2) (DRUMMOND & RAMBAUT 2007).Alignments of COI sequences, including outgroup, were imprted into BEAUti v.1.8.2.Substitution models followed jModelTest results and relaxed uncorrelated log-normal clocks were used.Four independent runs of 10,000,000 were conducted using Yule tree prior and default other priors and operator settings.The resulting files were combined in LogCombiner v.1.8.2 after removing the first 10% as burn-in.Outputs of the runs were examined for convergence and effective sample size in Tracer v. 1.6.A maximum clade credibility tree was obtained with TreeAnnotator v.1.8.2.

Results
The shells of the Bulgarian stygobiont taxa presents Fig. 2.
The test of XIA et al. (2003) revealed no saturation in any of the datasets.Topologies of the trees obtained from BI and ML analyses were identical.
In total we obtained 25 sequences of the COI gene (552 bp, GenBank Accession numbers MF179879-MF179903).All analysed COI sequences from Bulgarian cave snails belonged to one, 83% bootstrap-supported mitochondrial lineage (Fig. 3), with nucleotide diversity ð = 0.085.This lineage differed from the sister lineage, grouping the other Sadlerianinae, by 16.45%.The divergence time estimates between these two lineages ranged between 7±1 Mya, and 6.75±0.8Mya, according to the Bayesian (BEAST) and penalized likelihood (r8s) methods, respectively, showing that these methods give similar results.
The sequences of the Bulgarian cave snails lineage were grouped into six main mitochondrial clades (Fig. 3), with p-distances varying from 0.078 to 0.118 between them (Table 3).All those clades were very highly supported (bootstrap support 97-100%), but their mutual relationships remained unresolved.Considering the low bootstrap values, as well as the short branches, one large polytomy groups all the six clades.The first clade included sequences of Balkanica yankovi and Balkanospeum schniebsae.These two species were represented by one haplotype each, differing by p-distance only 0.0044.Two earlier described species also formed two different haplotypes in the second clade: genus Devetakia.Devetakia mandrica and D. krushunica differed by p-distance = 0.010.The third clade was also divided into two haplotypes (p-distance = 0.007), consisting of sequences from two populations of Cavernisa zaschevi.The sequences of Pontobelgrandiella belong to the fourth clade.The remaining two clades were more closely related (Fig. 3), p-distance = 0.078.The sequences of Stoyanovia stoyanovi were identical and come from one population.Sequences of Devetakiola devetakium (Table 1) formed one clade, differ by p-distance = 0.022, although they come from two populations from caves separated by approximately 170 km.
In addition, for each of the three nuclear genes, H3 (341 bp, GenBank Accession numbers MF179904-MF179921), 18S (318 bp, GenBank Accession numbers MF179922-MF179939) and 28S (707 bp, GenBank Accession numbers MF179940-MF179957), we analysed 18 sequences.The topologies of the trees, reflecting relationships between the clades, are different and poorly supported (Fig. 4), confirming the polytomy inferred with COI, most probably a hard one.However, the    support for each of the clades was high, confirming their distinctness.The p-distances between these clades for all three nuclear markers varied from 0.008 to 0.028 (Table 3).
The best resolution of the tree was given by the H3 marker (Fig. 4).Sequences of this gene formed six, well supported, main clades identical to those for COI.It was similar for 28S, were only Devetakia and Pontobelgrandiella did not form distinct clades.The 18S marker provided the least information, with sequences only forming Balkanica/Balkanospeum, Stoyanovia and Pontobelgrandiella as distinct clades.
We also analysed COI and nuclear genes sequences together (Fig. 5), since a partition homogeneity test indicated no conflicting phylogenetic signals between the datasets (p = 0.985).On this tree, cave populations also form six main distinct clades.Devetakiola and Stoyanovia clades are clearly distinct and in the outgroup position with respect to the remaining four, whose relationships remain unresolved, due to the hard polytomy.For large number of potential outgroup taxa, only COI and 18S sequences were available in GenBank (for comparison we sequenced also three additional taxa: Belgrandiella zermanica and two Pontobelgrandiella species).Thus, we analysed our COI and 18S sequences together, pooled with 12 outgroup taxa (Fig. 6).In the obtained tree, the whole cave lineage as well as its six main clades were confirmed and best supported.
The results of morphological studies of some of snails from localities presented in this study are presented elsewhere (GEORGIEV et al. 2017).

Monophyly and phylogeny of the Bulgarian stygobiont Truncatelloidea
All the studied Bulgarian cave truncatelloid gastropods belonged to the same, monophyletic, well supported lineage.It is noteworthy that the lineage does not include any taxa inhabiting surface habitats (with possible exception of Pontobelgrandiella: see below), all are stygobiont.The lineage belongs to the Hydrobiidae, subfamily Sadlerianinae.Thus Devetakiola and Stoyanovia, originally described as Bythiospeum, may not belong to the genus Bythiospeum Bourguignat, 1882.Moreover, "real" Bythiospeum (B.francomontanum, B. saxigenum), as a representative of the Moitessieriidae, formed a phylogenetically distant clade (Fig. 3), and p-distance between "real" Bythiospeum and Stoyanovia stoyanovi/Devetakiola devetakium was 0.187.
Balkanica yankovi is molecularly very close to Balkanospeum schniebsae, but the soft part mor-phology and, especially, the shells are very distinct (Fig. 2; GEORGIEV et al. 2017).The species distinctness of the two nominal species of Devetakia remains doubtful, since neither morphological nor molecular differences have been found between those inhabitants of two separate (although not distant one from another) caves.Moreover, in all the nuclear loci Devetakia do not differ, and in the COI tree D. krushunica is paraphyletic.
Within the main cave lineage, there are six highly supported clades, not less distinct than the other hydrobiid genera in the tree (Fig. 3).Considering the low values of support for the nodes connecting these six clades, one large polytomy should be assumed.Thus, the history of the cladogenesis cannot be reconstructed.Considering the fact that addition of more loci does not resolve this polytomy, and the bootstrap supports are always low, this polytomy seems to be a hard one, reflecting not too low number of characters considered, but rather parallel simultaneous speciation (in the sense of geological time measurement) of the six clades.There is no geographic pattern in the distribution of the clades: some of them have been found at one locality (Balkanica, Balkanospeum), or two localities (Devetakia) close to each other, while the other clades are distributed over wider areas that overlap one another, although it has to be stressed that our knowledge on their distribution remains poor.

Origin and history of the cave lineage
As stated above, two alternative models of the speciation of cave fauna have been proposed: "climatic-relict" and "adaptive-shift" (HOWARTH 1973;HOLSINGER 2000;RIVERA et al. 2002;JUAN et al. 2010).The first of these has been proposed for continental ecosystems of the temperate zone (HOLSINGER 1988(HOLSINGER , 2000;;PECK & FINSTON 1993).However, DANIELOPOL & ROUCH (2012) questioned this model, since in rapidly deteroriated environment there would be no time for evolutionary change: the invasion of the subterranean habitats must occur earlier, while the extinction of epigean populations may be rapid.
Perhaps somewhat modified model should be proposed: "climate-driven separation promoting speciation".Some specimens of an epigean population penetrate caves and/or other subterranean habitats, colonizing them gradually.Underground specimens gradually become better adapted to underground environment, and more sensitive for unstable conditions at the surface.Thus, when climate become severe, underground specimens abandon the zone close to the entrance of a cave, thus the zone of possible gene exchange disappear.The epigean population(s) may survive or not.
Our results do not allow to clearly conclude which of the possible mechanisms is responsible for observed differentiation of the Bulgarian stygobiont lineage.Presented data suggest, however, that "climatic-relict" model is more likely to occur.All six subterranean clades (Balkanica/Balkanospeum, Devetakia, Cavernisa, Pontobelgrandiella, Devetakiola, Stoyanovia) originated relatively long ago, and with one possible exception no parapatric surface population were discovered.Moreover, divergence of these clades took place in the same time, which corresponds to climatic and geological events.However, the case of Pontobelgrandiella makes the pattern more complicated.
16 populations of the Bulgarian Pontobelgrandiella have been studied (RYSIEWSKA et al. 2016), eleven COI haplotypes have been found.Divergence was small, and no intrapopulation polymorphism (with one exception) has been found.The haplotypes have formed four clades, p-distances between them have been small (0.007-0.015), but PCA on seven shell measurement confirmed, in part, the distinctness of those clades.The genetic distance between the clades have indicated the time of divergence from 0.38-0.43Mya to 0.81-0.92Mya, thus Calabrian and Middle Pleistocene.Pontobelgrandiella occurs in caves, but could also be found in springs close to the caves, as unnumerous specimens, but still alive, contrary to the representatives of all the other stygobiont Bulgarian clades.In fact, we do not know if those individuals are really epigean ones, or obligatory stygobionts only accidentally flooded from subterranean waters.Remarkably low genetic divergence between the population of Pontobelgrandiella was explained as a result of quite recent invasion of the subterranean habitats (RYSIEWSKA et al. 2016).However, our data indicate that Pontobelgrandiella clade is as old as the other five, undoubtedly stygobiont clades.Thus, we could hypothesize contrary that Pontobelgrandiella may rather recently invade epigean habitats.On the other hand, the observed lack of polymorphism may be a result of bottlenecks during subsequent local extinctions and/or drastic reductions of population size, and may be also a result of strong selection.
However, this low divergence coupled with relatively wide distribution at several localities may reflect high level of gene flow between the populations of Pontobelgrandiella, which may reflect adaptation of Pontobelgrandiella especially to the interstitial habitats, connecting other subterranean habitats; some analogies with the sea gastropods with larvae (ZHOU et al. 2016) may be drawn.In some cases, Pontobelgrandiella occurred in sympatry with undoubtedly obligatory subterranean species (e. g.Devetakia).Thus, one can speculate that Pontobelgrandiella may be ancestor of the other stygobiont clades, and this would support "adaptive-shift" concept, or "climate-driven separation promoting speciation".
The latter must be applied for the stygobiotic Bulgarian Bythinella spp.OSIKOWSKI et al. (2015) reported the representatives of this genus in four caves.Molecular analyses indicated that only one of these populations seemed to be clearly stygobiontic, being clearly phylogenetically distinct from the surface populations.The time of divergence of this clade (Clade V, formed by three haplotypes from the Vodni Pech Cave in NW Bulgaria) was estimated at 7.25±1.8Mya, which is strikingly similar to the estimates of separation of stygobiont subclades in the present paper.This is in congruence with the hypothesis that, about 7 Mya, the climatic changes resulted in gene flow restrictions between surface and cave populations, followed by the extinction of the surface populations.Later, the epigean populations were re-established by immigrants from surrounding epigean populations.The remaining sampled Bulgarian Bythinella cave populations did not form distinct clades, however, and their times of divergence from closely related surface relatives are much shorter than 7 Mya (OSIKOWSKI et al. 2015).
In Europe, the hot spots of underground diversity belong to Dinaric Karst (six caves).Among the subterranean aquatic invertebrates, Crustacea dominate, gastropods are the second group in species number (CULVER & SKET 2011).Subterranean species could be found in many phylogenetic lineages (VON RINTELEN et al. 2012). PAGE et al. (2008) studied wide-scale genetic diversity in subterranean (anchialine) shrimps Stygiocaris and Typhlatya, assigning the deepest cladogenesis to the time of Tethys.They found significant geographic intraspecific biodiversity within many groundwater species.Their results confirmed the hypothesis that common events or dispersal barriers, like sea-level changes, gorges cutting through to impermeable limestone layers, or hypersaline groundwater have isolated many populations into distinct geographic groups.As stated above, high genetic diversity and narrow geographic ranges are typical of stygobiont species.ZAKŠEK et al. (2009) studied interpopulation diversity in the cave shrimp Troglocaris anophthalmus (Kollar, 1848), known as stygobiont untypically widely distributed over almost the entire Dinaric Karst.They sequenced three loci in 232 specimens from the entire range, across more than 500 km.Their molecular data did not confirm such a wide range, unusual for subterranean animals, since four or five cryptic species were found.However, one of them -southern Adriatic phylogroup -presented an unusual pattern of recent dispersal across 300 kilometers of hydrographically frag-mented karst terrain.In the paper on widely distributed and ubiquitous groundwater amphipod Niphargus rhenorhodanensis, LEFÉBURE et al. (2007) found poor dispersal as main evolutionary factor in groundwater.All the studied genes supported the existence of numerous cryptic and mostly allopatric units within N. rhenorhodanensis, indicating that its apparently wide distribution range is an artifact generated by cryptic diversity.On the other hand, BUHAY and CRANDALL (2005) found in north american subterranean freshwater crayfishes extensive gene flow and surprisingly large evolutionary effective population sizes.In the paper on subterranean amphipod Niphargus virei (Chevreux, 1896) LEFÉBURE et al. (2006) found cryptic lineages 13 million years old; the diversity was shaped by extremely low dispersal abilities coupled with morphostasis, caused by the selection in extreme conditions of subsurface habitats.
Similar picture present some troglobionts.In the case of freshwater isopod Asellus aquaticus, two independent events of shifting to subterranean habitats occurred in different parts of Europe: Romania and Slovenia (KONEC et al. 2015).Within the clade formed by the Pyrenean coleopteran tribe Trechini (FAILLE et al. 2010) multiple invasions of the cave habitats has been postulated, in the midlate Miocene (ca. 10 Mya), when the climate in the Pyrenees area was, generally, warmer and wetter than today (Messinian Salinity Crisis, with desert at the place of Recent Mediterranean) confirming "climatic-relict" hypothesis.Nesticella cave spiders from China, an ancient group (originated in the middle Miocene), inhabited caves by multiple colonizations, as late as in Pleistocene, when the dramatic climate changes could have caused a burst of divergence in Nesticella (ZHANG & LI 2013).Subterranean beetles of the tribe Leptodirini inhabiting the mountains of the Western Mediterranean probably developed all modifications to subterranean life in the Early-Mid Oligocene-Miocene, with no evidence of extinct, less modified surface lineages (RIBERA et al. 2010).
The divergence time of the Bulgarian cave gastropod lineage was estimated at 7±1.0 Mya or 6.75±0.8Mya by the Bayesian and penalized likelihood techniques, respectively, which reflects the time of definitive genetic isolation from the hypothetic ancestral surface species.This gave a rate of 2.35 or 2.44±0.26% per Mya, which is only slightly higher than the estimated rate of 1.83±0.21%per Mya for COI provided by WILKE (2003).A similar estimate, of 1.62% per Mya, was obtained by HERSHLER and LIU (2008).The dating suggests the Messinian Salinity Crisis.During this period the Mediterranean Sea desiccated, which resulted in higher temperatures and aridity of the climate.This may have resulted in interruption of contact between subterranean and epigean populations, since the former tended to escape down the caves: climate-driven separation concept may be the most adequate.In the studied area, the neighborhood of Paratethys waters made the climate milder, but still the humidity must have been lower and temperatures much higher than before this period (MURPHY et al. 2009).

Fig. 3 .
Fig. 3.The maximum-likelihood phylogram for the COI gene.Bootstrap support and Bayesian posterior probabilities were shown.For geographical distribution of the clades see Fig. 1.Sequences obtained in present work are marked in bold.

Table 1
Sampling localities of the studied populations with their geographical coordinates.

Table 2
(SWOFFORD 2002)00)(2008)ic analyses with their GenBank accession numbers and references between these two, with correction according toFALNIOWSKI et al. (2008), was fixed at 5.96 Mya.The likelihoods for trees with and without the molecular clock assumption for a Likelihood Ratio Test (LRT)(NEI & KUMAR 2000)were calculated with PAUP(SWOFFORD 2002).The Relative Rate Test (RRT) time

Table 3 p
-distances between main clades of Bulgarian cave Hydrobiidae for the COI (below diagonal) and nuclear (above diagonal) genes