Geology correlates with gut microbial community composition in the Mountainsnails (Oreohelicidae: Oreohelix)

Sampling

We collected a total of 69 snail samples and 27 soil samples in July 2020, when snails were known to be active at the selected sampling locations. We selected sites roughly between Lewiston and Riggins, ID based on geological characteristics, using geologic maps and previous studies (Kauffman et al., 2014; Linscott, Recla & Parent, 2023) to characterize sites as either being on or off calcareous rock (i.e., CaCO₃ rock). Of the sites we chose, five were on and four were off CaCO₃ rock (Fig. 1).

At each site, we collected three surface soil samples and a minimum of three and up to ten live snails (mean number of samples per site = 7.66) opportunistically and placed each individual into sterile vials. We transported the samples to the University of Idaho where soil samples were stored at −80 °C and snails at room temperature. The snail species and subspecies sampled are morphologically and geographically distinct, and their identities were confirmed in prior work using the same populations (Linscott, Recla & Parent, 2023). Snails in our study, except for O. jugalis (site LUC01), include those in the recently diverged O. strigosa species complex (Linscott et al., 2020). A pruned tree based on Linscott et al. (2020) is presented in Fig. 1 to show the phylogenetic relationships among species/subspecies in the present study. None of the sampled populations occur in sympatry save O. jugalis and O. idahoensis which can sometimes be found in sympatry at the geologic boundary between volcanic and limestone populations.

Sample processing

We sacrificed live snails within two days of collection, removed them from their shells, and gut tissue (intestine from digestive gland to anal pore) was dissected, rinsed with sterile water, and stored in 95% ethanol at −20 °C until extraction. Our extractions of snail (whole gut tissue) and soil (approximately 100 mg of bulk soil per sample) DNA were conducted using the DNeasy Powersoil Pro Kit (Qiagen) using the manufacturer’s instructions. We eluted extractions into 50 µl of C6 solution and stored them at −20 °C until library preparation.

Library preparation and sequencing

We amplified the V4 hypervariable region of the 16S rRNA gene and barcoded it using the 515F/806R primer pair and protocol from the Earth Microbiome Protocol for 2-step PCR (Thompson et al., 2017). Positive and negative controls were used at each step and progress was checked using a Qubit Fluorometer (ThermoFisher Scientific) and gel electrophoresis. We utilized qubit scores and gel band intensity to ensure equal amounts of DNA were included in the libraries. All the libraries were sequenced on a 300 bp paired-end Illumina Miseq platform PE300 (Illumina Corporation, San Diego, CA, USA) lane by the Genomics and Bioinformatics Resources Core at University of Idaho.

Data analyses

We processed reads for further analyses using QIIME 2 2024.10 (Bolyen et al., 2019). We used the q2-demux plugin to quality filter the raw sequences, followed by denoising with DADA2 (via q2-dada2; Callahan et al., 2016). We trimmed the sequence length during DADA2 based on visual inspection of sequence quality using Phred scores while maintaining overlap between forward and reverse sequences. These thresholds were chosen as 230 bp for forward reads and 130 bp for reverse reads. Using extraction and PCR blanks, contaminants were then identified by prevalence using the Qiime2 plugins for DECONTAM (decontam-identify and decontam-remove) (Davis et al., 2018). We aligned the amplicon sequence variants (ASVs) using MAFFT (via q2-alignment; Katoh et al., 2002) and used the aligned sequences to construct a phylogeny with fasttree2 (via q2-phylogeny; Price, Dehal & Arkin, 2010). Samples below 5,000 reads were dropped (based on rarefaction curves; see Fig. S1) as well as ASVs with less than 10 reads in any sample to remove possible contaminants. This left 93 total samples (68 snails, 25 soils) with 14,184 ASVs and 3,075,768 total reads.

For alpha-diversity metrics (observed ASVs, Faith’s Phylogenetic Diversity (Faith, 1992), and Shannon’s diversity index (Shannon, 1948)), beta diversity metrics, Robust Aitchison distance (Martino et al., 2019) (using QIIME2 version 2023.5), and Principal Coordinate Analysis (PCoA) we used the plugins deicode and q2-diversity. To assign taxonomy to ASVs, we used the q2-feature-classifier (Bokulich et al., 2018) classify-sklearn using the pre-trained classifier Silva 138 99% OTU database (Robeson et al., 2021; Quast et al., 2013) trained on the 515F/806R region of the 16S gene. We used PERMDISP and Adonis (PERMANOVA) (Anderson, 2001; Oksanen et al., 2023) to compute beta-diversity differences between samples using q2-diversity beta-group-significance and q2-diversity adonis. We included presence of CaCO₃-rich substrates at the sampling location (Geology) and whether samples were soil or snail derived (Sample Type) in our assessment of overall variation in beta diversity. Furthermore, within soil samples we again assessed the variation in beta diversity for geology. For snail samples, we assessed beta diversity variation for geology and snail taxa. We also assessed the beta diversity variation explained by geology within two localities (PITS and RR) that had a single species and sites on and off limestone to mitigate the effect of host taxonomy.

