Recurrent founder-event speciation across the Mediterranean likely shaped the species diversity and geographic distribution of the freshwater snail genus Mercuria Boeters, 1971 (Caenogastropoda: Hydrobiidae)

Dispersal is known to play an important role in shaping the diversity and geographic range of freshwater gas-tropods. Here, we used phylogenetic methods to test for the influence of dispersal and other biogeographic processes (such as vicariance) on the speciation and distribution patterns of Mercuria Boeters, 1971, a snail genus widely distributed in the western Palaearctic. The 25 extant species traditionally thought to comprise the genus, which were described mainly on the basis of morphology, have been recorded from lowland waters in both the Mediterranean and Atlantic river basins of Europe and North Africa. Using molecular phylogenies based on three gene fragments (COI, 16S rRNA and 28S rRNA) from 209 individuals, four molecular species delimitation methods and a shell characterization, we identified 14 putative species in our dataset, nine of which correspond to species classified by traditional taxonomy. Furthermore, biogeographical modelling favoured a scenario in which recurrent founder-event speciation since the late Miocene is the most probable process explaining the species diversity and distribution of the Mediterranean clades, whereas episodes of postglacial northward colonization from Iberian refugia by the species M. tachoensis may explain the current presence of the genus in Atlantic lowlands. The dispersal events inferred for Mercuria , probably promoted by multiple factors such as the changing connectivity of drainage basins driven by climate change or better access for avian dispersal vectors in lowlands, may explain the rare case among hydrobiids of a species-rich genus containing individual species with a large distribution area .


Introduction
Knowledge of past distribution patterns and the processes that changed them over time is important for understanding the geographic range and species richness of present-day groups (Wiens and Donoghue, 2004).Several model-based methods developed in recent years have been used to infer the biogeographic patterns and processes (i.e., vicariance, dispersal or local extinction) of freshwater gastropod faunas throughout the western Palaearctic.Indeed, freshwater gastropods are suitable organisms for biogeographic studies given the dispersal limitation and habitat specialization of most of the species (Strong et al., 2008).Such studies have shown that changes in the climate and drainage systems of Eurasia during the Neogene period may have shaped the distribution of past freshwater gastropod species (Aksenova et al., 2018;Neubauer et al., 2015).The Pleistocene biogeography of European lacustrine and spring snails suggests a series of dispersal and vicariant events that were likely promoted by the glacial-interglacial cyclicity experienced by the region during this period (Anistratenko and Anistratenko, 2020;Benke et al., 2009;Georgopoulou et al., 2016).Although these findings revealed the processes that may have affected the past and present geographic distribution of Eurasian freshwater gastropods, we have very limited knowledge about the diversity and biogeography of the genera with a circum-Mediterranean distribution occurring in both Europe and North Africa.Due to the intense tectonic activity in the Mediterranean affecting the connectivity of drainage basins during the Miocene and Pliocene (Cohen, 1980;Rosenbaum et al., 2002), we expect that these events had a substantial impact on the ancestral geographic ranges of aquatic taxa.Though, emerging evidence on the evolution of Mediterranean freshwater gastropods suggests that observed biogeographic patterns are due to long-distance dispersal across drainage basins or the Mediterranean Sea rather than range fragmentation or vicariance (Boulaassafer et al., 2020;Neubauer et al., 2016;Sands et al., 2019).However, more research is needed to evaluate the processes responsible for such an intercontinental distribution.
The biogeographic origin of most genera within the highly diverse family Hydrobiidae Stimpson, 1865 has been poorly investigated as the great majority consist of a small number of springsnail species with narrow distributions (Miller et al., 2018).However, a few genera represent ideal candidate groups in which to study the biogeography of freshwater gastropods at a larger scale across the Mediterranean.These few genera are sufficiently diverse and widely distributed to be able to infer the potential imprints of Mediterranean geo-climatic changes on their biogeographic histories and speciation.For example, evolutionary studies of the widespread hydrobiid genera Pseudamnicola Paulucci, 1878 and Ecrobia Stimpson, 1865 have provided evidence for the relatively recent, long-distance dispersal of lowland populations across the Mediterranean or adjacent regions (Haase et al., 2010;Delicado et al., 2014;Boulaassafer et al., 2020).In another study, old vicariant events have been proposed to have split the Iberian and Moroccan clades of the genus Corrosella Boeters, 1970(Boulaassafer et al., 2021).However, for many hydrobiid groups, how species reached their current distribution is still unknown.
Mercuria Boeters, 1971(subfamily Mercuriinae Boeters & Falkner, 2017) is among one of the least explored hydrobiid genera, yet its relatively high species richness, wide distribution across the Mediterranean and Atlantic and potential parallel evolutionary history with the co-occurring hydrobiid genus Pseudamnicola make it of special interest for biogeographic studies.The current taxonomy of Mercuria comprises a relatively large number of species (25 extant and 5 extinct, according to MolluscaBase, 2021), with some species showing high levels of intraspecific variability of morphological characters, making the species delimitation based on morphology problematic.Recent studies have attributed the high level of intraspecific variation observed in shell and penis characters of some species to parasitism or ontogeny (Boulaassafer et al., 2018;Holyoak et al., 2017).In a few cases, some of the morphological species are hardly distinguishable from others and could be categorized as cryptic species (e.g., M. tachoensis and M. bayonnensis).All of these aspects, along with missing type material and/or unspecific descriptions of type localities of the original specimens (Boeters and Falkner, 2017), suggest the need for a re-evaluation of the currently recognized species of Mercuria using molecular phylogenetic methods.This would provide a more robust taxonomic classification of Mercuria species, which is required to determine the precise distribution range of the species, with the final aim of inferring more accurate biogeographic patterns and histories.
Previous taxonomic works have reported the presence of Mercuria within the Mediterranean in the Iberian and Italian peninsulas, southern France, some western Mediterranean islands and north-western Africa (Boeters and Rolán, 1988;Boeters and Falkner, 2017;Giusti et al., 1995;Glöer et al., 2010).In the Atlantic lowland regions, a few species have been reported from the southern British Islands, northern continental Europe, the Iberian Peninsula, Morocco and Macaronesia (Boeters and Falkner, 2017;Glöer et al., 2015;Kerney, 1992Kerney, , 1999;;Marquet and Lenaerts, 2008).A large proportion of the Mediterranean range of Mercuria also comprises part of the distribution of Pseudamnicola species (subfamily Pseudamnicolinae Radoman, 1977), with species of both genera living together, even within the same locality, in the Iberian Peninsula, North Africa, southern France, Malta and the Balearic Islands (Boeters, 1988;Boulaassafer et al., 2020;Delicado et al., 2014;Glöer et al., 2010).Competitive behaviour between co-inhabiting populations of the two genera has not been documented; however, its possibility cannot be excluded.In any case, the observed high degree of cooccurrence suggests that these species exhibit similar ecological preferences and could potentially present similar biogeographic histories.
For Pseudamnicola, a series of long-distance dispersal events by birds across the Mediterranean, starting from the late Miocene, has been proposed as the primary process explaining current species distributions (Boulaassafer et al., 2020;Delicado et al., 2015;Szarowska et al., 2016).This may represent a feasible biogeographic scenario for Mercuria given the common habitat preferences of both genera.An alternative hypothesis for the current distribution patterns of Mercuria is that these are a result of old vicariant events associated with the connection/fragmentation of landmasses in the Mediterranean during the Oligocene and Miocene (Rosenbaum et al., 2002).These geological episodes have been proposed to promote speciation in other hydrobiid genera such as Corrosella (Boulaassafer et al., 2021) and Salenthydrobia Wilke, 2003(Wilke, 2003).Sympatric speciation, although less plausible, may be another biogeographic process that could have played a role in the distribution of Mercuria.Sympatric speciation has been very rarely suggested in the hydrobiid literature (Falniowski, 2018), probably because there is insufficient knowledge about the ecological preferences and life history traits of hydrobiid populations.In any case, to our knowledge, there are no reported cases of two Mercuria species living together within the same locality.
The pipeline of this study has been designed to infer the biogeographic processes (i.e., dispersal, vicariance or sympatric speciation) that could potentially explain two unusual characteristics of Mercuria compared with other Hydrobiidae taxa: (1) its relatively high species richness (mainly in the Mediterranean region) and (2) its wide geographic distribution along Mediterranean and Atlantic coastal territories.For this purpose, we first used mitochondrial and nuclear DNA data sets to infer the phylogenetic position of the study populations, which include topotypes or 'near topotypes', and other records of ten nominal taxa of Mercuria.We then used various molecular species delimitation methods to classify specimens into species groups.Some groups could be straightforwardly identified with currently recognized species as they include individuals from other studies that were previously delimited on the basis of morphology.Finally, using the species scheme suggested in the previous step, we inferred a coalescent-based species tree that was then used for subsequent inference of the biogeography of the studied Mercuria species.The contrasting patterns of species diversity, distribution and biogeographic history found between the Mediterranean and Atlantic groups is information that can be taken into consideration in decision-making for conservation purposes.

