Does the genetic structure of spring snail Bythinella (Caenogastropoda, Truncatelloidea) in Bulgaria reflect geological history?

Abstract Bythinella is a minute dioecious caenogastropod that inhabits springs in central and southern Europe. In the Balkans, previous studies have addressed its morphological and genetic differentiation within Greece and Romania while the Bulgarian species have remained poorly known. The aim of the present paper has been to expand the knowledge on the subject in Bulgaria. Shell morphology and anatomy of the reproductive organs were examined, and a fragment of the mitochondrial cytochrome oxidase subunit I (COI) gene and the nuclear ribosomal Internal Transcribed Spacer 1 (ITS-1) were sequenced from 15 populations. Additional sequences from eight previously studied populations were included in our analyses. Phylogenetic analyses revealed five main mitochondrial DNA clades, which were partly confirmed by analyses of the ITS-1 sequences. The genetic differentiation between the clades was found to be in the range p=2.4-11.8%. Most of the populations belonged to clade I, representing Bythinella hansboetersi, and were distributed in SW Bulgaria. Clades II and III inhabit areas adjacent to clade I and were most closely related with the latter clade. Much more distinct were clade V, found at one locality in NW Bulgaria, and clade IV, found at one locality in SE Bulgaria, close to the sea. Four populations were found in caves, but only one of these represented a distinct clade. Considering the observed pattern of interpopulation differentiation of Bythinella in Bulgaria, we can suppose that isolation between clades I, II and III may have been caused by glaciations during the Pleistocene. The time of isolation between the above three clades and clade IV coincides with the Messinian Salinity Crisis, and the time of isolation between the clade V and the other four most probably reflects the isolation of the Rhodopes from western Balkan Mts by the seawater of the Dacic Basin.


2011).
In Bythinella, the shell as well as the penis exhibit substantial variability, often even within populations (Falniowski 1987), including all the character states considered by the authors cited above. Last but not least, the species-level taxonomy of Bythinella is heavily flawed by common assumption that different localities should harbour different species. Molecular studies by Falniowski et al. (2009a) on five populations from Bulgaria revealed low inter-population genetic distances, leading to the conclusion that all of them belonged to the same species. This result strongly suggests that it is necessary to incorporate the study of molecular markers in critical revisions of the Bythinella species in Bulgaria.
The palaeogeographic history of central and southern Europe has significantly influenced the distribution of fauna and flora in this region. An ecological event with a large impact on biodiversity in present day Bulgaria and Romania was the flooding of the Dacic Basin with seawater. This basin separated the Carpathians from central Bulgaria (Suc et al. 2011) from about 8 until about 1.8 Million years ago (Mya). Previous research has indicated that the genetic divergence of Bythinella is much higher in neighbouring Romania than in Bulgaria, suggesting that the Dacic Basin may have caused the extinction of this genus in a vast parts of Bulgaria (Falniowski et al. 2009b. However, this scenario was suggested based on limited sampling in Bulgaria and, therefore, requires further study. Among many nominal species of Bythinella described in Bulgaria on the basis of morphological characters, some cave taxa have been identified. B. markovi was reported from the Gargina Dupka Cave Georgiev 2009, 2011;Georgiev and Glöer 2013) and B. gloeeri from Lepenitsa Cave in the Rhodopes (Georgiev 2009). These two caves are situated about 80 km apart from each other as the crow flies in different ridges of the Rhodopes as parts of different river catchments. These two species differ substantially in their morphology. Current phylogeographic studies confirmed theoretical assumptions that cave animal taxa are often cryptic and possess highly restricted geographical distributions despite potential gene flow from surface populations (Juan et al. 2010). To test this hypothesis on Bulgarian Bythinella, confirmation of its distinctness, previously described on basis of morphological characters, is necessary using molecular techniques. Since low genetic differentiation between the surface Bythinella population in Bulgaria has been reported (Falniowski et al. 2009a), more data on the phylogenetic differentiation of potentially distinct cave species is needed.
The aim of our study has been to improve the knowledge of Bythinella distribution in Bulgaria through extended sampling and to answer the following questions: 1) Is the low genetic divergence previously reported for Bulgarian Bythinella a fact or the result of poor sampling? 2) Are the biogeographical patterns and phylogenetic relationships of Bythinella correlated with the geological history of the region? 3) Have the cave populations of Bythinella been isolated for a long time from the ones in surface water reservoirs? 4) Do molecular data support the opinion of Glöer and Georgiev (2011) about Bulgaria as a hot-spot of diversity of Bythinella? To answer these questions, both morphological (shell) and molecular (COI and ITS-1 genetic markers) characters were examined.

