Isolation and endemism in subterranean aquatic snails: unexpected case of Montenegrospeum bogici (Pešić et Glöer, 2012) (Gastropoda: Truncatelloidea: Hydrobiidae)

The subterranean aquatic snails may serve as a model of endemism and isolation vs. migration in subterranean habitats. The aim of the present paper is to verify the hypothesis that subterranean aquatic snails can migrate through diverse subterranean habitats, applying four molecular markers as well as a RAPD technique and shell morphometry. They were used to estimate the differences and gene flow between populations of the hydrobiid subterranean aquatic species Montenegrospeum bogici, collected in the Dinaric karst region. Three molecularly distinct taxonomic units were distinguished. The mOTU B was found at single locality, mOTU C at two, but the mOTU A at ten localities, scattered along 236 km distance, at two of them in sympatry with either mOTU B or C. Within mOTU A, the estimated levels of the gene flow were high. The pairwise measures of genetic differentiation were statistically significantly associated with geographic distances between the populations. In general, neither the infinite-island model of interpopulation differentiation, expected for isolated populations, nor the stepping-stone one, but rather the isolation-by-distance model explained the observed pattern. Our results suggest that interstitial habitats provide ways of migration for the stygobiont M. bogici, as has been already suggested for other subterranean gastropods.


Introduction
Of the approximately 20,000 worldwide described species of subterranean animals , stygobionts and troglobionts, there are over 350 described species of stygobiont (obligate subterranean-dwelling aquatic) gastropods, with 97% of them, and over 50 genera representing the superfamily Truncatelloidea (''Hydrobiidae s. lato'' in Bernasconi & Riedel, 1994;Culver, 2012). They are worldwide distributed in all kind of brackish and freshwater habitats, and grouping small or minute (shell height 1-3 mm), dioecious gill-breathing snails, whose phylogenetic relationships are still far from being clear, and whose systematics remains poorly understood. Sket et al. (2004) reported 169 stygobiont gastropod species inhabiting the Balkan Peninsula; 148 of them, representing 41 genera and 6 families, were found in the Western Balkans, 8 (4 truncatelloid genera) in Greece, and 13 (11 truncatelloid genera) in Bulgaria. Among the subterranean gastropods, aquatic species are much more numerous than the terrestrial ones-for example, in the USA stygobionts' species diversity is nearly six times higher than that of troglobionts (Hobbs, 2012;Gladstone et al., 2021). However, despite the obvious existence of many still undescribed species in these not easily accessible and explorable habitats, the number of those already known is overestimated (e.g. Osikowski et al., 2018). As pointed out by Falniowski (2018), using intraspecies variable shells as the only known structures, coupled with the widespread belief among many taxonomists in geographic isolation and unavoidable, immediate speciation of the cave animals has resulted in descriptions of new species (nearly) in each cave or other subterranean habitat (e.g. Reischütz et al., 2008;Georgiev, 2013;Glöer et al., 2015;Grego et al., 2017). The incredibly long list (more than 60 species and many subspecies) of the nominal species of Bythiospeum Bourguignat, 1882 presented by Glöer (2002), may serve as a good example. Even after including detailed, thorough morphological (e.g. shell morphology) and molecular (e.g. allozymes, mitochondrial data) information, the correct species identification in the Truncatelloidea still pose quite a challenge even for experienced taxonomic specialists (Falniowski, unpublished data). Usually one can find varying quantities of empty, often more or less damaged shells, found at the surface, at groundwater-fed springs. Among such shells only a few living individuals were ever found, washed out into springs at times of high flow, especially during the spring season (Haase, 1995;Richling et al., 2016). At a few localities, more living snails can sometimes also be found in springs and streams at the surface. There are also subterranean snails whose shells can only be found deep inside caves and are often incomplete, to say nothing of the extreme scarcity of living specimens. Collection of stygobiont Truncatelloidea is, therefore, difficult, especially in adequate numbers for any study of population genetic structure-subterranean populations are usually not dense, rich in individuals (e.g. Mammola et al., 2021). In fact, nothing but shell morphology is known for the vast majority of the subterranean gastropods. Even the soft parts, if accessible, are not usually informative enough to resolve taxonomic questions, since the animals are very tiny compared to other gastropods, and miniaturisation, coupled with convergent adaptations to freshwater habitats (osmoregulation, internal fertilisation, eggs in capsules, etc.) has resulted in simplification and unification of their anatomy (e.g. Fretter & Graham, 1994;Culver, 2012;Falniowski, 2018). Such incompleteness of the data, coupled with assumptions that they must be endemic, has resulted in strongly biased estimates of biodiversity, to say nothing about systematics, evolution, and ecology of these snails. Molecular data-DNA sequences-are helpful, but there are still only a few studies applying them (e. g. Grego et al., 2019;Hofman et al., 2019). Another limitation concerns the low number of applied molecular markers-in fact often COI only-which is a result of common problems with PCR amplification or lack of polymorphism. Many of the known stygobiont species primarily inhabit the interstitial habitats formed by the alluvial gravel (frequently not dependent on karst); thus, caves and spring are for them a secondary or occasional habitat, so their distribution ranges may be more widespread than originally supposed (Richling et al., 2016).
A relatively high species diversity of the subterranean snails could be expected, as the group is well preadapted to such an environment, as being nocturnal (or at least avoiding direct sunlight), opportunistic, ecologically (especially feeding) generalist animals, with well-developed mechano-sensory and/or chemosensory traits (Trajano & Cobolli, 2012). They are also well adapted to food-limited habitats, since their sizerelated energy requirements are low, being four times lower than in fishes and salamanders, and more than twice lower than in arthropod detritivores (Poulson, 2012). Thus, among the subterranean aquatic invertebrates, Crustacea dominate, but gastropods are the second group in species diversity (Culver & Sket, 2000). Animals adapted to subterranean, especially cave environments are often thought to be highly geographically isolated because of their limited dispersal ability, resulting from limited physiological tolerances and, especially in the case of snails, physical limitations of their locomotion (Purchon, 1977;Trueman, 1983). Nevertheless, the stygobiont gastropods (whose interstitial representatives hardly differ from the ones inhabiting caves) are known nearly only from taxonomic descriptions, mostly based on their shell morphology alone.
For example, Sbordoni et al. (2012) in their review of the structural studies of subterranean (meta) populations do not present any examples concerning the Gastropoda, as is also the case with Mammola et al. (2021). Simply there are no data about the Gastropoda from shallow hypogean habitats of many kinds, in similarity with hyporheic and deeper interstitial habitats along rivers and streams (Culver & Pipan, 2014). Apart from the descriptions of new species and faunistic lists, when considering the information on the stygobiont gastropod (meta)population genetic structure, robustly inferred levels of endemism, gene flow, etc. (perhaps because of the limited available material), the data are scarce. Thus, neither species boundaries nor phylogenetic relationships with the snails inhabiting surface environments are well understood (Falniowski, 2018). There are a few studies concerning the stygobiont gastropods, nearly exclusively of the Truncatelloidea. Falniowski et al. (1998Falniowski et al. ( , 1999 studied gene flow and metapopulation genetic structure (applying allozyme data) in Bythinella Moquin-Tandon, 1856 inhabiting springs. Bichain et al. (2007) also studied Bythinella in two caves (Padirac and Folatière) in France. Cytochrome oxidase (COI) sequences indicated several subterranean lineages, still belonging to the epigean species, thus, confirming that colonizations of hypogean habitats by epigean individuals are not rare events in the genus, as already suggested earlier (Boeters, 1979;Giusti & Pezzoli, 1982;Bole & Velkovrh, 1986;Bernasconi, 2000). Bichain et al. (2007) also found in sympatry or parapatry with those epigean lineages also distinct, exclusively hypogean phylogenetic lineages, thus, indicating multiple invasions of the caves and/or radiation already in caves. It has to be noted that Bythinella is, in general, not a stygobiont, but rather a stygophile, which means inhabiting both subterranean and epigean waters. In the western Balkans, Bythinella rather avoids the cave environment (JG personal observation).
COI sequences were used to estimate the time of divergence between Heleobia dobrogica (Grossu et Negrea, 1989) from the Movile Cave in Dobrogea, Romania and the epigean H. dalmatica (Radoman, 1974). The time 2.172 ± 0.171 Mya coincided with the period marking the beginning of the fall in temperature and precipitation that predated the Pleistocene, when H. dobrogica found a safe shelter within a warm cave (Falniowski et al., 2008). Falniowski & Sarbu (2015) described new species of Iglica A. J. Wagner, 1928 andDaphniola Radoman, 1973 from the Melissotripa Cave in Greece, inferring their molecular relationships (COI, 18S) with epigean species. Osikowski et al. (2015) studied the genetic and morphological differences between Bulgarian surface-and cave-dwelling populations of the genus Bythinella; while Rysiewska et al. (2016) Hershler & Holsinger (1990) described the zoogeography of the stygobiont North American Truncatelloidea, demonstrating that in some regions of the USA and Mexico, some stygobiont hydrobioids originated from their ancestors invading cave habitats directly from the sea during the late Cretaceous. Similar invasions from the sea have been also postulated for Dinarides' cave fauna (e.g. Sket, 2012), through the brackish/freshwater coastal sediments, especially during the Messinian Salinity Crisis (Boutin & Coineau, 2000). Contrary to common reports on the hydrobioid stygobiont species endemic to very restricted areaoften just one cave-the hydrobiid cave snail Fontigens tartarea Hubricht, 1963 is known from dozens of caves in West Virginia (North America), its geographic range reaching 200 km and showing patchy distribution, sporadic occurrence, and high variability of the shell (Hershler & Holsinger, 1990;Culver, 2012). In caves, one can also find some widely distributed generalist species, also inhabiting surface waters-such as the pulmonate gastropod Ancylus fluviatilis O. F. Müller, 1774 (reported by e.g. Sket, 2012; it has also been found in several interstitial pumped materials during our research), Radomaniola Szarowska, 2006, and the isopod crustacean Asellus aquaticus (Linnaeus, 1758) (reported by : Verovnik et al., 2003;Verovnik, 2012). Molecular studies on the stygobiotic genus Bythiospeum Bourguignat, 1882 revealed the relatively late post-glacial colonisation of large regions north of the Alps by one genetically uniform species, while two other species of this genus had a very limited distribution colonised from their post-glacial refugia (Richling et al., 2016). This pattern also suggested that the river alluvium and gravels are effective ways of dispersal for some species of Bythiospeum, but which are not so effective for some other species, not inhabiting gravel alluvia, which are rare in the karst Balkan regions.
In general, the central question for understanding the biogeography of the subterranean animals, for years known mostly from caves, is that of the relative role of surface and subterranean dispersal (Holsinger, 2005). In the case of obligatory subterranean aquatic species assuming isolation in caves is not obvious. In fact, there are many subterranean habitats which are not caves, but which also harbour eyeless, depigmented animals , 2014, especially in unconsolidated sediments bordering and underlying streams and rivers (hyporheic zone). They are parts of the interstitial habitat, neither rare nor discontinuous, thus, making possible migration between caves (Lamoreaux, 2004;Dole-Olivier, 2011). There are some studies confirming this possibility (e.g. Buhay & Crandall, 2005), but they are still few and are mostly devoted to the Crustacea (e.g. Lefébure et al., 2007;Eme et al., 2013), rather than the Gastropoda. The finds of live Paladilhiopsis oshanovae (Pintér, 1968) from the alluvial springs of the Danube (G. Majoros, pers. com.) clearly confirm alluvial gravel deposits as a habitat of this species; thus, the ways of its dispersal become obvious.
Pešić & Glöer (2013) Grego et Glöer, 2018, whose description was based on the slight differences in the shell morphology and geographic distribution, but which was not confirmed by the molecular differences (Grego et al., 2018). The aim of the present paper is to check the genetic structure (applying mtDNA cytochrome oxidase subunit I (COI), three nuclear markers and RAPD analysis) of the populations of Montenegrospeum inhabiting diverse subterranean habitats, including caves and interstitial ones, to estimate the levels of gene flow, to confirm the hypothesis that these snails can migrate through them and forms a genetically uniform metapopulation through a relatively long distance as a result of a high level of gene flow.

Material and methods
The snails were collected in 2018 and 2019, 27 specimens from 11 localities (Table 1), distributed in Croatia, Bosnia and Herzegovina, and Montenegro (Figs. 1 and 2). They were either collected by hand and sieve in caves and springs, or with a pump applying the Bou-Rouch technique (Bou & Rouch, 1967), to sample interstitial fauna below the bottom of streams, at a depth of about 50 cm. The tube was inserted in the streambed five times, and 20 L were pumped each time. Samples were sieved through a 500 lm sieve and fixed in 80% analytically pure ethanol, replaced twice, and later sorted. Next, the snails were put in fresh 80% analytically pure ethanol and kept at -20°C temperature in a refrigerator.
The shells were photographed with a Canon EOS 50D digital camera, under a Nikon SMZ18 microscope with dark field. Morphometric parameters of the shell (Table 2) were measured by one person using a Nikon DS-5 digital camera and ImageJ image analysis software (Rueden et al., 2017). The linear measurements were then logarithmically transformed; for angular measurements, the arcsine transformation was applied. Principal component analysis (PCA), based on the matrix of correlation, was computed, applying a descriptive, non-stochastic approach. The original observations were projected into PC space, to show relationships between the specimens, without any classification given a priori. For PCA analysis ClustVis 2.0 web tool (Metsalu & Vilo, 2015) was used.
DNA was extracted from whole specimens including their crushed shell; tissues were hydrated in TE buffer (3 9 10 min); then total genomic DNA was extracted with the SHERLOCK extraction kit (A&A  A fragment of mitochondrial cytochrome oxidase subunit I (COI) and three nuclear fragments (18S ribosomal RNA-18S, 28S ribosomal RNA-28S, histone H3-H3) were sequenced, as commonly used markers in phylogenetic studies of Truncatelloidea. Some other markers (nuclear internal transcription spacers ITS-1 and ITS-2, elongation factor-EF as well as mitochondrial 16S ribosomal RNA-16S) did not amplify despite numerous attempts. Details of PCR conditions, primers used, and sequencing techniques were given in Szarowska et al. (2016a). The Sanger sequencing was performed in Genomed Company in Warsaw, Poland. Sequences were initially aligned in the MUSCLE (Edgar, 2004) programme in MEGA 7 (Kumar et al., 2016). The correctness of the alignment was checked in BIOEDIT 7.2.5 (Hall, 1999), this programme was also used to translate sequences and check for reading frame and stop codons. Uncorrected p-distances were calculated in MEGA 7. The estimation of the proportion of invariant sites and the saturation test (Xia, 2000;Xia et al., 2003) were performed using DAMBE (Xia, 2018). In the phylogenetic analysis, additional sequences from GenBank were used (Table 3). The phylogenetic analysis was performed applying two approaches: Bayesian Inference (BI) and Maximum Likelihood (ML). In the BI analysis, the GTR ? I ? C model of nucleotide substitution was applied. Model was selected using MRMODELTEST 2.3 (Nylander, 2004). The Bayesian analyses were run using MrBayes v. 3.2.3 (Ronquist et al., 2012) with defaults of most priors. Two simultaneous analyses were performed, each with 10,000,000 generations, with one cold chain and three heated chains, starting from random trees and sampling the trees every 1000 generations. The first 25% of the trees were discarded as burn-in. The analyses were summarised as a 50% majority-rule tree. Convergence was checked in Tracer v. 1.7 (Rambaut et al., 2018). FigTree v.1.4.4 (Rambaut, 2010) was used to visualise the trees. The Maximum Likelihood analysis was conducted in RAxML v. 8.2.12 (Stamatakis, 2014) using the ''RAxML-HPC v.8'' on XSEDE (8.2.12) tool via the CIPRES Science Gateway (Miller et al., 2010). We applied the GTR model whose parameters were  (Stamatakis, 2014). Rapid bootstrap was calculated. Two species delimitation methods were performed: Poisson Tree Processes (PTP) (Zhang et al., 2013) and Automatic Barcode Gap Discovery (ABGD). The PTP approach was run using the web server https://species.h-its.org/ptp/, with 100,000 MCMC generations, 100 thinning, and 0.1 burn-in. We used RAxML output phylogenetic tree. The ABGD approach is using the web server (https://bioinfo.mnhn.fr/abi/public/abgd/abgdweb. html) and the default parameters. To infer haplotype network of the marker, we used a median-joining calculation implemented in NETWORK4 (Bandelt et al., 1999). The parameters describing a genetic variation and other statistics were calculated with DnaSP software (Librado & Rozas, 2009) and Arlequin v. 3.5 (Excoffier & Lischer, 2010). Arlequin was also used to perform Tajima's D test of neutrality (Tajima, 1989a(Tajima, , 1989b(Tajima, , 1993. The geographical distances between the localities were calculated with  (Ersts, 2020). Mantel tests were performed with NTSYSpc (Rohlf, 1998).
In the RAPD analysis, run to obtain another source of information about genetic diversity, five primers were used (Table 4), based on the data from Kulsantiwong et al. (2013). The RAPD reaction mix contained: 10 ll of 1 9 reaction buffer, 1 unit of Taq DNA polymerase (Thermo Fisher Scientific), 2.5 mM of potassium chloride, 0.2 lM of primer, 0.2 mM of dNTPs mix, and 10 ng of DNA0. The amplification conditions were as follows: the first step at 94°C for 5 min, followed by 41 cycles of 1 min each at 94°C, 1 min 30 s at 35°C, and then 2 min at 72°C. The final elongation step was for 5 min at 72°C. The amplification products were separated in 2% agarose gel with ethidium bromide in TBE buffer. The banding patterns were visualised under UV light. The bands were scored visually, and those of similar molecular size (for the same primer) were assumed to be homologous. Each band was treated as an independent locus with two alleles: the presence or absence of the band. The PopGen 32 software package (Yeh & Yang, 2000) was used to estimate the genetic variation. The following parameters were calculated: mean number of alleles per locus, effective number of alleles per locus, percentage of polymorphic loci and Nei's (1973) gene diversity, and the Shannon index. This same software was used to calculate the Nei's genetic distances between populations.

Results
Both techniques of species delimitation resulted in distinction of three mOTUs. The shells of the sequenced specimens ( Fig. 3A-Y) were slightly variable, only their size differed between the specimens. Neither the paratype of Montenegrospeum sketi (Fig. 3A), nor the specimens of molecularly distinct operational taxonomic unit (mOTU-see below) B and C differed morphologically from the other. Moreover, PCA analysis illustrates no variability and  Beran et al. (2015) no morphological distinctiveness of the three clades as well as for populations of the Montenegrospeum (Fig. 4).
The tests for COI and H3 revealed no saturation. We obtained 25 new sequences of each marker: COI (451 bp, GenBank Accession numbers OK001173-OK001197), 18S (401 bp, GenBank Accession numbers OK001198-OK001222), 28S (512 bp, GenBank Accession numbers OK001281-OK001305), and H3 (307 bp, GenBank Accession numbers OK033490-OK03351). For three nuclear fragments, all sequences were identical, so all analyses were made for the COI fragment. In the analysis of the COI, we included also all the sequences of Montenegrospeum available in the GenBank as well as for four sequences of the genera previously identified as the closest of Montenegrospeum for rooting the trees (Table 3).
We obtained 18 haplotypes for COI, with the haplotype diversity (Hd) 0.955 ± 0.024. In total, 44 sites were variable and nucleotide diversity per site (Pi) was 0.0276 ± 0.007.
The maximum likelihood (ML) phylogram (Fig. 5) presented the same topology as the consensus tree based on the Bayesian inference (BI). Reconstructing a shallow phylogeny of closely related taxa, it showed high values of bootstrap support or Bayesian probability, respectively. COI tree topology consisted of three molecular groups-mOTUs (confirmed by PTP and ABGD analysis), considering the identity of the sequences in the other loci they rather do not represent distinct species, although the values of pairwise p-distances for COI (Table 5) are within the interspecies level. All but four sequences of the snails, inhabiting ten localities (the type locality of M. sketi included: Fig. 1), formed mOTU A. Within this group, the genetic variation was very restricted, despite its wide (for stygobionts) geographic distribution, reaching 236 km (Table 5); three haplotypes were found in two (H9 and H15) or three (H11) populations (Table 1) mOTU B was found as a single specimen at one locality, and mOTU C as three specimens at two localities (Figs. 1 and 5). The above pattern of genetic differentiation was further supported by the results of AMOVA (Table 6). As much as 86.40% of the variance occurred between the mOTUs, confirming low genetic diversity within vs high between the mOTUs. Within the mOTUs, among the populations, only 7.36% of the variance was noted, and was not much more than 6.25% within the populations.
The results of Tajima's D test of neutrality were statistically insignificant for the mOTU A. For the other two clades (B and C) they were statistically significant, their values were 0.0000, rejecting selection (although in the mOTU B, with one sequence, the value was meaningless). In the clade A, the result was nearly significant (p = 0.8890), and D = 1.1291. Tajima's D [ 0 reflects, in general, scarce rare alleles (in our case substitutions), balancing selection and sudden population contraction. With one polymorphic locus only, this result should be interpreted with caution. The results of Fu's F test were statistically insignificant. The F ST pairwise values between the mOTUs were statistically not significant.
The Mantel test, checking the association between the matrices of pairwise p-distance and geographic distance, thus, testing the isolation-by-distance model ( Table 5,   populations 2 and 4, as well as 10 and 4, not statistically significant. Nm, an estimator derived from Fisher's F-statistics, may reflect the levels of gene flow, thus, migration in a subdivided population. In our populations, the values of Nm (Table 7) were, in general, high, with exceptions of the pairs of populations belonging to different mOTUs, or the pairs including the snails from the localities 10 and 5, Fig. 3 Shells of the sequenced specimens of Montenegrospeum: A-U-molecularly distinct operational unit (mOTU) A: A 2H70, locality 1, B 2D09, locality 2, C 2D10, locality 2, D 2D13, locality 2, E 1T57, locality 3, F 2D02, locality 4, G 2D04, locality 4, H 2J40, locality 5, I 2J41, locality 5, J 2J42, locality 5, K 2J37, locality 6, L 2J38, locality 6, M 2J39, locality 6, N 2C43, locality 8, O 2C44, locality 8, P 2C45, locality 8, Q 2C15, locality 9, R 1T19, locality 10, S 2C18, locality 10, T 2C40, locality 10, U 2C19, locality 10; V mOTU B, 2C28, locality 5; W-Y mOTU C: W 2C41, locality 7, X-Y locality 10 (2G20 and 2G23, respectively) inhabited by the representatives of two mOTUs, in sympatry. The network of 18 found haplotypes (Fig. 7) confirms the distinctiveness of the mOTU C, but not necessarily that of the mOTU B. Within the most diversified mOTU A, there are 15 haplotypes, and the pattern more or less resembles the so-called ''star phylogeny''.
In the RAPD analysis, each primer amplified from six to twelve fragments, with an average of eight, and among them, 38 bands (90%) were polymorphic (Table 4). Basic parameters of diversity are given in  Table 8. Nei's (1973) gene diversity ranged from 0.103 to 0.243 (mean 0.138, and the Shannon index ranged from 0.035 to 0.366, (mean 0.159)). The values were rather low, with the exception of the snails from the localities 7 and 10, inhabited each by two mOTUs. The Mantel test of association between the matrices of geographic distance and Nei's genetic distances (Table 9, Fig. 8) with 9999 permutations, resulted in p [random Z B observed Z] = 0.0189. There were not found to be any population-specific RAPD fragments, defined as RAPD bands found exclusively in one population in over 90% of individuals. Most of the detected RAPD bands were distributed across most populations (Fig. 9). In general, UPGMA clustering ( Fig. 10) showed groupings of specimens in accordance with the localities from where they came, and within six main clusters. The first of them groups all individuals from the localities 1, 3, and 4. The second groups individuals from the locality 7, belonging to mOTU A. The third groups all the individuals from the localities 5 and 6, belonging to mOTU A. The fourth consists of specimens from the locality 2, while the fifth consists of all the individuals of mOTU C (however, high genetic distance within mOTU C may suggest the existence of two divergent groups). The sixth clusters the specimens from the localities 8, 9, 10 and a single individual of mOTU B from the locality 5 (Fig. 10).
In the karst areas of the Balkans the gravel substratum is not common, with interstitial waters usually filling small fractures in the rocks, making the collection of the meiofauna problematic, especially when using the Bou-Rouch pump. Coupled with low densities of such subterranean snails' populations, most of the pumping does not result in the collection of live snails. We collected live specimens of Montenegrospeum at two localities: 2 and 10. Locality 2 was geographically most close (geographic distance 34 km) to the locality 4 (Table 4), which was a cave, and genetically, the populations were very close (p = 0.007; Table 4), confirming the possible migration from one locality to another. The locality 10 Table 5 The p-distances between Montenegrospeum populations (below the diagonal), p-distance within population (diagonal, bold, italic) and geographic distances (in km) between population (above diagonal) Statistical significance *P \ 0.001 harboured two mOTUs: A and C. The mOTU C was found also in a spring at locality 7, which is 66 km far from the locality 10, and p = 0.066 (Table 4), and Nm equalled 0.000. For the representatives of the mOTU A, the localities closest to 10 were 9 (24 km), 11 (25 km), and 8 (35 km; Table 4), all being springs. The representatives of mOTU A from the locality 10 shared haplotypes with either locality 8 or 9.

Discussion
Neither morphology nor molecular data (Grego et al., 2018) confirmed the species distinctiveness of Montenegrospeum sketi. Thus, M. sketi becomes a younger synonym of B. bogici. Such unjustified species descriptions are not rare in truncatelloid gastropods (e. g. Falniowski, 2018;Osikowski et al., 2018). Again, for the reasons outlined by e.g. Falniowski (1987Falniowski ( , 2018, the morphostatic mode of evolution, as defined by Davis (1992), and common cryptic species, morphologically based species identification in the Truncatelloidea is often impossible. On the other hand, in our Montenegrospeum molecular data, we uncovered three molecularly distinct units; their distinctiveness is within the range of COI interspecies differentiation in the Truncatelloidea (e.g. Bichain et al., 2007;Falniowski & Szarowska, 2011), although not confirmed by the other loci, thus, not satisfying the criteria of species delimitation (Fišer et al., 2018).   Perhaps mOTU C represents a distinct species, which is not reflected in its morphology-thus, there may be one more example of a cryptic species, not as rare as it had been thought before the common application of molecular techniques in taxonomy (e. g. Bickford et al., 2006;Pfenninger & Schwenk, 2007;Macher et al., 2016;Razkin et al., 2016Razkin et al., , 2017. Cryptic species seem especially common in subterranean habitats (e.g. Culver, 2012). In our case, the species distinctiveness of mOTU C is somewhat confirmed by the sympatric occurrence with mOTU A. Thus, instead of one stygobiont species, rather widely geographically dispersed, we may observe one widely distributed species (235 km of linear distance is neither characteristic nor common for subterranean species: e.g. Culver & Pipan 2009), and two possible cryptic species whose ranges are restricted to one or two localities known so far. Such a pattern has already been found for some stygobiont species (e.g. Eme et al., 2013;Gorički & Trontelj, 2006). Within the mOTU A, which includes the topotypes of M. bogici and, thus, represents M. bogici s. stricto, the interpopulation differentiation is low, lower than in, for example, the crenobiont/stygophile (inhabiting springs and facultatively subterranean habitats) Bythinella Moquin-Tandon, 1856 (Falniowski et al., 1998(Falniowski et al., , 2009Falniowski & Szarowska, 2011;Osikowski et al., 2015;Szarowska et al., 2016b). Low values of p-distance and low levels of intrapopulation polymorphism may suggest rather high levels of gene flow, and either sudden population contraction (suggested by Tajima's test, although the results have been marginally significant only), or a ''star phylogeny'' pattern visible in the haplotype network, which could be interpreted as an expected signature for a species that has expanded its range rather recently, from a small or modest number of founders (Avise, 2000). As was stressed by Endler (1977), adaptive differences among populations can be maintained by natural selection even under high levels of gene flow, as well as selection may mimic effects of gene flow. This similarity may be also due to numerous stochastic factors, area effects, etc. Thus, another possible background of low variation may be selection, which is rather obvious in COI, a part of the respiratory chain. However, such selection should act uniformly in such a rather homogenous (at least concerning respiration: in a general low oxygenated) habitat as a subterranean one, but the mOTUs B and C, living in the same habitats (even in sympatry with A) are molecularly markedly distinct, which makes selection a much less probable factor shaping low molecular diversity. If Nm [ 1, the allele frequencies in the subpopulations remain homogenised (Wright, 1931(Wright, , 1969. If Nm \ 1 but is still positive, an equilibrium based on the rate of mutation, migration, and genetic drift will be established. In our populations, the values of Nm were, in general, high, with exceptions of the pairs of populations belonging to different mOTUs, or the pairs including the snails from the localities 10 and 5, inhabited by the representatives of two mOTUs, in sympatry. The infinite-island model of interpopulation differentiation (Wright, 1978), expected for isolated populations, does not characterise the observed pattern. Together with the results of the Mantel test, it suggests isolation-by-distance (although slightly marked), rather than the stepping-stone model. RAPD confirms the same conclusion. Another three available, nuclear loci were not varied, which only confirms close relationships/high levels of gene flow of the studied individuals. The low subsurface dispersal hypothesis (Lefébure et al., 2006), most popular for many years, states that after a vicariant event, invasion of subterranean habitat (group of caves) followed by extinction of the surface ancestor, there is little or no dispersal between the subterranean microhabitats. Inhabitants of subterranean habitats were often traditionally understood as being either living fossils, relic species (last survivors of ancient radiation), or relict species (geographically separated from related species) and showing no progressive evolution, being isolated an extremely space-limited habitat. This already abandoned concept seems to be the background of the popular practice among some taxonomists, many of them rather speleologists than biologists-to describe distinct species or at least subspecies for nearly each cave (two species described from two caves at the same village may serve as an example), with no understanding of either mechanisms of speciation or of potentially available dispersal pathways. Molecular markers, especially DNA partial sequences, have identified numerous cryptic subterranean species (e.g. Trontejl et al., 2009) but have also rejected the species distinctiveness of many nominal species described following the low subsurface dispersal hypothesis (e.g. Osikowski et al., 2017Osikowski et al., , 2018. However, Christman et al. (2005) in their study of the troglobiotic (terrestrial) fauna in eastern North America reported as many as 45% (211 of 467 species) to be single cave endemics.
Such high levels of endemism seem not so characteristic for the aquatic, stygobiont fauna. Ward & Palmer (1994) in their overview of dispersal of meiofauna, applicable to stygobionts in general, stressed the connections of subterranean habitats, such as caves and springs, by alluvial aquifers, with their sand and gravel deposited by flowing waters. Williams (2008) stressed that in almost all caves, the rock around is fractured, forming small solution tubes that allow subsurface connections between caves; in karst areas, epikarst and associated vertical downward percolation of water is more or less continuous. Their opinion was confirmed by phylogenetic analysis of four species of Proasellus Dudich, 1925 (belonging to the Crustacea: Isopoda) by Eme et al. (2013). There are some more studies confirming this possibility, which are mostly devoted to the Crustacea (e.g. Lefébure et al., 2007;Villacorta et al., 2008;Malard et al., 2009;Eme et al., 2013), but also to anchialine  (Gonzalez et al., 2017) which are the prevailing group in the stygobiont fauna (Sket, 2012). To our knowledge, there are no such studies on any stygobiont gastropod species. Our study, although based on a few specimens-as already stressed stygobiont gastropods do not form dense populations and are not easily collectable-presents the same picture: dispersal through the subterranean interstitial waters which they also inhabit, leading to low levels of isolation and the resultant endemism.
There may be various hydrological routes of migration of the stygobiont fauna. For example, in the Dinarids, there are large flat depressions (karst poljes), some with neither aerial inflow nor aerial outflow of water (closed karst poljes) and having insoluble polje-floor sediments. They are periodically flooded. During such a flood, aquatic cave animals are frequently washed out of their subterranean habitat, and some of them may reach caves at other parts of the polje, or even reach a subterranean connection to another polje. Such a route of migration was suggested by Zakšek et al. (2009) for the shrimp Troglocaris Dormitzer, 1853. However, Troglocaris is a swimmer, and for creeping snails, such a means of dispersal seems less probable, although some may possibly passively wash out and disperse.
The relatively wide geographical distribution of Montenegrospeum bogici, over 230 km, with one or two cryptic species at one/two populations, presents a pattern already found for some other, non-gastropod more motile stygobionts. The phylogenetic analysis of five morphospecies of the crustacean Proasellus, widely distributed (more than 200 km) along the hyporheic corridors of rivers, resulted in the discovery of ten new cryptic species, seven known from a single locality and being peripheral isolates (Eme et al., 2013)-such peripheral isolates may be mOTUs B and C. The most widely known and familiar stygobiont animal Proteus anguinus Laurenti, 1768, the European Cave Salamander, has a wide geographical range along the western part of former Yugoslavia, although the range is discontinuous (Sket, 1997). Gorički & Trontelj (2006) sequenced several regions on its mtDNA and distinguished six groups of its populations, with genetic divergence between them higher than between most species of salamanders within a genus. The greatest geographic extent of a single group does not extend 200 km, compared with 500 km for all the nominal species. All the distinguished six species were cryptic. The dispersal pattern of Montenegrospeum also strongly resembles the one described for Bythiospeum presented by Richling et al. (2016).
The relatively broad geographic range and high levels of gene flow in M. bogici observed over areas without alluvial connection need some possible explanation. The hydrogeological map of the species range (Fig. 11) suggests that M. bogici may inhabit not only the alluvial gravels and the spring heads, but should also be capable of spreading through the cave streams and lakes. Additionally, the species may use the intermittent connections between the aquifers of river drainage basins. Such an intermittent connection between the Drin (Zeta) and Trebišnjica River Basins is relatively well understood through Nikšičko Polje (Gornjepolski Vir) to Nikšičko Vrelo near Bileća (divergence an in Fig. 11). The species also passes the second intermittent divergence line between the Trebišnjica and Neretva River Basins, where it exists between the Dabarsko and Fatničko Poljes. Although this hydrological connection has been previously proven to exist when the karst groundwater exceeds 400 m a.s.l. elevation, the identical stygobiont malacofauna of both Poljes, with Plagigeyeria reischuetzorum Grego, 2020 and Travunijana gloeri Grego, 2020, now adds further evidence to confirm the existence of this connection. Furthermore, the stygofauna of both basins may be in contact through the springs in the delta of the Lower Neretva River in Hutovo Blato/Deranjsko Jezero and also near Metković (springs Bajovci, Sjekoše, Glušći, Bijeli Vir, Mlinišče). The third, hydrogeologically least understood, unremittent connection possibly exists between the Neretva and Cetina Basins. Most likely it is formed by occasional groundwater flow around Imotsko Polje, from Prološko Blato to the spring of Ričina in Buško Blato, or may be directly to the spring Ruda Beguša. The high levels of gene flow presented above indirectly prove the existence of this groundwater connection, at least intermittently. The finds of Montenegrospeum shells at the nearby spring Peć Mlini (Tihaljina, BiH) and in the spring Lukavac (Veliki Prolog, HR) indicate the wider distribution of the genus in areas between the hitherto-proven localities.
Concluding our results, in general, neither the infinite-island model of interpopulation differentiation, which would be expected for isolated populations, nor the stepping-stone one, but rather the isolation-by-distance model explained the observed Estavelles; dark green diamond: Submarine or brackish karst spring; Grey fields: permeable (conductive) karst; White field: not permeable (non-conductive) non-karstic bedrock; Yellow fields: Permeable (conductive) alluvial gravel; Orange fields: Alluvial sands with limited permeability (conductivity) pattern. Our results suggest that interstitial habitats provide ways of migration for the stygobiont M. bogici, as has been already suggested for other subterranean gastropods.