Research Article
Print
Research Article
Does the genetic structure of spring snail Bythinella (Caenogastropoda, Truncatelloidea) in Bulgaria reflect geological history?
expand article infoArtur Osikowski, Dilian Georgiev§, Sebastian Hofman|, Andrzej Falniowski
‡ Jagiellonian University, Department of Comparative Anatomy, Kraków, Poland
§ Plovdiv University, Plovdiv, Bulgaria
| Institute of Zoology, Kraków, Poland
¶ Jagiellonian University, Cracow, Poland
Open Access

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 B. 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.

Keywords

Gastropoda, phylogeography, Balkans, Messinian Salinity Crisis, Dacic Basin

Introduction

The genus Bythinella Moquin-Tandon, 1856 consists of minute (2–4 mm shell height), dioecious, oviparous freshwater snails inhabiting springs and subterranean waters across Europe and western Asia. In Europe, their range extends from the Iberian Peninsula to the Ukraine, and from southern Poland and Germany, to Sicily and Crete. These snails also occur on numerous Mediterranean islands and in the western part of Turkey (Falniowski et al. 2012; Benke et al. 2011).

Early studies on Bythinella taxonomy mainly concentrated on the morphology of the shell and, later, the external morphology and anatomy of the body (e.g., Boeters 1973, 1998; Radoman 1976; Falniowski 1987). It has been demonstrated, however, that morphology alone cannot be used for unequivocal species delimitation due to the limited number of taxonomically useful characters and their large variability (Falniowski 1987, 1992; Falniowski et al. 2009a; Bichain et al. 2007a, b). Recent molecular studies (Bichain et al. 2007a; Haase et al. 2007; Benke et al. 2011; Falniowski et al. 2009a, b, 2012; Falniowski and Szarowska 2011, 2012) have led to revised species delimitations. For Bythinella an allopatric mode of speciation and distribution of the species has generally been postulated. This caused authors to overestimate the importance and effectiveness of geographical isolation. This partly triggered numerous descriptions of new species, which were based solely on the occurrence of snails in springs not previously studied. However, the real amounts of gene flow between local populations of Bythinella is high (Falniowski et al. 1998, 1999, Falniowski and Szarowska 2011), and there have been even sympatric occurrences of two species of Bythinella in the same spring (Radoman 1976, Falniowski et al. 2009b, Falniowski and Szarowska 2011).