Snail sampling and fixation
Bythinella snails were collected from 15 sites across Bulgaria (Fig. 1, Table 1). In four of these sites the snails were found in caves. Six of the studied populations were from the type localities of the nominal species.
Snails were collected by hand or with a sieve. Individuals for the morphological study were fixed in 4% formaldehyde and stored in 80% ethanol, while individuals for molecular analyses were washed in 80% ethanol and left to stand in it for about 12 hours. The ethanol was then changed twice during 24 hours and, after a few days, 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 the SHERLOCK extracting kit (A&A Biotechnology) and dissolved in 20 µl TE buffer. PCR was performed in the reaction mixture of 50 µl total volume using the following primers: LCOI490 (Folmer et al. 1994) and COR722b (Wilke and Davis 2000) for the COI gene, and two Bythinella-specific primers ITS1D and ITS1R for the ITS-1 (Bichain et al. 2007). The PCR conditions were as follows. COI -initial denaturation step of 4 min at 94 °C, followed by 35 cycles at 94 °C for of 1 min, 55 °C for 1 min, 72 °C for 2 min, and a final extension of 4 min at 72 °C; ITS-1 -initial denaturation step of 4 min at 94 °C, followed by 25 cycles at 94 °C for 30 s, 60 °C for 30 s, 72 °C for 30 s, and a final extension of 5 min at 72 °C. Ten µl of the PCR product was run on 1% agarose gel to check for quality. PCR products were purified using Clean-Up columns (A&A Biotechnology). The purified PCR products were 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.

Morphological studies
Snails were dissected under a NIKON SMZ-U stereo-microscope with a NIKON drawing apparatus, and a CANON EOS 50D digital camera was used to photograph the shells.
Sequences obtained from Bythinella specimens in the present work were used in a phylogenetic analysis with other sequences obtained from GenBank (Table 1). The data were analysed using Bayesian inference (BI) and the maximum likelihood (ML) approach. We applied the GTR + I + Γ model because over-parameterization seems to be less dangerous for BI analyses than under-parameterization (Huelsenbeck and Rannala 2004). For ML analyses, GTR + I + Γ is the only nucleotide substitution model implemented in RAxML.
The Bayesian analyses were run with MrBayes ver. 3.2.3 (Ronquist et al. 2012) using default priors. Two simultaneous analyses were performed, each lasting 10,000,000 generations, with one cold chain and three heated chains, starting from random trees and sampling trees every 1000 generations. The first 25% trees were discarded as burnin. The analyses were summarised on a 50% majority-rule tree.
A maximum likelihood (ML) approach was conducted in RAxML v8.0.24 (Stamatakis 2014). One thousand searches were initiated with starting trees obtained through the randomized stepwise addition maximum parsimony method. The tree with the highest likelihood score was considered as the best representation of the phylogeny. Bootstrap support was calculated with 1000 replicates and summarized onto the best ML tree. RAxML analyses were performed using free computational resource CIPRES Science Gateway (Miller et al. 2010).
To test the molecular clock, COI data were used. Two hydrobiids, Peringiaulvae Pennant, 1777 and Salenthydrobia ferreri Wilke, 2003 (AF478401, AF478410) were used as outgroups. The divergence time between these two species (5.96 Mya) was used to calibrate the molecular clock, with correction according to Falniowski et al. (2008), since the isolation started with the beginning, not the end of the Messinian Salinity Crisis. The likelihoods for trees with and without the molecular clock assumption for a Likelihood Ratio Test (LRT) (Nei and Kumar 2000) were calculated with PAUP. The Relative Rate Test (RRT) (Tajima 1993) was performed in MEGA. As Tajima's RRTs and the LRT test rejected the equal evolutionary rate throughout the tree for Pseudamnicola, time estimates were calculated using a non-parametric rate smoothing (NPRS) analysis with the recommended Powell algorithm, in r8s v.1.7 for Linux (Sanderson 1997(Sanderson , 2003.

