The Eastern Mediterranean land snails Levantina display a disjunct distribution spanning the Middle East (Levant), Cyprus, few locations along the Aegean Turkish coast between Bodrum and Datça and on the islands of Rhodes, Karpathos and a few surrounding islets (Dodecanese). These land snails are strictly bound to limestone; shell variability is noticeable with a pair of umbilicate and non-umbilicate forms parapatrically distributed in the Levant and along the Aegean Turkish coast; they overlap on the Dodecanese islands. We sequenced fragments of two mitochondrial genes (Cytochrome Oxidase I and 16S rRNA) from the historical Levantina materials available at the Museums of Hamburg and Berlin. The aim of the study is to explain the current distribution of Levantina in the Eastern Mediterranean in light of an earlier hypothesis suggesting anthropochory due to the movements of Crusaders across the area. The deeper nodes in our phylogeny indicate that Levantina reached the Dodecanese from continental Turkey during the Pliocene exploiting continuity of landmasses. In five circumstances the same haplotype co-occurs on two different islands; one haplotype is shared between one island (Rhodes) and the Levant. We suggest that the movements of Crusaders likely explain the current distribution of haplotypes. In particular, the Knights Hospitaller of St. John occupied Cyprus, the Dodecanese and the facing Turkish coasts for more than two centuries (1306–1522) after they withdrew from Jerusalem in 1187 and from the Levant in 1291. Snails could have been introduced as an item of food or transported with other material including limestone used for building.

phylogeography, paleogeography, mitochondrial DNA

In a study that appeared in the journal “Nature” on the 6^th of April 1882, just a few days before his death, Charles Darwin returned to his lifelong fascination with mechanisms of passive long-distance dispersal in mollusks (Darwin 1882). Ever since Darwin’s pioneer work, land snails occurring on islands have puzzled evolutionary biologists; it is indeed challenging to explain how these fragile and slow moving creatures could travel long distances across unsuitable ecological areas if not by passive transport. The advent of molecular tools in biogeography has proved Darwin’s intuition right. Recently, various studies convincingly demonstrated that passive dispersal by birds and humans is the most likely explanation for the current distribution of insular genetic lineages in land snails of the genera Balea and Cepaea (Gittenberger et al. 2006; Grindon and Davison 2013).

In the Eastern Mediterranean area, Levantina shows notable morphological shell variation and a geographic distribution that has been baffling biogeographers for decades (Glaubrecht 1993a, 1993b, 1995 and reference therein). This helicid genus includes rather large edible land snails strictly bound to limestone and displays a disjunct range, represented by two species each in the Middle East (Levant), on Cyprus, in a few locations along the Aegean Turkish coast between Bodrum and Datça and on Rhodes, Karpathos, and a few surrounding islets of the Dodecanese (see Glaubrecht 1993a for further details). Shell variability resides in the lack or presence of an umbilical opening in the adults (Glaubrecht 1993a, 1993b, 1995). In the Levant the two species differentiated by their shell form are distributed parapatrically (the non-umbilicate L. caesareana in the north and the umbilicate L. hierosolyma in the south) (Heller 1979). Heller (1979) identified a narrow hybrid zone between the two species. A similar pattern of parapatric distribution (but with no areas of hybridization identified thus far) is observed over a much-reduced geographic scale along the Turkish coasts. On the Dodecanese islands the umbilicate spiriplana is replaced in places by the non-umbilicate malziana, but this is not always the case. On Rhodes, the malziana type is widespread while the spiriplana type is restricted to a few places, among them the fortress of the Knights Hospitaller of St. John on the northern tip of the island (Glaubrecht 1993a), while the two species are distributed parapatrically on Karpathos and form mixed populations on Symi and surrounding islets (Glaubrecht 1993a).

