Recent integrative investigations of the terrestrial ellobiid genus, Zospeum, have revealed significant findings concerning its Alpine-Dinaric evolution and taxonomy. Due to the expected discrepancy between the useful, but limited, 1970s’ classification system based on shell data and the results of recent genetic analyses in the latest investigation, a revision of the entire radiation was undertaken, and a new classification system was devised by the present authors in an earlier paper. Concurrent to this work, molecular sequences from two Austrian caves were published independently of our revision by another research group. By incorporating these genetic data within our phylogenetic framework here, we show that the Austrian individuals are genetically most similar to Zospeum amoenum and consequently, classify them within that species. We additionally reveal two new genetic lineages from the largely under-sampled southern extension of Zospeum’s known distributional range. The first lineage, deriving from the region of Dubrovnik, Croatia, is a potential candidate for genetically clarifying Zospeum troglobalcanicum. The second lineage derives from the municipality of Tomislavgrad, Bosnia-Herzegovina and is herein, described a new species: Zospeum simplex Inäbnit, Jochum & Neubert, sp. nov.

Dinarides, microsnails, molecular phylogenetics, shell variability, subterranean ecology, troglobitic microsnails

The carychiid genus, Zospeum, consists of tiny (0.9–2.6 mm), troglobitic snail species that are distributed in two disjunct areas: a western zone, comprising the western Pyrenees and the Cantabrian mountains of Spain and France (Jochum et al. 2015a, 2019) and an eastern zone, encompassing the southeastern Alps and Dinarides of northeastern Italy, southern Austria, Slovenia, Croatia, Bosnia-Herzegovina and Montenegro (see Inäbnit et al. 2019). This work addresses the species rich, eastern radiation of Zospeum.

Until recently, the eastern radiation of Zospeum was largely classified using a scheme devised by Bole (1974), based solely on shell morphology. More recent studies (Weigand et al. 2011; Weigand et al. 2013; Jochum et al. 2015b), however, found Bole’s (1974) scheme, though effective for its time, now incongruent with genetic data, leading to a thorough revision by Inäbnit et al. (2019). They subdivided the eastern Zospeum radiation into 25 species that could be divided genetically into five species groups: the Z. spelaeum group (northeastern Italy, Slovenia, north-western Croatia; five species), the Z. alpestre group (Slovenian Alps and adjacent regions in Italy and Austria; four species), the Z. obesum group (southwestern Slovenia and adjacent Croatia; two species), the Z. pretneri group (Croatia, more or less close to the Adriatic coast; four species), and the Z. frauenfeldii group (southern Slovenia, northwestern Croatia, northwestern Bosnia-Herzegovina; five species); five species could not be assigned to any of the five groups due to lack of molecular data.

One of the issues raised in Inäbnit et al. (2019) is that Zospeum’s eastern distribution has been unevenly sampled throughout its history. Most studies covered almost only Slovenian (e.g., Frauenfeld 1854, 1856; Freyer 1855; Bole 1974; Weigand et al. 2013), Italian (Pezzoli 1992 and papers cited therein) and northwestern Croatian populations (Slapnik and Ozimec 2004; Inäbnit et al. 2019). The consequence of this sampling disparity is that we have very limited records from southern Croatia, Bosnia-Herzegovina and Montenegro (see Inäbnit et al. 2019: fig. 1a), none of which include genetic data. In fact, the only species described from the southern half of the Zospeum’s distribution range is Zospeum troglobalcanicum Absolon 1916. Shells that obviously belong to different species exist in museum collections (see Inäbnit et al. 2019: fig. 10W-Z; Gittenberger 1975), but genetic data from these southernmost populations is still lacking for a contemporary, integrative taxonomic assessment. In the current study, we add new sequences from 12 specimens, collected in southern Croatia and Bosnia-Herzegovina to the existing genetic dataset.

Approximately the same time as the revision by Inäbnit et al. (2019) was published, Kruckenhauser et al. (2019) published the results of a small barcoding study of specimens from Austria (for locations see Fig. 1). Due to this unfortunate overlap, their results could not be incorporated into the classification system proposed by Inäbnit et al. (2019). We have however, included these results in our work here.

Figure 1. 

Map showing the distribution of the Zospeum pretneri group and the Zospeum alpestre group (except Z. isselianum). Austrian specimens from Kruckenhauser et al. (2019) are labelled as “Z. cf. amoenum”.

Material is housed in the following collections:

AJC Adrienne Jochum Collection, Kelkheim, Germany;

