Phylogeography of the land snail genus Circassina (Gastropoda: Hygromiidae) implies multiple Pleistocene refugia in the western Caucasus region

Molecular Phylogenetics and Evolution 93 (2015) 129–142

Contents lists available at ScienceDirect

Molecular Phylogenetics and Evolution journal homepage: www.elsevier.com/locate/ympev

Phylogeography of the land snail genus Circassina (Gastropoda: Hygromiidae) implies multiple Pleistocene refugia in the western Caucasus region q Marco T. Neiber, Bernhard Hausdorf ⇑ Centre for Natural History, Zoological Museum, University of Hamburg, Martin-Luther-King-Platz 3, 20146 Hamburg, Germany

a r t i c l e

i n f o

Article history: Received 4 March 2015 Revised 16 July 2015 Accepted 18 July 2015 Available online 26 July 2015 Keywords: Caucasus Colchis Dart apparatus Historical biogeography Phylogeography Pleistocene refugia

a b s t r a c t The phylogeny and historical biogeography of the Caucasian land snail genus Circassina was reconstructed using multilocus amplified fragment length polymorphism (AFLP) data and mitochondrial DNA sequences. Diversification within the group started with a divergence of populations from the western Lesser Caucasus from those of the Greater Caucasus during the late Miocene. Distinct AFLP clusters and major mitochondrial clades separated by long internal branches lend evidence to the hypothesis of separate glacial refuges in the Lesser and Greater Caucasus during the Pleistocene. High genetic distances across low geographic distances and admixture analysis revealed a phylogeographic boundary running through the Colchis lowlands, which may have been established and maintained in part by repeated transgressions of the Black Sea during the Pleistocene and Holocene. Localities in Ciscaucasia were probably colonised through long-distance dispersal across the main ridge of the Greater Caucasus. The phylogeny implies multiple independent losses of accessory genital organs, i.e. dart sac and mucus glands, within Circassina. None of the anatomically defined (sub-) species distinguished so far is monophyletic and there is gene flow between the two main population groups across the Colchis lowlands. Thus, we propose to classify these population groups as subspecies of a single species. Ó 2015 Elsevier Inc. All rights reserved.

1. Introduction The Caucasus region is ranked among the 25 most important global biodiversity hotspots (Myers et al., 2000; Zazanashvili et al., 2004). The importance of the western Caucasus, especially the Colchis region, as a glacial refugium where, among others, Neogene relict species survived as well as a centre of ongoing radiation has increasingly been appreciated (Hewitt, 2000; Pokryszko et al., 2011; Tarkhnishvili et al., 2012; Nakhutsrishvili, 2013; Tarkhnishvili, 2014; Walther et al., 2014). There is dissent, however, on whether there existed a single continuous forest refugium at the eastern Black Sea coast as implied by some studies (van Andel and Tzedakis, 1996; Kikvidze and Ohsawa, 2001; Pokryszko et al., 2011; Tarkhnishvili et al., 2012; Volkova et al., 2013; Wielstra et al., 2013) or multiple forest refugia (Tarkhnishvili et al., 2000; Mumladze et al., 2013; Tarkhnishvili, 2014).

q

This paper was edited by the Associate Editor Jan Strugnell.

⇑ Corresponding author.

E-mail address: [email protected] (B. Hausdorf). http://dx.doi.org/10.1016/j.ympev.2015.07.012 1055-7903/Ó 2015 Elsevier Inc. All rights reserved.

Circassina Hesse, 1921 (Gastropoda: Hygromiidae) is a group of land snails endemic to the western and central Caucasus region and the eastern Pontus. It includes only two species, C. frutis (Pfeiffer, 1859) and C. lasistana Hausdorf, 2001. Abchasohela Hudec & Lezhawa, 1971, which was formerly classified as a subgenus of Circassina (Hausdorf, 2001), proved to be more closely related to other hygromiid groups than to Circassina (Walther et al., in press). The Circassina frutis complex is exceptional among helicoid land snails because it is polymorphic with regard to the dart apparatus, i.e. the snails possess either a complete dart apparatus with a dart sac plus an accessory sac and mucus glands (C. frutis circassica (Mousson, 1863)), only mucus glands (C. frutis frutis) or none of these accessory genital organs (C. frutis veselyi (Frankenberger, 1919)). Schileyko (1978) classified these morphotypes as subspecies, whereas Giusti and Manganelli (1987) questioned whether they are related at all. Hausdorf (2001) followed Schileyko (1978), but questioned whether the morphotypes are actually genetic entities, because their ranges interdigitate in western Georgia. Furthermore, Hausdorf (2001) separated C. lasistana from C. frutis because of the much shorter spermatophore-forming flagellum.

130

M.T. Neiber, B. Hausdorf / Molecular Phylogenetics and Evolution 93 (2015) 129–142

A dart apparatus with a calcareous dart and mucus glands is an apomorphy of a large group of helicoid land snails (Nordsieck, 1987; Hausdorf, 1998). The mucus glands produce an allohormone-like substance that is transmitted into the mating partner during courtship by the dart. This substance causes temporal changes in the female part of the genitalia of the mating partner, which delays the intake of the spermatophore into the bursa of the bursa copulatrix and allows more sperm to leave the spermatophore before it is digested in this organ, so that the fertilisation success of the dart shooter is increased (Koene and Chase, 1998; Chase and Blanchard, 2006; Chase, 2007). Since the fertilisation success of individuals that possess these organs is increased, sexual selection should favour the fixation of these traits. However, injuries or an increased infection risk of the recipient as well as a prolongation of the mating implied by the delayed intake of the spermatophore, during which the individuals are under risk of desiccation or predation, may have negative fitness effects. Thus, natural selection might counteract sexual selection. The population structure of land snails usually consists of local demes that often go extinct and may be founded again, potentially resulting in populations in which the lack of the dart sac or the complete dart apparatus is fixed. The purpose of this study was twofold. First, we studied the phylogeography of Circassina based on mitochondrial sequences and amplified fragment length polymorphism (AFLP) markers to infer whether the Colchis region can be viewed as a homogeneous Pleistocene forest refugium or whether there were multiple separate refugia in this region. Second, we tested for associations between environmental variables and the structure of the dart apparatus to infer potential reasons for the unusual polymorphism of this accessory genital organ in Circassina. We also reassessed the

taxonomy of Circassina, so far resting upon the anatomy of the genitalia, based on the results of our molecular phylogenetic and admixture analyses. 2. Material and methods 2.1. Sampling Mitochondrial cytochrome c oxidase subunit I gene (cox1) and 16S rDNA fragments as well as amplified fragment length polymorphism (AFLP) data of C. frutis frutis, C. frutis circassica, C. frutis veselyi and C. lasistana covering most of the ranges of these taxa in Georgia, the Russian part of the Caucasus region and north-eastern Turkey were analysed (Fig. 1; Supplementary Table S1). Individuals of Caucasigena eichwaldi (Pfeiffer, 1846), Fruticocampylaea narzanensis (Krynicki, 1836), and Kokotschashvilia holotricha (Boettger, 1884) were included as outgroups in the analyses. Circassina specimens were identified according to the diagnoses given in Hausdorf (2001), i.e. dart sac, accessory sac and mucus glands missing (C. frutis veselyi), mucus glands present, but dart and accessory sacs missing (C. frutis frutis), dart sac plus accessory sacs as well as mucus glands present, but epiphallus shorter than flagellum (C. frutis circassica) or epiphallus longer than flagellum (C. lasistana) (Fig. 2). Data on sampling sites, voucher numbers and the classification of each specimen used in this study are compiled in Supplementary Table S1. 2.2. DNA amplification and sequencing

at

Samples of foot muscle tissue were stored in 100% isopropanol 20 °C. DNA was extracted following a slightly modified version

D L1

Rus

D L2

Grea

ter C

auca

sus

D

sia

L3 D L4

L6

C

D

L5

L10 L19 L11 D L24 D D L12 L20 D D L17 D D L21 D L13 D L16 D L22 D L18 B D D B L14 L15 B L23 D B L36 L25 D L26 L27 C B L38 L46 L44 L28 L47 L35 L37 L43 A A C C L31 A B A L29 B L39 L34 C L48 A A L45 L33 C C A L51 L32 A L52 A A A A A A L49 L50 L42 L41 L40 A L54 A L53 A L55

Black Sea

Colchis Lowlands

Geo

C

rgia

L7 C L8 D L9

C

L30

Les

ser

A A L60

0

A

L59

Cau

cas

L58

A

L57 A L56

Turkey

us

Armenia

Az

erb

aija

n

200 km

Fig. 1. Sampling sites (see Supplementary Table S1) and distribution of mitochondrial haplotype clades (for letters in the circles; see Fig. 2) and of morphotypes (red: Circassina frutis fruits; blue: C. frutis circassica; green: C. frutis veselyi; yellow: C. lasistana) in the Caucasus region and north-eastern Turkey. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

131

M.T. Neiber, B. Hausdorf / Molecular Phylogenetics and Evolution 93 (2015) 129–142

accessory sac

vagina

dart sac atrium

penis

C. frutis circassica C. lasistana

C. frutis frutis

C. frutis veselyi

Circassina

