Micro-endemic species of snails and amphipods show population genetic structure across very small geographic ranges

Understanding variation in population genetic structure, even across small distances and for species with extremely limited ranges, is critical for conservation planning and the development of effective management strategies for imperiled species. Organisms that occupy the same geographic extent can maintain different population structures, ranging from highly diverged to panmictic. Such differences can result from differences in biological characteristics such as dispersal ability or demographic history. We used microsatellite loci to evaluate population genetic structure and variation of four desert spring invertebrates having high to low dispersal ability: the lung snail Physa acuta, two species of gilled snails (Juturnia kosteri and Pyrgulopsis roswellensis; family Hydrobiidae) and the amphipod Gammarus desperatus. The study location represents entire species ranges for the micro-endemic hydrobiids and G. desperatus, while P. acuta is ubiquitous throughout much of North America. We found little evidence of significant population genetic structure for P. acuta and J. kosteri, but much more for P. roswellensis and G. desperatus. Our results demonstrate differences in habitat preference and/or dispersal ability between the species. This information provides insight into how gene flow shapes varying population genetic structure between species across small spatial scales (<100 km2). Most importantly, our results suggest that conservation agencies should not consider these micro-endemic species to be composed of single populations, but rather, that management plans for such species should account for population genetic variation across the species’ ranges.