In Bulgaria, several Bythinella species have been described solely on the basis of morphological characters, including only shell and penis: B. ravnogorica, Glöer and Georgiev 2009, B. rhodopensis, Glöer and Georgiev 2009, B. srednogorica, Glöer and Georgiev 2009, B. markovi, Glöer and Georgiev 2009, B. walkeri, Glöer and Georgiev 2009, (Glöer and Georgiev 2009), B. rilaensis, Glöer and Georgiev 2011, B. slaveyae, Glöer and Georgiev 2011, B. angelovi, Glöer and Georgiev 2011 (Glöer and Georgiev 2011), B. cf. opaca, von Gallenstein 1848 (Georgiev and Stoycheva 2008; later described as B. srednogorica), B. gloeeri, Georgiev 2009 (Georgiev 2009), B. hansboetersi, Glöer and Pešić 2006 (Glöer and Pešić 2006) and B. stoychevae, Georgiev 2011 (Georgiev 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, 2012). 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 (Glöer and 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.

Materials and methods

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.

Figure 1.

Sampling sites used in the present study (red dots) and in phylogenetic analyses (blue dots: Falniowski et al. 2009a; green dots: Falniowski et al 2012). Compare with Table 1. The dotted line indicates an area that was searched but where no Bythinella sites were found.

The sampling localities with their geographical coordinates, and the haplotypes for COI and ITS genes detected in each locality. Sequences from GenBank are also included. Compare with Figure 1.

ID Taxon Site Coordinates COI haplotypes ITS haplotypes
B1 B. sp. Bulgaria - Bezbog Peak, Pirin Mts. 41°45'20"N; 23°32'39"E HB1×3 HB1×3
B2 B. sp. Bulgaria - Ribarits village, Stara Planina Mts. 42°49'33"N; 24°22'23"E HB2A×3, HB2B HB2A,HB2B
B3 B. sp. Bulgaria - Leshnishki Waterfall, Belasitsa Mts. 41°22'12"N; 23°11'12"E HB3A, HB3B×3 HB3
B4 B. sp. Bulgaria - Panagyurski Kolonii, Sredna Gora Mts. 42°34'20"N; 24°12'34"E HB2A×3 HB2B×2
B5 B. sp. Bulgaria - Nestinarka beach, Strandzha Mts. 42°09'16"N; 27°51'21"E HB5×6 -
B6 B. sp. Bulgaria - N of Pnitsite area, Stara Planina Mts. 42°39'44"N; 24°58'40"E HB6×5 HB2A×2
B7 B. gloeeri Bulgaria - Lepenitza cave, West Rhodopes Mts. 41°57'14"N; 24°00'43"E HB7A, HB7B×5 HB7×5
B8 B. stoychevae Bulgaria - Manuilova Dupka Cave, Rhodopes Mts. 41°42'53"N; 23°46'58"E HB8×5 HB8×2
B9 B. ravnogorica Bulgaria - Ravnogor village, West Rhodopes Mts. 41°57'00"N; 24°21'53"E HB9A, HB9B×2 -
B10 B. sp. Bulgaria - Vodni Pech Cave 43°30'11"N; 22°46'55"E HB10A×2, HB10B, HB10C×3 HB10A×2, HB10B, HB10C×2
B11 B. angelovi Bulgaria - Koprivskitsa town, Sredna Gora Mts. 42°08'14"N; 24°21'55"E HB11 HB2A
B12 B. rhodopensis Bulgaria - Modarskata Cave, West Rhodopes Mts. 41°52'40"N; 24°33'40"E HB12×3 HB2A×2
B13 B. rilaensis Bulgaria - Rila Mts., near Belovo 42°08'15"N; 23°58'00"E HB7B×2 HB7×4
B14 B. dierkingi Bulgaria - Ravnogor village, West Rhodopes Mts. 41°56'59"N; 24°22'02"E HB11×3 HB2A×4
B15 B. sp. Bulgaria - Koprivshtitsa town, Sredna Gora Mts. 42°38'15"N; 24°21'57"E HB11×4 HB2A×2, HB2B
Falniowski et al. 2009b (FJ545011FJ545131)
1 B. viseuiana Romania - Vişeu River Valley 47°52'14"N; 24°11'23"E HR1A×3, HR1B×6 HR1×2
2 B. molcsanyi Romania - Igniş Mts., western slope of Firiza Lake 47°43'02"N; 23°36'29"E HR2A×2, HR2B×5 HR2A, HR2B
3 B. molcsanyi Romania - Igniş Mts., upstream of locality 2 47°45'58"N; 23°38'32"E HR2A×3, HR3A×5, HR3B×2 HR2B×2
4 B. molcsanyi Romania - Igniş Mts., near Izvoare Resort 47°45'14"N; 23°42'28"E HR3A×2, HR4A×4, HR4B×3 HR4A×2, HR4B
5 B. molcsanyi Romania - Igniş Mts., Izvoare Resort 47°44'51"N; 23°43'03"E HR2B×12, HR4B×8 HR5A, HR5B, HR5C,HR5D
6 B. radomani Romania - Bihor Mts., close to Vârtop Pass 46°31'25"N; 22°37'25"E HR6A×5, HR6B×3 HR6A×2, HR6B, HR6C, HR6D
7 B. radomani Romania - Bihor Mts., Iarba Rea village 46°25'35"N; 22°46'29"E HR7×7 HR7×3
8 B. dacica Romania - Retezat Mts., La Beci, Buta river valley 45°18'26"N; 22°56'12"E HR8A×6, HR8B×4, HR8C×3 HR8A, HR8B
9 B. dacica Romania - Retezat Mts., Râu Şes valley 45°19'25"N; 22°40'51"E HR9A×3, HR9B×4 HR9A, HR9B×2, HR9C
10 B. dacica Romania - Cerna Valley 45°00'33"N; 22°32'40"E HR10A×2, HR10B×4 HR10A, HR10B, HR10C
11 B. dacica Romania - Cerna Valley, 3.5 km up from locality 10 45°02'10"N; 22°34'06"E HR10A×6, HR11×6 HR10B×2
12 B. calimanica Romania - Cӑlimani Mts. 46°57'10"N; 25°04'07"E HR12A×3, HR12B×4, HR12C×6 -
Falniowski et al. 2009a (GQ152518GQ152544)
13 B. hansboetersi Bulgaria - Smoljan town, below Smoljanske Lake 41°37'01"N; 24°40'31"E HBU13A×2, HBU13B×2, HBU13C HB2A, H13A, HBU13B
14 B. hansboetersi Bulgaria - Smoljan town, near Amzovo 41°33'42"N; 24°41'41"E HB2A×2, HBU14A×3, HBU14B×2 HBU13B, HBU14
15 B. hansboetersi Bulgaria - Anton town, Bolovan Hill 42°44'48"N; 24°16'51"E HB2A×3, HBU15A×3, HBU15B, HBU15C HB2B×5, HBU15
16 B. hansboetersi Bulgaria - Mugla village 41°37'43"N; 24°31'08"E HBU16A×4, HBU16B×3 HBU16A×2, HBU16B
Falniowski et al 2012 (JQ639859JQ639883)
17 B. hansboetersi Bulgaria - Stara Planina, spring of Cherni Osam 42°43'21"N; 24°46'47"E HBU15B -
18 B. srednogorica Bulgaria - Sredna Gora Mts., S. of Streltcha town 42°27'16"N; 24°20'27"E HB2A×2, HBU18 -
19 B. rhodopensis Bulgaria - West Rhodopes Mts., S. of Lilkovo village 41°52'39"N; 24°33'21"E HBU19A, HBU19B, HBU19C -
20 B. slaveyae Bulgaria - Belasits Mts., S. of Belasitsa village 41°21'07"N; 23°09'19"E HB2A×2, HBU20 -
21 B. nonveilleri Serbia - Rtanj Mt., Vrmd a spring 43°42'00"N; 21°49'00"E HSE21A×2, HSE21B -
22 B. pesterica Serbia - Pester Plateau, Djerekare village 43°00'00"N; 20°08'00"E HSE22A, HSE22B, HSE22C -
23 B. taraensis Montenegro - canyon of the river Tara, stream Ljevok 42°59'29"N; 19°25'53"E HMO23A×2, HMO23B -
24 B. luteola Montenegro - National Park Biogradska Gora 42°53'31"N; 19°36'16"E HMO24 -
25 B. dispersa Montenegro - spring in Petnjik village 42°49'35"N; 19°54'10"E HMO25A, HMO25B, HMO25C -

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.

Data analysis

Sequences were edited in Bioedit 7.1.3.0 (Hall 1999) and aligned with the ClustalW program in MEGA 6 (Tamura et al. 2013). Single ITS sequences were assembled and aligned using CodonCodeAlligner 4.2.7 (CodonCodeCorporation, Dedham, MA). Basic sequence statistics, including haplotype polymorphism and nucleotide divergence, were calculated in DnaSP 5.10 (Librado and Rozas 2009). The saturation test of Xia et al. (2003) was performed using DAMBE (Xia 2013).

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 burn-in. 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 infer haplotype networks of the markers used, a median-joining calculation was implemented in NETWORK 4.6.1.1 (Bandelt et al. 1999).

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, 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.

Figure 2.

Shells of Bythinella. A locality B1 B–C locality B7 D–F locality B5 G–Q locality B10; bar equals 1 mm.

We obtained 58 new sequences of COI (552 bp, GenBank Accession numbers KT381098KT381155) and 36 new sequences of ITS-1 (234–264 bp, GenBank Accession numbers KT381156KT381191). 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, 2012), twelve sites in Romania (22 haplotypes: Falniowski et al. 2009b) and five sites in Montenegro (six haplotypes) and Serbia (five haplotypes: Falniowski et al. 2012). 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 34). 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. hansboetersi, B. ravnogorica, B. rhodopensis, B. rilaensis, B. slaveyae, and B. srednogorica, in fact based mostly on their locations.

