Disparate dispersal limitation in Geomalacus slugs unveiled by the shape and slope of the genetic–spatial distance relationship

Long-term dispersal ability is a key species’ trait constraining species ranges and thus large-scale biodiversity patterns. Here we infer the long-term dispersal abilities of three Geomalacus (Gastropoda, Pulmonata) species from their range-wide genetic–spatial distance relationships. This approach follows recent advances in statistical modelling of the analogous pattern at the community level: the distance decay in assemblage similarity. While linear relationships are expected for species with high long-term dispersal abilities, asymptotic relationships are expected for those with more restricted mobility. We evaluated three functional forms (linear, negative exponential and power-law) for the relationship between genetic distance (computed from mitochondrial cox1 sequences, n = 701) and spatial distance. Range fragmentation at present time and at the Last Glacial Maximum was also estimated based on the projection of climatic niches. The power-law function best fit the relationship between genetic and spatial distances, suggesting strong dispersal limitation and long-term population isolation in all three species. However, the differences in slope and explained variance pointed to disparities in dispersal ability among these weak dispersers. Phylogeographic patterns of Geomalacus species are thus largely driven by the same major process (i.e. dispersal limitation), operating at different strengths. This strong dispersal limitation results in geographic clustering of genetic diversity that makes these species highly vulnerable to genetic erosion due to climate change