We migrated all output files to R using the package qiime2R (Bisanz, 2018). For soil and snail sample types, we used Mantel and partial Mantel tests with the vegan package in R (Oksanen et al., 2023) to test whether there was a correlation between microbiome composition using Robust Aitchison distance, geographic distance between samples, and/or presence of calcium-rich rocks (see Table 1 for details). Differential abundance analyses were done between sample types (snail vs soil) and between samples on/off CaCO₃ within sample types, using the full dataset as well as the subset of only RR and PITS sites. We performed analysis of compositions of microbiomes with bias correction (ANCOM-BC; Lin & Peddada, 2020) at various taxonomic levels (family level reported).

Effect size and power

We calculated effect size and power for the alpha diversity comparison between soil and snail samples and between snails on and off calcium-rich substrates. We calculated Cohen’s d for effect size using the command cohen.d in the effsize package (Torchiano, 2020) in R. Cohen’s d was used to calculate power using the pwr.t2n.test command in the pwr package (Champely, 2020).

Compositional differences between sample types

Our results showed that snail gut microbial compositions were significantly different from soil samples using Adonis (Oksanen et al., 2023) on Robust Aitchison distances (see Table 1). Using PERMdisp, we showed that there is also a difference in dispersion between sample types (F = 4.91, p = 0.028, permutations = 999). Figure 2 shows the separation of the two sample types, and the increased dispersion of snail samples compared to soil samples.

Compositional differences within sample types

Within sample types, we found that soil sample differences were significantly correlated with geographic distance but not significantly affected by CaCO₃ rock presence (Table 1). In snails, the reverse was true. There was no effect of geography (defined as spatial distance between sites), but CaCO₃ rock presence was found to be a significant factor for distinguishing microbiomes (Table 1). We also found that, when taking geographic distance into account using partial Mantel tests, the effect of CaCO₃ rock on microbiome dissimilarity remained significant (Table 1). Additionally, species identity was found to be a significant factor for snail microbiome differentiation, and geology (based on presence of CaCO₃-rich rock) remained a significant effect when accounting for species effect (Table 1). To further account for the effect of taxonomy, we assessed the compositional differences within a subset of microbiome samples from the same snail species. Namely, two locations, PITS containing O. strigosa and RR containing O. intersum, have the same taxa at sites on and off limestone. Using only these, we found that the effect of being on/off limestone is still significantly correlated with difference in microbiome composition (see Table 1).

Alpha diversity

Effect size and the power to detect difference was high for Shannon’s diversity between sample types but relatively low for all other alpha diversity measures (see all values in Table S1). With this in mind, we found no significant differences in observed richness (Kruskal–Wallace, H = 2.47, p = 0.12) and Faith’s phylogenetic diversity (Kruskal–Wallace, H = 0.54, p = 0.46) between soil and snail sample types. However, Shannon’s index showed significantly higher diversity in soil compared to the snail samples (Kruskal–Wallace, H = 12.87, p < 0.001—boxplots in Fig. S2). Within sample types, we found that alpha diversity was not significantly affected by CaCO₃ presence in either soil or snail samples.

Taxonomic results

We used taxonomic analyses to assess differences between samples based on microbe identity. Large differences were clear between snail and soil samples based on taxonomy (Fig. 3). For example, soil samples had significantly less of the phylum Bacteroidota and significantly more Crenarchaeota, an Archaeal phylum. In addition to having more Bacteroidota in the snail gut samples, there was also much more Deinococcota. At the family level, there was a lot of diversity, with only a few taxa being over 5% of any group of samples. The most abundant of these were Nitrososphaeraceae in soils and Spirosomaceae, Sphingomonadaceae, and Trueperaceae in snails.

To find core microbes, we looked for microbes that appeared in most snail samples. Only one ASV, family Trueperaceae (phylum Deinococcota), was found to be present in every snail sample. This ASV was also found in many of the soil samples at very low (<1%) relative abundance. In fact, nearly 85% of all samples contained this ASV. There were four other ASVs that were found in around 85% of Oreohelix snail samples. These were from the following families: Spirosomaceae, Sphingobacteriaceae, Sphingomonadaceae, and Nocardioidaceae. Spirosomaceae was the most abundant family found in the snail samples and was hardly found at all in soil samples. Sphingobacteriaceae was found at low relative abundance in snails, but even lower relative abundance in soils. Sphingomonadaceae was found in higher relative abundance in snails but was also present in most soil samples. Nocardioidaceae was found at similar relative abundances in snail and soil samples.

There were some snail samples that were dominated by outlier taxa. The main instance of this was bacteria in the family Francisellaceae, which dominated (over 20% of sample total—sometimes over 50%) the communities of samples from three disparate sites and was either absent or at much lower relative abundance in the rest of the samples.

Differential abundance