Based on the above geographic occurrence of shell forms and on field observations, two alternative scenarios have been hypothesized to explain such a peculiar insular distribution (Glaubrecht 1993a, 1993b). Following the paleogeographic history of the area, the Dodecanese could have been colonized during one of the multiple changes in sea level that created transitory land bridges between the islands and the continent from the late Miocene to the middle Pliocene. Similar scenarios have been validated on molecular grounds in a variety of scarcely vagile groups, including land snails (Ketmaier et al. 2006; Jesse et al. 2011; Kotsakiozi et al. 2012 and references therein). An alternative scenario argues for an historical, anthropogenic introduction of Levantina spiriplana on the Dodecanese islands. This hypothesis relies on the observation that most of the locations where the taxon is found – and in particular where the two shell forms co-occur – are actually historical human settlements, i.e. the castles of the Knights Hospitaller of St. John on Rhodes and in Bodrum (Turkey). It was suggested that the movements of Crusaders in the area could explain this distribution pattern (Glaubrecht 1993a, 1993b). Following their withdrawal from the Levant in 1291 after the fall of Acre (Israel), Crusaders conquered Rhodes in 1306 and kept the area under their control until the fall of the Rhodes Grand Master’s Palace in 1522 (for the history of the crusaders, see e.g. Mayer 2005; Murray 2006). For more than two centuries the Dodecanese – and in particular the island of Rhodes with the headquarter of the Knights Hospitaller of St. John – became a key area of occupation and activities, with frequent passages of ships sailing between the islands and the Turkish mainland, transporting Crusaders, the goods they traded and materials used for building.

Here we tested these two alternative scenarios by using sequences of two mitochondrial DNA (mtDNA) genes and the Levantina material available in the collections of the Zoological Museum Hamburg and the Museum für Naturkunde Berlin (Germany); this material was collected about two decades ago in a first attempt to verify the above biogeographical hypotheses (Glaubrecht 1993a, 1993b, 1995). The predictions for the above alternative hypotheses in terms of patterns of genetic structuring are the following. The paleogeographic scenario would translate into discrete amount of genetic divergence and the identification of allopatric haplogroups. In contrast, if the Crusades scenario were true, given the short (in evolutionary terms) time elapsed since the alleged translocation event(s) we would find very little (if any) genetic divergence between currently allopatric haplotypes.

The study is based on material collected by MG in April 1989 and May 1990 (Glaubrecht 1993a, 1993b); samples were preserved in 70%–95% ethanol and vouchers deposited in the Zoological Museum Hamburg (ZMH) with some additional material in the Museum für Naturkunde Berlin (ZMB). Details on sampling localities are given in Table 1. Fragments of the mitochondrial Cytochrome Oxidase I (COI) and 16S rRNA (16S) genes were amplified by Polymerase Chain Reaction (PCR) and directly sequenced as detailed in Ketmaier et al. (2006, 2010); given the age of the samples for about 50% of them we were unable to amplify directly the whole PCR fragments due to DNA degradation. For these samples, we used a nested PCR approach with internal primers designed anew for this study. Precautions were taken to avoid any cross-contamination of samples as detailed in Ketmaier et al. (2006). The first round of PCR amplification was followed by a second round of PCR cycling, using as template the undiluted PCR products from the first PCR amplification. PCR primers used for this second round were designed on the basis of Levantina sequences obtained in the first (not nested) PCR rounds. The primer pair for the nested PCR of the COI gene is (COI-Lev-for: 5’- TTGTAACWGCYCATGCATTTG- 3’ and COI-Lev-rev: 5’- AACTWATTCCAGGAGATCGYA- 3’). The primer pair for the nested PCR of the 16S gene is (16S-Lev-for: 5’- CCCTGACTGTGCAAAGGTAGC- 3’ and 16S-Lev-rev: 5’- GGCCCTAATCCAACATCGAGGTC- 3’). Nested PCR cycling conditions were as in Ketmaier et al. (2006, 2010).

Sequences were edited and aligned in SEQUENCHER 4.1 (Gene Code Corporation, Ann Arbor, MI, USA); the alignment was further checked by eye. We included a selection of Codringtonia sequences from (Kotsakiozi et al. 2012) for the molecular clock analyses and Cornu aspersum to root the phylogenetic searches (Table 1). Phylogenetic searches were carried out by Maximum Likelihood (ML) (run online on the Mobyle portal with the phyml option) and by Bayesian analyses (MrBAYES 3.2) (Ronquist et al. 2011). For both analyses the most appropriate model of sequence evolution was selected using JMODELTEST 2 (Darriba et al. 2012). JMODELTEST 2 returned the HKY+I+G as the best model fitting the concatenated dataset (proportion of invariable sites I=0.575; Gamma distribution shape parameter = 0.576). These settings were then adopted in the Maximum Likelihood, Bayesian and Beast searches.