Research 1230 associated with past environmental change is a challenging but essential task to enhance our understanding of biodiversity patterns (Pearse et al. 2018).Different methodological approaches have been proposed to achieve this goal, ranging from predictive models unveiling incomplete range filling (Svenning and Skov 2004), correlative studies with taxon's traits indicative of dispersal limitation (Gómez-Rodríguez and Baselga 2018) or intraspecific genetic divergence within geographic ranges (Petit et al. 2003).In fact, inference about processes of range expansion has been an integral part of phylogeography since its inception (Avise 2000, Hewitt 2000), even if the estimation of specific range expansion rates has been more elusive (Clobert et al. 2012).Towards this goal, the combination of large-scale genetic assessments with alternative ecological approaches opens exciting research avenues that can foster our understanding of diversity patterns and associated processes (Marske et al. 2013).Such an integration of methodological frameworks across disparate subdisciplines is both at the origin (Leitão et al. 2020) and future of our understanding of ecological change across large spatial, temporal and taxonomic scales (Rapacciuolo and Blois 2019).
At local scales, many complex demographic and genetic approaches can inform about contemporary dispersal among populations (Broquet andPetit 2009, Cayuela et al. 2018).At larger spatiotemporal scales, historical dispersal processes and long-term range expansion can be inferred from geographic patterns of putatively neutral molecular markers, such as mitochondrial DNA, mtDNA (Avise 2000), on the basis that, under no external restriction, the geographic distribution of mtDNA lineages broadens over time (Neigel and Avise 1993).This key idea has been formalised in a breakthrough conceptual synthesis of community ecology (Vellend 2010(Vellend , 2016)): the spatial distribution of neutral genetic variants within a species would be mostly driven by migration and drift, while other processes such as niche filtering would constrain spatial ranges at the species level and above.Thus, the geographic distribution of intraspecific genetic diversity depends on the species' dispersal ability to expand the lineages' ranges within ecologically suitable conditions.Following this rationale it is possible to infer that a species has high dispersal ability when weak geographic genetic structure has been observed among its populations (Schiebelhut and Dawson 2018), however, if dispersal ability is low, most genetic variants would be restricted to subsets of the species range, hence leading to strong geographic genetic structure (Bálint et al. 2011).As a result, long-term dispersal ability could be modelled using the relationship between mtDNA genetic distance and spatial distance, analogously to how the decay in assemblage similarity with spatial distance (i.e.distance decay patterns) informs about species dispersal at the community level (Chust et al. 2016).
Given the conceptual similarities between models for the genetic-spatial distance relationship and the distance decay pattern at the community level, recent advances in statistical modelling of decay in assemblage similarity with spatial distance (Nekola and McGill 2014) can be borrowed.
A refinement of models describing the relationship between genetic differentiation at neutral loci and geographic spatial distance would allow, at the least, a more accurate characterization of geographic patterns, a critical step in macroecological studies (McGill 2019) that allows the inductive generation of testable hypotheses (Currie 2019).One of the most relevant features of the genetic-spatial distance relationship would be its functional form, in the same way the functional form of distance decay patterns is used to infer the processes behind differentiation at the community level (Nekola and McGill 2014).Notably, at the genetic level, Hutchinson and Templeton (1999) also pointed to the use of the functional form in their seminal paper focused on contemporary dispersal processes at local scales.Two basic predictions can be made: a continuous linear relationship is expected when there is no dispersal limitation while an asymptotic relationship would be expected in weak dispersers because gene flow due to dispersal is only expected at the shortest distances while gene drift and mutation processes would be prevalent where the statistical relationship flattens out (Hutchison and Templeton 1999).While these predictions were elaborated in the context of population genetics and, thus, mostly targeting contemporary dispersal processes (Monaghan et al. 2002, van Strien et al. 2015), larger spatio-temporal patterns can be also described under this framework.Moreover, analogously to how the slope of the distance decay of assemblage similarity has been shown to be a reasonable surrogate of dispersal ability at the community level (Qian 2009, Saito et al. 2015, Gómez-Rodríguez and Baselga 2018), the slope of the genetic-spatial distance relationship could inform about species' long-term dispersal ability (Peterson and Denno 1998, Lester et al. 2007, Chust et al. 2016).Thus, a major aim of this paper is to use the tools developed for the assessment of distance decay patterns to characterise the shape and quantify the slope of the genetic-spatial distance relationship across a species' distribution in order to assess its long-term dispersal ability after historic range contractions, at least in relative terms.
Terrestrial molluscs offer a paradigmatic example of how low mobility and strong dependence on climatic conditions (e.g.humidity) shape biodiversity and phylogeographic patterns (Pfenninger and Posada 2002, Hausdorf 2006, Gómez-Rodríguez et al. 2019).This taxonomic group is particularly interesting due to the existence of deep evolutionary lineages within species with no morphological variability, that have been interpreted either as cryptic diversity (Dépraz et al. 2009, Razkin et al. 2017) or infraspecific variability (Thomaz et al. 1996, Watanabe and Chiba 2001, Pinceel et al. 2005).In the case of Geomalacus, a genus endemic to the Iberian Peninsula and Ireland (Reich et al. 2015), previous contributions have studied its phylogeography, speciation or habitat requirements using novel methodological approaches (Patrão et al. 2015, Cunha et al. 2017).However, these studies followed an outdated taxonomy based on limited morphological evidence that recognized four species (Castillejo et al. 1994).A recent comprehensive taxonomic review, based on > 1500 1231 specimens, has identified only three Geomalacus species based on distinct genital characters that are suggestive of reproductive isolation (Castillejo et al. 2019).
Here we 1) assess long-term dispersal ability of Geomalacus species based on the shape and slope of the relationship between genetic and spatial distances using range-wide mtDNA sequence data and 2) predict potential current and past distributions of each species from their realized niche (Nogués-Bravo 2009), in order to understand how the interplay between dispersal limitation and climatic constraints may have shaped their current distribution in the Iberian Peninsula.We expect an asymptotic relationship (negative exponential or power-law) to better fit the data given the low mobility of these taxa (Barker 2001) and the highly structured genetic variation reported in one of the species (G.maculosus, Reich et al. 2015), which probably results from the Iberian Peninsula having harboured multiple glacial refugia at very small scales (Gómez andLunt 2007, Schmitt 2007).Projections of the species' climatic niches at the LGM are therefore used to test the assumption of multiple glacial refugia.These projections will show whether the past distribution of Geomalacus was likely to be fragmented or not and, therefore, if range expansion has likely occurred from one or multiple refugia.Finally, to ensure our analyses of the spatial distribution of clades involve discrete evolutionary lineages, we validate the monophyly of the species revised in Castillejo et al. (2019) as a crucial previous step, by contrasting the morphological and evolutionary delimitation of species.

Field sampling, morphological identification and DNA sequences
Geomalacus specimens (n = 144) for this study were collected in a large-scale assemblage-focused survey of terrestrial molluscs from 20 localities (Fig. 1) in the Iberian Peninsula (Gómez-Rodríguez et al. 2019).Details on the dataset and survey can be found in Supplementary material Appendix 1.Additionally, four individuals from two other localities (Serra de Caldeirao [FAR] and Fuente-Dé/Picos de Europa [PIC-W]) were also considered.Specimens were identified based on external morphology and genitalia anatomy indicative of reproductive isolation, following a recent taxonomic revision by Castillejo et al. (2019).Short morphological diagnosis for each species can be found in Supplementary material Appendix 2.
DNA was extracted and the barcode region (cox1-5′) successfully amplified for 132 specimens.Details on the molecular lab methodology can be found in Gómez-Rodríguez et al. (2019).In brief, a piece of foot tissue (≈ 5-10 mg) was used for genomic DNA extraction with DNeasy Blood & Tissue Kit