Results
Selected shells of Bythinella from some of the studied localities are presented in Fig.  2A-Q. It is visible that the variability at one locality (B10: Fig. 2G-Q) is equivalent to the variation observed amongst all populations.
We obtained 58 new sequences of COI (552 bp, GenBank Accession numbers KT381098-KT381155) and 36 new sequences of ITS-1 (234-264 bp, GenBank Accession numbers KT381156-KT381191). For COI the saturation test of Xia et al. (2003) revealed no saturation. Seventeen COI haplotypes (haplotype diversity Hd = 0.932) and ten ITS-1 haplotypes (Hd = 0.837) were identified. For phylogenetic analyses, additional Bythinella sequences available in GenBank (Table 1) were included, including those from eight sites in Bulgaria (15 haplotypes: Falniowski et al. 2009a, twelve sites in Romania (22 haplotypes: Falniowski et al. 2009b) and five sites in Montenegro (six haplotypes) and Serbia (five haplotypes: . The topologies of the resulting ML and BI phylograms were identical. Sequences of Bythinella viridis were used as outgroup in all analyses to root the trees. In the COI trees five main clades could be distinguished for the Bulgarian populations (Figs 3-4). Clade I included the largest number of haplotypes covering an area from the Rhodopes Mts through the Maritsa Valley to the Stara Planina and Sredna Gora Mts. This clade is characterized by a low sequence divergence (p-distance within this group = 0.008, Table 2). The relationships between the haplotypes of this clade are depicted in a haplotype network in Fig. 5. Most haplotypes from this clade belonged to snails inhabiting surface waters while two of them represent cave populations (Fig.  3). This clade represents several nominal species: B. angelovi, B. dierkingi, B. gloeeri, B. Clades II and III were most closely related to Clade I differing by intercladal pdistances of 0.024 and 0.026 and inferred divergence times of 0.82 and 1.89 Mya, respectively (Table 2). Clade II contained haplotypes from two sites in the south-west of Bulgaria: the Pirin Mts and a cave in the Rhodopes Mts, from which B. stoychevae  has been described. Only one haplotype formed Clade III representing B. gloeeri. All other sequences of this nominal species belonged to Clade I, however. The haplotype from the easternmost site (B5, in the Strandzha Mts) formed Clade IV. It differed from clades I to III by genetic distances of 7.8 to 9.0% (inferred di-  Table 2). This clade is situated between the two reference clades from Serbia and Montenegro. The most divergent clade was Clade V (inferred divergence time 7.25 Mya), formed by three haplotypes from the Vodni Pech Cave in north-western Bulgaria.
The reference sequences formed three distinct clades (Fig. 3). First of them represent one population from eastern Serbia, second the populations from western Serbia and Montenegro. Haplotypes from Romania formed another distinct lineage. The level of divergence between haplotypes within this lineage has been much larger than within the other clades.
Unfortunately, due to technical problems, ITS-1 sequences were not available for samples from Clade IV and for the reference populations from Serbia and Montenegro. The ITS-1 tree (Fig. 6) confirmed the distinctiveness of the Bulgarian Bythinella from the Romanian ones. Three clades (A, B, C) could be distinguished. In correspondence with the COI tree, Clade C, containing sequences from the population from the Vodni Pech Cave, was found to be the most divergent (similarly as COI haplotypes from this population). Clade A comprised about half of all haplotypes of COI Clade I plus the COI Clade III. Clade B comprised all samples from COI Clade II and the rest of the samples from COI Clade I.
Shells from different clades (Fig. 2) were found to differ in morphology, but no clade-specific shell characters were found. Similar remarks concern both female reproductive organs, and the penes: no clade-specific character states were found.

Discussion
Species delimitation in the genus Bythinella remains unclear. New species descriptions were initially based mainly on shell morphology, and on the locality not studied so far. Even for Radoman (1976), whose experience and extensive studies on the truncatelloidean anatomy were a basis for the new taxonomy proposed by him, the shell characters alone were the only basis for species-level taxonomy. He even stated that there could not be any differences in soft part morphology and anatomy between the congeneric species. Later, anatomy, especially of the reproductive system was considered (Boeters 1973(Boeters , 1998Falniowski 1987). In many cases, the number of species recognized from different parts in Europe is probably overestimated and the characters traditionally used to delimit species should be re-evaluated (Bichain et al. 2007b). In general, in Bythinella it is impossible to distinguish particular species without molecular data Falniowski et al. 2009b, c;Falniowski and Szarowska 2011). Analysis of the COI genetic distances between recognized clades/ lineages could be considered an efficient tool for the rapid assessment of biodiversity in Bythinella (Bichain et al. 2007a). Bichain et al. (2007a), after COI analysis, proposed a K2P value of 1.5% as the species threshold for European Bythinella. However, such a threshold value may be biased as well. So delimitation of Bythinella species needs to be backed up by additional data (e.g., more nuclear genes and mtDNA fragments, as well as detailed morphological studies). The history of any DNA fragment not necessarily reflects the history of speciation (Avise 2000), so multilocus analyses are necessary. The five molecularly distinct clades we have found may represent five distinct species. Especially the amount of mitochondrial differentiation of the COI clade V: 8.4-12.6 %, is within the range characteristic of the species level, and this divergence was confirmed by the ITS-1 as well. As could be clearly seen in the trees (Figs 3 and 6), there is rampant incongruence between phylogenetic patterns and current species-level taxonomy. Only B. gloeeri (Clade III) and B. stoychevae (Clade II) have been found to be molecularly distinct from B. hansboetersi (Clade I). However, all morphological character states given in the descriptions of the Bulgarian species are variable even within a population in Bythinella (Falniowski 1987, Mazan 2000, Mazan and Szarowska 2000a. Most populations examined here occurred in the area from the Rhodopes Mts to the Sredna Gora Mts (Clades I, II, III). Only two other isolated populations were found in eastern Bulgaria, in the Strandzha Mts (Clade IV) and in the north-western part of this country (Clade V). Despite an extensive search, no members of this genus were found in the rest of Bulgaria (Fig. 1). However, there are several areas where no Bythinella occur, since these snails are sensitive to environmental conditions, such as high water calcium content and low temperature. They occasionally occur in spring outlets and creeks, and also in caves or groundwaters Pezzoli 1977, 1980). Water conditions may be one of the most important factors influencing the occurrence of Bythinella. On the one hand, regional and global environmental changes may have relatively small effects on this spring snail, since springs can buffer such changes. On the other hand, Bythinella is more resistant than had been expected for a long time (e.g., Szarowska 2000) and, since springs are ephemeral habitats that are certainly not long-lasting, there must be unexpectedly high gene flow between them to colonize/ recolonize them (e.g., Falniowski et al. 1998Falniowski et al. , 1999. Low, infraspecies-level diversity characterized Clade I, including most of the studied populations distributed across central and western Bulgaria. The representatives of Clades II, III and V either migrated from the west, or survived there from earlier time. It seems possible that both clades III and V survived glaciations inside the caves. Clade IV most probably originated in the present Asia Minor. The closest sequences to clade V come from the B. turca haplotype from the Egirdir Lake in Turkey (p distance = 0.055).
Mitochondrial interpopulation differentiation of Bulgarian Bythinella (p distance = 0.03) is much smaller than in neighbouring countries (p distance = 0.06-0.08). The greater genetic differences between Bythinella populations in Romania have been compared with the surprisingly low differentiation amongst the Bulgarian ones by Falniowski et al. (2009a,b). However, this analysis was only based on a small number of Bulgarian populations. The more detailed sampling in the present work confirmed this phenomenon. Moreover, lower interpopulation differentiation than in the Romanian Bythinella was also demonstrated for Greek populations (Falniowski et al. 2011) and throughout the East Balkans .
Within Bulgaria, some geological events could explain the low divergence in Clade I. The Dacic Basin, a vast water body that separated the Carpathians from the recent central Bulgaria before and just after the peak of Messinian Salinity Crisis (5.60-5.46 Mya) (Popov et al. 2004(Popov et al. , 2006Popescu et al. 2009;Suc et al. 2011), was a part of the Paratethys, connected with the Pannonian Basin in the west, the Euxinian Basin in the east, and directly with the present Aegean Sea in the south. Although its waterfilled area eventually decreased in size, it was still present until the middle Pleistocene, about 1.8 Mya. The Dacic Basin most probably separated the ancestors of the two large clades, about 8 Mya. Later, in the Pleistocene, the unstable fluviolacustrine system in south-western Bulgaria and northern Greece, with glaciers present in the Pirin and Rila Mts (Zagorchev 2007), probably formed effective, temporary barriers for Bythinella, and may have caused its extinction in most of Bulgaria. Considering the data known so far, the small differences among the Bulgarian populations representing Clade I may reflect the short history of Bythinella in the area, which was most probably recolonised from the south, certainly not from the north, no earlier than in the late Pleistocene. Benke et al. (2011) revealed that genetic diversity of Bythinella in Europe is not distributed equally, and identified five "hotspots": Massif Central and Pyrenees, western Alps and northern Apennines, eastern Alps, western Carpathians and eastern Carpathians. The authors of the present paper discovered another Bythinella hotspot in central Greece (Szarowska et al. 2015). Thus, all the hotspots occur in mountain areas, which strongly suggests, that this type of landscape is especially favourable for Bythinella.
Moreover, almost all these hotspots are in places that were previously identified as Bythinella Pleistocene glacial refugia (Benke et al. 2009;Falniowski et al. 2009b), so high differentiation level in Romania may be the result of glaciations. During the Pleistocene these areas were probably a set of small areas of a nunatak character, with a mild climate suitable for Bythinella survival (Falniowski et al. 2009). Habitat fragmentation and subsequent periods of isolation in such shelters must have promoted speciation and could explain high differentiation level. It is widely accepted that, during glacial periods, the Pontic-Mediterranean refugium included territory in present-day Romania (e.g., Falniowski et al. 2009b). It seems that there is no trace of such refugium in Bulgaria.
Caves are relatively stable long-lasting environments and individual ones often have an island character with no subterranean connections to any others. In some cases, particular caves can be characterized by endemic taxa with long, independent, evolutionary histories (e.g., Falniowski et al. 2008, Juan et al. 2010 for references) that differ strongly from their sister taxa occurring outside caves. Bythinella inhabits both surface and underground waters, providing thus opportunity to compare populations from those two kind of habitats. In Bulgaria, only clade V was formed by haplotypes clearly distinct from the remaining Bythinella populations, which may reflect their troglobiontic character, and longevity of isolation, approaching the Pliocene. The collected individuals were, in fact, found not inside the cave, but at its entrance, in the water running from the cave, but this is normal way of collecting of several troglobiontic gastropods. In this population three haplotypes were found, in COI as well as in ITS-1, which is not common in troglobiontic animals (Juan et al. 2010), whose populations are usually monomorphic. This polymorphism may also confirm the longevity of this population, or inhabiting the cave by more than one species. Clade III was formed by a single COI haplotype from Lepenitza cave, locality B7, and represents B. gloeeri. This haplotype, originated probably in the Calabrian, Pleistocene, also presents a distinct, probably troglobiont lineage. In both cases -COI clades V and III, the "climatic relict" hypothesis, as proposed by Howarth 1973(Holsinger 1988Peck andFinston 1993, Rivera et al. 2002). Accordin to this hypothesis, after the colonization of subterranean habitats there still takes place gene flow between the subterranean and surface population, but later, in strict allopatry, the subterranean population still speciates, and the surface population becomes extinct as a result of climatic changes, like glaciations or growing aridity. This seems a typical pattern for temperate climate. It has to noted, however, that at the same cave there was found another haplotype, belonging to the COI clade I. This is one more example of sympatric occurrence of more than one Bythinella. Such a situation is not unlikely and has been previously reported for Bythinella (Radoman 1976;Falniowski et al. 2009b, Falniowski andSzarowska 2011). This haplotype is close to the one inhabiting surface waters at locality B14, situated close to this cave. Similarly, the COI haplotype from troglobiont population at locality B12 (Modarskata Cave) was close to the one from the surface population 13, situated closely to population B12; the same concerns the troglobiont population B8 from Manuilova Dupka Cave compared with surface population B1. Low divergence between those troglobiont populations and their surface relatives may reflect either the early phase of the "climatic relict"-model processes, or "adaptive shift"-model: adaptive evolution of the lineages invading subterranean habitats, coupled with survival of the ancestral population at the surface.
Considering the observed pattern of interpopulation differentiation of Bythinella in Bulgaria, the facts listed above, and the divergence time estimates, we could suppose that the isolation between clades I, II and III (0.82 Mya and 1.89 Mya, respectively) may have been caused by subsequent glaciations during the Pleistocene. The time of isolation between the above three clades and clade IV from SE Bulgaria (4.39 Mya) coincides with the Messinian Salinity Crisis. Later, the low level of the present Black Sea promoted migration of the representatives of this clade from Asia Minor to Europe. The distinctness of clade V, found at NW Bulgaria, most probably reflects the isolation of the Rhodopes from the western Balkan Mts by the Dacic Basin (7.25-1.8 Mya).