The robustness of the ML hypothesis was tested by 1,000 bootstrap replicates; MrBAYES was run two times independently for 2,000,000 generations with a sampling frequency of 100 generations. We ran one cold and three heated Markov chains and two independent runs. To establish if the Markov chains had reached stationarity, we plotted the likelihood scores of the sampled trees against generation time. Trees generated before stationarity were discarded as burn-in (first 10% of the sampled trees) and posterior probability values for each node were calculated on the basis of the remaining 90% of sampled trees. We applied coalescence as implemented in the BEAST 1.7.2 package (Drummond and Rambaut 2007) to estimate divergence times in million years (Myr) for the supported clades found by the phylogenetic analyses. BEAST was used to estimate node ages to the most common recent ancestor (TMCRA) of the splits and substitution rates using an uncorrelated lognormal relaxed clock with a Yule or ‘pure birth’ prior process to model speciation. The output of each independent run was visualized using TRACER v1.4. Samples from both independent runs were then pooled after removing the first 10% as burn-in using LogCombiner 1.4.8. After an optimization step during which parameters were calculated to reach an optimum performance and achieve a reasonable effective sampling size (ESS, number of independent samples of the posterior distribution) for the parameters of interest, we carried out two independent runs of 30 million generations each, using a Yule tree prior and the default options for all other prior and operator settings. The age of the basal split in the geographically and phylogenetically closely related genus Codringtonia, estimated at 4.4 Myr in Kotsakiozi et al. (2012), was used to calibrate the tree with a normal distribution for the prior. Finally, we conducted a Bayesian binary MCMC dispersal-vicariance analysis in RASP (Yu et al. 2011). Analyses were run with the maximum units of areas allowed in ancestral nodes equal to three; other parameters were kept at the default settings.

The final alignment including all samples amplified by nested PCR was 482 base pair (bp) long, with 281 bp for COI and 201 for 16S and defined a total of 33 unique haplotypes (GenBank accession numbers KR080942-KR081055; additional sequences used for the phylogenetic searches are JQ239955, JQ239934, JQ239967, JQ239977, JQ240123, JQ240103, JQ240134, JQ240145, HQ203051 and EU912763). These haplotypes are robustly clustered in a monophyletic clade in all phylogenetic searches (Fig. 1). The two Assyriella taxa, viz. guttata and mardinensis, and the Levantina clade were found not to be reciprocally monophyletic and the placement of Assyriella is not supported statistically.

Figure 1. 

Evolutionary relationships in Levantina. Numbers at nodes are statistical support for the ML and Bayesian searches (first and second value above branches). Numbers below branches are age estimates in millions of years with the 95% highest posterior density (HPD) credibility interval in parentheses. Age estimates in bold are discussed in details in the text. Haplotype numbering is as in Table 1. The distribution of each haplotype and the relative shell shape is summarized in the column to the right of the haplotype identifiers (K = Karpathos Is.; R = Rhodes Is.; S = Symi Is.; N = Nimos Is.; CT = Continental Turkey; IS = Israel). Pictures illustrate how shell variability (closed or open umbilicus; squares and circles, respectively) is distributed in Levantina and Assyriella.