1232
(Qiagen, Germany) and LCO/HCO primers (Folmer et al. 1994) were used for amplification of all specimens.For 20 specimens representing different cox1-5′ evolutionary lineages in a preliminary analysis, a 1878-base pair fragment of the nuclear 18S rRNA gene was also amplified with primers 18S5′, 18S5.0rw,18Sai, 18Sbi, and 18sa2.0 and 18S3′I following Shull et al. (2001).Purification and sequencing in both directions with an ABI 3730xl DNA Analyzer were done by StabVida (Portugal).The GenBank accession numbers for these sequences can be found in Supplementary material Appendix 3.
Additionally, Geomalacus cox1-5′ sequences were also downloaded from GenBank in November 2017.The following species were represented, according to GenBank identification: 299 sequences of G. maculosus, 104 sequences of G. anguiformis, 99 sequences of G. malagensis and 67 sequences of G. oliveirae; totalling 569 sequences.Of these, 31 sequences were from Ireland, the only region outside the Iberian Peninsula were Geomalacus is found.Additionally, 78 sequences of a 527-base pair fragment of nuclear 18S rRNA were also downloaded from GenBank.These 18S sequences had been amplified for specimens that also had a cox1-5′ sequence in GenBank.The locality of the sequences was identified in the original sources (Patrão 2013, Reich et al. 2015) except for 18 sequences belonging to G. maculosus, that were lacking that information.The inclusion of GenBank sequences in this study aims to increase its spatial coverage and genetic diversity.The full list of downloaded sequences is provided in Supplementary material Appendix 3.

Phylogenetic analysis and haplotype networks
Sequences were aligned with transAlign (Bininda-Emonds 2005) and identical sequences of cox1-5′ (i.e. the same haplotype) were collapsed into a single sequence with a custom perl script ensuring that, whenever possible, the 'reference haplotypes' included in phylogenetic analyses belonged to specimens with an 18S sequence or sequenced in this study to ensure their identification by an expert taxonomist (Prof.J. Castillejo).In three cases, the same haplotype was represented by a specimen studied here plus a specimen from GenBank that also had an 18S sequence.In that case both records were included in the reference haplotype dataset in order to account for both morphological and nuclear information.
A total of 175 'reference haplotypes', including an outgroup sequence of Arion ater (GenBank HQ660058), were included in a phylogenetic tree of the cox1-5′ run in Beast 1.8 (Drummond and Rambaut 2007) in Cipres (Miller et al. 2010).Previously, jmodeltest (Posada 2008) with three substitution schemes had been used to identify the most likely evolutionary model for cox1-5′ and, independently, for 18S.Model priors were set in BEAUti v.1.8as a GTR + G + I and uncorrelated relaxed clock with mcmc = 1 000 000 000 and log every 100 000.The analysis of the mcmc with Tracer (Rambaut et al. 2018) showed high values of Effective Sample Size (all parameters with EES > 7750) and 5000 (50%) trees were discarded in TreeAnnotator to build a maximum clade credibility tree with median heights.A species tree, based on the cox1-5′ and 18S markers, was also built to assess deeper evolutionary relationships at the species level, using the StarBeast template in BEAUti v2.4 (Bouckaert et al. 2014) under a relaxed lognormal model with GTR + G + I for cox1-5′ and GTR + G for 18S (mcmc = 1 000 000 000; log = 5000).Preliminary cox1-5′ trees, based on specimens morphologically studied here, corroborated the monophyly of the three Geomalacus species revised in Castillejo et al. (2019).Preliminary trees also suggested that very closely related or identical GenBank sequences had been given a different Geomalacus species name or followed the outdated nomenclature/identifications in Castillejo et al. (1994).These GenBank sequences were renamed accordingly prior to the StarBeast species tree (Supplementary material Appendix 3).DensiTree v2.2 (Bouckaert 2010) was used to represent the species trees.Additionally, a supplementary maximum likelihood phylogenetic analysis, performed with RAxML on the same dataset, is provided in Supplementary material Appendix 4.
To visualize the relationships among haplotypes in their spatial context, haplotype networks were built with command haplonet() in 'pegas' using an infinite site model (Paradis 2010).To aide visual representation, nearby sites (< 40 km apart) were grouped together, yielding a total of 32 localities (Supplementary material Appendix 5).Additionally, to minimize the overrepresentation of the two G. maculosus haplotypes from Ireland, which were focus of Reich et al.'s (2015) study, only three and one sequences were included, depending on their relative abundance.

