New genetic data reveals a new species of Zospeum in Bosnia (Gastropoda, Ellobioidea, Carychiinae)

Abstract 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 Zospeumamoenum 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 Zospeumtroglobalcanicum. The second lineage derives from the municipality of Tomislavgrad, Bosnia-Herzegovina and is herein, described a new species: Zospeumsimplex Inäbnit, Jochum & Neubert, sp. nov.


Materials and methods
Material is housed in the following collections: 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 . Each specimen was inserted into a 0.2-ml PCR-tube and dried at room temperature. Eight μl ddH 2 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.
The PCR-solution included the following admixture: 2 μl template, 12.5 μl Go-Taq (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 H3-R 5'-ATATCCTT(AGGGCAT(AG)AT(AG)GTG-3' Colgan et al. (1998) 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;  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.

Phylogenetic trees
Both the ML and the BI trees (see Fig. 2 for the latter) are more or less identical. The specimens sequenced in this study clustered with Z. pretneri, Z. tholussum, and Z. manitaense. In both trees they form a badly supported monophyletic group that splits again into two groups in accordance with their geographical distribution (see Fig. 1) and could be separated at the species level: the two specimens from the region of Dubrovnik, Croatia (Špilja Jezero; referred to as Z. aff. troglobalcanicum), and the remaining specimens from Bosnia-Herzegovina (Jama u kamenolomu, Vranjača, Jama Dobravljevac; described as Z. simplex sp. nov. herein). The latter group is not supported in either tree but recovered in both. An additional specimen (NMBE 568054, Špilja Dahna), from which we were only able to amplify H3, didn't cluster with any species within the Z. pretneri group. The two groups were also recovered, though here with high node support, in the additional ML tree (Suplementary tree 1) calculated for the Z. pretneri group. The Austrian specimens from Kruckenhauser et al. (2019) form a strongly supported monophyletic group within Z. amoenum.

Divergences
For most markers, intraspecific divergences among the species in the Z. pretneri group are clearly smaller than the interspecific divergences (Table 3). This indicates that these species comprise separate lineages, especially the specimens classified as Z. aff. troglobalcanicum and those collected in Bosnia (henceforth referred to as Z. simplex sp. nov.), which were not included in previous genetic studies (see Inäbnit et al. 2019).
Zospeum amoenum shows a high intraspecific divergence when compared to other members of the Z. alpestre group (see Table 4), though other species (such as Z. aff. troglobalcanicum, see Table 3) show similarly high intraspecific divergence. When the Austrian populations from Kruckenhauser et al. (2019) are aligned within Z. amoe- Table 3. The number of base substitutions per site from averaging over all sequence pairs within (withinspecies divergences) and between (between-species divergences) species within the Z. pretneri group. Results shown for each marker separately. Between-species distances are listed below the black, empty boxes, the Standard errors above.

Species
No. of sequences than that amidst the other species within the Z. alpestre group, but still higher than the within-group divergence in both Z. amoenum and the Austrian specimens.

Automatic Barcode Gap Discovery (ABGD)
The ABGD run on the Z. pretneri-group COI alignment yielded two different possible subdivision schemes: one where the alignment was subdivided into five groups (five groups scheme; prior maximal distance P = 7.74e -03 ; barcode gap distance: 0.043) and a second where the alignment was subdivided into seven groups (seven groups scheme; prior maximal distance P = 4.64e -03 ; barcode gap distance: 0.003). Both subdivision schemes considered the previously published sequences of Z. pretneri, Z. tholussum, and Z. manitaense as separate groups. The five-group scheme separated the individuals sequenced in this study into a Croatian group (Špilja Jezero) and a Bosnian group (Jama Dobravljevac, Jama u kamenolomu, Vranjača), while the seven-group scheme separated those individuals into two Croatian groups (one for each of the two specimens from Špilja Jezero) and two Bosnian groups (1: specimens from Jama u kamenolomu; 2: specimens from Jama Dobravljevac and Vranjača).
Differing from Z. pretneri and Z. tholussum by its broader shell and the differentiated parietal shield; differs from Z. manitaense by the absence of a visible parietalis in the aperture; barely differs from Z. aff. troglobalcanicum morphologically, on average with reduced shell broadness and a slightly deeper suture (see Remarks).
Etymology. Named simplex (= simple, unsophisticated) due to the lack of any form of shell sculpture or lamellae.
Remarks. Difficult to separate from Z. troglobalcanicum without genetic data (which is not uncommon in Zospeum; see Inäbnit et al. 2019). Both species have a nondescript shell without prominent shell sculpture or lamellae within the aperture. Absolon's (1916) description of Z. troglobalcanicum consisted out of a photograph depicting multiple specimens haphazardly clustered together in various positions and a legend that established the name and type locality. The lack of a written characterisation of the species in the original description and the fact that the specimens in the photograph weren't depicted in any standardised position makes a characterisation of the species fairly challenging (putative syntype specimen, collected by K. Absolon from the type locality, was only located very recently by AJ in Vienna (NHMW Mol.Coll.Edlauer 32.749) and couldn't be studied yet). From the photograph in Absolon (1916), the species can be characterised as similar to Z. manitaense in shell shape, without any visible lamella in the aperture and with a comparatively large parietal shield. The larger parietal shield might serve as a distinguishing character between Z. simplex and Z. troglobalcanicum, though the illustration of a topotypic specimen in Bole (1974; fig. 3h) might indicate that this character is variable within the population. The two specimens we preliminarily assigned to Z. troglobalcanicum (Fig. 3, NMBE 568052;Inäbnit et al. 2019: fig. 7u) only have a small parietal shield. As of now, the shell height:shell width ratio seems to be the most effective way of separating the two specimens from Z. simplex (Z. simplex: generally higher than 1.3 (one exception); Z. aff. troglobalcanicum: below 1.3), but that might just be due to the low sample sizes. Investigation of the inner aspects of the shells will be presented in a later work.

Discussion
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. material1), 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. material1, 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 withingroup 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).