1528, L49 1730, L44 19/0.72 1732, L46 1738, L45 39/0.5 2205, L57 65/1 1512, L51 11/0.56 1735, L48 76/0.99 1729, L51 1531, L47 27/79 99/1 1529, L48 2195, L52 2198, L53 84/1 2199, L59 Clade A 2207, L60 51/0.52 2400, L58 67/0.98 49/0.6 2208, L55 95/1 2200, L56 1524, L43 57/0.82 1525, L33 100/1 1743, L43 100/1 1569, L50 88/0.99 2197, L54 100/1 2196, L54 1542, L42 100/1 1571, L40 1733, L41 66/0.8 1523, L37 98/1 1568, L38 1541, L39 72/1 1527, L36 48/0.73 Clade B 1567, L14 100/1 89/1 1530, L15 1526, L27 1510, L31 64/0.93 94/1 1516, L29 1570, L30 71/0.73 1519, L8 99/1 1741, L7 62/0.95 1540, L35 Clade C - 1517, L34 100/1 1539, L28 96/1 1739, L34 100/1 1511, L32 1538, L28 100/1 1731, L6 100/1 1742, L4 1737, L5 99/1 1533, L18 1513, L17 92/1 94/1 1508, L22 97/1 100/1 1740, L22 97/1 1515, L21 1514, L20 97/1 1661, L2 1665, L1 98/1 1593, L11 68/0.99 1520, L11 1518, L10 Clade D 63/0.97 1522, L9 1521, L13 46/0.8 1592, L13 95/1 1566, L12 87/1 1591, L13 1736, L3 95/1 1537, L16 98/1 1672, L16 78/0.99 1509, L24 85/0.98 90/1 1534, L23 79/1 1535, L26 1536, L25 60/0.84 1734, L19 F. narzanensis 1719 C. eichwaldi 1783 K. holotricha 1616 46/0.87

mucus glands

76/1

100/1

0.05

Outgroups

81/0.99

Fig. 2. Bayesian 50% majority consensus tree based on mitochondrial cox1 and 16S rDNA sequences. Numbers at the nodes refer to bootstrap support values (left) from the maximum likelihood analysis and posterior probability values (right) obtained from the Bayesian analysis. Pie charts at the tips indicate individual proportions of parentage in the two clusters delimited in the STRUCTURE analysis of AFLP data, the numbers refer to vouchers as detailed in Supplementary Table S1 and are coloured according to the (sub-) species sensu Hausdorf (2001) (see inset with the schematic representations of the distal genitalia). The L-numbers refer to localities (see Supplementary Table S1). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

132

M.T. Neiber, B. Hausdorf / Molecular Phylogenetics and Evolution 93 (2015) 129–142

of the protocol described by Sokolov (2000) as described in Scheel and Hausdorf (2012). Fragments of cox1 and the 16S rDNA were amplified by polymerase chain reaction (PCR) using the primer pairs LCO1490 plus HCO2198 (Folmer et al., 1994) and 16Scs1 plus 16Scs2 (Chiba, 1999), respectively (Supplementary Table S2). Amplifications were performed in 25 ll volumes containing 2 ll 10  amplification buffer B (biolabproducts), 4 ll MgCl2 (25 mM, biolabproducts), 1 ll dNTP mix (5 mM each, biolabproducts), 1 ll of each primer (10 lM), 0.2 ll Crystal Taq DNA polymerase (biolabproducts), 1 ll template DNA and 14.8 ll ddH2O under the following reaction conditions: an initial denaturation step at 94 °C for 2 min, 35 PCR cycles (94 °C for 30 s, 50 °C for 30 s, 72 °C for 30 s) and a final extension step at 72 °C for 5 min. Prior to sequencing, PCR products were enzymatically cleaned up by adding 0.65 ll FastAP thermosensitive alkaline phosphatase (Thermo Fisher Scientific) and 0.35 ll exonuclease I (Thermo Fisher Scientific) to a 5 ll volume of undiluted PCR mixture followed by an incubation step at 37 °C for 15 min. The enzymes were inactivated at 85 °C for 15 min. Both strands of the amplified products were sequenced at Macrogen Europe Laboratory (Amsterdam, The Netherlands).

consensus tree with average branch lengths and posterior probabilities as node support values were calculated in MrBayes. Heuristic Maximum likelihood (ML) analyses were performed with the Garli web service (Bazinet et al., 2014; Zwickl, 2006) allowing for independent estimation of parameters for individual partitions using best-fit models as suggested by Partitionfinder (GTR + I + G for first and second codon positions and TIM + G for third codon positions of cox1, TIM + I + G for 16S rDNA). Confidence values were computed by bootstrapping with 1000 replications. Bootstrap support values from the ML analysis and posterior probabilities from the BI analysis were mapped on the BI 50% majority consensus tree with Sumtrees 3.3.1, which is part of the Dendropy 3.8.0 package (Sukumaran and Holder, 2010). Bootstrap (BS) support values P70% and posterior probabilities (PP) P 0.95 were regarded as significant statistical support. GenBank accession numbers are listed in the Supplementary Table S1. Alignment files and phylogenetic trees are available from TreeBASE (study accession number 17711). 2.6. AFLP data scoring

2.3. Sequence analysis ChromasPro 1.7.1 (Technelysium) was used to assemble forward and reverse sequence reads. The sequences were aligned with MAFFT 7 (Katoh and Standley, 2013) using the Q-INS-i algorithm, the 1PAM/j = 2 option for the scoring matrix for nucleotide sequences and otherwise default settings. Ambiguously aligned positions in the 16S rDNA data set were removed with Gblocks (Castresana, 2000) allowing for smaller final blocks, gap positions within the final blocks and less strict flanking positions. To test if the cox1 and 16S rDNA sequences were evolving neutrally, we computed Tajima’s D and Fu’s Fs statistic (Fu, 1997; Tajima, 1989) in DnaSP version 5.10.01 (Librado and Rozas, 2009). 2.4. AFLP data AFLP data were generated using the protocol of Vos et al. (1995) with slight modifications as specified by Scheel and Hausdorf (2012). The used adapters and primers are compiled in Supplementary Table S2. 2.5. Phylogenetic analyses of mitochondrial DNA sequence data The concatenated data set of the cox1 and 16S rDNA fragments was initially divided into four partitions, the first, second and third codon positions of cox1 and 16S rDNA. An exhaustive search with Partitionfinder 1.1.1 (Lanfear et al., 2012) was conducted to select an appropriate partitioning scheme and evolutionary models based on the Akaike information criterion allowing for separate estimation of branch lengths for each partition. As best-fit partitioning scheme, the Partitionfinder analysis suggested to divide the data set into three partitions, first plus second codon positions of cox1, third codon positions of cox1 and 16S rDNA. Bayesian inference (BI) was performed using MrBayes 3.2.2 (Ronquist et al., 2012). Metropolis-coupled Monte Carlo Markov chain (MC3) searches in MrBayes were run with four chains in two separate runs for 20,000,000 generations with default priors, trees sampled every 1000 generations and separate estimation of parameters for individual partitions under default heating using best-fit models available in MrBayes as suggested by Partitionfinder (GTR + I + G for first and second codon positions and GTR + G for third codon positions of cox1, GTR + I + G for 16S rDNA). Diagnostic tools in MrBayes were used to ensure that the MC3 searches had reached stationarity and convergence. The first 2,000,000 generations were discarded as burn-in. A majority rule

Detections and size calculations of AFLP bands were performed with Peak Scanner 1.0 (Applied Biosystems) using default parameters except for allowing a light peak smoothing. The add-on package RawGeno (Arrigo et al., 2009) for the statistical software suite R (R Development Core Team, 2013) was used for automated binning and scoring of AFLP bands. Scoring range was set to 50–400 base pairs (bp), minimum intensity to 100 relative fluorescence units (rfu), minimum bin width to 1 bp, and maximum bin width to 1.5 bp. Holland et al. (2008) found that stringent settings for the reproducibility of bins, although increasing data quality, did not give the best phylogenetic resolution, as too many useful characters were discarded. Taking these results into account, we used somewhat relaxed conditions by only eliminating bins with a reproducibility of less than 80% (the default setting in RawGeno). A replicate reproducibility rate of 94% (n = 15) was calculated as the mean percentage of matching character states between replicates of the same individual (Pompanon et al., 2005). 2.7. AFLP data analyses To assess the amount of admixture among individuals, an admixture analysis of the AFLP data was performed with STRUCTURE version 2.3.4 (Pritchard et al., 2000; Falush et al., 2007). In order to determine the number (K) of genetically homogeneous groups of individuals, we carried out 20 runs with 800,000 iterations after a burn-in period of 200,000 iterations for K = 1 to K = 10 in STRUCTURE choosing the admixture model and allowing for correlated allele frequencies between populations. The resulting estimates of the posterior probabilities L(K) of the data for each K and each run were used to compute the ad hoc quantity DK proposed by Evanno et al. (2005). The value of K, for which DK reached its global maximum was assumed to represent an appropriate estimate of the number of clusters. The run with the highest L(K) values among the replicates was used for graphical display with Distruct version 1.1 (Rosenberg, 2004) and further analyses. We used the R package HIest (Fitzpatrick, 2012) to test whether individuals identified as hybrids in the STRUCTURE analysis (i.e. ancestry of less than 90% in one of the clusters detected) can be assigned to classes, i.e. F1, F2 or backcrosses to parental lineages, or represent hybrids that are at least several generations old. To determine if a classification is credible the log-likelihood of the best fit class was compared to the second best and accepted if it was two units greater and within two units of the maximum log-likelihood (Fitzpatrick, 2012). Only AFLP loci whose

M.T. Neiber, B. Hausdorf / Molecular Phylogenetics and Evolution 93 (2015) 129–142