Statistical fitting of the genetic versus spatial distance relationship
To assess whether genetic distance was related to Euclidean spatial distance within each species and, more importantly, to characterize the shape of such relationship as an indicator of long-term dispersal ability, three different GLM models were evaluated based on AIC.The three GLMs included genetic distance (i.e.proportion of sites that differ between each pair of sequences) as the response variable and spatial distance as the predictor variable.Each model was fitted with a different functional form: linear (Gaussian error and identity link), exponential (Gaussian error and log link) and power-law (Gaussian error, log link and log-transformed spatial distance).Exponential and power-law models were fitted using the decay.model()function in 'betapart' (Baselga and Orme 2012).To assess whether model slopes differed among Geomalacus species, we bootstrapped the parameter (1000 replicates) using the function boot.coefs.decay() in 'betapart'.Significance was estimated from pairwise comparisons of the bootstrapped distributions of slopes, computing p values as the proportion of times in which the bootstrapped parameters were larger (or smaller) in one of species than in the other.Sequences from Ireland were excluded from these analyses, as they would be an outlier case given the long distance between Ireland and the Iberian Peninsula and because 1233 the presence of Geomalacus in Ireland seems to be the result of accidental anthropogenic introduction rather than natural range expansion (Reich et al. 2015).

Climatic niche and potential distribution
Climatic niches and their similarities among Geomalacus species were estimated from their current distributions using environmental niche models in package 'humboldt' in R (Brown and Carnaval 2019).This software allows the characterization of environmental niches under the assumption that most species' contemporary distributions may be in nonequilibrium state, as it is expected in species with low dispersal ability.Niche similarity among species was assessed with function humboldt.doitall()for each pair of species and represented in relation to the accessible environmental space for each species.Similarity values range between 0 and 1: a value of 1 equals identical niches and a value 0 is completely different niches.A buffered minimum-convex-polygon (buffer distance = 50 km) was applied to 10 000 background points, allowing non-analogous environments in the comparison among species' niches.The variables used to describe the climatic niche were annual mean temperature [Bio1], maximum temperature [Bio5], minimum temperature [Bio6], annual precipitation [Bio12], precipitation of wettest quarter [Bio16] and precipitation of driest quarter [Bio17]) at 30 arc sec resolution from CHELSA (Karger et al. 2017).These variables represent standard measures of yearly average and extremes in temperature and precipitation, as commonly used in biogeographic studies.In this particular case, temperature extremes and water availability are also key ecophysiological constraints for terrestrial molluscs (Nicolai and Ansart 2017).CHELSA is a climate data set for the earth land surface based on a quasi-mechanistical statistical downscaling global reanalysis and global circulation model output.
Species distribution models at present time and at the Last Glacial Maximum (LGM) were also built using CHELSA variables.Predictions for species distributions at past time periods were used to identify regions of climatic stability where a species may have persisted over time (Alvarado-Serrano and Knowles 2014), independently of its dispersal ability.Following recommendation in Varela et al. (2015), the potential distribution range of each morphological species at the LGM (≈ 20 ka ago) was projected considering four different models of paleoclimate: NCAR-CCSM4, MIROC-ESM, MPI-ESM-P and IPSL-CM5A-LR.Bioclimatic congruence (Morales-Barbero and Vega-Álvarez 2019) among paleoclimatic models for the study area was also assessed with Pearson correlation.Two different modelling approaches were followed using default settings: bioclim and maxent in package 'dismo' (Hijmans et al. 2011).All records in our dataset were used for model calibration using current climate data and, in the case of maxent, 1000 background points were randomly generated in the Iberian Peninsula.After calibration based on current species distribution and climate data, the species' niches were independently projected both to the present and LGM climate conditions.