Levantina hierosolyma, L. caesareana and the two L. spiriplana forms from the Dodecanese traditionally considered as subspecies (i.e. spiriplana spiriplana and spiriplana malziana) (Glaubrecht 1993a, 1993b) were not retrieved as separate clades. This claims for a revision of the systematics of the group, which is outside the scopes of this study and would require a denser sampling in terms of taxa, geographic coverage and molecular markers. Such a larger study would also allow alternative hypotheses such as incomplete lineage sorting and introgression to be tested in details. We found strong support for three main haplogroups, in spite of the short mtDNA fragments we were able to amplify and sequence. Haplotypes 12, 16, 31 and 33 are from non-umbilicate shelled populations distributed in Turkey, Israel as well as on Symi and their last common ancestor is not older than about 3 Myr old. A large group not older than 3 Myr clusters populations with umbilicate shells distributed over a broad geographic area from Israel, Turkey and the Dodecanese (Karpathos, Rhodes and Symi Islands); haplotypes 23 and 27 are present on Rhodes and Symi Islands and on Rhodes and in Israel, respectively. The remaining haplotypes are gathered in a clade with non-umbilicate shells only; this cluster is limited to the Dodecanese and dates back to about 2.6 Myr ago. Within this clade, in three circumstances the same haplotype is shared between two islands. The RASP analysis (Fig. 2) postulates an origin of the Levantina / Asyriella clade in Turkey (node I; occurrence of this range 92.2%; 80% marginal probability) followed by two dispersal events to the Dodecanese (node II) and the Levant (node IV), respectively. Node V grouping umbilicate populations implies a back dispersal from the Levant to the island of Rhodes followed by multiple dispersal events across Rhodes, Karpathos and Symi. The distribution of haplotypes 19, 27 and 30 supports a close association of umbilicate forms in the Levant, continental Turkey and Rhodes where the umbilicate form (spiriplana) exclusively occurs in the vicinity of the Crusaders’ fortresses (node VII; 82% marginal probability). The lineage with a covered umbilicus (malziana) originated on Rhodes (node VIII; occurrence of this range 93.7%; 79% marginal probability) and progressively spread and diversified across Rhodes, Karpathos and Symi. Haplotypes found on Symi and on the neighboring islet of Nimos (haplotypes 13, 14, 24 and 15) reached the islands via two independent dispersal events from Rhodes (nodes XI and XIII).

Figure 2. 

Historical biogeography in Levantina. On the left is the cladogram (as in Fig. 1 but pruned of the outgroup taxa) summarizing the Bayesian dispersal – vicariance analysis. The distribution of each haplotype and the relative shell shape is summarized in the column to the right of the haplotype identifiers (K = Karpathos Is.; R = Rhodes Is.; S = Symi Is.; N = Nimos Is.; CT = Continental Turkey; IS = Israel). Pie charts and numbers next to them indicate marginal probabilities of alternative ancestral ranges; colors identify the different geographic areas considered and match those in Fig. 1. Roman numbers identify events discussed in the text. On the right is the schematic of the proposed biogeographic history of Levantina. Arrows indicate the direction of the dispersal events inferred by the Bayesian dispersal – vicariance analysis and discussed in the text; roman numbers are the same as in the cladogram shown on the left. The bottom left panel details events within the umbilicate clade (circles), the bottom right panel those within the insular non-umbilicate clade (squares).

We could not retrieve two reciprocally monophyletic clades for Levantina and Assyriella. The validity of the distinction between Assyriella and Levantina was already questioned (Glaubrecht 1995); this is an issue in need of further work based on a taxonomically more exhaustive sampling. Our results do not fully embrace neither the paleogeographic nor the Crusaders hypothesis but rather suggest that both played a role in shaping the observed mtDNA diversity in Levantina. The umbilicate and non-umbilicate shell forms, although clearly distinguishable morphologically (Glaubrecht 1993b, 1995), do not cluster in two reciprocally monophyletic lineages. #The umbilicate populations are restricted to a single and highly supported monophyletic clade whereas the non-umbilicate shell is displaced as a paraphyletic trait in the tree of Fig. 1. It is also worth emphasizing that in our phylogeny only the non-umbilicate shell form is associated to both basal and terminal branches.