frequencies in the parental taxa differed absolutely by at least 0.333 were used in the analysis. A non-metric multidimensional scaling (Kruskal, 1964) of the Jaccard distances based on the AFLP data set was performed with the program package PRABCLUS (Hennig and Hausdorf, 2013), an add-on package for the software suite R (R Development Core Team, 2013). A neighbour-net (Bryant and Moulton, 2004) was constructed with SplitsTree4 version 4.13.1 (Huson and Bryant, 2006) using Jaccard distances based on the AFLP data set containing all Circassina individuals. Multilocus genomic data such as AFLPs usually contain information on the genealogy of the group under study as well as information on the amount of gene flow between populations. The latter effect may severely affect the reconstruction of the phylogeny from these data using tree-building algorithms. Individuals with high levels of parentage in more than one genetic lineage or cluster are usually resolved in intermediate positions between the parental lineages or even basally and typically have a decreasing effect on the statistical support values for the clades that contain the parental lineages (Seehausen, 2004; Reeves and Richards, 2007). To reduce the effects of gene flow on the reconstruction of the phylogeny of Circassina with AFLP data, we only included individuals with a parentage of more than 90% parentage in one cluster based on individual assignments in the STRUCTURE analysis for the DK maximum. The phylogenetic analysis was performed with a restriction site model (lset coding = noabsencesites) as implemented in MrBayes. MC3 searches were run with four chains in two separate runs for 50,000,000 generations with default priors, trees sampled every 5000 generations. Diagnostic tools in MrBayes were used to ensure that the MC3 searches had reached stationarity and convergence. The first 5,000,000 generations were discarded as burn-in. A majority rule consensus tree with average branch lengths and posterior probabilities as node support values were calculated in MrBayes. Alignment files and phylogenetic trees are available from TreeBASE (study accession number 17711). 2.8. Inference of the time frame We dated the divergences within Circassina using the Bayesian algorithm implemented in Beast 2.1.3 (Bouckaert et al., 2014) based on the concatenated cox1 and 16S rDNA data set with the same partitioning scheme and nucleotide substitution models as in the MrBayes analysis described above. Molecular clock tests were performed for each data partition independently in MEGA 6 (Tamura et al., 2013) using the ML tree from the Garli analysis for comparisons. The null hypotheses of equal evolutionary rates throughout the tree were rejected at a 0.05 significance level for each partition assuming the same models of nucleotide substitution as in the MrBayes analysis. Thus, a linked uncorrelated relaxed lognormal molecular clock was assumed for the Beast analysis. As tree prior the Yule speciation model was chosen. The Beast analysis was performed with ten independent runs for 100,000,000 generations each with a sampling frequency set to 50,000. Tracer v1.6 was used to assess convergence of runs, which were combined with LogCombiner 2.1.4 included in the Beast package discarding 10% of generations as burn-in. A maximum clade credibility tree with median node heights was calculated with TreeAnnotator 2.1.3. Dated fossils assignable to the clade of Caucasian Hygromiidae that includes Circassina are extremely scarce. Steklov (1966) described two species from the Caucasus region, which were assigned to Caucasigena by Schileyko (1978). One of these species, C. fortangensis was discovered in Maeotian deposits (8.3–6.05 Ma before present, see Gradstein et al., 2012) of the Fortanga river valley in Chechnya (Russia). Considering the

133

difficulty to assign fossil shells of Caucasian hygromiids to a specific genus, we deemed it more appropriate to consider the minimum age of the fossil described by Steklov (1966) of 6.05 Ma before present as a minimum age for the root of the tree in our analysis rather than assigning it to a specific genus. We further assumed that the radiation to which Circassina belongs originated on the Greater Caucasus land mass since the highest diversity of the group is found in that area (Walther et al., in press). Taking this into account and the possibility that the age of the root may be underestimated by using the minimum age of the fossil of 6.05 Ma before present as calibration, if it was actually representing an extant genus, we applied a lognormal distributed prior for the root of the tree choosing the mode of the distribution so that its 95% quantile (standard deviation at its default value) corresponded to the Eocene/Oligocene boundary (33.9 Ma before present), the approximate emergence of the land mass of the present day Greater Caucasus from the closing Tethys Sea (Popov et al., 2004; Vincent et al., 2007). We would like to stress that the datings obtained in this manner, ought to be regarded as rough estimates of the time frame rather than accurate point estimates. 2.9. Isolation by distance vs. genetic clustering In order to assess whether the phylogeographic structure of Circassina can be better explained by isolation by distance or genetic clustering we used the distanced-based redundancy analysis (dbRDA) implemented in DISTLM v.5 (McArdle and Anderson, 2001; Anderson, 2004). Pairwise Jaccard distances based on the AFLP data were calculated with SplitsTree4 and GTR + G distances (shape parameter: 0.157) based on the mitochondrial cox1 and 16S rDNA data with PAUP⁄ version 4.0 (Swofford, 2002). Geographical latitude and longitude were used as spatial variables. The genetic cluster affiliation of each individual was modelled as the proportion of parentage obtained from the STRUCTURE analysis based on the AFLP data. The affiliation of each individual to a main mitochondrial clade (as distinguished in the phylogenetic analyses) was modelled with a set of n 1 dummy variables, where n was the number of main mitochondrial clades. The affiliation to the nth clade was excluded as redundant information. The amount of the variation explained by isolation by distance was then inferred from the genetic distance matrices with dbRDA using the spatial variables as predictor variables. To control for possible effects of genetic clustering, analyses were also performed with genetic cluster or mitochondrial clade affiliations as covariables. The other way around, the amount of variation explained by genetic clustering was inferred from the genetic distance matrices using genetic cluster/mitochondrial clade affiliations as predictor variables. To control for possible effects of isolation by distance, analyses were repeated with spatial variables as covariables. Pseudo-p values were estimated with 9999 random permutations of the predictor variables. We also plotted geographical distances versus genetic distances (Jaccard distances for AFLP data and GTR + G (shape parameter: 0.2195) for mitochondrial sequences) and used a Mantel test with 9999 random permutations as implemented in the add-on package ade4 (Dray and Dufour, 2007) for the software suite R (R Development Core Team, 2013) to test for correlation between geographic and genetic distances applying a slight jittering (if necessary) to the matrices to ensure full rank. 2.10. Genetic barriers and long-distance dispersal events We used the method described by Sauer et al. (2013) to identify barriers to gene flow and long-distance dispersal events. Pairwise geographic distances were calculated in Mathematica 9.0.1 (Wolfram Research, Champaign) based on the geographic

134

M.T. Neiber, B. Hausdorf / Molecular Phylogenetics and Evolution 93 (2015) 129–142