Material examined
A total of 79 populations of Mercuria were genetically analysed in this study (Fig. 1), some of which had been previously assigned to ten nominal species.These populations cover most of the geographic range of the genus (see Introduction).Species of Mercuria usually live in habitats in the lower courses of rivers (Boeters and Falkner, 2017).Snails from the majority of the localities in Europe were newly collected for this study by handpicking from stones or from sediment recovered in sieves after washing the stones.For the genetic analyses, specimens were either frozen at − 80 • C or preserved in absolute ethanol and then conserved at − 20 • C. Vouchers and DNA samples were deposited in the Tissues and DNA Collection at the Museo Nacional de Ciencias Naturales (MNCN), Madrid, Spain.The samples from North Africa, Madeira, Italy and Malta were previously collected and are part of the University of Giessen Systematics and Biodiversity (UGSB) Collection (Diehl et al., 2018) in Germany.Some of the analysed specimens are paratypes, topotypes or 'near-topotypes' of nine known Mercuria species, which helped us to assign clades to valid taxa.However, the taxonomic coverage of our study is uncertain as the taxonomy of most of the valid species needs to be reassessed using integrative approaches.Information on studied specimens is provided in the supplementary material (Supplementary Material, Table S1).
Phylogenetic relationships of Mercuria were assessed for individual gene fragments and the concatenated (COI, 16S and 28S) dataset under Bayesian inference (BI) and maximum likelihood (ML).The BI analysis was conducted using Markov chain Monte Carlo (MCMC) sampling in MrBayes v.3.2.7a (Ronquist et al., 2012) under the best-fit substitution model for each gene partition as selected by jModelTest v.2.1.6(Darriba et al., 2012) according to the corrected Akaike information criterion (AICc) (Akaike, 1974).The selected models were: HKY (Hasegawa et al., 1985) + Γ for COI, K80 (Kimura, 1980) + Γ for 16S and TrN (Tamura and Nei, 1993) + I + Γ for 28S.MCMC sampling was run for 20 × 10 6 generations, four parallel chains and saving a tree every 1,000 generations.After assessing chain convergence by ensuring a standard deviation of the split frequency < 0.01, the first 10% of the sampled trees were discarded as burn-in.The robustness of the inferred tree was quantified using Bayesian posterior probabilities (BPP).The ML analyses were computed with RAxML-NG v.1.0.2 (Kozlov et al., 2019), applying the best-fit substitution model for each partition and 10 random starting trees.Branch supports were assessed by heuristic bootstrapping (BS) with a stopping threshold of 0.03 and later quantified using the transfer bootstrap expectation (TBE; Lemoine et al., 2018) metric.Final trees and branch supports were visualized in FigTree v.1.4.3 (Rambaut, 2012).