The biogeographic reconstruction presented in Fig. 2 supports multiple long-distance and over-sea dispersal events; some of these events are phylogenetically very recent (Figs 1 and 2). The scenario we detail here requires particular caution in its interpretation for the following reasons. First, inferences are based on mtDNA only, a single genetic locus and hence do not give any indication on the geographic origin of the rest of the genome. Second, the uncertainty of the tree topology is not taken into account. Nonetheless, in our opinion it represents a likely scenario given the phylogenetic and geographic information available and the ecology of the group. The origin of the whole Levantina clade, with the inclusion of one of the “Assyriella” taxa, should be sought in continental Turkey (node I; Fig. 2). Nodes V and VIII are roughly coeval (Fig. 2; 3.03 – 2.64 Myr, ages in bold in Fig. 1) and group only insular non-umbilicate populations (malziana) on one hand and almost exclusively insular umbilicate populations (spiriplana) on the other. The exceptions to insularity in the umbilicate clade will be discussed later on; we anticipate that these are, in our opinion, due to anthropogenic translocations. The above time estimates perfectly overlap with the last connection of Karpathos and Rhodes to the mainland, which dates back to 3.5 – 2.8 Myr ago (Beerli et al. 1996); this suggests that Levantina reached the Dodecanese exploiting continuity of landmasses. This, however, did not happen via a single colonization event but rather through a twofold process, which brought on Rhodes two immigrant lineages of independent origin. The reconstruction of Fig. 2 identifies a first dispersal event from continental Turkey to Rhodes (nodes II), which resulted in the insular non-umbilicate clade widespread there (malziana). A second wave of colonization arrived on Rhodes following an earlier dispersal from Turkey to the Levant (nodes IV and V); this clade had the latter as the most likely area of origin and marked the appearance of the umbilicate shell type (spiriplana). Subsequently, both lineages colonized Karpathos and Symi (nodes VI-XIII); the non-umbilicate lineage realized more across-island movements (nodes VIII-XIII) than the umbilicate lineage (nodes VI-VII). How and when Levantina dispersed from Turkey to the Levant and why the umbilicate shell type is mostly insular are questions that remain to be answered. Recent evidence suggests that shell shape diversification in Mediterranean land snails might be due to the interplay between historical and ecological factors (Fiorentino et al. 2013). The positioning of haplotypes from the Levant at the base and within the otherwise poorly differentiated umbilicate clade suggests that the occurrences of this form in the Dodecanese and in southwestern Turkey is probably the result of (anthropogenic?) introductions. Such a shallow level of differentiation in the umbilicate clade is not mirrored in the non-umbilicate clade.

Even though the hypothesis of the non-umbilicate clade being of insular origin is appealing from an evolutionary perspective, we should not overlook the fact that this result could be an artifact due our limited sampling. In particular, we cannot completely rule out the hypothesis that the non-umbilicate clade originated in the Levant and subsequently reached the Dodecanese. We indeed identified a non-umbilicate clade grouping haplotypes from both the Levant and the Dodecanese (haplotypes 12, 16, 31, and 33). The non-umbilicate samples from the Levant are genetically close to a few samples from Symi. The other non-umbilicate individuals from that island are spread in the upper part of the tree of Fig. 1, scattered among haplotypes found on other islands. In order to test this hypothesis adequately we would have needed a better survey of the genetic variation in Levantina from the Levant that is, on the contrary, too limited to address the issue.