coordinates and used along with Jaccard distances for AFLP data and GTR + G distances (shape parameter: 0.2195) for mitochondrial sequences, respectively, to determine pairs of individuals that belong to both the quintile of the highest genetic distances and the quintile of lowest geographic distances (indicating barriers to gene flow) as well as pairs of individuals that belong to both the quintile of lowest genetic distances and the quintile of largest geographic distances (indicating long-distance dispersal events). Results were visualized by plotting lines between pairs of individuals on a map of the Caucasus region. 2.11. Associations between traits of the dart apparatus and environment variables We used an approach based on logistic regressions to model the probability of presence/absence of the dart sac and presence/absence of mucus glands given the environmental conditions of the sampling localities implemented in Sambada (Stucki, 2014). Although Sambada was written to work with genomic data, the implemented logistic regression algorithm is appropriate for any nominal variables (presence/absence data) and measurement variables as independent variables (McDonald, 2014). Geographical coordinates and altitudinal data were taken at the sampling sites with a GPS device. For each sampling site, 19 biologically meaningful climatic variables (BIOCLIM data set; Nix, 1986) that represent annual trends, seasonality and extreme or limiting environmental factors as well as minimal, mean and maximal monthly temperatures and monthly precipitation data were acquired from the WorldClim database (http://www.worldclim.org; 1 km2 resolution; Hijmans et al., 2005; Supplementary Table S1). Hence, including geographical latitude and longitude, a total of 70 environmental variables were available for each locality. When testing for associations between individual morphotypes and their local environment, it is important to correct for spatial population structure to avoid bias deriving from the overlap between the spatial component in the genetic and environmental data (Meirmans, 2012). Therefore, we included in addition to the environmental variables the proportions of parentage from the STRUCTURE analyses as predictor variables in the Sambada analyses. We only considered associations between morphotype and environment that were not also associated with genotype (i.e. with the proportions of parentage from the STRUCTURE analyses). This approach decreases the risk of erroneously identifying an association of a morphotype with an environmental parameter as significant, when this association may actually be the consequence of population structure. The statistical significance threshold in Sambada was set to 0.01 before applying Bonferroni correction. Global and local indicators of spatial autocorrelation provided in Sambada were computed with a weighting scheme based on the 20 nearest neighbours using 1000 permutations to determine pseudo-p values as described in Stucki (2014). 3. Results 3.1. Mitochondrial phylogeny The aligned cox1 sequences had a length of 655 base pairs (bp). The 16S rDNA alignment had a length of 925 bp, of which 849 remained after exclusion of ambiguously aligned positions with Gblocks. Tajima’s D and Fu’s Fs test for cox1 and 16S rDNA sequences did not reject the null hypothesis of neutral evolution (p > 0.1). Both, BI and ML analyses resulted in mostly well-resolved phylogenies recovering Circassina as a monophyletic group with maximal statistical support (Fig. 2). Four major mitochondrial lineages, hereafter referred to as clades A–D (Fig. 2),

received significant statistical support (BS: 72–100, PP: 1.00), however none of the morphologically defined taxa C. frutis frutis, C. frutis circassica, C. frutis veselyi and C. lasistana (briefly referred to as frutis, circassica, veselyi or lasistana morphotypes) were recovered as monophyletic groups in our analyses. Clade D was consistently placed in a significantly supported sister group relationship with the clade comprised of the other major mitochondrial lineages A–C (Fig. 2), albeit the relationships between the latter were not resolved. For the geographic distribution of morphotypes and mitochondrial lineages, see Fig. 1. Clade A included individuals of the frutis morphotype from the Pontus Mountains in north-eastern Turkey (L55–L60) and adjacent south-western Georgia (L33, L40–L51) as well as individuals of the lasistana morphotype from the Turkish province Artvin (L52–L54) and a single individual from the western-most part of Georgia from a population (L51) otherwise composed exclusively of individuals of frutis morphotype (Fig. 2, Supplementary Table S1). Individuals of the lasistana morphotype (with a complete dart apparatus, Fig. 2) did not form a monophyletic group, but rather were scattered among individuals of the frutis morphotype (lacking dart and accessory sac). Clade B comprised exclusively individuals of the frutis morphotype from the Colchis lowlands and adjacent southern slopes of the Greater and northern slopes of the Lesser Caucasus in western Georgia (L14–L15, L27, L36–L39; Fig. 2, Supplementary Table S1). Individuals in clade C were predominantly of the veselyi morphotype with a completely reduced dart apparatus except for two individuals of the circassica morphotype with a complete dart apparatus (Fig. 2) from the Lesser Caucasus between Borjomi and Bakuriani (L31–L32). Most individuals with a veselyi morphotype in this clade were from localities in the Lesser Caucasus in southern to central Georgia or the southern slopes of the Greater Caucasus north of Tbilisi (L7–L8, L28–L30, L34–L35). Only one individual of veselyi morphotype in clade C came from a locality further to the west (L6). This lineage represented the sister group to the remaining individuals in this clade (Fig. 1). Clade D included individuals of circassica and veselyi morphotypes from the Greater Caucasus in Georgia (L3–L5, L9–L13, L16– L26) as well as the Ciscaucasian plain in the surroundings of Stavropol (L1–L2). Although individuals were mostly grouped according to morphotype within this clade, C. frutis circassica and C. frutis veselyi did not form monophyletic units, but rather sub-clades with individuals of circassica morphotype interspersed between clades of veselyi morphotype (Fig. 2). Within clade D, two individuals of veselyi morphotype from the vicinity of Mestia (L4–L5) were placed in a sister group relationship with the remaining individuals in this clade. The next lineage to branch off grouped individuals of circassica morphotype distributed on the southern slopes of the Greater Caucasus (L17–L18, L20–L22). Individuals of veselyi morphotype from the Ciscaucasian plain (L1–L2) were placed into a sub-clade with individuals of the same morphotype from north-western Georgia (L10–L11). The last lineage to branch off in clade D included individuals of circassica morphotype from central Georgia (L19, L23–L26) and individuals mostly of veselyi morphotype from northern-central to north-western Georgia (L3, L9, L12–L13, L16). Only one individual of circassica morphotype was included in the latter group (L16). Interestingly, this individual was found at one of only two sampling sites where two different morphotypes occurred sympatrically (Fig. 2, Supplementary Table S1). 3.2. Analyses of AFLP data After a distinct increase from K = 1 to K = 2, mean L(K) increased only slowly with increasing K (Supplementary Fig. S1A). DK showed a maximum at K = 2 (Supplementary Fig. S1B). Estimated

135

L58 L59 L60

L57

L54 L55 L56

L53

L52

L51

L49

L50

L48

L45 L46 L47

L44

L43

L35 L36

L37 L38 L39 L40 L41 L42

L34

L28

L29 L30 L31 L32 L33

L24

L25 L26 L27

L23

L22

L20 L21

L14 L15 L16 L17 L18 L19

L13

L05 L06 L07 L08 L09 L10 L11 L12

L04

L02

L03

L01

M.T. Neiber, B. Hausdorf / Molecular Phylogenetics and Evolution 93 (2015) 129–142

Fig. 3. Results of the admixture analysis with the AFLP data of Circassina with STRUCTURE for K = 2. Individuals are grouped according to geographic origin (locality numbers L01–L60, see Supplementary Table S1).

ancestries of individuals for K = 2 are shown in Fig. 3 (for the run with the highest L(K) value). Individuals were grouped into geographically coherent clusters: a ‘south-western’ cluster including individuals of frutis and lasistana morphotypes from Turkey and westernmost Georgia and a ‘northern’ to ‘north-eastern’ cluster including individuals of circassica or veselyi morphotypes and to a lesser extent also individuals of frutis morphotype from the Greater Caucasus, the Ciscaucasian plain and the Lesser Caucasus in central/southern Georgia (Fig. 4). There is a narrow zone with populations with mixed ancestry between the ‘south-western’ and ‘northern’ clusters in the Colchis lowlands. The admixture analysis also indicated extensive gene flow between the ‘northern’ and ‘south-western’ cluster along an east/west transect through the Georgian Lesser Caucasus (L37–L42; Fig. 4).

The majority of individuals that were grouped in the mitochondrial clade A had a high proportion of ancestry in the ‘south-western’ cluster, while those from mitochondrial clades C and D usually had a high proportion of ancestry in the ‘northern’ cluster (Fig. 4). Almost all individuals that have mitochondrial clade B haplotypes showed high levels of admixture, which is in line with the occurrence of these individuals in the centre of the contact zone between the ‘northern’ and ‘south-western’ clusters (Fig. 4). The HIest analysis detected a single credible F2 hybrid (L37, individual 1523, see Supplementary Table S1) among the 67 individuals that had an ancestry of less than 90% in one of the two clusters detected in the STRUCTURE analysis. No credible F1 hybrids or backcrosses to parental taxa were detected.

L1

Rus

L2

Grea

ter C

auca

sia

L3

sus

L4 L5

L6

Black Sea

L19 L20 L21

Colchis Lowlands L46

L44

L52

L54

L58 L60

0

L59

L49

Geo

L14

L27

L37 L39

L50

L55

L26

L31 L33

L42

L53

L41

L40

L7

rgia

L8 L9

L28

L35

L29

L34

L45

L51

L13

L15

L43

L48

L12

L16 L18

L38

L11

L17

L22 L23 L36 L25

L47

L10 L24

L30

L32

Les

ser

Cau

cas

L57 L56

Turkey

us

Armenia

Az

erb

aija

n

200 km

Fig. 4. Admixture of the Circassina populations in the Caucasus region and north-eastern Turkey calculated based on the AFLP data using STRUCTURE. Proportions of parentages in the pie charts indicate mean assignments for the population (see text) to the ‘northern’ cluster (blue) and the ‘south-western’ cluster (red). Data for sampling sites L01–L60 are given in Supplementary Table S1. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

M.T. Neiber, B. Hausdorf / Molecular Phylogenetics and Evolution 93 (2015) 129–142

0.6

0.8

136

0.2 0.0

1530, B 1568, B 1523, B

1527, B

1526, B

1542

−0.2

dimension 2

0.4

1972

1567, B

1541, B

−0.4

2406

2002 2001

−0.6

2405

3.3. Time frame of the diversification of Circassina

1999

−0.4

−0.2

0.0

cluster, individuals of circassica morphotype from the southern slopes of the Greater Caucasus included in the analysis (L19, L21, L23–L24, L26) formed a monophyletic group (albeit not significantly supported; PP: 0.81). This group was joined with individuals belonging to the veselyi morphotype from the Greater Caucasus (L04–L05, L16) and a maximally supported clade including individuals of veselyi morphotype from localities in the Greater Caucasus (L07–L13) and the Ciscaucasian Plain (L01–L03) and a single individual of frutis morphotype from a locality east of Kutaisi (L14). Individuals of veselyi morphotype included in the analysis from localities in the Lesser Caucasus (L28–L29) were recovered as the sister group to all other specimens in the ‘northern’ cluster (Fig. 6).

0.2

0.4

dimension 1 Fig. 5. First two dimensions of a non-metric multidimensional scaling based on Jaccard distances between AFLP data of 159 individuals of Circassina. Dots, ‘southwestern’ cluster; squares, ‘northern’ cluster. Individuals are coloured according to the (sub-) species sensu Hausdorf (2001) (see Fig. 2). The numbers refer to vouchers as detailed in Supplementary Table S1 and the letter B indicates individuals belonging to the mitochondrial clade B (see Fig. 2). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

A non-metric multidimensional scaling (Fig. 5) and the neighbour-net (Supplementary Fig. S2) constructed using Jaccard distances based on the AFLP data support a division of Circassina into two main clusters corresponding to the ‘northern’ and ‘south-western’ cluster that were inferred in the admixture analysis. The ‘south-western’ cluster included individuals of frutis and lasistana morphotypes (with few exceptions corresponding to populations with mitochondrial haplotypes from clade A), the ‘northern’ cluster individuals of veselyi, circassica and frutis morphotypes (corresponding to populations with mitochondrial haplotypes from clades B–D). Most outliers in the non-metric multidimensional scaling (Fig. 5) were admixed individuals belonging to the frutis morphotype, and often had mitochondrial clade B haplotypes (Fig. 2). The individual identified as a F2 hybrid (1523) in the analysis with HIest was included in the ‘south-western’ cluster in the neighbour-net (Supplementary Fig. S2) as well as the non-metric multidimensional scaling (Fig. 5) although in both analyses in a somewhat marginal position. The Bayesian analysis of the AFLP data that included only individuals with at least 90% parentage in one of the clusters inferred from the admixture analysis (Fig. 6) maximally supported the monophyly of Circassina. Likewise, the monophyly of the ‘south-western’ cluster (including individuals of frutis and lasistana morphotypes) was maximally supported in this analysis. Nested within a clade of individuals with frutis morphotypes from Georgia (L33, L43–L46, L48–L51), a clade with individuals belonging to the lasistana morphotype (L52–L54) was recovered as a maximally supported monophyletic group with the exception of one individual from Gonio fortress in western Georgia (L51) that, albeit possessing a dart sac, was intermediate between typical C. frutis frutis and C. lasistana with regard to the epiphallus/flagellum ratio. Individuals of frutis morphotype from Turkey (L55, L57–L58) were recovered as a monophyletic group (PP: 1.0) that was placed in a sister group relationship to all remaining individuals included in the ‘south-western’ cluster (Fig. 6, Supplementary Table S1). The monophyly of the ‘northern’ cluster was also significantly supported (PP: 0.98) and a sister group relationship with the ‘south-western’ cluster was recovered (Fig. 6). In the ‘northern’

The maximum clade credibility tree obtained from the Beast analysis recovered the same major clades (A–D) with maximal support (Fig. 7) as the BI with MrBayes (Fig. 2). Given the calibration imposed on the root, the split between clade D and clade A–C was dated to an age of 6.4 Ma before present (95% highest posterior density interval (HPD) 3.8–10.7) in the late Miocene, while the most recent common ancestor of clade A–C was dated to an age of 4.7 Ma before present (95% HPD 1.9–4.5) in the early Pliocene. Lineages within clade D started to diverge 4.3 Ma before present (95% HPD 2.3–7.6) in the first half of the Pliocene, while lineages within clades A–C started to diverge slightly earlier in the early Pliocence (clade A: 2.9 Ma, 95% HPD 1.4–5.5; clade B: 2.3 Ma, 95% HPD 0.9–4.4; clade C 3.0 Ma, 95% HPD 1.4–5.4). 3.4. Isolation by distance vs. genetic clustering The results of the dbRDA of Jaccard distance (AFLP data) and GTR + G distance (mitochondrial sequences) matrices as response variables and the geographic coordinates of the sampling sites as predictor variables were significant (p < 0.01), but the proportion of variation of the genetic distances explained by isolation by distance was low for the AFLP data (14%) and moderately high for the mitochondrial data (37.8%; Table 1). Incorporating information on AFLP cluster affiliation or affiliation to one of the major mitochondrial clades as covariables indicated that a very low proportion of the variation of the genetic distances was explained by isolation by distance, both for the AFLP (5.3%) and mitochondrial (0.7%) data sets. Contrarily, the results of the dbRDA of AFLP cluster affiliations or affiliations to major mitochondrial clades as predictor variables and the genetic distance matrices as response variables, respectively, showed that genetic clustering explained a low to moderately high proportion of the variation of the Jaccard distances (14.7%) and a very high proportion of the variation of the GTR + G distances (91.7%) (Table 1). Incorporating information on the geographical origin of specimens as covariables showed that a low proportion of the variation of the Jaccard distances (6.1%) was explained by genetic clustering in the AFLP data set, while a high proportion of the variation of the GTR + G distances was explained by genetic clustering (54.6%) in the mitochondrial data set (Table 1). Taken together, both effects, isolation by distance as well as genetic clustering explained a significant proportion of variation of the genetic distances, but the effect of genetic clustering was generally higher. The Mantel tests between geographic distances and Jaccard distances based on AFLP data as well as geographic distances and GTR + G distances between the mitochondrial sequences were both significant (p < 0.001). The observed positive correlations between geographic distances and Jaccard distances (r = 0.454) and geographic distances and GTR + G distances (r = 0.285) were moderate or weak, respectively (see also Supplementary Fig. S3 for scatter plots of geographic distances versus genetic distances).

137

M.T. Neiber, B. Hausdorf / Molecular Phylogenetics and Evolution 93 (2015) 129–142

2208, A, L55

0.05

1

2401, L58

1

2403, L57

1

2404, L57

1

2406, L57 1524, A, L43

0.96

1743, A, L43 1525, A, L33

0.97 0.62 0.98

0.97

1

1983, L43 2067, L43 1981, L43

1

0.7

1982, L43 2004, L49 2021, L50

0.84

1730, A, L44

0.5

1732, A, L46

1

2038, L44

0.77

2037, L44

0.96

2039, L44

0.98

2044, L46

1

2045, L46 2010, L48

0.93

2011, L48

0.99

0.71

2012, L48

1

2013, L48 2195, A, L52 2198, A, L53

0.63

2196, A, L54

1

2391, L53

1

2393, L53

0.72

2392, L53

0.99

1

2394, L53 2395, L53

0.99

1

2397, L53

0.96 0.71

'south-western' cluster

1738, A, L45

0.54 0.88

2396, L53

0.96

2399, L53 2398, L53 2033, L51

1 1

2035, L51 1992, L51

0.95

1993, L51

0.77

2036, L51 1

1512, A, L51

0.89

1515, D, L21

1

1729, A, L51 1734, D, L19 2048, L19

1 0.97

0.95

2049, L19 2094, L24 1973, L23

1 0.88

0.81

2093, L24

1

2095, L24 1535, D, L26

0.73

1672, D, L16 1737, D, L5

1 1

0.75 0.98 1

2074, L5 2088, L4 2089, L4

0.97

1593, D, L11

0.81

1520, D, L11

0.99

1

1567, B, L14 1519, C, L8

0.99

1741. C, L7 2009, L8

1 1

0.68

1522, D, L9 1521, D, L13

0,84

1566, D, L12 1591, D, L13

0.98

1

'northern' cluster

2090, L4 1518, D, L10

0.97

1592, D, L13

0.7

1736, D, L3

1

2017, L13

1

1

0.98

1742, D, L4 2087, L4 2073, L5

1

2018, L13 2019, L13

0.99

2071, L3

1

2072, L3 1661, D, L2

1

2024, L2

0.96

1665, D, L1

1 1

2026, L1

1

1 1

2027, L1 2029, L1 2028, L1 2025, L2 1538, C, L28

1

1 0.98

1 0.99

1975, L29 1976, L29 1979, L28 Fruticocampylaea narzanensis 1719 Kokotschashvilia holotricha 1616 Caucasigena eichwaldi 1783

Outgroups

Fig. 6. Bayesian 50% majority rule tree based on multilocus AFLP data. Included in the analysis were individuals that had at least 90% ancestry in the ‘northern’ cluster (blue bar) or ‘south-western’ cluster (red bar) as delimited in the STRUCTURE analysis. The numbers at the nodes indicate posterior probability values. The numbers at the tips refer to vouchers as detailed in Supplementary Table S1 and are coloured according to the (sub-)species sensu Hausdorf (2001) as in Fig. 2. The letters A–D indicate the main mitochondrial haplotype group an individual was assigned to (see Fig. 2). The L-numbers refer to localities (see Supplementary Table S1). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

138

M.T. Neiber, B. Hausdorf / Molecular Phylogenetics and Evolution 93 (2015) 129–142 0.98

Outgroups

4.4

1

Clade A

2.9

7.4 1 4.7 1

Clade B

2.3 0.5 4.2 1

1

6.4

3.0

Clade C

1

Clade D

4.3 0.5

Miocene

Ma

11

10

9

Pliocene

8

7

6

5

4

Pleistocene

3

2

1

0

Fig. 7. Bayesian chronogram for major mitochondrial lineages in Circassina estimated by Beast from cox1 and 16S rDNA sequences. Scale bars at nodes indicate the 95% credibility intervals of node ages. The numbers above the scale bars indicate the posterior probability value for the nodes and numbers below the scale bars the estimated age of the node in Ma before present. For the composition of clades A–D, see Fig. 2.

Table 1 Results of distance-based redundancy analyses (dbRDA) indicating the proportion of variation of the Jaccard distances (AFLP data) and GTR + G distances (mitochondrial data) explained by isolation by distance (modelled by geographic coordinates) or genetic clustering (modelled as proportions of parentage from the STRUCTURE analysis at K = 2 for the AFLP data or affiliation to the major mitochondrial clades). Pseudo-F

p

Proportion of variation explained

A. AFLP data 1. Isolation-by-distance Covariables None Proportions of parentage

12.667 5.058

0.000 0.007

0.140 0.053

2. Genetic clustering Covariables None Geographic coordinates

13.410 5.877

0.000 0.000

0.147 0.061

B. Mitochondrial sequences 1. Isolation-by-distance Covariables None Major mitochondrial linages

20.64 2.92

0.000 0.019

0.378 0.007

2. Genetic clustering Covariables None Major mitochondrial linages

246.793 15.379

0.000 0.000

0.917 0.546

3.5. Barriers to gene flow and long-distance dispersal Pairs of individuals that show high genetic distances across short geographic distances are concentrated in the Colchis lowlands and the adjacent southern slopes of the Greater Caucasus and the northern slopes of the Lesser Caucasus, both with regard to Jaccard distances obtained from the AFLP data set and GTR + G distances obtained from cox1 and 16S rDNA sequences (Fig. 8). This finding was in agreement with the presence of highly admixed individuals in this area (Fig. 3). An additional barrier across the

Kura valley in central Georgia is indicated by the mitochondrial data alone (Fig. 8B). Pairs of individuals that showed low genetic distances, congruently for Jaccard distances from AFLP data and GTR + G distances from mitochondrial sequences, across large geographic distances uncovered two routes along which potential long-distance dispersal events may have occurred. The first route was indicated by connections between individuals mainly from the southern slopes of the Greater Caucasus and individuals from the Ciscaucasian Plain across the main crest of the Greater Caucasus (Fig. 8). The second route was indicated by connections between individuals from south-western Georgia with individuals from the surroundings of Trabzon in north-eastern Turkey (Fig. 8). 3.6. Associations between traits of the dart apparatus and environmental variables The Sambada analysis detected four environmental variables that were significantly associated with the presence/absence of the dart sac, while the presence/absence of the mucus glands was, in addition to environmental variables, also significantly associated with genetic clustering. Therefore, we regarded associations between the presence/absence of the mucus glands and the environmental variables as possible false positives, while associations between the presence/absence of the dart sac could not be explained by genetic clustering. The presence/absence of the dart sac was significantly associated (at a confidence level of 1.43  10 4 after Bonferroni correction) with maximum temperature in March and May and with precipitation of the driest month, which was usually May, and precipitation seasonality. 4. Discussion 4.1. Origin and diversification of Circassina Circassina belongs to a group of hygromiids that can be traced back to a single colonisation of the region of the present day

M.T. Neiber, B. Hausdorf / Molecular Phylogenetics and Evolution 93 (2015) 129–142

A

0

B

200 km

0

C

0

139

200 km

D

200 km

0

200 km

Fig. 8. Visualization of past/present barriers to gene flow and potential long-distance dispersal events. Potential barriers to gene flow are visualized by plotting lines between the pairs of individuals that belong to both, the quintile of highest genetic distances and the quintile of lowest geographic distances. A: Using Jaccard distances based on AFLP data, and B: using GTR + G distances for mitochondrial cox1 and 16S rDNA sequences. Potential long-distance dispersal events are visualized by plotting lines between the pairs of individuals that belong to both, the quintile of lowest genetic distances and the quintile of highest geographic distances. C: Using Jaccard distances based on AFLP data, and D: using GTR + G distances for mitochondrial cox1 and 16S rDNA sequences.

Greater Caucasus during the Miocene (Walther et al., in press), which then was an island in the Paratethys Sea (Popov et al., 2004). The restriction of most genera of this endemic radiation to the Greater Caucasus indicates that also Circassina probably originated in the Greater Caucasus. The tree based on mitochondrial sequences (Figs. 1 and 2) as well as the non-metric multidimensional scaling, the tree and the network based on AFLP markers (Figs. 4–6, Supplementary Fig. S2) show an early divergence of Circassina along a line running through the Colchis lowlands in the late Miocene (6.4 Ma before present; see Fig. 7). The boundary running through the Colchis lowlands is also visualized by high genetic distances between individuals across short geographic distances with regard to both mitochondrial sequences and AFLP markers (Fig. 8A and B). Divergence of population groups across the Colchis lowlands indicates that this area did not act as a homogeneous refuge as implied in some studies (van Andel and Tzedakis, 1996; Kikvidze and Ohsawa, 2001; Volkova et al., 2013; Wielstra et al., 2013). The observed clinal pattern across the Colchis Lowlands cannot be explained by a primary cline because our data disfavour a clear-cut pattern of isolation by distance as shown by the results of the STRUCTURE analysis (Figs. 3 and 4, Supplementary Table S1), the results of the dbRDA (Table 1) and Mantel tests (Supplementary Fig. S3) as well as the results of the non-metric

multidimensional scaling (Fig. 5), but clearly support the presence of two genetic clusters that came secondarily into contact in this region. One mechanism that may have contributed to the establishment and maintenance of the boundary through the Colchis lowlands are marine transgressions. The Colchis lowlands have been strongly influenced by sea level changes of the Black Sea during the Pleistocene and Holocene, with some transgressions surpassing present day sea levels (Muratov et al., 1974; Connor et al., 2007; Yanina, 2014). During these transgressions the westernmost parts of the Colchis lowlands were partially flooded or paludified rendering them uninhabitable for Circassina, while during Black Sea regressions Circassina may have re-colonised the area from the adjacent slopes of the Greater and Lesser Caucasus, which represented separate refuges, resulting in the present day observed pattern of a contact zone between the two population groups (Fig. 8, Supplementary Table S1). The contact zone may remain narrow as a result of limited migration into the zone or partial pre- or postzygotic barriers as can be indirectly inferred from the rarity of young hybrids (only one credible F2 hybrid was detected in the HIest analysis, see also Supplementary Fig. S4). Within the ‘northern’ cluster individuals from the Lesser Caucasus (L28–L29) were resolved as the sister group of individuals from the Greater Caucasus indicating an additional refuge in the Lesser Caucasus (Fig. 6, Supplementary Table S1) corresponding to a part of the

140

M.T. Neiber, B. Hausdorf / Molecular Phylogenetics and Evolution 93 (2015) 129–142

range of mitochondrial clade C haplotypes. This refuge is approximately concordant with a refuge proposed for the Caucasian salamander Mertensiella caucasica (Waga, 1876) by Tarkhnishvili et al. (2000), who found a phylogeographic break between populations from the Pontus Mountains and the western Lesser Caucasus and populations from the north-eastern Lesser Caucasus in this species. There are long internal branches between the most recent common ancestor of Circassina and the diversification within the main mitochondrial haplotype groups, which occurred during the Pliocene and Pleistocene (Fig. 2). This possibly indicates bottlenecks during these periods that may be related to gradual global cooling or the Pleistocene glacials. A similar, but even more extreme pattern has been found in the Caucasian Fruticocampylaea species (Walther et al., in press). Several studies identified the Greater Caucasus chain, which rises to 5642 m above sea level, as a major barrier hampering colonisation of Ciscaucasia by mammals from Transcaucasian refugia (Orth et al., 2002; Seddon et al., 2002; Dubey et al., 2006). In contrast, the mainly passively dispersing snail Circassina, which is known from several areas north of the main ridge of the Greater Caucasus (Supplementary Table S1; Hausdorf, 2001), and cold-adapted steppe vipers (Zinenko et al., 2015) crossed this barrier several times and lack deep phylogeographic splits corresponding to the main ridge. Low genetic distances, both for AFLP and mitochondrial data, between Circassina specimens from south of the main crest of the Greater Caucasus and specimens from Ciscaucasia indicate that long-distance dispersal played a role in shaping the current distribution of Circassina (Fig. 8). We also found some evidence for long-distance dispersal between localities in the Lesser Caucasus of Georgia and the westernmost localities in Turkey (Fig. 8). 4.2. Dart apparatus evolution in Circassina If the dart sac and the mucus glands were lost only once in Circassina, we would expect that the specimens without dart sac and without mucus glands respectively form monophyletic groups. However, individuals with a dart apparatus with dart sac and mucus glands were nested in clades with individuals lacking a dart sac and partially also the mucus glands (Fig. 6). The number of losses cannot be easily inferred from the tree because the accessory genital organs might have been regained or lost by interbreeding with individuals possessing or lacking dart sac and the mucus glands. In any case, it is obvious that there were several changes of these character states. We tested for environmental correlates of genital morphotype (presence/absence of the dart sac, presence/absence of mucus glands) to infer conditions that favour the loss of the accessory genital organs. We found significant associations between the presence/absence of the dart sac and precipitation during the driest month, which was usually May, and the maximum temperatures in March and May. May falls in the mating season of Circassina so that these correlations might actually point to a causal relationship. Differing effects of dart shooting on the duration of the mating behaviour of helicoids have been reported. While Adamo and Chase (1988, 1990) found that dart shooting decreased courtship duration in Cornu aspersum (but see Chase, 2007), Jeppesen (1976), Lind (1976) observed a decrease in courtship duration for Helix pomatia when no dart shooting behaviour was performed. Adamo and Chase (1990) also reported a retarding effect on locomotion by the mucus transmitted via the dart. The mucus causes a reconfiguration of tracts of the bursa copulatrix complex in Cornu aspersum, facilitating the uptake of the spermatophore into the bursa tract diverticulum by letting this tract become more accessible, while closing off the tract leading to the bursa. Thus, more sperm of a successful dart shooter is able to

escape digestion in the bursa and fertilisation success is increased (Koene and Chase, 1998; Chase and Blanchard, 2006). Assuming a similar function of the mucus in Circassina, in which a bursa tract diverticulum is lacking (as in all hygromiids), the mucus’ influence on the bursa copulatrix tract could cause retardation in spermatophore uptake and thus prolong mating duration. Similarly, performing dart shooting behaviour as such could also result in a prolongation of mating behaviour. The correlation of a lack of a dart sac with drier and cooler conditions during the mating season could then be viewed to indicate that natural selection for reducing the desiccation risk by shortening the mating may have contributed to the loss of the dart apparatus in populations of Circassina in drier areas. In moister areas the desiccation risk is smaller so that the fitness gain resulting from the increased fertilisation success of a successful dart shooter, which may be associated with a prolongation of the mating, is not overcome by natural selection for shorter mating. 4.3. Taxonomy of Circassina Our phylogenetic analyses based on mitochondrial loci and on AFLP markers showed that none of the anatomically defined (sub-) species distinguished by Schileyko (1978) and Hausdorf (2001) represent reciprocally monophyletic groups and that the presence or absence of structures of the dart apparatus are inadequate characters to delimit taxa in Circassina. Taking into account that the admixture analysis based on the AFLP markers indicated not complete reproductive isolation, we classify all Circassina populations in only a single, albeit genetically deeply structured and phenotypically polymorphic species. This species has to bear the oldest available name for the group, C. frutis (Pfeiffer, 1859). We propose to classify the ‘northern’ and ‘south-western’ clusters as delimited in the admixture analysis, the non-metric multidimensional scaling and the phylogenetic analysis of the AFLP data (Figs. 3–7), that represent the most divergent, but geographically coherent clusters, as subspecies. For the decision, which of the two groups ought to be named C. f. frutis, the type locality of this taxon is decisive. The lectotype of C. frutis designated by Hausdorf (2001) was collected in Qulevi (=Redut-Kale) near the Black Sea coast at the mouth of Khobistsqali River in Georgia. Unfortunately, this locality is situated in the contact zone between the two genetic entities. In accordance with the closer proximity to localities assigned to the ‘northern’ cluster, the admixture analysis inferred an ancestry of 53.3–60.3% in the ‘northern’ cluster and 39.7–46.7% ancestry in the ‘south-western’ cluster for individuals from locality L36 and a mitochondrial haplotype belonging to clade B, which is present in both clusters (Figs. 1, 2 and 6, Supplementary Table S1). Individuals from the vicinity of the type locality (L36) belong to the morphotype without dart sac, but with mucus glands, which is the predominant morphotype of the ‘south-western’ cluster. Because there is no objectively correct solution in this case, we propose applying the name C. frutis frutis to the ‘south-western’ cluster to conserve its usage for the cluster, in which most specimens lack a dart sac, but possess mucus glands, the morphotype for which this name was used so far (Hausdorf, 2001). Circassina lasistana differs from the surrounding C. frutis frutis populations in possessing a dart sac and a much shorter flagellum (Hausdorf, 2001). The flagellum forms the tail of the spermatophore. Dart sac and spermatophore (and, hence, flagellum) are under sexual selection (Koene and Schulenburg, 2005) and are often used to delineate species of land snails. In the AFLP tree (Fig. 6), all C. lasistana individuals except one from Gonio fortress form a well-supported monophyletic group (PP: 1.0), whereas they are scattered among C. frutis frutis in the mitochondrial gene tree (Fig. 2). This discrepant pattern may be the result of reconstructing

M.T. Neiber, B. Hausdorf / Molecular Phylogenetics and Evolution 93 (2015) 129–142

the phylogeny from a single marker, i.e. mitochondrial sequences, versus the summary consensus of multiple AFLP loci. Another explanation of the discrepancy between nuclear AFLP markers and mitochondrial sequences might be that at least a part of the nuclear markers became sorted more rapidly than the presumably neutral mitochondrial DNA because of sexual selection. Another, not mutually exclusive explanation would be that the polyphyly with regard to mitochondrial DNA is the result of gene flow between C. lasistana and typical C. frutis frutis without dart sac. Actually, individuals that are intermediate between C. frutis sensu Hausdorf (2001) and C. lasistana with respect to the epiphallus/flagellum ratio (1.0–1.1, n = 3) have been found in Gonio fortress in western Georgia (L51), in one of which a dart sac was present. Thus, we place C. lasistana in the synonymy of C. frutis frutis, although the distinct clustering of the C. lasistana populations in the AFLP network (Supplementary Fig. S2) and tree (Fig. 6) along with morphological differentiation could be indicative of an ongoing differentiation process triggered by sexual selection that may lead to speciation, if not merely reflecting relative geographical proximity of samples. The oldest available name for the ‘northern’ cluster is C. frutis circassica (Mousson, 1863), the lectotype of which was collected in Nakalakevi in Georgia (Hausdorf, 2001). Individuals from that locality (L22) were assigned to the ‘northern’ cluster in the analyses based on AFLP data (Figs. 3–6, Supplementary Table S1) and possessed mitochondrial haplotypes (clade D) that exclusively occur in the ‘northern’ cluster (Figs. 1 and 2, Supplementary Table S1). Since C. frutis veselyi (Frankenberger, 1919) described from Ananuri in Georgia also belongs to the ‘northern’ cluster, we place this name in the synonymy of C. frutis circassica. Acknowledgments We would like to thank F. Walther (Hamburg) for help with data and maps, P. Kijashko (St. Petersburg), L. Mumladze (Tbilisi) and F. Walther for their support during fieldwork in the Caucasus, J. Sauer (Bielefeld) for laboratory help and the Volkswagen foundation for funding the project ‘‘Biogeography of the land molluscs of the Caucasus region’’. Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.ympev.2015.07. 012. References Adamo, S.A., Chase, R., 1988. Courtship and copulation in the terrestrial snail Helix aspersa. Can. J. Zool. 66, 1446–1453. Adamo, S.A., Chase, R., 1990. The ‘‘love dart’’ of the snail Helix aspersa injects a pheromone that decreases courtship duration. J. Exp. Zool. 255, 80–87. Anderson, M.J., 2004. DISTLM v. 5: a FORTRAN computer program to calculate a distance-based multivariate analysis for a linear model. Department of Statistics, University of Auckland, Auckland. Arrigo, N., Tuszynski, J.W., Ehrich, D., Gerdes, T., Alvarez, N., 2009. Evaluating the impact of scoring parameters on the structure of intra-specific genetic variation using RawGeno, an R package for automating AFLP scoring. BMC Bioinformatics 10, 33. Bazinet, A.L., Zwickl, D.J., Cummings, M.P., 2014. A gateway for phylogenetic analysis powered by grid computing featuring GARLI 2.0. Syst. Biol. 63, 812– 818. Bouckaert, R., Heled, J., Kühnert, D., Vaughan, T., Wu, C.-H., Xie, D., Suchard, M.A., Rambaut, A., Drummond, A.J., 2014. Beast 2: a software platform for Bayesian evolutionary analysis. PLoS Comput. Biol. 10, e1003537. Bryant, D., Moulton, V., 2004. Neighbor-net: an agglomerative method for the construction of phylogenetic networks. Mol. Biol. Evol. 21, 255–265. Castresana, J., 2000. Selection of conserved blocks from multiple alignments for their use in phylogenetic analysis. Mol. Biol. Evol. 17, 540–552. Chase, R., 2007. The function of dart shooting in helicid snails. Am. Malacological Bull. 23, 183–189.

141

Chase, R., Blanchard, K.C., 2006. The snail’s love-dart delivers mucus to increase paternity. Proc. R. Soc. B 273, 1471–1575. Chiba, S., 1999. Accelerated evolution of land snails Mandarina in the oceanic Bonin Islands: evidence from mitochondrial DNA sequences. Evolution 53, 460–471. Connor, S.E., Thomas, I., Kvavadze, E.V., 2007. A 5600-yr history of changing vegetation, sea levels and human impacts from the Black Sea coast of Georgia. The Holocene 17, 25–36. Dray, S., Dufour, A.B., 2007. The ade4 package: implementing the duality diagram for ecologists. J. Stat. Softw. 22, 1–20. Dubey, S., Zaitsev, M., Cosson, J.F., Abdukadier, A., Vogel, P., 2006. Pliocene and Pleistocene diversification and multiple refugia in a Eurasian shrew (Crocidura suaveolens group). Mol. Phylogenet. Evol. 38, 635–647. Evanno, G., Regnaut, S., Goudet, J., 2005. Detecting the number of clusters of individuals using the software Structure: a simulation study. Mol. Ecol. 14, 2611–2620. Falush, D., Stephens, M., Pritchard, J.K., 2007. Inference of population structure using multilocus genotype data: dominant markers and null alleles. Mol. Ecol. Notes 7, 574–578. Fitzpatrick, B.M., 2012. Estimating ancestry and heterozygosity of hybrids using molecular markers. BMC Evol. Biol. 12, 131. Folmer, O., Black, M., Hoeh, W., Lutz, R., Vrijenhoek, R., 1994. DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates. Mol. Mar. Biol. Biotech. 3, 294–299. Fu, Y.-X., 1997. Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics 147, 915–925. Giusti, F., Manganelli, G., 1987. Notulae Malacologicae, XXXVI. On some Hygromiidae (Gastropoda: Helicoidea) living in Sardinia and in Corsica. (Studies on the Sardinian and Corsican Malacofauna VI). Bolletino Malacologico 23, 123–206. Gradstein, F., Ogg, J., Schmitz, M., Ogg, G., 2012. The Geological Time Scale 2012. Elsevier, Oxford, Amsterdam, Waltham. Hausdorf, B., 1998. Phylogeny of the Limacoidea sensu lato (Gastropoda: Stylommatophora). J. Molluscan Stud. 64, 35–66. Hausdorf, B., 2001. A systematic revision of Circassina from the western Caucasus region (Gastropoda: Hygromiidae). J. Molluscan Stud. 67, 425–446. Hennig, C., Hausdorf, B., 2013. The prabclus package version 2.2-4. Department of Statistical Science, University College London, London, . Hewitt, G.M., 2000. The genetic legacy of the Quaternary ice ages. Nature 405, 907– 913. Hijmans, R.J., Cameron, S.E., Parra, J.L., Jones, P.G., Jarvis, A., 2005. Very high resolution interpolated climate surfaces for global land areas. Int. J. Climatol. 25, 1965–1978. Holland, B.R., Clarke, A.C., Meudt, H.M., 2008. Optimizing automated AFLP scoring parameters to improve phylogenetic resolution. Syst. Biol. 57, 347–366. Huson, D.H., Bryant, D., 2006. Application of phylogenetic networks in evolutionary studies. Mol. Biol. Evol. 23, 254–267. Jeppesen, L.L., 1976. The control of mating behavior in Helix pomatia L. (Gastropoda: Pulmonata). Anim. Behav. 24, 275–290. Katoh, K., Standley, D.M., 2013. MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol. Biol. Evol. 32, 772– 780. Kikvidze, Z., Ohsawa, M., 2001. Richness of Colchic vegetation: comparison between refugia of south-western and East Asia. BMC Ecol. 1, 6. Koene, J.M., Chase, R., 1998. Changes in the reproductive system of the snail Helix aspersa caused by mucus from the love dart. J. Exp. Biol. 201, 2313–2319. Koene, J.M., Schulenburg, H., 2005. Shooting darts: co-evolution and counteradaptation in hermaphroditic snails. BMC Evol. Biol. 5, 25. Kruskal, J., 1964. Multidimensional scaling by optimizing goodness of fit to a nonmetric hypothesis. Psychometrika 29, 1–27. Lanfear, R., Calcott, B., Ho, S.Y.W., Guindon, S., 2012. Partitionfinder: combined selection of partitioning schemes and substitution models for phylogenetic analyses. Mol. Biol. Evol. 29, 1695–1701. Librado, P., Rozas, J., 2009. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25, 1451–1452. Lind, H., 1976. Causal and functional organization of the mating behaviour sequence in Helix pomatia (Pulmonata, Gastropoda). Behaviour 59, 162–202. McArdle, B.H., Anderson, M.J., 2001. Fitting multivariate models to community data: a comment on distance-based redundancy analysis. Ecology 82, 290–297. McDonald, J.H., 2014. Handbook of Biological Statistics, third ed. Sparky House Publishing, Baltimore. Meirmans, P.G., 2012. The trouble with isolation by distance. Mol. Ecol. 21, 2839– 2846. Muratov, V.M., Ostrovskiy, A.B., Fridenberg, E.O., 1974. Quaternary stratigraphy and palaeogeography on the Black Sea coast of Western Caucasus. Boreas 3, 49–60. Mumladze, L., Tarkhnishvili, D., Murtskhvaladze, M., 2013. Systematics and evolutionary history of large endemic snails from the Caucasus (Helix buchii and H. goderdziana) (Helicidae). Am. Malacological Bull. 31, 225–234. Myers, N., Mittermeier, R.A., Mittermeier, C.G., da Fonseca, G.A.B., Kent, J., 2000. Biodiversity hotspots for conservation priorities. Nature 403, 853–858. Nakhutsrishvili, G., 2013. The Vegetation of Georgia (South Caucasus). Springer, Berlin, Heidelberg. Nix, H.A., 1986. A biogeographic analysis of Australian elapid snakes. In: Longmore, R. (Ed.), Atlas of Elapid Snakes of Australia, Australian Flora and Fauna Series 7. Australian Government Publishing Service, Canberra, pp. 4–15.

142

M.T. Neiber, B. Hausdorf / Molecular Phylogenetics and Evolution 93 (2015) 129–142

Nordsieck, H., 1987. Revision des Systems der Helicoidea (Gastropoda: Stylommatophora). Archiv für Molluskenkunde 118, 8–50. Orth, A., Auffray, J.-C., Bonhomme, F., 2002. Two deeply divergent mitochondrial clades in the wild mouse Mus macedonicus reveal multiple glacial refuges south of Caucasus. Heredity 89, 353–357. Pokryszko, B.M., Cameron, R.A.D., Mumladze, L., Tarkhnishvili, D., 2011. Forest snail faunas from Georgian Transcaucasia: patterns of diversity in a Pleistocene refugium. Biol. J. Linn. Soc. 102, 239–250. Pompanon, F., Bonin, A., Bellemain, E., Taberlet, P., 2005. Genotyping errors: causes consequences and solutions. Nat. Rev. Genet. 6, 847–859. Popov, S.V., Rögl, F., Rozanov, A.Y., Steininger, F.F., Shcherba, I.G., Kovac, M. (Eds.), 2004. Lithological-Paleogeographic Maps of Paratethys. 10 Maps Late Eocene to Pliocene. Courier Forschungsinstut Senckenberg, vol. 250, pp. 1–46 (10 maps). Pritchard, J.K., Stephens, M., Donnelly, P., 2000. Inference of population structure using multilocus genotype data. Genetics 155, 945–959. R Development Core Team, 2013. R: A language and environment for statistical computing, version 3.0.1. R Foundation for Statistical Computing, Vienna. Available via: . Reeves, P.A., Richards, C.M., 2007. Distinguishing terminal monophyletic groups from reticulate taxa: performance of phenetic, tree-based, and network procedures. Syst. Biol. 56, 302–320. Ronquist, F., Teslenko, M., van der Mark, P., Ayres, D.L., Darling, A., Höhna, S., Larget, B., Liu, L., Suchard, M.A., Huelsenbeck, J.P., 2012. MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space. Syst. Biol. 61, 539–542. Rosenberg, N.A., 2004. Distruct: a program for the graphical display of population structure. Mol. Ecol. Notes 4, 137–138. Sauer, J., Oldeland, J., Hausdorf, B., 2013. Continuing fragmentation of a widespread species by geographical barriers as initial step in a land snail radiation on Crete. PLoS ONE 8, e62569. Scheel, B.M., Hausdorf, B., 2012. Survival and differentiation of subspecies of the land snail Charpentieria itala in mountain refuges in the Southern Alps. Mol. Ecol. 21, 3794–3808. Schileyko, A.A., 1978. Fauna SSSR. Molljuski. Tom. III, vyp. 6. Nazemnye molljuski nadsemejstva Helicoidea. Nauka, Leningrad. Seddon, J.M., Santucci, F., Reeve, N., Hewitt, G.M., 2002. Caucasus Mountains divide postulated postglacial colonization routes in the white-breasted hedgehog, Erinaceus concolor. J. Evol. Biol. 15, 463–467. Seehausen, O., 2004. Hybridization and adaptive radiation. Trends Ecol. Evol. 19, 198–207. Sokolov, E.P., 2000. An improved method for DNA isolation from mucopolysaccharide-rich molluscan tissues. J. Molluscan Stud. 66, 573–575. Steklov, A.A., 1966. Terrestrial Neogene Mollusks of Ciscaucasia and their Stratigraphic Importance. Nauka, Moscow. Stucki, S., 2014. Développement d’outils de géo-calcul haute performance pour l’identification de régions du génome potentiellement soumises à la sélection naturelle: analyse spatiale de la diversité de panels de polymorphismes nucléotidiques à haute densité (800k) chez Bos taurus et B. indicus en Ouganda. Ph.D thesis, École Polytechnique Fédérale de Lausanne, Lausanne. Sukumaran, J., Holder, M.T., 2010. Dendropy: a Python library for phylogenetic computing. Bioinformatics 26, 1569–1571.

Swofford, D.L., 2002. PAUP⁄. Phylogenetic Analysis Using Parsimony ⁄and other methods. Version 4.0b10. Sinauer Associates, Sunderland MA. Tajima, F., 1989. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics 123, 585–595. Tamura, K., Stecher, G., Peterson, D., Filipski, A., Kumar, S., 2013. MEGA6: molecular evolutionary genetics analysis version 6.0. Mol. Biol. Evol. 30, 2725–2729. Tarkhnishvili, D., 2014. Historical Biogeography of the Caucasus. Nova Science Publishers, New York. Tarkhnishvili, D., Gavashelishvili, A., Mumladze, L., 2012. Palaeoclimatic models help to understand current distribution of Caucasian forest species. Biol. J. Linn. Soc. 105, 231–248. Tarkhnishvili, D.N., Thorpe, R.S., Arntzen, J.W., 2000. Pre-Pleistocene refugia and differentiation between populations of the Caucasian salamander (Mertensiella caucasica). Mol. Phylogenet. Evol. 14, 414–422. van Andel, T.H., Tzedakis, P.C., 1996. Palaeolithic landscapes of Europe and environs 150,000–25,000 years ago: an overview. Quatern. Sci. Rev. 15, 481–500. Vincent, S.J., Morton, A.C., Carter, A., Gibbs, S., Barabadze, T.G., 2007. Oligocene uplift of the Western Greater Caucasus: an effect of initial Arabia–Eurasia collision. Terra Nova 19, 160–166. Volkova, P.A., Schanzer, I.A., Meschersky, I.V., 2013. Colour polymorphism in common primrose (Primula vulgaris Huds.): many colours–many species? Plant Syst. Evol. 299, 1075–1087. Vos, P., Hogers, R., Bleeker, M., Reijans, M., van de Lee, T., Hornes, M., Frijters, A., Pot, J., Peleman, J., Kuiper, M., Zabeau, M., 1995. AFLP: a new technique for DNA fingerprinting. Nucleic Acids Res. 23, 4407–4414. Walther, F., Kijashko, P.V., Harutyunova, L., Mumladze, L., Neiber, M.T., Hausdorf, B., 2014. Biogeography of the land snails of the Caucasus region. Tentacle 22, 3–5. Walther, F., Neiber, M.T., Hausdorf, B., 2015. Systematic revision and molecular phylogeny of the land snail genus Fruticocampylaea (Gastropoda: Hygromiidae) from the Caucasus region. Syst. Biodiversity (in press). Wielstra, B., Crnobrnja-Isailovic, J., Litvinchuk, S.N., Reijnen, B., Skidmore, A.K., Sotiropoulos, K., Toxopeus, A.G., Tzankov, N., Vukov, T., Arntzen, J.W., 2013. Tracing glacial refugia of Triturus newts based on mitochondrial DNA phylogeography and species distribution modelling. Frontiers Zoology 10, 13. Yanina, T.A., 2014. The Ponto-Caspian region: environmental consequences of climate change during the Late Pleistocene. Quatern. Int. 345, 88–99. Zazanashvili, N., Sanadiradze, G., Bukhnikashvili, A., Kandaurov, A., Tarkhnishvili, D., 2004. Caucasus. In: Mittermeier, R.A., Gil, P.G., Hoffmann, M., Pilgrim, J., Brooks, T., Mittermeier, C.G., Lamoreux, J., da Fonseca, G.A.B. (Eds.), Hotspots Revisited: Earth’s Biologically Richest and Most Endangered Terrestrial Ecoregions. CEMEX, Mexico City, pp. 148–153. Zinenko, O., Stümpel, N., Mazanaeva, L., Bakiev, A., Shiryaev, K., Pavlov, A., Kotenko, T., Kukushkin, O., Chikin, Y., Duisebayeva, T., Nilson, G., Orlov, N.L., Tuniyev, S., Ananjeva, N.B., Murphy, R.W., Joger, U., 2015. Mitochondrial phylogeny shows multiple independent ecological transitions and northern dispersion despite of Pleistocene glaciations in meadow and steppe vipers (Vipera ursinii and Vipera renardi). Mol. Phylogenet. Evol. 84, 85–100. Zwickl, D.J., 2006. Genetic algorithm approaches for the phylogenetic analysis of large biological sequence datasets under the maximum likelihood criterion. Ph.D thesis, University of Texas at Austin, Austin.