MCSMNH Malacological Collection of the Slovenian Museum of Natural History (former CSR SASA, MZBI & SMNH) Ljubljana, Slovenia;

NHMW Naturhistorisches Museum Wien, Wien, Austria;

NMBE Naturhistorisches Museum der Burgergemeinde Bern, Bern, Switzerland;

RSC Rajko Slapnik Collection, Kamnik, Slovenia;

SMF Senckenberg Forschungsinstitut und Naturmuseum, Frankfurt am Main, Germany.

In order to preserve the shell from dissolution during the extraction, our DNA extraction protocol was based on a method initially described in Schizas et al. (1997) and partially modified after Böttger-Schnack and Machida (2011). DNA extraction was conducted on 12 ethanol-preserved individuals (NMBE 568052-568063). Each specimen was inserted into a 0.2-ml PCR-tube and dried at room temperature. Eight μl ddH₂O and 2 μl 5× PCR-buffer (Promega 5× Colorless GoTaq Reaction Buffer) were added and the mixture was heated at 94 °C for 2 min. whereby 1.3 μl proteinase K solution (from the DNEasy Blood and tissue kit, Qiagen) were then added and the solution was homogenised and then incubated in a PCR-thermocycler at 55 °C for 15 min., afterwards at 70 °C for 10 min. The incubation was repeated once. Ten μl of Gene Releaser (Bioventures Inc.) were then added and the mixture was inserted into a thermocycler with the following protocol: 65 °C for 30 s, 8 °C for 30 s, 65 °C for 1.5 min., 97 °C for 3 min., 8 °C for 1 min., 65 °C for 3 min., 97 °C for 1 min., 65 °C for 1 min., 80 °C for 5 min. The mixture, including the intact shell, was centrifuged for 1 min. using a table centrifuge and the clear phase with the DNA was transferred to another 0.2 mL PCR-tube, where 15 μl of AE-Buffer (DNeasy Kit, Qiagen) was added. The shell was cleaned from the remains of the Gene Releaser chemicals by rinsing with 80% EtOH.

We used five markers, two mitochondrial (COI (658 bp), 16S (483 bp)) and three nuclear markers (H3 (330 bp), ITS2 (809 bp), 28S (590 bp)) with a total length of 2870 bp (for primers, see Table 1).

The PCR-solution included the following admixture: 2 μl template, 12.5 μl GoTaq (Promega) polymerase, 8.5 μl of nuclease-free water, and 1 μl of both forward and reverse primer (10 μmol) respectively. In cases where the PCR signal was judged too weak, the reaction was repeated using 3 μl template DNA, 3 μl of the previous PCR product, and 5.5 μl of nuclease-free water. The amount of GoTaq and primers remained the same. The amplification was conducted using the following cycling protocols: For COI, the admixture was first heated up to 95 °C for 1 min, followed by 30 cycles of 30 s (of denaturation at 95 °C for 30 s, annealing at 52 °C for 30 s, extension at 72 °C for 1 min), and a final extension at 72 °C for 3 min. For 16S, the protocol started with 2:30 min at 90 °C, followed by 10 cycles of 30 s at 92 °C, 30 s at 44 °C, and 40 s at 72 °C, followed again by 30 s at 92 °C, 40 s at 48 °C, and 40 s at 48 °C. The protocol for 28S started with 1 min at 96 °C, then went into 35 cycles of 30 s at 94 °C, 30 s at 50 °C, and 1 min at 72 °C, finishing with 10 min at 72 °C. The ITS2 protocol started with 1 min at 96 °C, followed by 35 cycles of 30 s at 94 °C, 30 s at 44 °C, and 1 min at 72 °C, ending with 10 min at 72 °C. For H3, the admixture was first heated up to 95 °C for 3 min, followed by 40 cycles of 45 s at 94 °C, 45 s at 50 °C, and 2 min at 72 °C, finishing with 10 min at 72 °C. The protocols for COI and H3 could be used for both markers. The PCR products were sequenced at the LGC Genomics GmbH (Berlin, Germany) using their standard protocol.

Sequences received from LGC were imported into the Geneious 5.4.7 software (Kearse et al. 2012). The forward and reverse sequences for each gene and individual were combined and edited. In addition to the sequences that were generated during this study, we used the sequences previously used and generated in Inäbnit et al. (2019), as well as those generated by Kruckenhauser et al. (2019). The name of some of the Spanish specimens were updated based on the results of Jochum et al. (2019). A total list of samples can be found in Table 2. For each marker, sequences were aligned in Geneious using the MAFFT multiple sequence alignment plugin version 1.3.6 (based on MAFFT v7.308; Katoh et al. 2002; Katoh and Standley 2013), allowing the program to choose the most appropriate algorithm. The sequence length of each alignment was standardised to the length mentioned above.

