RADIATION IN BYTHINELLA MOQUIN-TANDON, 1856 (MOLLUSCA: GASTROPODA: RISSOOIDEA) IN THE BALKANS

A BSTRACT: 61 sequences of cytochrome oxidase subunit 1 (COI), 570 bp long, from Central Europe (8), Slovenia (2), and from the Balkans (51: 6 from, Romania, 15 from Greece, 15 from Bulgaria, 6 from Serbia, and 9 from Montenegro), 33 of them new, together with 61 sequences of the ribosomal 18S, 450 bp long, all of them new, were analysed to infer the pattern of radiation of Bythinella in the Balkans. Thirty two nominal taxa of Bythinella (22 nominal species: B . austriaca , B . calimanica , B . charpentieri , B. compressa, B . dacica , B . dispersa , B . grossui , B . hansboetersi , B . luteola , B . micherdzinskii , B . molcsanyi , B . nonveilleri , B . pannonica, B . pesterica , B . radomani , B . rhodopensis , B . robiciana , B . schmidti , B . slaveyae , B . srednogorica , B . taraensis and B . viseuiana ; one nominal subspecies: B . austriaca ehrmanni ; and nine Greek species not yet described) were included, represented mostly by paratypes or at least topotypes, collected at 31 Balkan localities. The phylogeny, inferred on the combined data set with the ML approach, showed two large clades, although they were weakly supported. One of them comprised the Romanian and Montenegro populations, and one Serbian population, the other (less genetically diversified) consisted of one Serbian and all the Bulgarian/Greek populations. The origin of those two clades was dated, with the external data, to be more than 4.341±0753 MYA old, thus its origin was assigned to isolation by the Dacic Basin (part of Paratethys). All the Bulgarian populations presumably belong to one species, which may be assigned to the recent recolonisation of this territory from the south.