Rhodes, the largest island of the Dodecanese and the closest to the mainland, is identified as the first colonized by the two Levantina lineages (nodes V and VIII). The non-umbilicate form is widespread across the island while umbilicate-shelled individuals are restricted to the fortress of the Knights Hospitaller of St. John in the northern part of island (city of Rhodes) (Glaubrecht 1993b, 1995). This has been taken as an indication that humans introduced this shell type on the island (Glaubrecht 1993b, 1995). Our molecular time estimates reject this hypothesis, at least for the basal haplotypes involved in the event (haplotypes 11 and 22). In our phylogeny, however, we found five cases of haplotypes shared between islands or between islands and distant locations. These distributions are difficult to justify on natural bases; the lack of any genetic differentiation indicates very recent dispersal in spite of the intervening marine barriers. Within the umbilicate clade, node VII suggests a recent dispersal event across Rhodes, continental Turkey and Levant; the same haplotype 27 co-occurs on Rhodes and in the Levant. Similarly, haplotype 23 is present on Rhodes and Symi. Haplotypes 1, 2, 3 and 5 (all non-umbilicate) are shared between two islands (alternatively Karpathos and Rhodes or Karpathos and Symi); most interestingly, on Rhodes these haplotypes are confined to human settlements (Table 1). How could we explain the allopatric occurrence of the same (or slightly diverging) haplotypes currently separated by geographic barriers theoretically insurmountable by land snails? Over-sea passive dispersal through i.e. hitchhiking on birds and/or surviving the passage through birds’ gut has been documented in land snails (Gittenberger et al. 2006; Miura et al. 2012). It should also considered that large helicids are edible and they have been found associated with human settlements in many Mediterranean archeological sites (Grindon and Davison 2013); furthermore, the limestone they dwell upon had been used for a long time as a building material and transported in large blocks along historical trading routes across the Mediterranean Sea (Fiorentino et al. 2008). The geographic distributions of haplotypes 27 and 30 (Rhodes, Israel and continental Turkey) and that of the closely related haplotypes 17 (Karpathos) and 29 (Israel) are particular striking and advocates for bringing the Knights Hospitaller of St. John back into the play. After the rising power of Islam expelled in 1291 the Knights from Jerusalem (where haplotypes 27 and 29 are found), they conquered Rhodes (where we found haplotype 27), the neighboring islands (i.e. Karpathos with haplotype 17) and the coast of Anatolia (haplotype 30). The Knights kept the Dodecanese and the Anatolian ports facing it under their control for 200 years before being defeated in 1522 by Sultan Suleiman the Magnificent and forced to withdraw to Malta (Mayer 2005; Murray 2006). During this period of time, they built anew or fortified with huge walls the already existing castles. In particular, they fortified the city of Rhodes with the Palace of the Grand Master where one of the two only umbilicate (spiriplana) populations of the island are found (Glaubrecht 1993a, 1993b); haplotype 27 (co-occurring on the island and Israel) is indeed carried by land snails collected on the rampart of one of the city gates (D’Amboise gate; Table 1). It is thus not difficult to envision these historical events as being responsible for the unusual geographic distribution of some of the mtDNA haplotypes we identified in the study (Figs 1 and 2). It would be interesting to expand this study to include samples of Levantina from Cyprus; the island was never connected to the Anatolian mainland but served as the first stronghold or retreat of the Crusaders after they had to leave the Levant following the fall of Acre in 1291.

We are aware that the scenario presented in here – although fascinating – is not the only possible one. Due to the sub-optimal quality of most of the samples at our disposal, we were able to sequence short gene fragments. This implies that we could have easily missed out on rare genetic variants. In addition, Örstan (2004) already suggested that three disjunctive records of L. spiriplana in western Turkey, just north of the region considered of this present study, could be due to introduction by humans during the Ionian period, perhaps on ballast rocks. Similarly, Welter-Schultes (1998) suggested that some Albinaria species that are found aestivating on rocks in Crete might have been carried to places outside their natural ranges on rocks used for construction or as ballast in ships. More recently, the same author (Welter-Schultes 2008) provided direct evidence that land snails have been carried on ships for more than 3,000 years in the Mediterranean area by describing shells discovered in the underwater archaeological excavations of a Late Bronze Age (3,300 years BP) shipwreck at a Southern Turkey location. Däumer et al (2012 and references therein) revealed a complex scenario in the invasive land snail Theba pisana pisana suggesting that primarily human activities rather than natural processes have shaped (and still are) the distribution of the taxon. The authors also suggested that different lineages identified on genetic bases only could have different adaptive and invasive potentials, unveiling a complex scenario where different forces at different levels (from the ecological to the genomic one) could come into play.

The data presented in here, along with the similar evidence existing for the area mentioned in the previous paragraph, suggest that two different layers of complexity (natural colonization vs. historical human activities) should be considered when addressing puzzling distributions in an area interested by intense human activities since historical times. Also, this study represents a starting point for further investigations based on a more extensive sampling in terms of geographic and taxon coverage as well as molecular markers.

We wish to thank Binia De Cahsan for producing the shell drawings in Figs 1 and 2; Annabell Szymanski and Nora Mass have helped with retrieving the museum materials. Martin Haase and Bernhard Hausdorf provided valuable criticisms.