Topologies were estimated using two different phylogenetic methods: Maximum Likelihood (ML) and Bayesian Inference (BI). The five markers were set as partitions in both of these methods, using a distinct model for the third codon in protein-coding genes (COI, H3). The maximum likelihood (ML) topology was estimated using the RAxML 7.2.8 (Stamatakis 2014) plugin of Geneious with the GTR gamma nucleotide model and 1000 bootstrap replicates. An additional ML tree was calculated for the Z. pretneri group (with Z. robustum NMBE 548777 as an outgroup) without H3 and 28S.

The Bayesian tree was reconstructed with MrBayes 3.2.6 (Huelsenbeck and Ronquist 2001) using the substitution models suggested by PartitionFinder (Lanfear et al. 2016, Lanfear et al. 2012, Guindon et al. 2010), a Markov Chain Monte Carlo (MCMC) chain length of 10000000 generations, a subsampling frequency of every 4000 generations, the first 100000 generations were discarded as burn-in, four heated chains and a chain temperature parameter of 0.2. Calculations were performed on the UBELIx (http://www.id.unibe.ch/hpc), the HPC cluster at the University of Bern.

The single gene alignments of COI, 16S, and ITS2 were imported into MEGA X 10.1.7 (Kumar et al. 2018) and the various sequences grouped into species. The average evolutionary divergence between sequence pairs within species (subsequently referred to as within-species divergence) was estimated where possible (only for species with more than one sequence present) using the Maximum Composite Likelihood model (Tamura et al. 2004) on standard settings. The Maximum Composite Likelihood model was also used to estimate the average evolutionary divergence between sequence pairs between species (subsequently referred to as between-species divergence). The focus of the analyses lay on the Z. pretneri group (as defined by Inäbnit et al. 2019; all markers) and the Z. alpestre group (only COI, with the Austrian specimens from Kruckenhauser et al. 2019) classified as separate species or included in Z. amoenum.

Additionally, an Automatic Barcode Gap Discovery (ABGD; Puillandre et al. 2011; https://bioinfo.mnhn.fr/abi/public/abgd/abgdweb.html) analysis was performed on the COI alignments of the Z. pretneri group and of the Z. alpestre group using the default settings (Pmin = 0.001, Pmax = 0.1, Steps = 10, X = 1.5, Nb bins = 20, distance = Jukes-Cantor).

A map (Fig. 1) was constructed using the Natural Earth dataset in QGIS 3.16.3. Most locality data was taken from Inäbnit et al. (2019), and the coordinates for the Austrian sites were taken from Kruckenhauser et al. (2019). Locality data of the specimens sequenced in this study were provided by the various collectors.

The phylogenetic tree reconstructions (Fig. 2) agree mostly with those figured in Inäbnit et al. (2019). The main difference is that the node support values within the Z. pretneri group and in that of Z. amoenum are now fairly low and the topology is different. This can be explained by the high number of new specimens that sometimes are only represented by one marker (especially in Z. amoenum). It should also be noted that our current trees resolve Z. robustum, for which we didn’t have any new specimens, with a significant node support as a monophyletic group (node support was not significant in Inäbnit et al. 2019, but the classification as an independent species could be justified via species delimitation methods). Since its position was not resolved with significant node support in either tree, the specimen from Tonkovića špilja is not included in Z. robustum in this tree, as was the case in Inäbnit et al. (2019). Due to lack of additional material, the classification within Z. robustum remains unchanged in this work.

The 12 Zospeum individuals from Bosnia-Herzegovina and Croatia, are the first to be molecularly assessed from the greatly understudied, southern extension of Zospeum’s distribution. Within the phylogenetic trees (Fig. 2, Suppl. material 1), these specimens form a monophyletic group with a deep split between the two specimens from Croatia and nine of the ten specimens from Bosnia-Herzegovina (the remaining specimen from Špilja Dahna is only represented by a sequence of the conservative histone H3 gene, which doesn’t usually resolve to species level).While recovered in all phylogenetic trees calculated for this work, this arrangement only has high node support values in the Suppl. material 1, which was calculated without the conservative H3 and 28S nuclear markers. This result might indicate that conservative markers may have a destabilising effect on species level phylogeny within this group. Both ABGD schemes support the separation of the Croatian and Bosnia-Herzegovina individuals from each other at species level, though the seven-group scheme further subdivided the specimens from both geographical regions. We prefer to use the five-group scheme for the following reasons here: a) The barcode gap of the seven-group scheme is much lower (0.003) than the barcode gap (0.032) that was detected in the Carychiidae alignment in Weigand et al. (2011), while the barcode gap in the five-group scheme was slightly higher (0.043) than in Weigand et al. (2011); b) both individuals from Croatia (considered separate groups in the seven-group scheme) derive from the same cave and are unambiguously recovered as monophyletic and closely related in all trees, making their status as separate taxa unlikely. The divergence analysis further corroborates the results of the ABGD five-group scheme whereby the between-group divergence between the Croatian and the Bosnian groups (see Table 3) was within the general range of interspecific divergence within the Z. pretneri group. We thus, propose separating the individuals sequenced in this study into two species:

• A species encompassing all ten specimens from Bosnia-Herzegovina. This species is described as Z. simplex sp. nov. above. Since we do not have enough molecular and morphological data for the individual from Špilja Dahna, we cannot confidently place it within Z. simplex right now. However, due to its close geo graphical proximity (less than 1 km) to one of the caves with genetically identified specimens (Jama u kamenolomu), we expect it could well be assignable to Z. simplex as no external morphological inconsistencies separate it from other Z. simplex specimens in our study.
• A species comprising two specimens from Špilja Jezero in the region of Dubrovnik. This locality is fairly close (around 22 km) to the type locality (Benetina pećina) of Z. troglobalcanicum Absolon, 1916. The sequenced specimens do not show any major external morphological differences from the specimen identified as Z. troglobalcanicum (as figured in Bole 1974: fig. 3h) and from those imaged in Inäbnit et al. 2019: fig. 7u), though the adult specimen clearly has a smaller parietal shield than the specimens figured in Absolon (1916). We propose tentatively classifying those specimens within Z. troglobalcanicum until genetic material from the type locality can clarify its status and the morphological investigation of the singular syntype (NHMW Mol.Coll.Edlauer 32.749) of this species can be taxonomically and nomenclaturally clarified in a separate work.

Even if it is not as large as the between-group divergence of other species pairs within the Z. alpestre group, our divergence analysis revealed that the between-group divergence between Z. amoenum and the two Austrian populations is greater than the within-group divergence of either lineage. Our analysis also found that the within-group divergence in Z. amoenum is only slightly increased if the Austrian populations are included within this species. These results agree with the tree reconstruction published in Kruckenhauser et al. (2019), which resolved the Austrian population as the sister group of Z. amoenum. Our trees, as mentioned above, lack the resolution to separate the Austrian populations from Z. amoenum and can thus, not confirm this conclusion. The ABGD scheme for the Z. alpestre group recovers the Austrian population as a separate group from Z. amoenum and splits the latter species into three groups. The barcode gap in this scheme is, however, much lower (0.016) than the one proposed for Carychiidae in Weigand et al. (2011), which was used for species classification within the Z. alpestre group before (e.g., in Weigand et al. 2013). We are thus, reluctant to draw conclusions regarding Z. amoenum and the Austrian specimens from the ABGD scheme. It may indicate some large intraspecific genetic variability within Z. amoenum (with the possibility of the presence of several species) that might coincide with the large morphological variation found in this species (Inäbnit et al. 2019), which would need to be addressed in a separate study with better sampling.

Zospeum amoenum described in Inäbnit et al. (2019) bears either a small parietalis that does not expand within the shell or it is lacking completely. Kruckenhauser et al. (2019) did not figure a specimen in which the configuration of the parietalis within the last whorl could be seen, but Gittenberger (1982) figured one specimen from the Hafnerhöhle (one of the two caves sampled by Kruckenhauser et al. 2019), where the parietalis was exposed. The parietalis of this specimen is slightly broadened three quarters of a whorl into the shell and seems to decrease expansion again further into the shell. Though the syntype of Z. amoenum (see Inäbnit et al. 2019: fig. 6L) shows a similar configuration of the parietalis, it is not congruent with the description of this structure in Z. amoenum assessed in Inäbnit et al. (2019).

Our study suggests that a final species assignment for the two Austrian populations is not possible until further supporting information becomes available. Until then, we classify these two Austrian populations as Z. amoenum, avoiding the now outdated classification of these populations with Z. isselianum (as was done in Kruckenhauser et al. 2019).

We thank Estée Bochud for her help in data retrieval. We also thank Jana Valentinčič for her assistance in the field and Axel Schoenhoffer for providing us with samples for analysis.