INTRODUCTION
The characterization of population genetic structure and dispersal patterns offers fundamental information on the ecology and evolution of species and provides vital knowledge for the conservation of threatened organisms (Cegelski et al. 2003;Walters and Schwartz 2020). Species that occupy similar geographic ranges may exhibit varying patterns of population structure determined by interactions of the species' abilities to move between habitat patches, the intervening matrix between suitable habitat patches, and the spatial distribution of such patches. Differences in landscape genetics will, in turn, vary with biological characteristics of species including body size, vagility, habitat requirements, and life history (Hamrick and Godt 1996;Waples 1998). For many species, the movement of individuals between patches is infrequent when suitable habitat is separated by, or when the suitable patches are positioned in, a resistant landscape matrix. Under such circumstances, high degrees of isolation cause local populations to become more susceptible to demographic and stochastic effects that tend to reduce genetic diversity, increase inbreeding, and eventually lead to local extirpation (Frankham 2005). At the same time, such population isolation can also lead to local adaptation to environmental conditions (Kawecki and Ebert 2004).
While dispersal ability may be an important factor determining linkages between populations, direct observations of such dispersal using telemetry or mark-and-recapture techniques are difficult for freshwater invertebrates (Hughes 2007). Molecular genetic tools provide an alternative method for inferring the relative influences of landscape features and dispersal capability on connectivity of populations. Gene flow, which reduces genetic differentiation between populations, tends to be positively correlated with dispersal ability (Slatkin 1987;Bohonak 1999); as such, quantifying genetic similarity between populations provides a "long-term" average of gene flow (Hughes 2007;Zickovich and Bohonak 2007). For example, highly vagile species tend to be genetically more homogenous throughout the landscape than more-sessile species (Avise et al. 1987;Bohonak 1999), allowing for the inference of dispersal ability based on genetic structure.
Comparative genetic studies are particularly powerful for determining the influence of landscape patterns on population genetic structure (Hickerson and Cunningham 2005). If species with different phylogenetic and life histories demonstrate similar population genetic structure, then one may infer that the landscape has played a central role in shaping the structure of the community. Sympatric taxa may have different responses to landscape fragmentation, which often reflects variability in lifehistory and ecological traits . For example, gene flow is reduced for sedentary species or habitat specialists, compared to species with greater dispersal capacities or broader ecological niches (Douglas et al. 2006;Callens et al. 2011). The influence of the landscape on the distribution of species can be determined by comparing the genetic structure of multiple species occupying the same geographic range with varying life histories.
In arid regions, aquatic ecosystems have been influenced by historic climate change and the patchiness of freshwater habitats embedded in a dry matrix. This often results in isolated populations separated by significant stretches of unsuitable terrestrial environments. Among these habitats are desert springs possessing communities dominated by non-vagile taxa (Stanislawczyk et al. 2018). Spring taxa are predicted to display high levels of genetic differentiation in arid regions (e.g., Marten et al. 2006). This pattern has been observed in multiple spring systems (Lefébure et al. 2007;Murphy et al. 2010) and it helps explain why highly localized taxa are common, with some species restricted to single springs (Ponder et al. 1989;Ponder and Colgan 2002;Witt et al. 2006;Murphy et al. 2009;Seidel et al. 2009;Adams et al. 2018). Such so-called micro-endemic species (Ochoa-Ochoa et al. 2011) have been described from desert springs throughout the world (Ponder and Colgan 2002;Witt et al. 2006;Murphy et al. 2013). Geographic isolation among desert springs may also influence demographic processes such as population bottlenecks or expansions, and gene flow (Harpending et al. 1998).
The Chihuahuan Desert of the southwestern United States and northern Mexico is one of the most biologically rich and diverse deserts in North America (Olson and Dinerstein 1998); its ecosystems harbor high numbers of endemics that are often evolutionarily unique and adapted to harsh environments (Metcalf and Smartt 1997;Kodric-Brown and Brown 2007). These high levels of biodiversity and endemism are especially evident in the desert's aquatic ecosystems (Abell et al. 2000). This region was part of the Western Interior Seaway of the mid-late Cretaceous and early Paleogene periods (~95 mya; Bousfield 1958), and the formation of artesian wells and multiple spring systems likely arose as a result of sea-level declines (Brune 1981). Investigations of phylogenetic structure of amphipods and hydrobiid gastropods suggest geographic isolation of springs has had a major influence on patterns of diversity within these taxa (Liu et al. 2003;Gervasio et al. 2004;Seidel et al. 2009;Adams et al. 2018;Walters et al. 2020).
In this study, we compared population genetic structure of four invertebrate species with overlapping ranges in the Chihuahuan Desert: the pulmonate gastropod Physa acuta, the prosobranch gastropods Juturnia kosteri and Pyrgulopsis roswellensis, and the amphipod Gammarus desperatus. Physa acuta is a globally widespread freshwater snail native to North America (Dillon et al. 2002;Lydeard et al. 2016) which now occurs on all continents except Antarctica (Wethington and Lydeard 2007;Bousset et al. 2014). The hydrobiid snails Juturnia kosteri and Pyrgulopsis roswellensis, as well as Gammarus desperatus, are micro-endemics restricted to springs and sinkholes located on Bitter Lake National Wildlife Refuge (BLNWR), Chaves County, New Mexico, USA. All three species are on the US Endangered Species List. Interestingly, previous research indicates that Gammarus on the refuge may represent multiple cryptic taxa (Adams et al. 2018); however, all populations are managed by the United States Fish and Wildlife Service as G. desperatus. Although all three species are locally abundant, their small geographic range provides the opportunity to investigate the response of these organisms to the landscape and therefore, estimate within-and among-population genetic variation across the entire geographic range of each species. The comparison of four sympatric taxa-three with narrow ranges compared to the more widespread, but less abundant, P. acutaprovides insight about landscape patterns that influence genetic structure at small spatial scales.
Our goal was to determine the relative impacts of landscape structure and dispersal ability on population genetic structure by using a comparative, multispecies approach to examine genetic variation for obligate, aquatic species in a highly managed system with natural and human-made barriers to dispersal. We used microsatellite DNA to (1) characterize genetic diversity within populations and genetic divergence among populations for each species and (2) contrast genetic structure among co-distributed species. These organisms represent the dominant macroinvertebrates at BLNWR (Stanislawczyk et al. 2018), but typically occupy different microhabitats throughout the refuge and likely possess different dispersal abilities. Previous genetic studies utilizing nuclear and mitochondrial DNA (mtDNA) gene sequences found evidence of strong differentiation among populations of hydrobiids (Morningstar et al. 2014) and gammarid amphipods (Seidel et al. 2009;Adams et al. 2018); conversely, physid snails showed little variation in mtDNA (K. Inoue, unpublished data). We employed genetic markers with higher rates of evolution to test the hypothesis that population genetic structure and gene flow would vary with life histories. Interestingly, little knowledge is available about dispersal of both freshwater amphipods and gastropods; however, passive dispersal may occur via birds (Ponder et al. 1989;Ponder et al. 1994;Rachalewski et al. 2013). Because it occurs throughout southeastern New Mexico, we predicted that the more widespread P. acuta would exhibit higher levels of genetic connectivity and less population differentiation compared to the three micro-endemic species. Furthermore, P. acuta is a pulmonate ("lunged") snail and can tolerate extended periods of drying (Collas et al. 2014), which may provide a mechanistic explanation for its widespread distribution. Alternatively, Gammarus desperatus and the hydrobiids are concentrated near spring sources, with density declining downstream (Noel 1954;Taylor 1987;Hershler 1994;Lang 2005). Therefore, we predicted that these three species would display low dispersal abilities and show significant population differentiation. Because the hydrobiid snails possess opercula that can decrease desiccation, we predict that they will be intermediate in dispersal ability and gene flow, while the amphipod will show the greatest degree of isolation. Multi-species comparisons of population genetic structure are valuable for understanding ecological and evolutionary processes that shape variation across landscapes. In the cases of imperiled species, it is essential to understand how processes such as gene flow and genetic drift affect diversity to inform conservation and management actions (Toro and Cabalero 2005;Allendorf et al. 2012). This is especially true in desert springs as these systems are threatened by groundwater withdrawal and reduction of spring connectivity (Rice and Emery 2003).