Phylogenetic relationships and spatial distribution of genetic diversity (haplotype networks)
Geomalacus species delimited by Castillejo et al. (2019) show unequivocal correspondence with three clades of different age in the phylogenetic tree based on cox1-5′ sequences for specimens studied here (Fig. 2).However, the G. squammantinus and G. anguiformis clades also include interspersed GenBank sequences attributed to different names in their original sources (see details in Supplementary material Appendix 3).Once these GenBank sequences are relabelled accordingly, the number of sequences is: G. anguiformis = 117 sequences; G. maculosus = 364 sequences and G. squammantinus = 220.The species tree (StarBeast template), based on cox1-5′ and 18S markers, shows G. squammantinus and G. maculosus to be sister species, with G. anguiformis being more distantly related.The same major lineages were also recovered using maximum likelihood phylogenetic analysis (Supplementary material Appendix 4).
Geomalacus anguiformis, the species with the smallest distribution range, was only found in the southwest of the Iberian Peninsula while Geomalacus maculosus occurs in the northwest of the Iberian Peninsula (Fig. 1a).Geomalacus squammantinus shows a fragmented distribution in a wide area from the southwest to the centre of the Iberian Peninsula.Haplotype networks showed strong spatial structure of genetic diversity within each species, especially in the case of G. squammantinus, with most haplotypes restricted to a single locality where closely related haplotypes also occur (Fig. 1).In few cases, several distantly related haplotypes also occurred in the same locality (e.g.G. anguiformis sequences in Gibraltar, GIB; or G. squammantinus sequences in Monchique, MON).

Genetic versus spatial distance assessment
The range of intraspecific genetic distances largely differed among species, with divergence reaching up to 19.5% within G. squammantinus, while only to 6.3% within G. anguiformis (Fig. 3a).Genetic distance was significantly related to spatial distance in the three species (Fig. 3), and the power-law model fitted the data better than either the exponential or linear models (Table 1).Among species, the fit was best in G. squammantinus (Pseudo-R 2 = 0.71) and worst in G. anguiformis (Pseudo-R 2 = 0.09).Slope values also suggested a stronger effect of spatial distance in G. squammantinus (b = 0.030) than in G. maculosus (b = 0.017) or G. anguiformis (b = 0.002).All slopes were significantly different among Geomalacus species (p < 0.001).A gap is observed in the genetic distances in the G. maculosus and G. anguiformis scatterplots 1234 1235 (y-axis in Fig. 3b-c), pointing to deeper spatial genetic structure in these species.

Climatic niche and potential distribution
Geomalacus maculosus and G. anguiformis show markedly allopatric distributions while G. squammantinus and G. anguiformis show partial sympatric distributions in the southwest of the Iberian Peninsula (Fig. 1a).Niche similarity between G. anguiformis and G. maculosus was lower (0.003) than between G. maculosus and G. squammantinus (0.112) or G. anguiformis and G. squammantinus (0.06).Notably, and despite the low value of niche similarity, the climatic niche of G. anguiformis was largely nested within the one of G. squammantinus (see Supplementary material Appendix 6 for a graphical representation of niche overlap).
Geomalacus anguiformis and G. maculosus show concordance between their observed distributions (i.e.recorded data) and the potential distributions predicted from climate envelope models, although the potential distribution of G. maculosus would not be limited to the Iberian Peninsula (Fig. 4).On the contrary, the observed distribution of G. squammantinus is more restricted than its potential distribution; note the unfragmented area of suitable climatic conditions Table 1.Comparison of three GLM models (linear, exponential and power-law) assessing the shape of the relationship between genetic and spatial distance for each Geomalacus species.Models are evaluated according to their AIC and the one with lowest value is shown in bold.The variance explained by each model (Pseudor-R 2 ) is also shown.1237 observed northwest of its current distribution range.Bioclim and maxent models predict similar current potential distributions (Fig. 4, Supplementary material Appendix 7).Bioclim projections at the Last Glacial Maximum show G. anguiformis to have been practically absent in the Iberian Peninsula, maybe with some suitable locations along the coast according to the NCAR-CCSM4 paleoclimate estimation (Fig. 4, Supplementary material Appendix 8).Although paleoclimatic estimations were largely congruent for the study area (Supplementary material Appendix 9), some differences were observed in projected distributions for G. maculosus and G. squammantinus.In G. maculosus, despite differences in extension, all models predict a contraction and fragmentation of its potential distribution at the LGM, pointing to the existence of multiple microrefugia for this species in the Iberian Peninsula.In G. squammantinus, the NCAR-CCSM4 projection shows a large area of continuous suitable habitat in the centre and southwest of the Iberian Peninsula (Fig. 4), with a northern limit at the LGM very similar to that of the current species distribution.The models based on the other paleoclimate estimations suggest a more restricted distribution in the south of the Iberian Peninsula (Supplementary material Appendix 8).In the case of maxent models, wider LGM distributional ranges are predicted (Supplementary material Appendix 7), as expected from limitations of maxent models to project distributions to novel climatic conditions (Merow et al. 2013).