Molecular species delimitation methods
To re-evaluate the morphology-based taxonomy of the Mercuria species included in our dataset, we used four molecular species delimitation methods that are based on genetic distance or Bayesian and maximum likelihood approaches.These methods have previously proven useful to delimit molluscan species, including hydrobiids, and the molecular-based results are often congruent with those based on morphology (e.g., Araujo et al., 2017;Delicado et al., 2019;Padula et al., 2016;Strong and Whelan, 2019).
The distance-based automatic gap discovery (ABGD) method (Puillandre et al., 2012), which is based on the detection of a barcode gap between inter-and intraspecific variability, was used as a threshold to delimitate species.The COI alignment was uploaded to the ABGD server (http://wwwabi.snv.jussieu.fr/public/abgd/abgdweb.html), and a Kimura 2-parameter (Kimura, 1980) distance matrix was generated with the parameters of an intraspecific divergence between 0.001 and 0.1 and a relative gap width of 1.0.
We also applied the Poisson Tree Processes (PTP) method (Zhang et al., 2013) in two variants: the Bayesian single-rate (bPTP) and the multi-rate PTP (mPTP; Kapli et al., 2017).In general, PTP delimits species based on the substitutions between and within species independent of time, but mPTP incorporates different rates of coalescence within clades and determines the number of species that best fits the data using the AIC.As the PTP method requires a non-ultrametric and bifurcated tree as the input file, we used the best topology based on the concatenated mitochondrial (COI and 16S) dataset estimated with RAxML after pruning identical haplotypes with the R v.4.1.2(R Development Core Team, 2021) package ape v.5.5 (Paradis et al., 2004).The bPTP and mPTP analyses were conducted in mPTP v.0.2.4 (downloaded from GitHub at https://github.com/Pas-Kapli/mptp;Kapli et al., 2017).Priors for bPTP were set to 50 × 10 6 MCMC generations (first million discarded as burn-in), saving a tree every 100,000 iterations, and a minimum branch length of 0.001.The same value of minimum branch length was used for the mPTP analysis.
The final delimitation method used was the single-threshold generalized mixed Yule-coalescent method (ST-GMYC) (Pons et al., 2006), which defines either coalescent or speciation processes of ultrametric trees (i.e., calibrated based on the molecular clock approach) by detecting a time threshold on their branching patterns.To do this, an ultrametric tree was first inferred using BEAST v.1.8.4 (Drummond et al., 2012) based on the concatenated mitochondrial dataset.Clock models and trees were linked, and jModelTest was used to select the best model for each partition.The analysis was conducted in relative divergence time (i.e., all branches evolving with a rate of one substitution per site per unit of time), and the prior for the tree was set to the Yule speciation process (Yule, 1925).The MCMC sampling was run for 100 × 10 6 generations, saving a tree every 2,000 iterations and discarding the first 10% as burn-in.Stationarity in the posterior distribution was ensured by having an effective sample size (ESS) above 200, as calculated in Tracer v.1.7.1 (Rambaut et al., 2018).The maximum clade credibility (MCC) tree was identified from the sampled trees using TreeAnnotator v.1.8.4.Subsequently, all zero branch lengths were removed using the R package ape.The ultrametric MCC tree with unique haplotypes pruned to single tips was used as the input, and the analysis was performed on the GMYC web server (https://species.h-its.org/gmyc/).
To assess the performance of the species delimitation methods, we calculated the match ratio as in Ahrens et al. (2016) according to the following formula: Match ratio = 2 × Nmatch/(Ndelimited + Nmorph), whereby Nmatch refers to the number of delimited species that exactly match a defined morphospecies, Ndelimited is the total number of species delimited by the method and Nmorph refers to the total number of morphospecies.Lower values of match ratio indicate an overestimation of the species delimitation method.

Time calibrated tree and biogeographic modelling
We inferred a time calibrated species tree of the combined (COI, 16S and 28S) dataset without outgroups under a Bayesian multispecies coalescent framework using the extension *BEAST (Heled and Drummond, 2009) of the package BEAST v.1.8.4.The 16S and 28S amplified fragments from the specimen of M. rolani were excluded from this analysis because the COI gene fragment could not be amplified and this missing information distorted the species coalescence-based inference by *BEAST.This extension of BEAST combines datasets from several gene fragments and multiple individuals per species, clustered according to a grouping file, to generate a species tree.The input grouping file was configured according to the species groups obtained by the species delimitation methods described above.Despite recent discoveries of Mercuria fossils (Esu and Girotti, 2015a, 2015b, 2019;Sesé et al., 2004), the evolutionary convergence of shell shape among hydrobiid genera (e.g., Bodon et al., 2001;Delicado et al., 2016) makes it difficult to accurately assign a fossil to a specific genus and species.Therefore, in order to calibrate the molecular phylogenies, we used an external COI substitution rate of 0.81 ± 0.24% Myr − 1 (the percentage of substitutions per lineage per million years) calculated in a previous study of the hydrobiid genus Corrosella (Delicado et al., 2013) and applied only for the COI partition.This rate represents an average of substitution rates calculated for other hydrobiid genera using different biogeographic events [e.g., 0.91% by Wilke (2003); 0.81% by Hershler and Liu (2008); 0.72% (mean rate) by Wilke et al. (2009)].The analyses were run for 120 × 10 6 iterations sampling every 2,000 generations.The substitution rates of J.P. Miller et al. the 16S and 28S partitions were estimated relative to the COI rate in the initial runs using a normal distribution through MCMC sampling.
The results of the analyses were evaluated for convergence (ESS > 200) using Tracer v.1.7.1.The MCC tree was selected using TreeAnnotator with a 10% burn-in and visualized in FigTree.Node ages and 95% highest posterior density (HPD) intervals were also displayed on the MCC tree.Bayes factor (BF) was used to assess the applicability of strict vs. relaxed clock models (i.e., substitution rate constancy among branches) by comparing their marginal likelihoods (mL).Marginal likelihoods were calculated in BEAST v.1.8.4.by stepping-stone sampling using 75 paths and 2,666,666 iterations.To place the tree into a geological context, the tree was visualized using the R packages strap (Bell and Lloyd, 2015) and phytools (Revell, 2012).
Ancestral areas were inferred with the R package BioGeoBEARS v.0.2 (Matzke, 2013a).Each putative species was first assigned to one or more of the eight freshwater ecoregions (Abell et al., 2008) encompassing the present distribution of Mercuria.The applied algorithm then estimates the ancestral area of distribution and tests possible models of evolutionary processes such as dispersal, vicariance and founder-event speciation.
The following models were tested: Dispersal-Extinction-Cladogenesis (DEC; Ree et al., 2005), Dispersal-Vicariance (DIVA; Ronquist, 1997) and BayArea (Landis et al., 2013).The likelihood versions of the last two models as implemented in BioGeoBEARS are known as DIVALIKE and BayAreaLIKE, respectively.These models have two free parameters estimated along the phylogenetic branches: d = dispersion or range extension and e = extinction or range contraction.Additionally, model improvement was tested using the founder effect parameter j (occurring at cladogenetic events) for each of the tested models, resulting in three additional models: DEC + J, DIVALIKE + J and BayAreaLIKE + J.The AICc was used to select the best model.The most likely ancestral range was calculated for each node in agreement with the best-fit model and the results were plotted on the MCC tree.

Genetic variation and sequence divergence
The three analysed gene fragments (COI, 16S and 28S) were amplified from a total of 183 specimens of Mercuria; for the specimen of M. rolani, only two of the gene fragments (16S and 28S) were amplified.Another 26 COI sequences of Mercuria species available from GenBank were included in the analyses, as well as the corresponding 16S and 28S sequences that were derived from the same genomic DNA as the COI sequences (i.e., from the same study).The concatenated alignment consisted of 209 sequences and was 2,088 base pairs (bp) in length.From these sequences, we observed 80 identical haplotypes in our search for unique haplotypes to include in the bPTP analysis.Average base frequencies for the complete dataset were 25.2% A, 20.3% C, 25.1% G and 29.4 % T. Of the sites, 1,579 (75.51%) were conserved, 512 (24.48%) were variable and 336 (16.07%) were parsimony-informative.

Mitochondrial dataset
An alignment of 1,165 bp was generated for the mitochondrial genes COI (658 bp) and 16S (507 bp).For each of the gene partitions, the BI method recovered 13 species-level clades within Mercuria (see Supplementary Material, Figures S1-S2) with high supports (BPP > 0.95).Eight of these clades individually represent seven formally recognized morphospecies and the species M. melitensis and M. rolani, which were collapsed into a single species-level clade.For the combined mitochondrial dataset, BI analyses recovered four well-supported (BPP > 0.95) groups of clades: (i) M. similis + Mercuria sp. 1 + Mercuria sp.2; (ii) M. balearica + M. tensiftensis + Mercuria sp.4; (iii) M. midarensis + M. targouasensis and (iv) Mercuria sp. 3 + Mercuria sp. 5.The remaining species grouped with one of the four aforementioned clade groups with low support (BPP < 0.95).The ML analyses recovered clade groups (i) and (ii) with high support (BS > 75); however, the rest of the species grouped with low support (BS < 75).

Nuclear dataset
The final alignment of the partial nuclear 28S sequences consisted of 923 bp.The BI analyses (see Supplementary Material, Fig. S3) yielded four supported species-level clades: Mercuria sp. 3, M. saharica, M. targouasensis and M. melitensis (BPP > 0.95).These clades were also recovered in the ML analyses with high support (BS > 75), except for M. saharica.The rest of the species grouped with low support (BS < 75).

Combined dataset
The combined mitochondrial and nuclear dataset (COI, 16S and 28S) consisted of an alignment with 2,088 bp.Despite the low support of most branches in the 28S phylogeny, phylogenetic analysis of the concatenated dataset using both BI and ML approaches recovered a total of 13 species-level clades with high support (BS > 75, BPP > 0.95) (Fig. 3).Of these clades, seven correspond to previously described morphospecies, one clusters M. melitensis and M. rolani and five are putative species not yet described. Both

Molecular species delimitation methods
The species delimitation methods applied to our dataset (209 sequences) differed in the number of putative species inferred (Fig. 3), ranging between 8 (ST-GMYC) and 21 (bPTP) species groups.Overall, ABGD presented the highest match ratio of all the methods (0.85), and mPTP presented the highest ratio among the tree-based methods (0.57) (Table 1).Only M. saharica was consistently recovered as an independent unit in all the methods.The other two putative new species, Mercuria sp. 1 and Mercuria sp. 2, were considered as a single group by all of the methods.Delimitation of the remaining species varied among the methods.
The ABGD method detected 12 species groups at a prior maximal distance of 0.007.The number of groups inferred by this method is nearly congruent with the number recovered by the phylogenetic analyses and by morphology-based taxonomic classifications ( 14), except for the clustering of two different species pairs: Mercuria sp. 1 and Mercuria sp. 2 and M. rolani and M. melitensis.Using ST-GMYC, the most likely model suggested eight groups [confidence interval (CI) = 2-41], of which four merged pairs of putative species suggested as independent entities by the other methods.Examples include the M. tensiftensis-Mercuria sp. 4 or M. rolani-M.melitensis groups.The bPTP analysis inferred 21 species groups, splitting the morphology-based species M. tachoensis, M. balearica, Mercuria sp. 5, M. melitensis, M. midarensis and M. targouasensis into 2, 2, 2, 2, 3 and 2 groups, respectively.The latter three species were likely oversplit considering the low support (BPP < 0.9) obtained for their suggested groups.The mPTP analysis delimited 14 putative species, which was congruent with the number of putative species identified by the phylogenetic and morphological analyses; however, it split the morphology-based species M. tachoensis into 2 groups, and clustered together the putative species Mercuria sp. 1 and Mercuria sp. 2.
Based on the information provided by the inferred phylogenies, the species delimitation methods and the observed morphological shell variation among the delimited groups (Fig. 2), we conclude that our dataset comprises 14 Mercuria species, of which five may represent new taxa.This species partition scheme was applied in the subsequent analyses of geographic range evolution.Despite being recovered as a single group by all species delimitation methods, we considered Mercuria sp. 1 and Mercuria sp. 2 to be distinct entities based on different shell morphotypes (and other anatomical characters not shown here).These differences reside mainly in shell shape (shell slender and narrower in Mercuria sp. 2), aperture size (generally larger in Mercuria sp. 1) and colour (whitish in Mercuria sp. 1 and yellowish in Mercuria sp. 2).

Divergence times and ancestral range estimation
Using the levels of evidence suggested by Kass and Raftery (1995), the Bayes factor estimation favoured the relaxed lognormal clock approach for the studied genus (BF = -0.92).Substitution rates and absolute ages were estimated for the 16S and 28S partitions of Mercuria based on a prior COI rate of 0.81 ± 0.24% (Delicado et al., 2013).From this, we estimated a substitution rate of 0.36 ± 0.09% for 16S and 0.07 ± 0.02% for 28S.The time-calibrated species tree (i.e., the MCC tree derived from *BEAST, which is depicted in Fig. 4) recovered six clades of Mercuria species and showed that the earliest divergence (i.e., the crown node) in Mercuria occurred around 7.5 Mya (HPD: 12.81-4.11Mya).Overall, phylogenetic relationships within the recovered clades  S1).Scale bar: expected change per site.resembled those inferred by the ML and BI analyses, except for the phylogenetic positions of M. tachoensis (closely related to the clade iii in the *BEAST analysis) and M. melitensis (within clade vi in the *BEAST analysis).
According to the BioGeoBEARS analysis, the best-fit model was DEC + J, with an AICc = 79.64 (Table 2).This model selection suggests that dispersal was the most commonly estimated scenario for Mercuria's biogeography.Along the time-calibrated *BEAST tree, BioGeoBEARS estimated ten dispersal events since the late Miocene: three anagenetic range expansions and seven range changes leading to speciation (i.e., cladogenetic (jump) dispersal or founder-event speciation).Results from this biogeographic model in combination with the time-calibrated *BEAST tree are depicted in Fig. 4.
The earliest ancestral lineage of Mercuria most likely diverged in the Iberian Peninsula.During the associated cladogenetic (jump) dispersal event, the ancestral lineage of the species in western Europe and the African M. tensiftensis (i.e., clades iiii) split from the ancestral lineage of the species in North Africa, Italy and Malta (i.e., clades ivvi).Nmatch = number of species delimited that exactly match a defined morphospecies; Ndelimited = total number of species delimited by the method; Nmorph = the total number of morphospecies.

Species diversity in the Mercuria dataset
Species richness of Mercuria has been based primarily on conchological and genital characters, specifically those related to the male penis and penial appendix and the female bursa copulatrix (e.g., Boeters and Falkner, 2017;Glöer et al., 2010Glöer et al., , 2015)).Based largely on these characters, 30 recognized species, 25 extant and 5 extinct, have been ascribed to this genus.However, a major problem associated with establishing the taxonomy of hydrobiids on the basis of such characters is the high level of convergence and similarities found in shell features (Bodon et al., 2001;Hershler and Ponder, 1998), which may conceal the presence of cryptic species (i.e., species that are not morphologically distinguishable from each other).On the contrary, high intraspecific variability in morphological characters has challenged the species assessment of several populations of Mercuria (Boulaassafer et al., 2018;Holyoak et al., 2017).An integrative approach that combines morphological, molecular and geographic data is therefore needed to enhance species identification and discovery in Mercuria.
For instance, the specimens from Arundel and Barking in the United Kingdom and Hoogvliet in The Netherlands, which were attributed to M. confusa (Kerney, 1999) and M. anatina (Boeters and Falkner, 2017), respectively, resolved as M. tachoensis in our molecular analyses.This finding suggests that these specimens were misidentified in the previous studies.Also, the specimens from Bidart in France, which were attributed to M. bayonnensis (Boeters and Falkner, 2017), are here identified as M. tachoensis.Though these findings highlight the need for a taxonomic revision of some Mercuria species, for robust species assignments and/or potential synonymizations, we first need to better characterize the morphology and genetics of populations living at the type localities of such species and then compare them with other collected samples.
Of the applied molecular delimitation methods, ABGD, with a match ratio value of 0.85, most closely reflected the number of putative Mercuria morphospecies, delimiting 12 of the 14 (nine previously recognized plus five new ones).The tree-based methods ST-GMYC and bPTP both presented a relatively low match ratio value and were less consistent with the number of morphology-based species.Inconsistencies between the number of morphospecies versus delimited species have been found in other studies.Several authors caution interpretation of results derived from these molecular delimitation methods for groups with low vagility and poor dispersal capacities, such as gastropods (e.g., Delicado et al., 2019;Razkin et al., 2017;Sauer and Hausdorf, 2012;Strong and Whelan, 2019), thereby contrasting the validity of the delimited species with the morphological data under an integrative taxonomic approach (Padial and De La Riva, 2010).
The main problem with using methods based on the generation of partitions (GMYC and PTP) is that they utilize information obtained through the modelling of speciation and coalescent processes.These processes can be affected by inter-and intraspecific variation (Ross et al., 2008;Zhang et al., 2013), leading to, in many cases, the delimitation of entities within species with a strong population structure (Sukumaran and Knowles, 2017).Many gastropod species have been shown to have a strong population structure (Delicado et al., 2019;Razkin et al., 2017), and according to the genetic distances observed in our study, most of the putative Mercuria morphospecies display some degree of population structure.Consistent with this, the bPTP method inferred 21 putative species of Mercuria, likely reflecting the population structure of species rather than speciation events.Indeed, only the species with little genetic structure, such as M. balearica and M. tensiftensis, were supported by all of the delimitation methods.
Despite its generally good performance, ABGD may underperform with respect to closely related species (Puillandre et al., 2012); therefore, it should not be used as a single source of evidence.For instance, a lower mean COI divergence (1.3%) was observed between the morphologically differentiated species Mercuria sp. 1 and Mercuria sp. 2 than between other hydrobiid species [e.g., 8% between Corrosella species; Delicado et al. (2012), or 3-5.5% between Hydrobia species; Wilke et al. (2000)].Despite this low level of genetic divergence, the morphological differences between Mercuria sp. 1 and Mercuria sp. 2 (see shell morphologies in Fig. 2) sufficiently support their delineation as two distinct species.A low level of genetic distance in conjunction with high morphological divergence has been detected between other molecularly closely related species of hydrobiids as a result of evolutionary processes such as adaptive radiation, dispersal to a new environment or recent speciation (Delicado et al., 2014;Falniowski et al., 2016;Stelbrink et al., 2020).In the case of these two Mercuria species, the inferred patterns could be explained by recent founder-event speciation (see Fig. 4).
Based primarily on the information provided by the phylogenetic inferences and shell observations, we consider that our dataset comprises 14 species of Mercuria (Fig. 3), of which five may represent new species.Using this approach, species known since the mid-nineteenth century, such as M. similis and M. saharica, were recovered as valid entities.Three of the species, M. tensiftensis, M. midarensis and M. targouasensis, correspond to the three clades inferred by Boulaassafer et al. (2018), though with different phylogenetic positions.In our phylogeny, sequenced populations of M. midarensis and M. targouasensis formed a well-supported monophyletic group, whereas the monophyly of the latter species was unsupported.Both results are in agreement with those of Boulaassafer et al. (2018).By contrast, M. tensiftensis showed a sister taxon relationship with the clade of M. bakeri and M. tingitana in that study, but with M. balearica in our study.Our analyses did not include the species M. bakeri and M. tingitana, which may account for the difference in the phylogenetic relationship of the aforementioned species.Also, as in Boulaassafer et al. (2018), neither of our phylogenetic inferences recovered the species M. targouasensis as a well-supported group, despite the inclusion of a greater number of species and gene fragments, which provided a better coverage of Mercuria's species diversity, genetic diversity and geographic distribution.
Uncorrected pairwise COI distances for Mercuria species ranged from 1.3% to 9.3%, with an average of 7.5%.This is similar to the situation found for the western Mediterranean species of Pseudamnicola, which have presented distances ranging from 0.5% to 10%, with an average of 6.7% (Delicado et al., 2015).By contrast, for Corrosella, Delicado et al. (2015) found an average genetic distance close to 8% (range: 5.3%-  (Hershler et al., 2003).
Highly variable genetic distances have also been observed in other Truncatelloidean groups, such as the spring snail genus Bythinella Moquin-Tandon, 1856, which occurs at relatively high elevations over glacial zones and ranged from 1.8% to 4.7% (Fehér et al., 2013), or the groundwater dwelling species of Bythiospeum Bourguignat, 1882, which ranged from 3.6% to 16.4% (Richling et al., 2017).

Time of divergence and biogeography of Mercuria: Evidence for multiple dispersal events
Our biogeographic model comparison selected the DEC model incorporating founder-event speciation (+J) as best fit for the delimited species of Mercuria (Table 2).This model has been used to predict longdistance dispersal (Matzke, 2013b), a process more concordant with the inferred age and modes of speciation within Mercuria than vicariance via range fragmentation.The latter process is not very plausible considering that, for instance, the formation of the western Mediterranean continental regions and islands (where most species of Mercuria occur today) resulted from the break-up of the continental Hercynian Belt ca. 25 Mya in the Oligocene (Rosenbaum et al., 2002), pre-dating the estimated divergence times in Mercuria.In our analyses, the selected biogeographic model (DEC + J) estimated founder-event speciation more frequently than ecoregion-wide sympatry, and vicariance (which may occur after plates break-up) was never predicted.Further supporting a dispersal scenario is the fact that no vicariant events were estimated by the second best-fitting biogeographic model (i.e., BayArea + J), which had an AIC value almost identical to that of the DEC model (see Table 2).Collectively, these findings suggest that the presence of Mercuria in the Mediterranean and Atlantic regions is likely a result of a series of dispersal and colonization events followed by subsequent isolation rather than range fragmentation, a biogeographic scenario that has been suggested for other Palaearctic freshwater molluscs (Boulaassafer et al., 2020;Neubauer et al., 2016;Pfenninger et al., 2010).
The DEC model inferred that the MRCA of the studied Mercuria species occurred in the Iberian Peninsula and diverged as early as the late Miocene (7.5 Mya;.The putative fossil record of Mercuria suggests the presence of the genus in northern Europe and the Italian Peninsula by the Oligocene and middle Miocene (Esu and Girotti, 2015a;Marquet and Lenaerts, 2008).Although these records are consistent with Mercuria's habitats and geographic distribution, their assignment to this genus is debatable.Similar shell features observed in extant species of Pseudamnicola and other non-hydrobiid genera such as Bithynia Leach, 1818 have contributed to genus misidentification in some Mercuria taxa (e.g., the taxonomic transfers of M. emilianus to Pseudamnicola or B. meridionalis to Mercuria; Boeters and Falkner, 2017).Therefore, the fossil species that have been previously classified as Mercuria should be interpreted with caution.
Although the earliest divergence of Mercuria is well supported in our phylogeny, subsequent ones are not, perhaps because of incomplete taxon sampling of the African species or the low resolution of the studied genetic fragments within this time frame.Regardless, our dataset that covers Mercuria's entire distribution range (except Corsica and Sardinia) indicates a potential jump-dispersal event of individuals from the ancestral populations of Mercuria from the Iberian Peninsula to North Africa during the late Miocene.During this period multiple sea-level oscillations resulted in Europe and North Africa being connected through land bridges between Iberia and Morocco (Strait of Gibraltar) and between Tunisia, Sicily and Italy (Strait of Sicily) (Duggen et al., 2003;Krijgsman et al., 1999;Rindi et al., 2019), which facilitated faunal interchange between the two continents when the sea level drop (Giusti et al., 1995;Sherpa et al., 2018).The subsequent changes of the sea level acted as a barrier to gene flow between populations of freshwater taxa living throughout the Mediterranean, including gastropods, and thus promoted speciation (Sherpa et al., 2018).
The Mediterranean species belonging to clades i, iii, iv and vi (see Fig. 4) diverged after the Miocene, suggesting that these clades were most likely influenced by overseas dispersal events via birds rather than the aforementioned palaeo land bridges.An example from our data that supports this hypothesis is the colonization of the Balearic Islands of Majorca and Minorca in the western Mediterranean.We found that dispersal events of M. similis and M. balearica from the Iberian Peninsula to these two islands occurred relatively recently.In fact, our results indicate that the entire clade (i.e., clade i with M. balearica + M. tensiftensis + Mercuria sp. 4) originated in the southern Iberian Peninsula and then dispersed towards the Balearic Islands and North Africa during the Pliocene.Colonization from the Iberian Peninsula towards the Balearic Islands has been reported for other freshwater groups (e.g., Abellán et al., 2009;Lázaro et al., 2009) and for land snails (Chueca et al., 2015(Chueca et al., , 2017;;Neiber et al., 2021).Interestingly though, the biogeographic origins of the Mercuria populations in Majorca differ from those inferred for Pseudamnicola, which probably colonized this island from Africa (Boulaassafer et al., 2020).
Since the late Miocene, Mercuria was estimated to have extended its distribution from the Iberian Peninsula to other Atlantic regions, probably through dispersal mechanisms other than bird vectors or land bridges.In particular, our findings indicate an anagenetic northward range expansion of populations of M. tachoensis from the Iberian Peninsula to north-western Europe.Past connections between river basins during Quaternary glacial cycles have been inferred along this pathway (Dias et al., 2014), which may have facilitated the dispersal of Mercuria in the Atlantic lowland drainages in Europe.This pattern of colonization of north-western European regions from southern glacial refugia has been proposed for other continental aquatic species living in running waters (Benke et al., 2009;Marrone et al., 2017); however, to our knowledge, it has been very rarely suggested for taxa within Hydrobiidae.Based on their phylogeographic analyses, Wilke and Davis (2000) proposed that individuals from the ancestral populations of two brackish-water hydrobiid species from northern Europe dispersed southwards during glacial periods, returning to northern waters under warmer conditions.The phylogenetic relationships of populations of M. tachoensis from different geographic regions (see Fig. 3) suggest that dispersal events occurred exclusively from southern to northern regions during the evolutionary history of the genus.
Several factors may account for the hypothesized dispersal scenario in Mercuria: 1) better access for birds at low elevations, which has been suggested to explain the inferred biogeographic patterns of the (lowland) hydrobiid genera Pseudamnicola and Ecrobia (Delicado et al., 2018;Haase et al., 2010;Vandendorpe et al., 2019), 2) greater connectivity of lower river courses (that evolved into larger rivers) during glacial times (Dias et al., 2014), which could have promoted range expansion in M. tachoensis along the Atlantic lowland drainages of Europe and in M. balearica in the Mediterranean and 3) greater ecological plasticity of the species, which may also account for the relatively high dispersal and jump values observed in our study [as concluded by Delicado et al. (2018) for Pseudamnicola].
Based on our phylogenetic evidence and morphological observations, we conclude the diversity of the Mediterranean clades of Mercuria is higher than previously thought, whereas in the Atlantic ones, it is lower.The numerous long-distance dispersal events that occurred across the Mediterranean in the Plio-Pleistocene likely account for an increase in cladogenetic speciation of the genus in this region, whereas the (anagenetic) northward range expansion of M. tachoensis, and the lack of co-occurrence with other species seems to explain the presence of Mercuria in the Atlantic lowland drainages of Europe.These findings support the general hypothesis that dispersal modes have important consequences on the evolution and species richness of different organismal groups (Ricklefs, 2004;Shurin and Allen, 2001).For freshwater organisms, variations in dispersal strategies have been attributed to differences in habitat types, elevational preferences or body sizes (Hof et al., 2008;Ribera et al., 2003;Shah et al., 2017;Shurin, 2000).Habitat types (e.g., headwater springs, lakes, subterranean waters, brackish waters) linked to factors such as elevation or latitude differ from each other in the extent of their geographic connectedness, natural environmental conditions and anthropogenic impact (Miller et al., 2018), which in turn, affect the survival and dispersal of populations (Hof et al., 2008).In our study, we observed that most of the analysed Mediterranean species of Mercuria are restricted to a single ecoregion, and an anagenetic range expansion was inferred only for M. tachoensis and M. balearica.The persistent low level of connectivity between Mediterranean river basins caused by sea-level lowstands during glacial periods (Dias et al., 2014) may have hindered short-distant dispersal between adjacent areas during the Quaternary, thereby contributing to the high levels of endemicity found among extant Mediterranean congeners.Given all of this evidence, we conclude that the tectonic and climatic history of the Mediterranean region affected potential dispersal pathways and limited the range expansion of the Mercuria species.However, their occurrence in lowlands facilitated overseas dispersal by avian vectors and, ultimately, speciation.
Our biogeographic reconstruction method has several limitations in estimating ancestral areas.For instance, given the low support at the basal branches, we were not able to accurately infer the ancestral areas from which the groups originated during their evolution.An example of this is the earliest ancestral area estimations within Mercuria: the low probabilities obtained in the phylogenetic reconstruction do not support the selected areas as those that can be reasonably related to the first dispersal events in the history of the genus.Also, the absence of many of the North African species of Mercuria, especially those from Algeria and Tunisia, might have introduced a bias in our analyses.Despite these caveats, multiple more recent speciation events promoted by dispersal were well supported across the inferred phylogenies.Therefore, in general, the biogeographic history of the Mercuria species studied here can be inferred with confidence and compared using our calibrated phylogeny.

Conclusions
Based on the molecular phylogenetic methods used here, the species diversity of Mercuria is higher than previously thought in the sampled areas of the Mediterranean but lower in the Atlantic ones.The relatively young age of the inferred species and the ancestral range reconstruction suggest that founder-event speciation since the late Miocene period played a major role in shaping the diversity and distribution of the Mediterranean clades, whereas episodes of postglacial northward colonization from Iberian refugia by the species M. tachoensis may explain the current presence of the genus in Atlantic lowlands.The greater species richness and higher levels of endemicity found in the Mediterranean freshwater habitats may have been a consequence of multiple factors such as palaeo-disconnected drainage basins driven by climate change; better access for avian dispersal vectors in lowlands, where Mercuria typically occurs; the high adaptability of some species of Mercuria to live in different habitat types and the role of southern regions as glacial refugia.The here inferred past and present distribution of Mercuria species, as well as their driving causes, can provide useful information for the conservation of freshwater species in the future.

Author statement
JPM, DD and MR conceived and designed the study; JPM and FG collected material; JPM performed identifications, which were verified by MR; JPM and FG conducted the laboratory work; JPM and DD curated the sequence data; JPM analysed the data, with supervision of DD; JPM produced the figures, with help of MR and DD; all authors contributed to the interpretation of data; JPM and DD drafted the initial manuscript and all authors contributed to later versions.All authors read and approved the final manuscript.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Map showing the locations of the Mercuria populations included in the molecular study.Bottom boxes represent enlarged sections of two regions of the Iberian Peninsula.Black lines on the general map contour freshwater ecoregions sensu Abel et al. (2008), which were labelled with the same letter scheme as in Fig.4.Code names for each of the studied population, species assignment and other information are provided in the Supplementary Material (TableS1).

Fig. 3 .
Fig. 3. Bayesian phylogenetic tree of Mercuria species based on the concatenated dataset (COI, 16S and 28S).Bootstrap values from ML analysis and Bayesian posterior probabilities are provided above the branches when greater than 75% and 0.95, respectively.ABGD, distance-based automatic gap discovery; GMYC, singlethreshold generalized mixed Yule-coalescent; bPTP, Bayesian approach of the Poisson Tree Processes; mPTP, multi-rate Poisson Tree Processes.bPTP bars are displayed with Bayesian support values.Species-level groups are indicated by numbers on the bars of the right column.Label abbreviations are explained in Supplementary Material (TableS1).Scale bar: expected change per site.
The next inferred sequence of dispersal events involved first a cladogenetic (jump) dispersal of the progenitors of M. tachoensis from southern to northern Iberia around 5.79 Mya (HPD: 9.80-3.11Mya) and then a range expansion along the M. tachoensis lineage to northern Europe.Another cladogenetic (jump) dispersal event occurred about 1.64 Mya (HPD: 3.24-0.42Mya) and resulted in the divergence of the Tunisian species Mercuria sp. 5 and the Italian species Mercuria sp. 3. The expansion of the most recent common ancestor (MRCA) of the Moroccan species M. targouasensis + M. midarensis and the Maltese species M. melitensis occurred around 2.43 Mya (HPD: 4.38-1.08Mya).The remaining Mercuria species, represented by clade i in the phylogenetic reconstruction (M.similis, Mercuria sp. 1 and Mercuria sp.2; see Fig. 4),

Fig. 4 .
Fig. 4. Biogeographic history of Mercuria species based on the lognormal uncorrelated relaxed clock approach and the combined (COI, 16S and 28S) dataset.An external COI substitution rate calculated from other hydrobiid substitution rates based on biogeographic events was used for the divergence time estimation.Blue bars represent 95% HPD (highest posterior density) intervals and black dots on nodes correspond to well-supported species clades (BPP > 0.95).Coloured branches indicate expansion via anagenesis (pink) or jump-dispersal (yellow).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Mercuria sp. 3 + Mercuria sp. 5 + M. saharica + M. melitensis + M. rolani and (v) M. midarensis + M. targouasensis.Although all species and clades were well supported, the relationships among them differed in the different analyses, and the basal phylogenetic relationships were not well supported.

Table 1
Match ratio (MR) of each of the species delimitation methods.

Table 2
Comparison of the models estimated in BioGeoBEARS for Mercuria species.The best-fit model selected was based on the AICc.