MATERIALS AND METHODS
BLNWR is located in the Pecos Valley of the northern Chihuahuan Desert near Roswell, Chaves County, New Mexico, USA. The 94.5 km 2 area is a primarily karst landscape of natural groundwater discharge at the downgradient end of the regional hydrological system known as the Roswell Artesian Basin (Land 2005;Land and Newton 2007;. This discharge, which amounts to roughly 37 million m 3 yr −1 , reaches the surface via springs, seeps, sinkhole lakes and extensive wetlands formed in gypsum bedrock on the west side of the Pecos River floodplain. These aquatic habitats contain many range-restricted taxa, including several endemic macroinvertebrates that are listed as endangered by the US Fish and Wildlife Service (Cole 1981; U.S. Fish and Wildlife Service 2005). The landscape of BLNWR has been highly modified by human activity. The refuge was established in 1937 as wintering and breeding grounds for migratory birds (U.S. Fish and Wildlife Service 1998) with management directed primarily at creating dikes and permanent ponds. Since 1994, management for waterfowl has focused on providing seasonal wetlands through moist soil management techniques which disconnect springs from seasonally flooded areas (U.S. Fish and Wildlife Service 1998). The water management on the refuge has resulted in fragmentation of aquatic ecosystems.
From 2014 to 2016, we sampled 16 sites throughout the refuge for all four species, attempting to obtain 24 individuals of each species per site; however, presence and abundance varied among species and across sites.
The sampled area includes the entire ranges of the hydrobiid snails and amphipods (Fig. 1). Each species was sampled throughout the refuge within a single sampling effort (i.e., all Gammarus and Physa specimens were collected in 2015, all hydrobiids in 2014) with the exception of Unit 5 sites (U5A, U5B, U5C) which were sampled in 2016. We sorted individuals while in the field; P. roswellensis and J. kosteri, were distinguished by operculum color using a hand lens or dissecting microscope. All specimens were stored in 95% ethanol until DNA was extracted using DNeasy kits Fig. 1 Map of Bitter Lake National Wildlife Refuge, Chaves County, New Mexico, USA. Points denote collection sites for Physa acuta, Juturnia kosteri, Pyrgulopsis roswellensis, and Gammarus desperatus with black lines indicating management units. (Qiagen, Inc.). We used published microsatellite primers and methods to genotype 11 loci for Physa (Ansah et al. 2014), 16 loci for J. kosteri and 20 loci for P. roswellensis (both described in Holste et al. 2016) to assess genetic diversity of these species.
For G. desperatus, we developed novel microsatellite primers using the Plant-Microbe Genomics Facility at Ohio State University to perform de novo pyrosequencing in 1/4th of a picotiterplate using a Roche 454 FLX Titanium Genome Sequencer which produced 55,692 reads. From sequencing outputs, we used MSATCOMMANDER v1.0.8 (Faircloth 2008) to identify putative microsatellite loci and PRIMER3 (Rozen and Skaletsky 2000) to design primers on flanking regions. Tail sequences of M13R or CAG were added to the 5'-end of either the forward or reverse primer. We identified 202 putative microsatellite markers and we were able to design primers for 120 microsatellite loci. Trinucleotides were the most abundant repeat type (46.7% of all loci), followed by dinucleotides (42.5%), and tetranucleotides (10.8%). Twenty loci (Table S1, Supplementary Information) were consistently amplified and used for population-level analyses. These loci were amplified in 10 µL polymerase chain reactions (PCR) using Gotaq Master Mix (Promega Corporation), 0.2 µM of universal fluorescently labeled primer (M13R or CAG) including the non-tailed primer, 0.04 µM of tailed primer, and 10 ng of DNA. For each locus, we identified optimum annealing temperature using the following PCR conditions: initial denaturing at 95°C for 2 min; followed by 35 cycles at 95°C for 30 s, annealing at 48-62°C (or 55-65°C as a higher temperature gradient, if needed) for 45 s, extension at 72°C for 1 min; and final extension at 72°C for 30 min.
For all species examined, fragment analyses were performed on an ABI Genetic Analyzer with LIZ600 size standard (Applied Biosystems, Inc.). We used PEAKSCANNER v1.0 (Applied Biosystems, Inc.) to score alleles and TANDEM v1.07 (Matschiner and Salzburger 2009) to assign integers to DNA fragment sizes. We used MICRO-CHECKER (van Oosterhout et al. 2004) to detect null alleles and estimate allele frequencies, and GENEPOP (Rousset 2008) to conduct exact tests of pairwise linkage disequilibrium and deviation from Hardy-Weinberg expectation (HWE) for each combination of locus and locality. Sequential-comparison Bonferroni correction (s-cB; Rice 1989) was applied to adjust p values to account for multiple simultaneous tests. For polymorphic loci, we used GENALEX v6.5 (Peakall and Smouse 2006) to estimate allelic richness (number of alleles per locus; N A ), observed and expected heterozygosities (H O and H E ), and mean number of private alleles (N P ) for each locus-by-population combination. Allelic richness (A R ) was estimated using a rarefaction procedure in the program FSTAT v.2.9.3 (Goudet 1995) to account for unequal numbers of populations and samples sizes (Kalinowski 2004).
To assess population genetic structure, we estimated pairwise genetic differences among populations. Genetic divergence among sites for each species was assessed by calculating pairwise F ST and R ST . While F ST is influenced by the level of heterozygosity, R ST takes into account differences in the number of repeats between microsatellite alleles (allele size). To choose the most suitable estimator, we implemented a randomization procedure of allele sizes in SPAGeDi v1.5 (Hardy and Vekemans 2002) to determine whether stepwise-like mutations contributed to genetic differentiation . We compared the observed R ST values with the distribution of R ST obtained with 10,000 allele size permutations (pR ST ). R ST would be expected to be significantly higher than pR ST when the migration rate is lower than the mutation rate, suggesting that populations had been isolated for a long time and/or populations exchanged migrants at a rate similar to or less than the mutation rate . Additionally, FreeNA was used to assess the degree of genetic differentiation among sites by calculating pairwise F ST values, employing the ENA (excluding null alleles) correction, which has been shown to effectively correct for the positive bias of F ST that may result from the presence of null alleles (Chapuis and Estoup 2007). To assess the effect of geographic distance on genetic structure, we examined the correlation of pairwise genetic distances and geographic distances between sites using Mantel tests. Surface water connections among a majority of sites are either lacking or ephemeral in nature; therefore, we measured Euclidean (overland) distances between pairs of sites using ARCGIS v10.2 (ESRI, Inc). Mantel tests were performed in GenAlEx with 9999 permutations. The use of Mantel tests for spatial analyses in ecological studies has been the subject of some controversy in the literature because statistical weaknesses of the test can lead to false positives (Legendre and Fortin 2010;Guillot and Rousset 2013;Legendre et al. 2015). However, it has been suggested that the use of Mantel tests for analysis of isolation-by-distance can be useful especially when results are interpreted conservatively (Diniz-Filho et al. 2013). Additionally, genetic variation in our system is not hierarchically clustered (as in a phylogeny) so that the assumptions of the Mantel test are not violated. A Principal Coordinate Analysis (PCoA) for each species was conducted in GenAlEx using a covariance matrix created from all polymorphic loci to visualize population genetic structure.
BAYESASS 3.0.4 (Wilson and Rannala 2003) was used to infer recent gene flow (over the last several generations) between pairs of sampling localities via a Bayesian MCMC approach that relies on multilocus gametic disequilibrium information to estimate asymmetrical migration rates (and their respective 95% credible sets). BAYESASS assumes low levels of migration (not to exceed one-third of the population) and does not assume that populations are in Hardy-Weinberg equilibrium. Each run was performed with 10 million generations, sampling every 1000 generations after a burn-in period of 10 million generations. Delta values (i.e., maximum parameter change per iteration) were adjusted to achieve optimal acceptance (i.e., 20-40%; Wilson and Rannala 2003). Average migration rates and 95% credible intervals (calculated as migration rate ± 1.96 standard deviation) were estimated.

RESULTS
Multi-locus genotypes were obtained from a total of 678 individuals across the four species. Across the 16 sampling locations on the refuge, each species was present at 6-12 sites.
Within-population diversity For P. acuta, a total of 123 individuals from seven populations were successfully genotyped for 11 polymorphic loci, with a range of two (locus Tetra3) to 26 (Tri3) alleles per locus. We found considerable evidence of null alleles (50.65% of all locality-by-locus pairs; Table  S2, Supplementary Information) and no evidence of linkage disequilibrium. A substantial proportion of locus-by-population combinations deviated from HWE, with 43 of 77 tests showing significant deviations after s-cB correction (Table S2). Of these, 29 showed heterozygote deficits. Some departures from HWE may be the result of the presence of null alleles (Table S2, Supplementary  Information); however, the consistent heterozygote deficit over localities and loci is likely explained by a high selfing rate for this species (Monsutti and Perrin 1999). Therefore, all loci were included in subsequent analyses. Average number of alleles scored per locus across all loci ranged from 3.09 at site SS to 6.09 at RH. Overall mean rarefied allelic richness standardized for N = 7 varied from 2.62 (BC-DS) to 4.63 (S-31). Average number of private alleles varied from 0.18 (SS, HM, and BC-LRC) to 0.91 (LFS and RH; Table 1). The highest and lowest mean levels of heterozygosity (H O = 0.54 and 0.35) were found at RH and S-31, respectively (Table 1).
For J. kosteri, a total of 269 individuals from 12 populations were genotyped for 16 polymorphic loci. We found evidence of null alleles (34.38% of all locality-by-locus pairs; Table S3, Supplementary Information), no evidence of linkage disequilibrium, and significant deviation from HWE for 24 of 192 locality-by-locus pairs after s-cB correction (Table S3). Of these, 21 showed a heterozygote deficit. Some departures from HWE appeared to be the result of the presence of null alleles (Table S3). Locus JKTetra1 deviated from HWE in eight of the 12 populations examined, with a high frequency of null alleles in nine populations; therefore, this locus was removed from subsequent analyses. All remaining loci were retained for analyses because violations of HWE assumptions were inconsistently distributed among loci and populations. We scored a total of 145 alleles across 15 loci, ranging from three (locus JKTetra12) to 19 (JKTri5) alleles at a locus. The mean allelic richness across all 15 loci varied from 2.73 at site SS to 5.20 at HM (Table 1). Overall mean rarefied allelic richness standardized for N = 7 ranged from 2.44 (SS) to 3.47 (HM). Average number of private alleles varied from 0 (LSF and USU) to 0.80 (HM; Table 1). The highest and lowest mean levels of heterozygosity (H O = 0.47 and 0.28) were found at sites HM and BC-DS, respectively (Table 1).
For P. roswellensis, a total of 135 individuals from six populations were genotyped for 20 polymorphic loci. We found evidence of null alleles (32.5% of all locality-by-locus pairs; Table S4, Supplementary Information), no evidence of linkage disequilibrium, and significant deviation from HWE for 31 of 120 localityby-locus pairs after s-cB correction (Table S4). Of these, 19 showed heterozygote deficits. Some departures from HWE appeared to be the result of the presence of null alleles (Table S4). Locus PRTri18 deviated from HWE across all populations, and therefore was removed from subsequent analyses. All remaining loci were retained for analyses because violations of assumptions were inconsistently distributed among loci and populations. We scored a total of 266 alleles across 19 loci, ranging from four (locus PRTetra04) to 27 (PRTetra10) alleles at a locus. The mean allelic richness across all 19 loci ranged from 3.96 (site USU) to 8.47 (SS-L; Table 1). Overall mean rarefied allelic richness standardized for N = 16 ranged from 4.06 (U5A) to 7.32 (SS-L). Average number of private alleles varied from 0.16 (SS) to 1.79 (SS-L; Table 1). The highest mean level of heterozygosity was found at U6 (H O = 0.66) and the lowest at SS-L and U5A (H O = 0.49; Table 1).
For G. desperatus, a total of 151 individuals from seven populations were genotyped for 20 polymorphic loci. Over these loci, we found evidence of null alleles (40.0% of all locality-bylocus pairs; Table S5, Supplementary Information), no evidence of linkage disequilibrium, and significant deviation from HWE for 56 of 140 locality-by-locus pairs after s-cB correction (Table S5). Of these, 37 showed a heterozygote deficit. Some departures from HWE appeared to be the result of the presence of null alleles (Table S5). The following loci were removed from subsequent analyses because of deviation from HWE and/ or high frequency of null alleles: GD11, GD13, GD14, GD16, GD17. The remaining 15 loci were retained for analyses because violations of assumptions were inconsistently distributed among loci and populations. We scored a total of 180 alleles across these loci, ranging from three (site GD20) to 20 (GD22) alleles at a locus. The mean allelic richness across all 15 loci ranged from 4.2 (S-31) to 6.53 (U7). Overall mean rarefied allelic richness standardized for N = 14 ranged from 3.65 (U6) to 5.67 (U7;

Among-population divergence
For P. acuta, all pairwise comparisons of F ST were significantly greater than zero but the range of values was small (0.04-0.12; Table 2). ENA-corrected estimates of F ST were slightly higher, ranging from 0.03 between sites BC-LRC and BC-DS to 0.23 between sites SS and HM (Table 2). Pairwise comparisons of R ST ranged from 0.01 between sites BC-LRC and RH to 0.76 between sites SS and HM. R ST was significantly greater than pR ST (P = 0.003) suggesting that F ST is likely to provide a biased estimate of gene flow parameters; therefore, a Mantel test was performed on R ST values. There was no significant pattern of isolation-by-distance (Mantel r = 0.0.02, P = 0.27, Fig. 2). The PCoA recovered a single cluster on the refuge, with overlap between all populations (Fig. 3a). The first and second principal coordinate axes explained 8.64% and 7.42% of the genetic variation, respectively (Fig. 3a). The estimated proportion of migrants was low (m < 0.07) among a majority of sites with no current migration inferred between most population pairs because migration rates included 0 in the 95% credible interval (Table 3). Asymmetric migration from BC-DS to directly downstream BC-LRC (m = 0.21, 95% CI: 0.11-0.31), as well as further downstream S-31 (m = 0.23, 95% CI: 0.15-0.31), HM (m = 0.21, 95% CI: 0.13-0.29), and RH (m = 0.11, 95% CI: 0.03-0.19), was evident, suggesting that BC-DS was a net producer of migrants entering southern sites (Table 3).
For J. kosteri, all pairwise comparisons of F ST were small (F ST = 0.03-0.14; Table 4), but significantly greater than zero. ENAcorrected estimates of F ST were higher, ranging from 0.03 between sites BC-F and BC-LRC to 0.24 between sites SS and U5A (Table 4). Pairwise comparisons of R ST ranged from 0.03 between sites LSF and BC-F to 0.45 between sites BC-F and U5A (Table 4). Again, R ST was significantly greater than pR ST (P = 0.000); a Mantel test of R ST values failed to show isolation-by-distance (Mantel r = 0.02, P = 0.22; Fig. 2).The same was true when U5 sites were removed from the Mantel test to control for temporal variation (Mantel r = −0.06, P = 0.11; Fig S1, Supplementary Information). The PCoA did not identify discrete clusters on the refuge, with overlap between all populations (Fig. 3b). The first and second principal coordinate axes explained 22.12 and 17.88% of the genetic variation, respectively (Fig. 3b). The first and second principal coordinate axes explained 7.10 and 13.94% when U5 sites were removed (Fig S2, Supplementary Information). Estimates of migration rates were low (m = 0.03) among a majority of sites, with no current migration inferred between most population pairs (Table 5). For sites with m > 0.03, asymmetric migration was observed as a north-to-south pattern of migration (with one exception), possibly indicating gene flow overland or through groundwater connections (  (Table 6). Pairwise F ST values ranged from 0.03 between sites U7 and USU to 0.21 between SS and U5A. ENA-corrected estimates of F ST were higher, ranging from 0.03 between sites USU and U7 to 0.32 between sites SS and U5A (Table 6). Pairwise comparisons of R ST ranged from 0.02 between sites USU and U7 to 0.42 between sites SS and U5A (Table 6). R ST was significantly greater than pR ST (P = 0.03). A strong pattern of isolation-by-distance was found for P. roswellensis (R ST Mantel r = 0.60; P = 0.003; Fig. 2). Such a pattern was even stronger when U5 sites were removed (Mantel r = 0.84 P = 0.03, Fig S3, Supplementary Information). The PCoA revealed three discrete clusters: Sago Spring sites (SS and SS-L), U5A, and the southernmost sites each formed a cluster (Fig. 3c). The first and second principal coordinate axes explained 37.28% and 24.86% of the genetic variation, respectively (Fig. 3c). A similar pattern was observed when U5 sites were removed; the first and second principal coordinate axes explained 17.80% and 25.17%, respectively ( Fig  S4, Supplementary Information). Estimates of migration rates were low (m = 0.04) among a majority of sites, with no current migration inferred between population pairs, except from site SS into the adjacent SS-L site (m = 0.18, 95% CI: 0.12-0.24) and U7 into the adjacent USU site (m = 0.23, 95% CI: 0.17-0.29; Table 7), consistent with downstream migration. Fig. 2 Empirical relationships between genetic distances and geographic distances between pairs of populations for Physa acuta, Juturnia kosteri, Pyrgulopsis roswellensis, and Gammarus. Genetic distance is represented by R ST for P. acuta, J. kosteri, and P. roswellensis and F ST for Gammarus. Significant (P < 0.05) regression lines provided. For G. desperatus, all pairwise comparisons of F ST were significantly greater than zero (Table 8). Pairwise F ST values ranged from 0.02 between sites BC-LRC and SS to 0.25 between S-31 and RH, U6 and USU, and U6 and RH. ENA-corrected estimates of F ST were higher, ranging from 0.02 between sites BC-LRC and SS to 0.32 between U6 and RH (Table 8). Pairwise comparisons of R ST ranged from 0.02 between sites BC-LRC and SS to 0.52 between sites SS and USU (Table 8). R ST was not significantly greater than pR ST (P = 0.20) suggesting that allele size is not informative for population differentiation; therefore, a Mantel test was performed on F ST . A strong pattern of isolation-by-distance was observed for Gammarus on the refuge (F ST, Mantel r = 0.61, P = 0.003; Fig. 2) when all sites were included. Because of previous research indicating potentially cryptic species of Gammarus in the southern portion of the refuge (Adams et al. 2018), we performed a PCoA without USU and RH sites. When these sites were removed, there was no significant pattern of isolation-by-distance (F ST , Mantel r = 0.65, P = 0.06; Fig. S1, Supplementary Information). The PCoA revealed three discrete clusters: the northernmost sites, USU, and RH (Fig. 3). The first and second principal coordinate axes explained 21.60% and 8.49% of the genetic variation, respectively (Fig. 3). When USU and RH were removed, the PCoA revealed three clusters: the northern sites (BC-LRC, SS, and S-31), U6, and U7, with some overlap between the clusters. The first and second principal coordinate axes explained 10.88% and 7.94% of the genetic variation, respectively ( Fig. S2; Supplementary Information). Estimates of migration rates were low     (m < 0.02) among sites, with no current migration between most population pairs (Table 9). However, the BC-LRC site was a net producer of migrants into SS (m = 0.26, 95% CI: 0.20-0.32) and S-31 (m = 0.23, 95% CI: 0.17-0.29; Table 9).

DISCUSSION
We tested the hypothesis that differences in inferred dispersal ability of four co-distributed species, based on size of range and morphological traits allowing survival in dry environments, were    Each value represents the fraction of individuals in a row population that is derived from a column population (per generation). Values in parentheses are ±95% confidence interval.
A.D. Walters et al. correlated with variation in genetic structure at a small spatial scale across a highly managed groundwater system. Our multi-species approach revealed that genetic patterns and levels of connectivity varied among these species, including those for which this region represents the entire geographic range, suggesting that the relevant spatial scale for genetic management is speciesdependent even within small areas (~100 km 2 ). While quantitative information regarding dispersal ability is lacking for most invertebrates (Benton and Bowler 2012), genetic approaches can be useful for investigating dispersal among populations (Cayuela et al. 2018). As predicted, the micro-endemic P. roswellensis and Gammarus exhibited patterns of isolation-by-distance; however, this pattern was not observed in the micro-endemic J. kosteri. Additionally, the pattern of isolation-by-distance pattern was not present in the more-widespread P. acuta, consistent with the inference that this species has a higher dispersal ability possibly due to the ability to exchange atmospheric gases using the pallial lung. The fact that the isolation-by-distance relationships for the micro-endemic P. roswellensis and Gammarus are not asymptotic suggests that dispersal limitation per se is not present at this spatial scale (Gómez-Rodriguez et al. 2020), but rather that habitats outside the boundaries of the refuge are unsuitable creating a hard "edge" for the ranges of the species. Whether this is due to the refuge boundaries coinciding with suitable aquatic habitat or to alteration of habitat with human development outside of the refuge is not clear. We expected that the hydrobiids would have similar population genetic structure, however P. roswellensis showed distinct population clusters on the refuge that were not seen with J. kosteri. While we hypothesized that Gammarus would show the greatest degree of isolation, once we accounted for the likely presence of multiple species, G. desperatus showed reduced genetic structure across the refuge. Overall, the observed differences in genetic patterns can likely be attributed to species-specific combinations of dispersal ability and habitat preference. The genetic diversity observed in P. acuta was higher than previously documented for a single site on the refuge (Ansah et al. 2014), which is likely a result of increased sampling effort in the current study. Interestingly, P. acuta on the refuge is characterized by the absence of isolation-by-distance along with admixture among spring sites; however, R ST values suggest that P. acuta populations are structured throughout the refuge. This pattern corresponds to an island model in which dispersal is limited between the sampled populations, even if not limited by distance. This pattern is consistent with the result that current migration went undetected between most P. acuta populations. Similar with our results, a study examining the population genetic structure of P. acuta failed to detect a pattern of isolation-by-distance at a larger geographic scale (Bousset et al. 2014). Additionally, all P. acuta populations were characterized by heterozygote deficits, which are consistent with high rates of self-fertilization. Many species within the genus are functionally self-fertile hermaphrodites, but preferentially cross-fertilize if potential mates are available (reviewed by Duncan 1975;Jarne and Charlesworth 1993). While P. acuta is distributed throughout the refuge, this species occurs at very low densities (Stanislawczyk et al. 2018), potentially reducing mate availability and thus, increasing the likelihood of self-fertilization.
While J. kosteri displayed a similar pattern to P. acuta, P. roswellensis showed a pattern of isolation-by-distance. This species also displayed discrete populations on the refuge. While information is lacking for direct measures of dispersal abilities of the hydrobiids, the strong positive correlation between geographic distance and genetic distance for P. roswellensis indicates that dispersal distances are restricted or that populations recently invaded the region where the isolation-by-distance pattern arose in nonequilibrium conditions (Puechmaille et al. 2011). The existence of distinct genetic groups in P. roswellensis could indicate either discrete populations or incomplete sampling across the refuge with data missing from intermediate populations.
Limited dispersal leads to geographic isolation which, when combined with small population size, can lead to increased risk of extinction due to inbreeding depression even when environmental conditions are suitable (Spielman et al. 2004). Population genetic and metapopulation theories suggest that isolated populations should have lower effective population sizes and less genetic variation (Slatkin 1987), although P. roswellensis has the highest levels of genetic diversity of the species examined in this study. This species may persist as a result of inherently high levels of genetic diversity, as is the case with other dispersal-limited taxa (Witt et al. 2008;Guzik et al. 2009).
It is important to note that while a majority of individuals were collected during the 2014 sampling effort, the U5 sites were collected in 2016 and this potentially adds temporal genetic variation for both J. kosteri and P. roswellensis. While the U5 sites are in the central part of the refuge, they provide the highest pairwise F ST and R ST values for both species, with residuals for U5 sites positive in regression plots. The temporal component of genetic structure is also visible in the PCoA analyses. The temporal pattern may be confounded with the spatial pattern of population genetic structure for these species in the middle section of the refuge; however, the relationship between geographic distance and genetic distance and clustering patterns were similar when temporal variation was negated by removing the U5 sites from analyses.
Our analyses suggest that Gammarus on the refuge consists of multiple populations, with the two southernmost sites (USU and RH) distinct from the northern springs. Similar patterns were shown for nuclear and mitochondrial gene sequences with divergence between sites on the refuge and the RH site, which has been identified as an undescribed species that likely diverged from G. desperatus~5.8 million years ago (Adams et al. 2018). While the USU population was not included in Adams et al. (2018), our microsatellite data suggest that it, too, is a distinct population. There is high genetic differentiation and extremely low gene flow between these southern populations (USU and RH) and populations in the northern portion of the refuge as seen by high F ST (0.20-0.25), low admixture amongst PCoA genetic clusters, fixed allelic differences and high frequencies of private alleles in the southern populations. The Rio Hondo itself is disconnected from surface water flows linking other sites on BLNWR. Historical surface water flow may have facilitated some seasonal gene flow among sites; however, extensive diversion of groundwater in the Pecos River drainage for irrigation has altered surface water flows and hence, the connectivity of this system.
The estimates of migration suggest that dispersal occurs in a north-to-south direction on the refuge, and that migration can occur between sites lacking surface flow connections. If so, migration may occur through groundwater connections. Lake Saint Francis (LSF) is a sinkhole with no surface water connections; however, it was a net producer of J. kosteri migrants for sites in the northern portion of the refuge. Additionally, the Bitter Creek site (BC-LRC) was a net producer of G. desperatus migrants for Sago Spring sites (SS and S-31). While the groundwater sources for Bitter Creek and Sago Spring seem to be distinct and independent (Lang 1998), the overland distance between these locations is less than 0.5 km, with a river distance between the two sites just over 2 km. While little information is available regarding dispersal of these species, passive dispersal may occur via birds for both snails (Ponder et al. 1989(Ponder et al. , 1994 and amphipods (Rachalewski et al. 2013), while groundwater connections have been shown to facilitate dispersal in amphipods (Harris et al. 2002). Estimates of migration rates based on allele frequencies assume populations are in long-term gene flow-drift equilibrium (Bohonak and Jenkins 2003). Although BAYESASS is relatively free from assumptions, this method may be sensitive to the presence of null alleles, low population structure, and limited sample size (Wilson and Rannala 2003). Additionally, our results should be treated with caution because the dynamic nature of springs and frequency of local extinction and recolonization events can lead to nonequilibrium conditions over short time periods (Ponder and Colgan 2002;Worthington-Wilmer et al. 2011); however, temporal monitoring across the refuge does not suggest recent local extinction events (Johnson et al. 2019). Nevertheless, our results support migration across the refuge that does not appear to be constrained to surface water connections; therefore, future research should include efforts at delineating both surface and subsurface flows and discharges for the refuge.
The distribution of habitat patches scattered throughout a hostile terrestrial matrix in arid landscapes provides ideal conditions for highly structured populations of non-vagile, obligate aquatic species (Ribera and Vogler 2000;Burridge et al. 2008;Hershler and Liu 2008). In arid spring systems, natural fragmentation has been shown to promote speciation (Murphy et al. 2015;Adams et al. 2018) and influence population structure (Murphy et al. 2010). Additionally, traits such as mobility can explain disparities in patterns of structure if lower dispersal capability results in lower gene flow and higher population divergence (Henle et al. 2004). While the water management activities on the refuge, including flooding and fragmentation, have likely influenced contemporary dispersal and gene flow among populations, the results of our study examining spatial patterns of co-distributed taxa suggest that population structure is highly dependent on characteristics of each species. Similarly, genetic patterns that vary with dispersal abilities have been observed for sympatric salmonids (Gomez-Uchida et al. 2009), aquatic insects (Phillipsen and Lytle 2012;Phillipsen et al. 2015) and amphibians , suggesting that speciesspecific dispersal characteristics can lead to major differences in population structure. Comparisons from multi-species studies is an important first step in understanding how landscape features and variation in traits interact to drive gene flow and population genetic structure in fragmented landscapes.
The combination of limited dispersal and a highly altered flow regime on the refuge in a naturally fragmented landscape has influenced the connectivity of populations. The micro-endemics consist of multiple populations, distributed over relatively limited geographic distances. Significant genetic structure has been observed in other micro-endemics (Ponder and Colgan 2002;Murphy et al. 2010). Further fragmentation of aquatic habitats could increase the effects of dispersal limitation with negative consequences for imperiled taxa, because geographic isolation can lead to increased risk of extinction (Spielman et al. 2004). Our multispecies comparison of population structure provides valuable insight into the landscape ecology of micro-endemic aquatic species in this region and other deserts around the world.

Conservation implications
Conservation of endangered species is dependent on knowledge of the genetic structure and diversity of individual populations (Frankham 2003). Genetic diversity is often spatially structured because suitable habitat is not always continuous or may vary in quality across a species' range. Additionally, differences in dispersal ability and environmental tolerances may lead to differences among co-distributed taxa. Anthropogenic activities, resulting in fragmentation and habitat loss, also play a role in shaping the genetic structure of species. At the population level, fragmentation can reduce suitable habitat, restricting gene flow and promoting genetic drift (Frankham 2005). Our analyses show incongruent patterns of population structure and varying degrees of genetic diversity among four invertebrate species, three of which are narrowly endemic and of conservation concern. Information regarding environmental and genetic features aids the development of species recovery strategies including the designation of management units. The management units on BLNWR are currently identified based on the characterization of threats, without regard for biological population groupings (U.S. Fish and Wildlife Service 2019). Our results provide vital information that identifies population genetic structure for these species. Incorporation of spatial processes in species recovery plans is critical for decisions regarding which habitats, and the quantities of these habitats, to restore (Huxel and Hastings 1999) so as to mitigate the effects of habitat loss and fragmentation (Lewis et al. 1996). Maintaining gene flow among populations and/or preserving the spatial component of genetic diversity is fundamental to ensuring the viability and the evolutionary potential of species in the face of anthropogenic challenges including ongoing climate change (Frankham et al. 2002;Lemer and Planes 2014). Habitat alteration in these spring systems through groundwater extraction and pollution is perhaps the greatest threat to these species and has resulted in the extirpation and extinction of endemics in desert spring systems throughout the world (Ponder et al. 1995;Sada and Vinyard 2002;Witt et al. 2006). Thus, even at spatial scales which humans perceive as quite small, conservation planning should account for multiple populations in fragmented habitats because such populations are likely to contain unique diversity within them.

Data archiving
The data have been deposited on dryad with the unique digital object identifier https://doi.org/10.5061/dryad.6hdr7sr2w.