Discussion
Geomalacus species in the Iberian Peninsula offer a paradigmatic example of how long-term dispersal ability may differ among very closely related species.Their strong spatial genetic structure can be interpreted as evidence for weak long-term dispersal ability from glacial refugia, this finding is further supported by the functional shape (i.e.power-law) of their genetic-spatial distance relationships.An asymptotic shape is in agreement with ongoing long-term isolation in moderately distant populations, where genetic differentiation due to drift would be predominant.However, while the three species can be considered poor dispersers, differences in the rate at which genetic distance increases with spatial distance (i.e.function slope) and goodness of fit (i.e.explained variance) also point to marked differences in dispersal ability among them.Our results have two important implications: 1) Geomalacus species are revealed to have different ages and evolutionary histories and 2) the varying strength of long-term dispersal limitation among the three species is identified as a key factor in structuring phylogeographic patterns.
In a wider context, our results exemplify how long-term dispersal limitation can be inferred from the asymptotic (i.e.power-law) relationship between genetic and spatial distance.The asymptotic relationship suggests that lineage expansion (after range fragmentation and contraction) has only occurred at the very short distances, where a steep positive relationship is observed.This relationship tends to flatten out at distances > 100 km (i.e.asymptotic region of the curve, indicative of negligible gene flow), thus pointing to a large number of relatively close populations evolving independently mostly due to mutation and genetic drift.Such genetic differentiation due to drift can be explained by long-term isolation of populations within geographically separate refugia (Bennett and Provan 2008), a likely scenario in the Iberian Peninsula (Gómez and Lunt 2007).Therefore, the asymptotic relationship between genetic and spatial distance in Iberian Geomalacus species can be interpreted as the genetic signature of historic processes (i.e.range contraction into glacial refugia), still detectable today due to the weak long-term dispersal ability of these species.
The parameterization of asymptotic functions describing the relationship between genetic and spatial distance not only allows the inference of dispersal limitation, but also quantifying its relevance and strength.Goodness of fit and slope can shed light on both the existence of additional explicative factors and the relative strength of long-term dispersal limitation, respectively.Here, long-term dispersal limitation is stronger in G. squammantinus (i.e.steeper slope and high explained variance) than in G. anguiformis (i.e.low explained variance suggests factors other than space explain its genetic structure).Although the interest of significant genetic-spatial distance relationships has been debated (Guillot et al. 2009, Diniz-Filho et al. 2013) and we acknowledge that our study may not explicitly model the process of range expansion, we want to stress the value of comprehensive and alternative characterizations of geographic genetic patterns, which may aide our understanding of long-term mechanisms structuring diversity (Marske et al. 2013).
Assessing long-term dispersal using this type of parameterization can complement approaches with similar goals, such as range filling estimation from climatic niche models (Svenning and Skov 2004).This would be particularly relevant in weak dispersers, for which realized niches used to estimate range filling may be largely different from fundamental niches due to strong dispersal limitation.Similarly, in species with historic range fragmentation due to multiple glacial refugia, weak long-term dispersal limitation may not be easily identified from comparing its current distribution with suitable climatic areas, as it would likely be the case of G. maculosus in this study.We acknowledge that our estimation of dispersal limitation could be considered coarse, but it also has the advantage of being feasible for multiple species for which georeferenced sequences are available in public databases (Gratton et al. 2017) and thus would potentially allow comparisons among varying organism groups of different geographic regions and ecological guilds (Bálint et al. 2011).Besides, long-term dispersal ability can be used as an indicator of contemporary dispersal ability, a species' attribute that is lacking for many species, mostly active invertebrate dispersers (Clobert et al. 2012), and thus, ameliorate the problems associated to assessing species' vulnerability to anthropogenic change (Franklin 2010).
1238 Historic range contractions and glacial refugia have been here inferred from paleodistributional models, a common approach in phylogeographic studies (Richards et al. 2007, Alvarado-Serrano andKnowles 2014).Taken altogether, our results suggest that the genetic structure of Geomalacus species in the Iberian Peninsula is the result of different evolutionary histories intertwined with different dispersal abilities: an old species with very limited dispersal ability that mostly occupies climatically-stable areas (G.squammantinus); a species that has been able to expand its distribution from several LGM refugia (G.maculosus) and a species that probably arrived to the Iberian peninsula after the LGM, with relatively higher dispersal ability and lack of clear spatial structure (G.anguiformis).Deep intraspecific lineages with strong geographic structure are common in southern Europe and have been associated to limited range expansion following range fragmentation in climatic refugia during glaciations (Hewitt 2001).In the case of terrestrial molluscs, large intraspecific variability and range fragmentation have been associated with glacial refuges and slow range expansions, mostly in the case of snails (Pfenninger and Posada 2002, Dépraz et al. 2008, Scheel and Hausdorf 2012).Thus, low dispersal ability in G. squammantinus would explain why its current distribution is smaller than its potential distribution, as inferred by niche models here and in a previous study (Patrão et al. 2015), and the strong phylogenetic structure which is usually associated to climatically stable areas (Carnaval et al. 2009).Geomalacus anguiformis appears to have arrived in the study area relatively recently given the likely lack of suitable climatic conditions at the LGM and its very recent diversification (based on mtDNA data), despite having diverged from the G. squammantinus-G.maculosus clade a long time ago.However, further exploration of its relationship with Letourneuxia, the sister north African genus (Cunha et al. 2017), would be needed to fully understand its likely origin outside the Iberian Peninsula.
This study exemplifies how an important biological trait (i.e.long-term dispersal ability) can be easily assessed from genetic data usually collated for phylogeographic studies.Although geographic patterns in mitochondrial DNA had been previously used to qualitatively assess dispersal ability (Lourie et al. 2005, Papadopoulou et al. 2009, Fraser et al. 2018), here we follow an statistical framework analogous to that of distance decay analysis of community similarity to provide quantitative estimates of long-term dispersal ability, at least in relative terms.We propose that combined interpretation of the shape, slope and fit of the relationship between spatial and genetic distance informs about range expansion rate after glacial retreat and, thus, how dispersal ability has shaped the genetic structure within a species in the long term.We show that phylogeographic patterns in Geomalacus species appear to be largely driven by the same major process (i.e.dispersal limitation), although operating at different strengths in these very closely related species.This strong dispersal limitation results in geographic clustering of genetic diversity, making these species highly vulnerable to genetic erosion due to climate change (Bálint et al. 2011, Pauls et al. 2013, Carvalho et al. 2017).Some lineages are likely to become extinct because they will not be able to shift their ranges as a response to climate change.Assessing dispersal limitation and spatial genetic structure is therefore critical to forecast the potential impact of climate change on diversity.