Differential abundance testing, using ANCOM-BC, confirmed the differences between sample types with many differentially abundant taxa. There were 93 microbial families that were differentially expressed between soil and snail samples (Table S2). Within snail samples, we found that there were two families, Sandaracinaceae and Hyphomonadaceae, that were differently abundant between samples on CaCO₃ rock and off. Both were found to be more abundant in the off-limestone samples. For the subset dataset using only RR and PITS sites, we also looked at taxa significantly different based on CaCO₃ proximity and found four significant families including Sandaracinaceae again as well as Chroococcidiopsaceae, Rubrobacteriaceae, and Pyrinomonadaceae. In soil samples, the only differently abundant families were unclassified bacteria in the phylum Chloroflexi and the family Amb-16S-1323 from the Hyphomicrobiales order. There was an increase in the Chloroflexi bacteria in samples on limestone and the Amb-16S-1323 group was more abundant off limestone.

Differences between host and environment

We found that the host snails’ microbiomes were significantly different from that of soil samples. There were differences in alpha diversity, dispersion, and composition between sample types. Out of the alpha diversity measures, only Shannon’s diversity had high detection power so unsurprisingly, only Shannon’s diversity between sample types was shown to be significantly different, with soil diversity being higher. The similarity of snail microbiome richness to soils, known to be one of the richest reservoirs of microbial diversity, mirrors the results from previous work that shows that Oreohelix has elevated richness compared to other land snail taxa (Chalifour & Li, 2021). The dominant taxa were different between these samples. While they shared the dominant Phyla Proteobacteria and Actinobacteria, soils contained more Thermoproteota (an Archaea), and snails contained more Bacteriodota and Deinococcota. The differences between the composition of snail samples and their environment were also reflected in the results from differential abundance. We found many differently abundant taxa between sample types. Recent studies have also shown differences between Oreohelix host and surrounding environment, with snail gut samples more closely resembling those from the plant phyllosphere (Chalifour & Li, 2021; O’Rorke et al., 2015). These studies also found that some of the same bacterial families, Sphingobacteriaceae, Sphingomonadaceae, and Trueperaceae, dominated the snail microbiomes from populations across the large native range (Chalifour, Elder & Li, 2023). Lastly, Francisellaceae is extremely abundant in a few snail samples, spanning multiple sites, but not found in over half of the samples. This family appears to become successful when present but is not ubiquitous. There are known endoparasites within this family (Duron & Gottlieb, 2020), but more work needs to be done to identify if there is a meaningful relationship between this microbe and the snail host. While we decided not to speculate on the function of specific microbes using the present dataset, future work will aim to identify the sources of these microbes from phyllosphere and detritus samples, and understand the function of core taxa (present in all samples and relatively abundant) in the snails, such as Trueperaceae.

Effects of CaCO₃ rock on microbial communities

We found that the presence of CaCO₃ rock had a small and significant effect on the compositional differences of snail gut microbiomes. In soils, however, we found no effect of CaCO₃ rock on compositional differences. We also found that neither sample type showed a difference in alpha diversity measures based on CaCO₃ rock presence, though our detection power was low for those tests. Furthermore, geographic distance correlated with microbial dissimilarity in soils, but not snail samples. This could be due to dispersal limitation of soil microbes or another factor that changes spatially but was not measured in this study. We found that taxonomy (for snails), geographic distance (for soils), and site (for both) were important factors driving differences between samples. The sampling design of our study limited the inferences we can make about the effects of host taxonomy and site variables. Host taxonomic differences could be caused by differences in diet, habitat use, or differences in vertical transmission while site differences could be due to an even larger variety of factors including vegetation, aridity, and population structure for available environmental microbes. Each of these factors could be valuable avenues for future work to understand how the environment is affecting microbial community assembly in the snail gut. Previous work in Bornean microsnails living on limestone outcrops (Hendriks et al., 2021) found a link between host diet diversity, host diversity, and microbiome diversity. Resource availability and distance to caves where the snails are from were found to have significant effects on the snail microbiome diversity. Further work in similar land snail systems is needed to identify and study key microbes involved in host-habitat relationships, which could enable meaningful comparisons in microbial taxonomy. The differences in microbiomes between snails on and off CaCO₃ rock could reflect functional needs of snails in their habitat or important host filtering. We would like to note that recent work by Chalifour, Elder & Li (2022) found that the process of snail tissue preservation can have a significant compositional impact on the microbiome, with relative increases in Enterobacteriaceae in preserved samples compared to fresh ones. Interestingly, they also show that after the initial preservation process, tissue stored in ethanol preserved a stable microbiome composition for decades. In the present study we did not test for the effects of preservation on Oreohelix microbiomes, as we processed the samples before (Chalifour, Elder & Li, 2022) was published. We also did not observe high relative abundances for Enterobacteriaceae and therefore cannot make any claims regarding these effects.

The finding that the gut microbiome of snails were significantly impacted by the presence of CaCO₃ rock but the soil microbial communities were not could point to a dosage effect driven by the snails’ rasping (feeding by scraping the substrate with their radula) of the CaCO₃ rock surface. For example, perhaps relatively little CaCO₃ is being eroded from weather into soils and we know that snails are actively rasping the rock surface during their waking hours and likely ingest more of the substrate than would be leached into the soils. Captive breeding, calcium amendments, and transcriptomics could be used to provide more insight into any functional changes or dosage effect that CaCO₃ rock might have on these microbial communities.