Figure 3.

The maximum-likelihood phylogram for COI gene. Haplotypes obtained in present work are indicated in bold. Arrows and the letter C indicate cave haplotypes.

Figure 4.

Geographical distribution of COI clades. Compare with Figure 3.

Figure 5.

The median-joining haplotype network of COI haplotypes for clades I, II and III. Sequences from Falniowski et al. 2009a and Falniowski et al. 2012 are also included. Arrows and the letter C indicate cave haplotypes.

Mean distances within clades (italics) and p-distances between main COI clades of Bythinella.

clade_I clade_II clade_III Serbia clade_IV Serbia/Mont. clade_V Romania
0.008
clade_II 0.024 0.002
clade_III 0.026 0.031 -
Serbia 0.063 0.070 0.054 0.002
clade_IV 0.089 0.090 0.078 0.089 -
Serbia/Mont. 0.109 0.112 0.096 0.119 0.122 0.023
clade_V 0.091 0.090 0.084 0.090 0.118 0.126 0.004
Romania 0.112 0.106 0.097 0.120 0.122 0.118 0.109 0.079

Clades II and III were most closely related to Clade I differing by intercladal p-distances 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 divergence time 4.39 Mya) (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.

Figure 6.

The maximum-likelihood phylogram for the ITS-1 gene. Haplotypes obtained in present work are shown in bold. The COI clades are also shown.

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, 1998; Falniowski 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 (Szarowska and Falniowski 2008; 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, b).

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 (Giusti and 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. 1998, 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 (Falniowski et al. 2012).

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, 2006; Popescu 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 water-filled 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 1988, 2000; Peck and Finston 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 and Szarowska 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).

Acknowledgements

The study was supported by a grant from the National Science Centre (2012/05/B/NZ8/00407) to Magdalena Szarowska. We wish to thank Dr Helen McCombie-Boundry for improving the English and anonymous reviewers and the editor for valuable comments on the text.

References

  • Avise JC (2000) Phylogeography. The history and formation of species. Harvard University Press, Cambridge, MA and London.
  • Bandelt HJ, Forster P, Röhl A (1999) Median–joining networks for inferring intraspecific phylogenies. Molecular Biology and Evolution 16: 37–48. doi: 10.1093/oxfordjournals.molbev.a026036
  • Benke M, Brändle M, Albrecht C, Wilke T (2009) Pleistocene phylogeography and phylogenetic concordance in cold-adapted spring snails (Bythinella spp.). Molecular Ecology 18: 890–903. doi: 10.1111/j.1365-294X.2008.04073.x
  • Benke M, Brändle M, Albrecht C, Wilke T (2011) Patterns of freshwater biodiversity in Europe: lessons from the spring snail genus Bythinella. Journal of Biogeography 38: 2021–2032. doi: 10.1111/j.1365-2699.2011.02527.x
  • Bichain JM, Gaubert P, Samadi S, Boisselier-Dubayle MC (2007a) A gleam in the dark: Phylogenetic species delimitation in the confusing spring-snail genus Bythinella Moquin-Tandon, 1856 (Gastropoda: Rissooidea: Amnicolidae). Molecular Phylogenetics and Evolution 45: 927–941. doi: 10.1016/j.ympev.2007.07.018
  • Bichain JM, Boisselier Dubayle MC, Bouchet P, Samadi S (2007b) Species delimitation in the genus Bythinella (Mollusca: Caenogastropoda: Rissooidea): A first attempt combining molecular and morphoetrical data. Malacologia 49: 293–311. doi: 10.4002/0076-2997-49.2.293
  • Boeters HD (1973) Die Gattung Bythinella and Gattung Marstoniopsis in Westeuropa, 1. Westeuropäische Hydrobiidae, 4 (Prosobranchia). Proceedings of Fourth European Malacological Congress. Malacologia 14: 271–285.
  • Boeters HD (1998) Mollusca: Gastropoda: Superfamilie Rissooidea. In: Schwoerbel J, Zwick P (Eds) Süsswasserfauna von Mitteleuropa. Begründet von A. Brauer, 5/1–2, Gustav Fischer Verlag, Jena – Lübeck – Ulm.
  • Falniowski A (1987) Hydrobioidea of Poland (Prosobranchia: Gastropoda). Folia Malacologica 1: 1–122.
  • Falniowski A (1992) Shell outer and inner structure and rissoacean phylogeny. III. Truncatella subcylindrica (Linnaeus, 1767) (Prosobranchia: Rissoacea: Truncatellidae). Malakologische Abhandlungen, Staatliches Museum für Tierkunde Dresden 16: 7–11.
  • Falniowski A, Horsak M, Szarowska M (2009a) Bythinella hansboetersi Glöer et Pešić, 2006 (Gastropoda: Rissooidea) in Bulgaria: its morphology, molecular distinctness, and phylogeography. Folia Malacologica 17: 11–20. doi: 10.2478/v10125-009-0002-3
  • Falniowski A, Mazan K, Szarowska M (1999) Homozygote excess and gene flow in the spring snail Bythinella (Gastropoda: Prosobranchia). Journal of Zoological Systematics and Evolutionary Research 37: 165–175. doi: 10.1111/j.1439-0469.1999.tb00980.x
  • Falniowski A, Szarowska M, Fiałkowski W, Mazan K (1998) Unusual geographic pattern of interpopulation variation in a spring snail Bythinella (Gastropoda, Prosobranchia). Journal of Natural History 32: 605–616. doi: 10.1080/00222939800770311
  • Falniowski A, Szarowska M, Sirbu I (2009b) Bythinella Moquin-Tandon, 1856 (Gastropoda: Rissooidea: Bythinellidae) in Romania: species richness in a glacial refugium. Journal of Natural History 43: 2955–2973. doi: 10.1080/00222930903359636
  • Falniowski A, Szarowska M, Sirbu I (2009c) Bythinella Moquin-Tandon, 1856 (Gastropoda: Rissooidea: Bythinellidae) in Romania: its morphology with description of four new species. Folia Malacologica 17: 21–36. doi: 10.2478/v10125-009-0003-2
  • Falniowski A, Szarowska M (2011) Radiation and phylogeography in a spring snail Bythinella (Mollusca: Gastropoda: Rissooidea) in continental Greece. Annales Zoologici Fennici 48: 67–90. doi: 10.5735/086.048.0201
  • Falniowski A, Szarowska M (2012) Sequence-based species delimitation in the Balkan Bythinella Moquin-Tandon, 1856 (Gastropoda: Rissooidea) with general mixed Yule coalescent model. Folia Malacologica 20: 111–120. doi: 10.2478/v10125-012-0017-z
  • Falniowski A, Szarowska M, Sirbu I, Hillebrand A, Baciu M (2008) Heleobia dobrogica (Grossu & Negrea, 1989) (Gastropoda: Rissooidea: Cochliopidae) and the estimated time of its isolation in a continental analogue of hydrothermal vents. Molluscan Research 28: 165–170.
  • Falniowski A, Szarowska M, Glöer P, Pešić V, Georgiev D, Horsák M, Sirbu I (2012) Radiation in Bythinella Moquin-Tandon, 1856 (Mollusca: Gastropoda: Rissooidea) in the Balkans. Folia Malacologica 20: 1–10. doi: 10.2478/v10125-012-0006-2
  • Folmer O, Black M, Hoeh W, Lutz R, Vrijenhoek RC (1994) DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates. Molecular Marine Biology and Biotechnology 3: 294–299.
  • Georgiev D (2009) Bythinella gloeeri n. sp. – a new cave inhabiting species from Bulgaria (Gastropoda: Rissooidea: Hydrobiidae). Acta Zoologica Bulgarica 61: 223–227.
  • Georgiev D (2011) New species of snails (Mollusca: Gastropoda: Rissooidea) from cave waters of Bulgaria. Buletin Shkencor, ser. Shkencat Natyrore [Universiteti i Shkodrës “Luigi Gurakuqi”] 61: 83–96.
  • Georgiev D, Glöer P (2013) Identification key of the Rissooidea (Mollusca: Gastropoda) from Bulgaria with a description of six new species and one new genus. North-Western Journal of Zoology 9: 103–112.
  • Georgiev D, Stoycheva S (2008) A record of Bythinella cf. opaca (Gallenstein, 1848) (Gastropoda: Prosobranchia: Hydrobiidae) in Bulgaria. Malacologica Bohemoslovaca 7: 51–54.
  • Giusti F, Pezzoli E (1977) Primo contributo alla revisione del genere Bythinella in Italia. Natura Bresciana, Annurario del Museo Civico di Storia Naturale di Brescia 14: 3–80.
  • Giusti F, Pezzoli E (1980) Gasteropodi, 2 (Gastropoda: Prosobranchia; Hydrobioidea, Pyrguloidea). In: Collana del progetto finalizzato “Promozione della qualità dell’ambiente,” AQ/1/47. Guide per il riconoscimento delle specie animali delle acque interne Italiane 8. Consiglio Nazionale delle Ricerche, Verona, 1–67.
  • Glöer P, Georgiev D (2009) New Rissooidea from Bulgaria (Gastropoda: Rissooidea). Mollusca 27: 123–136.
  • Glöer P, Georgiev D (2011) Bulgaria, a hot spot of biodiversity (Gastropoda: Rissooidea)? Journal of Conchology 40: 1–16.
  • Glöer P, Pešić V (2006) Bythinella hansboetersi n. sp., a new species from Bulgaria. Heldia 6: 11–15.
  • Haase M, Wilke T, Mildner P (2007) Identifying species of Bythinella (Caenogastropoda: Rissooidea): A plea for an integrative approach. Zootaxa 1563: 1–16.
  • Holsinger JR (1988) Troglobites: the evolution of cave-dwelling organisms. American Scientist 76: 146–153.
  • Holsinger JR (2000) Ecological derivation, colonisation and speciation. In: Wilkens H, Culver D, Humphreys W (Eds) Subterranean Ecosystems. Elsevier, Amsterdam, 399–415.
  • Howarth FG (1973) The cavernicolous fauna of Hawaiian lava tubes. I. Introduction. Pacific Insects 15: 139–151.
  • Huelsenbeck JP, Rannala B (2004) Frequents properties of Bayesian posterior probabilities of phylogenetic trees under simple and complex substitution models. Systematic Biology 53: 904–913. doi: 10.1080/10635150490522629
  • Juan C, Guzik MT, Jaume D, Cooper SJB (2010) Evolution in caves: Darwin’s ‘wrecks of ancient life’ in the molecular era. Molecular Ecology 19: 3865–3880. doi: 10.1111/j.1365-294X.2010.04759.x
  • Librado P, Rozas J (2009) DnaSP v5 a software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25: 1451–1452. doi: 10.1093/bioinformatics/btp187
  • Mazan K (2000) Morphologica and allozymic variation within and between populations of Bythinella Moquin-Tandon, 1855 (Gastropoda: Prosobranchia). I. Morphological characters. Folia Malacologica 8: 107–139. doi: 10.12657/folmal.008.007
  • Mazan K, Szarowska M (2000a) Morphologica and allozymic variation within and between populations of Bythinella Moquin-Tandon, 1855 (Gastropoda: Prosobranchia). II. Phenetic analysis. Folia Malacologica 8: 189–213. doi: 10.12657/folmal.008.014
  • Mazan K, Szarowska M (2000b) Morphologica and allozymic variation within and between populations of Bythinella Moquin-Tandon, 1855 (Gastropoda: Prosobranchia). III. Phylogenetic analysis. Folia Malacologica 8: 257–269. doi: 10.12657/folmal.008.022
  • Miller MA, Pfeiffer W, Schwartz T (2010) Creating the CIPRES Science Gateway for inference of large phylogenetic trees. Proceedings of the Gateway Computing Environments Workshop (GCE), 14 Nov., New Orleans, LA, 1–8. doi: 10.1109/gce.2010.5676129
  • Nei M, Kumar S (2000) Molecular Evolution and Phylogenetics. Oxford University press, Oxford, UK and New York.
  • Peck SB, Finston TL (1993) Galapagos islands troglobites: the questions of tropical troglobites, parapatric distributions with eyed-sister-species, and their origin by parapatric speciation. Memoires de Biospieliologie 20: 19–37.
  • Popescu S-M, Dalesme F, Jouannic G, Escarguel G, Head M, Melinte-Dobrinescu M-C, Süto-Szentai M, Bakrac K, Clauzon G, Suc JP (2009) Galeacysta etrusca complex: dinoflagellate cyst marker of Paratethyan influxes to the Mediterranean Sea before and after the peak of the Messinian Salinity Crisis. Palynology 33: 105–134. doi: 10.1080/01916122.2009.9989688
  • Popov SV, Rögl F, Rozanov AY, Steininger FF, Shcherba IG, Kovac M (2004) Lithological-paleogeographic maps of the Paratethys. 10 maps Late Eocene to Pliocene. Courier Forschungsinstitut Senckenberg 250: 1–46.
  • Popov SV, Shcherba IG, Ilyina LB, Nevesskaya LA, Paramonova NP, Khondkarian SO, Magyar I (2006) Late Miocene to Pliocene palaeogeography of the Paratethys and its relation to the Mediterranean. Palaeogeography, Palaeoclimatology, Palaeoecology 238: 91–106. doi: 10.1016/j.palaeo.2006.03.020
  • Radoman P (1976) Speciation within the family Bythinellidae on the Balkans and Asia Minor. Zeitschrift für Zoologische Systematik und Evolutionsforschung 14: 130–152. doi: 10.1111/j.1439-0469.1976.tb00522.x
  • Rivera MAJ, Howarth FG, Taiti S, Roderick GK (2002) Evolution in Hawaiian cave isopods (Oniscidea: Philosciidae): vicariant speciation or adaptive shifts? Molecular Phylogenetics and Evolution 25: 1–9. doi: 10.1016/S1055-7903(02)00353-6
  • Ronquist F, Teslenko M, van der Mark P, Ayres D, Darling A, Hohna S, Larget B, Liu L, Suchard MA, Huelsenbeck JP (2012) MrBayes 3.2: Efficient Bayesian phylogenetic inference and model choice across a large model space. Systematic Biology 61: 539–542. doi: 10.1093/sysbio/sys029
  • Sanderson MJ (1997) A nonparametric approach to estimating divergence times in the absence of rate constancy. Molecular Biology and Evolution 14: 1218–1231. doi: 10.1093/oxfordjournals.molbev.a025731
  • Sanderson MJ (2003) R8s: inferring absolute rates of molecular evolution, divergence times in the absence of a molecular clock. Bioinformatics 19: 301–302. doi: 10.1093/bioinformatics/19.2.301
  • Suc J-P, Do Couto D, Carmen Melinte-Dobrinescu M, Macaleţ R, Quillévéré F, Clauzon G, Csato I, Rubino J-L, Popescu S-M (2011) The Messinian Salinity Crisis in the Dacic Basin (SW Romania) and early Zanclean Mediterraneannerranea Paratethys high sea-level connection. Palaeogeography, Palaeoclimatology, Palaeoecology 310: 256–272. doi: 10.1016/j.palaeo.2011.07.018
  • Stamatakis A (2014) RAxML Version 8: A tool for Phylogenetic Analysis and Post-Analysis of Large Phylogenies. Bioinformatics 30: 1312–1313. doi: 10.1093/bioinformatics/btu033
  • Szarowska M (2000) Environmental stress and stability of Bythinella populations in South Poland (Gastropoda: Prosobranchia: Hydrobioidea). Malakologische Abhandlungen 20: 93–98.
  • Szarowska M, Falniowski A (2008) There is no philosopher’s stone: coup de grace for the morphology-based systematics in the rissooidean gastropods? 5th Congress of the European Maacological Societies, Ponta Delgada, 28.
  • Szarowska M, Osikowski A, Hofman S, Falniowski A (2015) Do diversity of the spring-inhabiting snail Bythinella (Gastropoda, Bythinellidae) on the Aegean Islands reflect geological history? Hydrobiologia. doi: 10.1007/s10750-015-2415-x
  • Tajima F (1993) Simple methods for testing molecular clock hypothesis. Genetics 135: 599–607.
  • Tamura K, Stecher G, Peterson D, Filipski A, Kumar S (2013) MEGA6: Molecular evolutionary genetics analysis version 6.0. Molecular Biology and Evolution 30: 2725–2729. doi: 10.1093/molbev/mst197
  • Wilke T, Davis GM (2000) Infraspecific mitochondrial sequence diversity in Hydrobia ulvae and Hydrobia ventrosa (Hydrobiidae: Rissoacea: Gastropoda): Do their different life histories affect biogeographic patterns and gene flow? Biological Journal of the Linnean Society 70: 89–105. doi: 10.1111/j.1095-8312.2000.tb00202.x
  • Xia X (2013) DAMBE: A comprehensive software package for data analysis in molecular biology and evolution. Molecular Biology and Evolution 30: 1720–1728. doi: 10.1093/molbev/mst064
  • Xia X, Xie Z, Salemi M, Chen L, Wang Y (2003) An index of substitution saturation and its application. Molecular Phylogenetics and Evolution 26: 1–7. doi: 10.1016/S1055-7903(02)00326-3
  • Zagorchev I (2007) Late Cenozoic development of the Strouma and Mesta fluviolacustrine systems, SW Bulgaria and northern Greece. Quaternary Science Reviews 26: 2783–2800. doi: 10.1016/j.quascirev.2007.07.017