Figure 1 .
Figure 1.Distribution map of Geomalacus species (a).Different shapes and main colours are used for each species: G. anguiformis (triangles/ green), G. squammantinus (squares/blue) and G. maculosus (circles/yellow).Localities with closely related haplotypes (genetic distance < 5%, see haplotype networks) have the same colour shade.Size is proportional to the number of different haplotypes in the locality.Haplotype networks for each species: Geomalacus anguiformis (b); G. squammantinus (c) and G. maculosus (d).Size is proportional to the number of haplotypes and colours correspond to different localities.Numbers represent the number of mutations and dashed lines are used when the number of mutations is > 5%.Note that haplotypes of the same locality have been connected regardless of the number of mutational steps (thin dashed line).* To aide visual representation, nearby sites (< 40 km apart) were grouped together (Supplementary material Appendix 5) and the number of sequences from Ireland reduced (see main text).

Figure 3 .
Figure 3. (a) Density plot of uncorrected pairwise genetic distances (i.e.proportion of different sites) within each Geomalacus species.(b-d) Power-law relationship between spatial distance (km) and uncorrected pairwise genetic distance among specimens of each Geomalacus species: G. anguiformis (b); G. maculosus (c) and G. squammantinus (d).The power-law function is plotted using command plot.decay() in 'betapart'.An outlier sequence (mac_ANC_KM102947) has been removed from G. maculosus scatterplot for visual purposes, but not from the analyses. 1236