INTRODUCTION
The genus Bythinella Moquin-Tandon, 1856 is distributed from west Europe (the Iberian Peninsula), to west Asia. These dioecious, oviparous snails with minute, dextral, ovoid shells inhabit freshwater springs, small brooks and groundwaters (GIUSTI & PEZZOLI 1977, FALNIOWSKI 1987. They are found in abundance on mosses and other aquatic plants, fallen leaves, dead wood, etc. The rich literature on Bythinella (RADOMAN 1976, 1983, GIUSTI & PEZZOLI 1977, FALNIOWSKI 1987, GLÖER 2002, SZAROWSKA & WILKE 2004) pertains mainly to western, southern and central Europe. The characters of either the shell or the soft parts do not allow for clear delineation of Bythinella species. Recently, BICHAIN et al. (2007a, b), HAASE et al. (2007), FALNIOWSKI et al. (2009b, c) and BENKE et al. (2009) applied molecular data and proved the distinctiveness of several species of the genus. The ranges of some Bythinella species studied molecularly so far (e.g. BICHAIN et al. 2007a, b, FALNIOWSKI et al. 2009b, c, BENKE et al. 2009, FAL-NIOWSKI & SZAROWSKA 2011) are restricted to a single spring, spring complex, or local watershed, while those of the other congeners stretch across one or more drainage divides. There were several species described and redescribed in central Europe (JUNGBLUTH & BOETERS 1977, FALNIOWSKI 1987, GLÖER 2002, SZAROWSKA & WILKE 2004, HAASE et al. 2007) and the former Yugoslavia (RADOMAN 1976, 1983, GLÖER 2008, GLÖER & PEŠIAE 2010. Recently we found six molecularly distinct species of this genus in Romania (FALNIOWSKI et al. 2009a, b) and ten molecularly distinct (9 of them not yet described) species in Greece (FALNIOWSKI & SZAROWSKA 2011). Several nominal species were recently reported from Bulgaria (GLÖER & PEŠIAE 2006, FALNIOWSKI et al. 2009a, GLÖER & GEORGIEV 2009) and Montenegro (GLÖER & PEŠIAE 2010). FALNIOWSKI et al. (2009b, c) and FALNIOWSKI & SZAROWSKA (2011) demonstrated that (1) the morphostatic mode of evolution as defined by DAVIS (1992) is common in Bythinella; (2) it is impossible to distinguish Bythinella species if there are no molecular data, although the morphology must be considered as well (BICHAIN et al. 2007b, HAASE et al. 2007. In this study we applied mitochondrial cytochrome oxidase subunit I (COI), and ribosomal 18S to infer the general pattern of radiation in this genus in the Balkan Peninsula. Most of the specimens sequenced were paratypes or at least topotypes of the nominal taxa; considering all the controversies concerning the species distinctness within Bythinella (GIUSTI & PEZZOLI 1977, FALNIOWSKI 1987, BICHAIN et al. 2007a, b, HAASE et al. 2007, FALNIOWSKI et al. 2009b, c, FALNIOWSKI & SZAROWSKA 2011, we leave unresolved the taxonomic status of those nominal taxa. The Balkans once harboured one of the main European glacial refugia (HEWITT 1996, VOGEL et al. 1999, PETIT et al. 2003, SCHMITT 2007, FALNIOW-SKI et al. 2009b). Hence, they provide an interesting area for a study of radiation. In Romania (FALNIOW-SKI et al. 2009a, b) we found greater genetic differences between Bythinella species than in Greece (FAL-NIOWSKI & SZAROWSKA 2011). However, the data on the genus inhabiting the zone between Romania and Bulgaria are scarce (FALNIOWSKI et al. 2009a). Thus the aim of the present paper was to answer the following questions: (1) Can the pattern of genetic differentiation, decreasing from north to south, be applied to all the Balkans? (2) Does the pattern of (presumably) vicariant speciation reflect any palaeogeographic conditions in the Balkans?

MATERIAL COLLECTION AND FIXATION
We collected snails, with a sieve or by hand, at 30 localities in the Balkans (Fig. 1, Table 1). Snails were washed twice in 80% ethanol and left to stand in it for ca. 12 hours, after which the ethanol was changed twice in 24 hours and finally, after a few days, the 80% solution was replaced with a 96% one. The material was stored at -20°C.

MOLECULAR TECHNIQUES APPLIED
Snails were hydrated in TE buffer (3 × 10 min.) and their DNA was extracted with the SHERLOCK extracting kit (A&A Biotechnology); the final product was dissolved in 20 µl of TE buffer. The PCR reaction was performed with the following primers: LCOI490 (5'-GGTCAACAAATCATAAAGATATTGG-3') and COR722b (5'-TAAACTTCAGGGTGACCAAAAAA TYA-3') for the COI gene (FOLMER et al. 1994) and SWAM18SF1 (5'-GAATGGCTCATTAAATCAGTCGA GGTTCCTTAGATGATCCAAATC-3'), SWAM18SR1 (5'-ATCCTCGTTAAAGGGTTTAAAGTGTACTCATT CCAATTACGGAGC-3') for the 18S gene (PALUMBI 1996). The PCR conditions were as follows: 1. For CO1: 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, after all cycles an additional elongation step of 4 min. at 72°C was performed; 2. For 18S: 4 min. at 94°C followed by 40 cycles of 45 sec. at 94°C, 45 sec. at 51°C, 2 min. at 72°C, after all cycles an additional elongation step of 4 min. at 72°C was performed. The total volume of each PCR reaction mixture was 50 µl. To check the quality of the PCR products, 10 µl of the PCR product was run on 1% agarose gel. The PCR product was purified using Clean-Up columns (A&A Biotechnology). The purified PCR product was dye-terminator amplified in both directions  using BigDye Terminator v3.1 (Applied Biosystems), following the manufacturer's protocol and with the primers described above. The sequencing reaction products were purified using Ex-Terminator Columns (A&A Biotechnology), and the sequences were read using an ABI Prism sequencer.

DATA ANALYSIS
The COI sequences were aligned by eye using Bio-Edit 5.0.0 (HALL 1999), and edited with MACCLADE 4.05 (MADDISON & MADDISON 2002). We examined mutational saturation for the COI dataset with saturation test of XIA et al. (XIA et al. 2003), performed with DAMBE 5.2.9 (XIA 2000. In 18S alignment was performed using CLUSTALX 1.82 (THOMPSON et al. 1997).
Phylogenies of COI and 18S haplotypes were inferred using maximum-likelihood (ML), minimum evolution (ME) and maximum parsimony.
The maximum-likelihood technique of phylogeny reconstruction (SWOFFORD et al. 1996, NEI & KUMAR 2000, FALNIOWSKI 2003) although criticised, is commonly used for molecular data and considered the most reliable one, so we decided to use the ML approach for each of the two data sets. For each maximum likelihood analysis, we tested different models of sequence evolution using Modeltest v3.06 (POSADA & CRANDALL 1998, POSADA 2003. Following the recommendations of POSADA & BUCKLEY (2004) and SOBER (2002), the best model for each dataset was chosen, using the Akaike Information Criterion (AKAIKE 1974). ML analyses were performed in PAUP*4.0b10 (SWOFFORD 2002), using an heuristic search strategy with stepwise addition of taxa, 10 random-sequence addition replicates, and tree-bisection-reconnection (TBR) branch swapping (SWOFFORD et al. 1996). Nodal support was esti-  (FELSENSTEIN 1985). Bootstrap values for ML trees were calculated using 1,000 bootstrap replicates, the "fast" heuristic search algorithm, and the same model parameters as used for each ML analysis. Minimum evolution (ME) and maximum parsimony (MP) were run on PAUP*; nodal support was estimated using the bootstrap approach (full heuristic search) with 1,000 replicates. PAUP was used to test the molecular clock hypothesis for COI with a likelihood ratio test (LRT) upon the estimated tree and best-fit model (FELSENSTEIN 1981, NEI & KUMAR 2000, POSADA 2003. In the phylogeny reconstruction for COI, five central European Bythinella species (B. austriaca, B. compressa, B. pannonica, B. robiciana and B. schmidti: Table 2), and B. micherdzinskii (Table 1) were used as outgroups. Next, the partition homogeneity test (FARRIS et al. 1995) was performed (1,000 replicates) with PAUP*, to check whether the two genes could be analysed together. The maximum likelihood heuristic search was then run for the combined molecular data. K2P (KIMURA 1980) distances for the COI data were calculated using PAUP*4.0b10 (SWOFFORD 2002). There is a rich literature concerning the use of K2P distances for COI data, thus we chose this distance for the sake of comparison of the levels of differentiation.

RESULTS
In total, we analysed 63 sequences of COI, 570 bp long. Two of the sequences represented the outgroup species (Bithynia tentaculata and Marstoniopsis insubrica: Table 2 The partition homogeneity test (p=0.732) showed that the two genes could be analysed together. For the combined data set the Akaike Information Criterion (AIC) with ModelTest selected the model HKY+I+Ã, with base frequencies: A = 0.3077, C = 0.2397, G = 0.1821, T = 0.2704; substitution rate matrix: ti/tv ratio: 5.2693, proportion of invariable sites: (I) = 0.7340, and Ã distribution with the shape parameter = 1.0465. Figure 2 shows the ML tree computed for the combined data, in which Bythinella pannonica from Slovakia forms a distinct Bythinella clade closest to the root. B. micherdzinski, B. austriaca and B. austriaca ehrmanni from Poland form another clade. The Balkan taxa form two clades that are distinct but weakly supported (59/64/89, depending on the technique): clade A (Fig. 2) consisting of populations 1-9 (Table 1), B. robiciana, and B. schmidti, although population 4, representing B. radomani (marked with asterisk in Fig. 2) falls between the clade of the Polish Bythinella and the other Balkan clade (Fig. 2: A'). The second large Balkan clade (B in Fig. 2) clusters all the other Balkan populations 10-31 (Table 1) (this bifurcation is marked with arrow in Fig. 2, the dotted line in Fig. 1 indicates the geographic border between the two large clades). Within clade A the differences between the nominal taxa, reflected in branch lengths (Fig. 2), as well as in K2P distances (within group mean: 0.1004), are greater than in clade B (within group mean: 0.0499). The pairwise K2P distances between the two clades are 0.09-0.1744, mean 0.1224, net mean 0.0473. FALNIOWSKI & SZAROWSKA (2011) demonstrated that in the Greek Bythinella there was much less interspecific differentiation than in the Romanian Bythinella, and the diversity decreased from north to south. In this study we observed the same general pattern throughout the East Balkans. The somewhat larger distances between the studied Romanian populations may contribute to this pattern. However, this cannot be the main explanation, as between the two sympatric species at locality 4 in Romania the distance was not smaller than between the allopatric Romanian species. Thus another explanation should be invoked.

DISCUSSION
The Romanian Pleistocene refugium, or rather a set of small territories with a mild climate, was situated farther to the north than the other European refugia. In fact, it was rather a group of nunatak-like habitats, the fauna of which was probably under a strong influence of unstable climatic conditions. Across the country there were several scattered glacial refugia, which, during the interglacials and especially in the postglacial, served as distribution centres (FALNIOWSKI et al. 2009b). The expansions and contractions of the glacier proceeded by way of climatic changes and were followed by the contractions and expansions of the nunatak refugia, thus Bythinella ranges in the studied territory may have promoted genetic differentiation within the genus (HEWITT 1989). This contrasts with the south Balkan (Greek) refugium, which was rather continuous geographically and characterised by relatively stable conditions. Perhaps those more stable conditions during the Pleistocene are reflected in smaller differences between the species of Bythinella in the south of the Balkans.
The two large clades in the inferred phylogeny are well-defined geographically. The dotted line in the map (Fig. 1) marks the border between these two clades. The first (clade A, localities 1-9) includes the six Romanian species (localities 1-5), the three Montenegro populations (localities 6 and 8-9), and one of two Serbian populations of Bythinella (locality 7). This group is closer to B. pannonica from Slovakia and Hungary, B. austriaca and B. micherdzinskii from Poland, B. compressa from Germany, as well as B. schmidti and B. robiciana from Slovenia, and to the outgroup. In fact, B. radomani from locality 4 (A' in Fig. 2) clusters between the Polish taxa and the second large clade. This may be explained, in part, by the location of that species in the northernmost part of Romania, thus not far from the Polish localities. One can speculate that B. radomani invaded the locality from the north/northwest.
The second clade (clade B) clusters all populations from Greece and Bulgaria together with one population (10) from Serbia. As stated above, the genetic diversity, as marked in COI, was lower within this clade. It is noteworthy that presumably only one species occurs across Bulgaria. As could be seen in the map (Fig. 1), there are, unfortunately, no samples close to the border marked on the map, thus its localisation is a very rough estimate only. Moreover, we cannot prove the real distinctness of those two clades, since their bootstrap supports are rather low (Fig. 2). Thus we can only hypothesise about this part of the history of Bythinella in the Balkans. There are no data on the fossil occurrence of Bythinella in the studied part of Europe. Fossil Bythinella were found in the Italian Pliocene/Pleistocene boundary (southern margin of the Alps, north of Bergamo: ESU & GIANOLLA 2009). In the Balkans -their southern part, at least -the history of this genus cannot be shorter. As we rejected the molecular clock hypothesis Radiation in Bythinella in the Balkans -and had not any precise point of callibration that would enable us to apply the techniques of relaxing molecular clock assumptions -we could not put a precise time frame on our tree. However, we can use, with some caution, the estimate of FALNIOWSKI & SZAROWSKA (2011), dating the origin of the SE clade (containing all the Greek and Bulgarian populations), thus the time of cladogenesis marked with the arrow in Fig. 2 certainly was before 4.341±0753 MYA (middle Pliocene). This coincides in time with the existence of the Dacic Basin, a vast water body that separated the present Carpathians from the present middle Bulgaria (POPOV et al. 2004, 2006, POPESCU et al. 2009). The Dacic Basin, 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, eventually decreasing in size, was still present till the latest Pliocene, 1.8 MYA. The Dacic Basin most probably separated the ancestors of the two large clades. Later, in the Pleistocene, the unstable fluvio-lacustrine system in SW. 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 the vast part of Bulgaria. Considering the data known so far, the small differences among the nine Bulgarian populations may reflect the short history of Bythinella in the area that most probably was recolonised from the south not earlier than in the late Pleistocene.