Medium sized Tegula eiseni, Tegula funebralis, and Tegula gallina adults 15-20mm in shell diameter (n = 50 per species) were collected in the summer of 2022 from the southern California site Bird Rock in San Diego County (32°48’N, 117°15’W; Figure 1 ). Ideally Tegula would be collected from multiple southern California locations, but as of 2022, sites that have been used for previous research (Gleason and Burton, 2013, 2016b) no longer support robust enough populations of T. funebralis for collection from the field (Sato, 2001). Within 12 hours of collection from Bird Rock, snails were transported to California State University, Sacramento (Sac State) in plastic Nalgene bottles ~one quarter filled with seawater from the collection site. Bottles were stored in a cooler on ice to prevent extreme temperature fluctuations during transport. At Sac State snails of all three species were moved to a flow-through recirculating saltwater aquarium system containing artificial seawater created using Instant Ocean Sea Salt mixed with non-sterilized DI water to 32 ppt. All tanks were filled approximately 80% with water, which allowed Tegula individuals of all species to position themselves either fully submerged underwater or towards the top portion of the tank out of water, in air. The aquarium system was set to 15°C (with an allowable offset of 2°C) and no temperature deviations above 17°C occurred. Snails of all three species were regularly fed dried green algae sheets ad libitum. Before the start of any experiments, snails (n = 50 per species) were kept at these temperature and food conditions for a common garden acclimation period of three weeks. Snails not used in the experiments for this current study were used for other projects in the Gleason Lab.

Following the common garden period T. eiseni, T. funebralis, and T. gallina individuals were exposed to either 1) control 15°C conditions (n = 3 per species) or 2) a 34°C heat stress over 5.5 hours to simulate a single low tide event, following the same 3°C increase per 30 min ramping protocol as described in Gleason and Burton, 2013 (n = 3 per species). Heat stress and control exposures were conducted in air to simulate low tide conditions Tegula individuals experience in the field. Each individual snail was placed in an empty, sterile, capped 50 mL falcon tube that was then placed in a temperature-controlled water bath to reach the desired experimental temperatures (either 15°C for controls or 34°C for heat stress). A foam tube holder was used to float each falcon tube at the top of the water line and to ensure that each tube was held underwater and thus exposed to the desired temperatures. At the conclusion of the control or heat stress exposures, snails were frozen in liquid nitrogen and kept at -80°C. Samples were collected immediately after the 5.5-hour temperature exposure because a heat stress event of this duration is sufficient to cause mortality for all Tegula species, albeit at different maximum temperatures (authors’ unpubl. data). Therefore, this study examines the potential microbiome responses correlated with these survival differences. For sample processing, the shell of each frozen Tegula individual was removed with an ethanol-cleaned woodworking vise, the remaining body was rinsed with ethanol to remove surface bacteria, the gonads and muscular foot were removed, and DNA was extracted from 250 mg of the remaining whole-body tissue using a Qiagen DNeasy PowerSoil Pro kit. This tissue processing ensured that 1) microbiome differences based on gonadal sex differences did not influence the data, and that 2) the microbiome living in/on organs such as the gill, stomach, heart, etc. were examined (as opposed to the exterior surface bacteria on the muscular foot). A previous study found no differences in the microbiome of individual organs in T. funebralis and T. gallina (Neu et al., 2019), and thus whole-body tissue samples were used here. Using whole body samples, as opposed to a single organ such as gill tissue that is in direct contact with the surrounding seawater, also minimizes the likelihood that the microbiome recovered from the Tegula hosts is heavily influenced by the ambient seawater microbiome (Li et al., 2022). No negative controls were used during DNA extraction because swabbing the air or an empty microfuge tube is not a true negative control for a whole-body tissue sample (Hornung et al., 2019); however, we were consistent when taking all samples, which is important for minimizing technical variation and is especially important when a true negative control is difficult to obtain (Vandeputte et al., 2017). All extracted DNA from 18 total samples was quantified with a UV-Vis Nanodrop spectrophotometer (range 51.2 – 329.9 ng/uL).

A 350 bp region of the bacterial V4 region of the 16S rRNA gene was amplified using indexed forward (5’-GTGCCAGCMGCCGCGGTAA-3’) and reverse (5’-GGACTACHVGGGTWTCTAAT-3’) primers (Kozich et al., 2013) using the following protocol: 94°C for 3 min, 35 cycles of 94°C for 45 sec, 50°C for 60 sec, 72°C for 90 sec, and finally 72°C for 10 min. PCR reactions were performed in 25 uL volumes and contained 12.5 uL 2x Phusion High-Fidelity DNA Polymerase Master Mix (NEB), 5 uL ddH2O, 1.25 uL of each primer (10 uM), and 5 uL of template DNA. Amplifications were verified on a 1% agarose gel with GelRed and pooled. A double-sided bead cleanup using Cytiva Sera-Mag SpeedBeads Carboxyl Magnetic Beads (Fisher Scientific) was carried out to remove primer-dimers and a low amount of off-target larger PCR products. Quality and concentration of the pooled library was checked using a Bioanalyzer (Agilent, Santa Clara, CA, USA) and NEB’s Library Quant Kit for Illumina following the manufacturer’s standard protocol. The library representing 18 total samples was then sequenced on an Illumina MiniSeq, yielding 300bp PE reads, using a mid-output reagent cartridge at the James Madison University Center for Genome and Metagenome Studies (CGEMS). Before loading, the library was combined with Illumina’s PhiX control (30:70 16S:PhiX) to ensure a high-quality run despite the low diversity of the 16S library.

All 16S rRNA amplicon sequences were processed using the QIIME2 bioinformatics pipeline as implemented in the Purple Line workflow of DNA Subway that is designed for analyzing microbiome metabarcoding data (Bolyen et al., 2019). Before bioinformatics analysis all primers and sequencing adapters were removed. DADA2 (Callahan et al., 2016) was used to trim low quality reads using the following parameters: truncLenF = 250; truncLenR = 231. After trimming to retain only high-quality reads, all samples were rarefied to a maximum depth of 33,529 sequences. We used a sampling depth of 3933 and the classifier Greengenes2 (515F/806R) to identify bacterial taxa in each sample. The Greengenes2 classifier was chosen due to its compatibility with NCBI taxonomy, optimization for 16S rRNA amplicon sequencing, and design to work seamlessly with QIIME 2 workflows and plugins. The plot_bar function of the phyloseq package was used to visualize the relative abundances of the 25 most common bacterial genera with abundance counts of 400 and above.

Alpha diversity was calculated in the metabarcoding Purple Line workflow of DNA Subway using Pielou’s Evenness (Pielou, 1966) and Faith’s Phylogenetic Diversity (Faith, 1992), following the approach of previous microbiome papers such as O’Connell et al., 2018. Pielou’s Evenness assesses, independent of sample size, the relative abundance distribution of species, with a higher value indicating more evenness (i.e., a relatively equal number of bacteria in each taxon). Faith’s Phylogenetic Diversity assesses how many phylogenetically diverse bacterial taxa are present, with a higher value indicating more richness (i.e., more evolutionarily distinct bacterial taxa). Pielou’s Evenness and Faith’s Phylogenetic Diversity are appropriate alpha diversity metrics for our dataset because, unlike alternative metrics such as Chao1 and Observed ASVs, these chosen metrics do not depend on singleton values that are removed as part of the denoising algorithm during DADA2 processing (Cassol et al., 2025). Moreover, given that our sample size per treatment group is relatively small, Faith’s Phylogenetic Diversity is an appropriate metric for this study because it requires a lower sample size to identify statistical significance compared to alternative alpha diversity metrics such as Shannon and Chao1 indices (Kers and Saccenti, 2022).

We used Bray-Curtis distances to assess beta diversity of samples. NMDS plots were created using the metaMDS function in the R package vegan (Oksanen et al., 2001). To determine whether species or treatment significantly affected the microbiome, we conducted a Permutational multivariate analysis of variance (PERMANOVA) using the adonis2 function in the R package vegan with 9,999 permutations. We used the betadisper function of the vegan package to assess PERMANOVA results for heterogeneity of variance. We used the linear discriminant analysis Effect Size (LEfSe) tool (Segata et al., 2011) as implemented in the R package yingtools2 version 0.0.1.174 (Taur, 2024) to identify specific microbial taxa that are enriched, or found significantly more often, in each Tegula species compared to the others. The R package ggplot2 was used to generate and customize all figures. Lastly, for each treatment (control vs. heat stress), bacterial genera that are found 1) only in a single species, 2) in two of the three species, and 3) in all three species were identified and visualized in Venn diagrams using the online tool available at https://bioinformatics.psb.ugent.be/webtools/Venn/.

A total of 348,662 reads were sequenced across all 18 samples (9 total control samples plus 9 total heat stress samples), and 284,649 reads were retained after all filtering steps. From these filtered reads 1406 total amplicon sequence variants (ASVs) were identified ( Figure 2 ). The dominant taxa identified in Tegula samples were the phylum Proteobacteria. Regarding alpha diversity, average (± SEM) Pielou’s Evenness values for T. eiseni, T. funebralis, and T. gallina were 0.63 ± 0.042, 0.07 ± 0.041, and 0.69 ± 0.027, respectively ( Figure 3A ). Average (± SEM) Faith’s Phylogenetic Diversity values for T. eiseni, T. funebralis, and T. gallina were 11.27 ± 0.59, 15.19 ± 1.09, and 14.19 ± 1.42, respectively ( Figure 3B ). There were no significant differences in alpha diversity across the three Tegula species using either Pielou’s Eveness or Faith’s Phylogenetic Diversity indices (Kruskal-Wallis [pairwise], p>0.05; Figure 3 ). Therefore, there are no significant differences in the relative abundance distribution of bacteria (Pielou’s Evenness) or in the amount of evolutionarily distinct bacterial taxa (Faith’s Phylogenetic Diversity) across the three Tegula species. PERMANOVA results indicate that the microbiome communities were significantly different across the three Tegula species (df = 2, F = 2.758, p = 0.0016, and betadisper p = 0.392). Control vs. heat stress treatments were not significantly different from each other (df = 1, F = 1.958, p = 0.0511, and betadisper p = 0.0568; Figure 4 ), although the interaction between species and treatment was significant (df = 2, F = 1.789, p = 0.0370).

Significantly enriched bacterial genera with higher relative abundance were identified in each of the three Tegula species ( Table 1 ). For T. eiseni, the genera Tunicatimonas, Halomicronema, and Mycobacterium were enriched (linear discriminant analysis [LDA] log scores = 2.89, 3.24, and 4.90, respectively; Figure 5 ). T. funebralis had four enriched bacterial genera: Polynucleobacter, Shewanella, Fluviicola, and Rubritalea (LDA log scores = 2.67, 3.59, 3.73, and 3.89, respectively). Lastly, the three bacterial genera Pelagicoccus, Polaribacter, and Lutimonas (LDA log scores = 2.12, 3.26, and 3.52, respectively) were significantly more abundant in T. gallina. These enriched bacteria are unlikely to be a sampling artifact of the ambient environment because they occur at relatively low frequencies ( Figure 2 ); in previous Tegula microbiome studies, most bacterial OTUs were not shared between T. funebralis and algal and water samples, and the OTUs that were shared were the most abundant ones (Neu et al., 2019).

As shown in Figure 6A , under control conditions T. eiseni had the least unique bacterial taxa (36 genera, 26.5% of the total genera identified; Supplementary Table 1 ), including three members of the Flammeovirgaceae family. In contrast, T. gallina had the most unique bacterial taxa (50 genera, 38.2% of the total genera identified). Unique taxa in T. gallina include two members of each of the following families: Alcaligenaceae, Cytophagaceae, Peptostreptococcaceae, and Rhodobacteraceae. Under heat stress conditions this pattern was different: T. gallina had the least unique bacterial taxa (20 genera, 15.4% of the total genera identified; Supplementary Table 2 ), including two members each of the Flavobacteriaceae and Verrumicrobiaceae families. Heat-stressed T. funebralis had the most unique bacterial taxa (61 genera, 32.3% of the total genera identified), including three members each of the families Moraxellaceae and Pirellulaceae ( Figure 6B ). There were three bacterial genera found only in T. eiseni under both control and heat stress conditions: Owenweeksia, Tunicatimonas, and Jannaschia. Five bacterial genera were only found in T. funebralis under both control and heat stress conditions: Clostridium, Aquicella, Croninitomix, Plantomycete, and Roseovarius. There were no taxa present only in T. gallina under both control and heat stress conditions.

Proteobacteria is the most common bacterial phyla across all three Tegula species – this is consistent with previous work in marine animals such as snails, limpets, red abalone, corals, copepods, fish, and barnacles (Bayer et al., 2013; Brown et al., 2020; Dorosz et al., 2016; Dudek et al., 2014; Givens et al., 2015; Neu et al., 2019; Ousley, 2023). However, at the genus level the microbial community composition significantly varies across Tegula species. Notably, these differences persisted after a common garden acclimation period in which all individuals were housed in the same flowthrough aquarium system, ate the same dried algae diet, and were exposed to the same submersion times in the same recirculating artificial seawater. Overall, these differences across species are consistent with previous work in Tegula (Chlorostoma) eiseni and Tegula (Chlorostoma) funebralis examined directly from the field from La Jolla, California (Neu et al., 2019). Our results also match findings in other marine species such as sponges and algae. De Castro-Fernández et al. (2023) reported that in four different demosponge species, samples clustered together in non-metric multidimensional scaling (NMDS) space according to species (De Castro-Fernández et al., 2023). Moreover, in an experimental design similar to this current study, the microbiome composition and structure significantly varied across three Fucus algae congeners occupying different heights of the intertidal zone and thus experiencing different levels of abiotic stress (Quigley et al., 2020). Notably, our results provide evidence for phylosymbiosis, in which the microbiomes of more closely related species are more similar to each other. As seen in Figure 4 , in multivariate NMDS space the microbiome of T. eiseni is distinct from the other two species T. funebralis and T. gallina, who are more closely related to each other (Hellberg, 1998). Such phylosymbiosis has also been observed when comparing the microbiomes of Chlorostoma (now Tegula) and Littorina intertidal snails (Neu et al., 2019) and of different tropical sponge species (Easson and Thacker, 2014).

Evolutionary codivergence of the microbiome with these Tegula species is one possible explanation for microbiome differentiation, but several other explanations, such as habitat filtering, should also be considered and investigated further (Mazel et al., 2018). Microbial communities are sensitive to environmental parameters such as pH and temperature; potential differences in the gut environment between T. eiseni, T. funebralis, and T. gallina could be contributing to their distinct microbiomes (Neu et al., 2019). Neu et al. also hypothesized that field-collected T. eiseni and T. funebralis could be differentially ingesting microbes through distinct dietary inputs (Neu et al., 2019). Moreover, the cirri microbiome differences in barnacles occupying low vs. high intertidal microhabitats was thought to be due to differential exposure to and/or time underwater (Brown et al., 2020). Comparing the microbiome of T. eiseni, T. funebralis, and T. gallina individuals collected directly from the field to those that have been in a lab common garden environment could help determine how much influence the distinct microhabitats of each species have on their microbiome. Similarly, explicitly manipulating diet and/or time underwater in a laboratory setting could also provide further insight into the mechanisms of microbial differentiation across species.

We did not observe significant microbiome differentiation between individuals exposed to control vs. heat stress conditions for 5.5 hours. The microbiome of Mytilus coruscus larvae exposed to elevated seawater temperature for 4 hours was also similar to the microbiome of control individuals (Zhu et al., 2024), perhaps indicating that a heat stress period longer than several hours is required to elicit clear differentiation of the microbiome. This lack of differentiation contrasts what has been observed in other intertidal marine mollusks. In Haliotis rufescens red abalone, another type of marine snail, control individuals had a distinct microbial composition compared to heat stressed individuals (Ousley, 2023). Similarly, in the intertidal Sydney rock oyster Saccostrea glomerata exposure to elevated heat wave temperatures significantly changed bacterial community composition (Scanes et al., 2023). Microbial communities also changed in response to elevated temperature in the mussel Mytilus galloprovincialis (Li et al., 2019). One potential cause for the lack of significant differentiation between control and heat stress Tegula microbiomes is high interindividual variability magnified by a small sample size (n = 3 per species per treatment). For example, as shown in Table 1 , the abundance of bacterial genera varied widely across individuals of the same species and treatment, especially among T. eiseni control samples with regards to Mycobacterium and among T. gallina heat stress samples with regards to Lutimonas and Polaribacter. Ultimately, the small sample size used in this study limits the statistical power and our ability to reliably detect meaningful biological difference across treatments. Thus, the fact that we did not identify significant differences in the microbiome between control and heat stress treatments could be an artifact of a small sample size.

In T. eiseni, the species with the lowest thermal tolerance, enriched bacteria could be contributing to infection and apoptotic cell death under heat stress conditions. Mycobacterium is much more abundant in T. eiseni under both control and heat stress conditions compared to T. funebralis and T. gallina. This is a pathogenic bacteria that is known to infect other marine mollusks (Davidovich et al., 2020), and enrichment of this taxa could suggest an increased risk of disease in T. eiseni, especially under heat stress conditions. Another taxa enriched in T. eiseni, Halomicronema, is a cyanobacteria that has also been found in marine sponges (Caroppo et al., 2012). Notably, some compounds produced by the Halomicronema genus can trigger apoptotic cell death (Mutalipassi et al., 2019). In other marine mollusks such as the oyster Crassostrea virginica, the density of apoptotic cells increased after exposure to high temperatures of 26 and 30°C (Rahman and Rahman, 2021). Rahman and Rahman (2021) hypothesized that a high level of reactive oxygen species (ROS) could be contributing to this increase in apoptosis (Nash et al., 2019; Nash and Rahman, 2019). Our results suggest components of the microbiome could also be contributing to the induction of apoptosis in heat sensitive marine mollusks, and should be investigated further. The last enriched bacterial genus in T. eiseni, Tunicatomonas, has also been isolated from sea anemones, but there is only one known species, and not much information is available on this taxon (Yoon et al., 2012).

Three bacterial genera were significantly enriched in T. funebralis and were more abundant in heat stress vs. control conditions: Rubritalea, Fluviicola, and Polynucleobacter. Rubritalea belongs to the class Verrucomicrobiae (Yoon et al., 2007) and bacteria of this genus also increase in abundance following high temperature exposure in corals and sponges. For example, in Turbinaria peltata corals repeated heatwaves led to increases of the beneficial bacteria Rubritalea tangerine that facilitates coral health and growth (Zhai et al., 2024). Zhai et al. (2024) hypothesized that such changes in the microbiome could represent adaptive stress responses to improve survival of the coral host under marine heat wave conditions. Moreover, heat wave exposed Crella incrustans sponges produced larvae with significantly more Rubritalea marina (Strano et al., 2023). How this bacterial genus functions during high temperatures, and the potential mechanisms that allow these bacteria to benefit marine invertebrate hosts, remain unknown and require further investigation. Fluviicola, the second genera enriched in T. funebralis, is a gram-negative rod-shaped bacterium in the class Flavobacteriia, family Cryomorphaceae. In this current study Fluviicola was more abundant in heat stress vs. control samples, which contrasts previous results in green-lipped mussels: in Perna canaliculus Fluviicola was more abundant in control vs. heat stress individuals (Ericson et al., 2024). Other studies in marine mollusks found that Fluviicola abundance was high in razor clams exposed to hyposalinity stress (Yang et al., 2024) and that Fluviicola was more abundant in younger vs. older Pacific oyster spat (Zhong et al., 2024). Lastly, bacteria in the genus Polynucleobacter were also significantly more abundant in T. funebralis. Polynucleobacter can be either free-living or symbiotic (Miklós et al., 2023) and are often associated with ciliates (Vannini et al., 2007) and hydra (Fraune et al., 2010; Fraune and Bosch, 2007). In the cold-water adapted hydra Hydra oligactis Polynucleobacter was positively correlated with hydra population size, although the mechanism for this effect requires further research (Miklós et al., 2023). In contrast to the current findings in this study, Polynucleobacter abundance was lower at higher temperatures in H. oligactis (Miklós et al., 2023). Overall, any hypotheses about the potential relationship between Polynucleobacter and T. funebralis should be interpreted with caution because the overall abundance of Polynucleobacter was relatively low ( Table 1 ).

In contrast to the other three bacterial genera that were significantly more abundant in T. funebralis, the genus Shewanella was more abundant in control vs. heat stress samples. These results match previous papers that have identified bacteria in this genus that are adapted to low temperatures (Kloska et al., 2020; Zhao et al., 2010). In general, Shewanella are gram-negative, aerobic and facultatively anaerobic γ-Proteobacteria (Garrity and Holt, 2001; Gauthier et al., 1995; MacDonell and Colwell, 1985; Satomi et al., 2003; Venkateswaran et al., 1999). Shewanella have been isolated from a variety of marine invertebrates including hydrocorals, the benthic marine “peanut worm” (Ivanova et al., 2004), the sea urchin Sterechinus neumayeri (González-Aravena et al., 2016), and the soft coral Alcyonium antarcticum (Webster and Bourne, 2007). Notably, exopolysaccharides (EPS) produced by Shewanella have been isolated from Antarctic sponges (Caruso et al., 2018). These EPS were more abundant at low temperatures, and they are thought to serve protective functions under low temperature stress conditions (Lo Giudice and Rizzo, 2022). Although this current study did not expose any individuals to low temperatures, the fact that Shewanella was significantly more abundant in T. funebralis, the only Tegula species examined in this study whose geographic range extends up to Vancouver Island on the western coast of North America, suggests that this bacterium could be contributing to T. funebralis’ ability to withstand the lower temperatures indicative of the northern parts of its range. For example, sites occupied by T. funebralis in northern California can reach absolute minimum temperatures of 8.8°C, which is roughly 5-6°C colder than the absolute minimum temperatures in southern California, where T. eiseni and T. gallina reside (Gleason and Burton, 2016b).

Lutimonas was found in relatively high abundance in T. gallina under both control and heat stress conditions. These results contrast a previous study that observed a decrease in Lutimonas under heat stress conditions in Pacific white shrimp Litopenaeus vannamei (Duan et al., 2021). Lutimonas is a nitrifying bacteria that degrades ammonium (Fu et al., 2009). In other marine mollusks such as the ark shell Scapharca subcrenata (Jiang et al., 2020) and in Daphnia (N. Nash et al., 2022), elevated temperatures increase ammonia excretion rates. High levels of ammonia are known to have an array of negative consequences for aquatic invertebrates (Zhang et al., 2023); thus, the relatively high abundance of Lutimonas in T. gallina, even under heat stress conditions, could suggest that the microbiome helps regulate ammonia levels to prevent toxic overaccumulation.

Our finding that the genus Polaribacter is enriched in T. gallina, including under heat stress conditions, contrasts a previous study in another marine mollusk Mytilus galloprovincialis that found Polaribacter was a dominant genus under control conditions, but its abundance decreased when water temperature was increased (Li et al., 2019). However, in other marine invertebrates, including the mussel Mytilus coruscus (Li et al., 2019) and the common yellow sponge M. acerate (De Castro-Fernández et al., 2023), abundance of Polaribacter was also high under heat stress conditions. Polaribacter has also been detected in seawater (Fukui et al., 2013; Yoon et al., 2006) and in diatom phytoplankton blooms (Xing et al., 2015). Overall, not much is known about this genus, and further research must be done before hypotheses regarding the ability of this genus to increase heat tolerance of T. gallina can be formed.

The third and final bacterial genus that was enriched in T. gallina is Pelagicoccus. Pelagicoccus, a member of the family Puniceicoccaceae, contains four different species (Feng et al., 2021; Yoon et al., 2007). These species have previously been isolated from sea grass (Yoon et al., 2007), seawater (Yoon et al., 2007), and marine sediment (Feng et al., 2021). Based on the genomic analysis performed by Feng et al., it is thought that all bacteria in this genus produce exopolysaccharides (EPS) to better cope with extreme stress conditions, including high temperatures, low nutrients, drought, salinity stress, and antimicrobial agents (Liu et al., 2013; Nichols et al., 2005; Wagh et al., 2022; Wang et al., 2019). For example, EPS are produced in hydrothermal vent communities exposed to extreme high temperatures (Nichols et al., 2005). Importantly, there is evidence in other species that EPS produced by bacteria have positive effects on their host: EPS produced by plant-associated plant-growth promoting rhizobacteria reduce the negative effects of heat shock and growth in plant hosts such as tomatoes (Morcillo and Manzanera, 2021). It is possible the Pelagicoccus observed in this study could similarly be contributing to host thermal tolerance in T. gallina, although current evidence is merely correlational and it is not currently known if the relatively low abundance of Pelagicoccus would be sufficient to produce beneficial effects. Ultimately, more direct functional experiments are needed to confirm whether this particular bacterium plays a role in T. gallina’s high thermal tolerance.

As with any study, there are limitations that should be considered. Most notably, we did not collect or sequence environmental samples from the field or the lab to characterize the surrounding microbial habitat (e.g., seawater and/or artificial seawater in the aquarium system) of the Tegula snails. Although we did surface sterilize each sample by thoroughly rinsing tissue samples with ethanol before extracting DNA (Brown et al., 2020; Neu et al., 2019), we nevertheless cannot determine whether the Tegula microbiome communities we characterized are subsets of the surrounding bacterial community in the artificial seawater and/or the algal diet of the snails. In particular, because Tegula snails were housed in non-sterile artificial seawater during the three-week common garden period, it is possible this water was an environmental source of bacteria that contributed to the Tegula microbiomes analyzed in this study. We do, however, note that a similar study conducted in the Gleason Lab’s same aquarium system found that the fecal microbiome of juvenile red abalone Haliotis rufescens was significantly different from the microbiome of the surrounding aquarium water (Ousley, 2023). Regarding the potential effects of feeding snails dried green algae sheets, a previous study of Haliotis rufescens red abalone marine snails fed red vs. brown vs. green algae found significant differences in the microbiomes of snails fed the distinct diets (Guo, 2017). Thus, the choice of diet in our current study could have affected the microbiome composition of the Tegula snails. Specifically, Guo (2017) found that bacteria in the Mollicutes and Alphaproteobacteria classes and in the Mycoplasma genus were more abundant in H. rufescens fed green algae diets compared to red and brown algae diets. Thus, the abundance of these bacteria may also be enriched in our current study that used a green algae diet.

Additional limitations of study design include the fact that this study used a relatively small number of samples (n = 3 per species per treatment). As noted above, this small sample size could have amplified the high degree of interindividual variability and contributed to the lack of significant variation detected between control and heat stress microbiome samples. Although previous similar microbiome studies in marine snails have used similar sample sizes (Y.-J. Zhu et al., 2021), we do acknowledge that our results may not be representative of each Tegula species’ full range of biological responses to high temperature conditions. Similarly, this work investigates the microbiome from a single field site in San Diego, California. Neu et al., 2019 found that for T. funebralis, Bray-Curtis dissimilarities significantly distinguish snails from north vs. south of Point Conception along the California coast (although the southern site only had a sample size of two). Thus, our findings likely cannot be generalized to T. funebralis individuals from further north in California. Previous studies using RNA sequencing and ddRAD sequencing have identified gene expression and gene sequence differences across geographically distinct T. funebralis populations in northern and southern California that experience unique climates and differ in thermal tolerance (Gleason and Burton, 2015, 2016a). To determine if there are also microbiome differences within a single Tegula species across different geographic regions, future work could sample and subsequently compare the microbiome of T. eiseni and T. gallina from multiple additional sites in southern California (e.g., Abalone Cove in Los Angeles County, Aliso Beach in Orange County, La Jolla in San Diego County), and T. funebralis from northern (e.g., Slide Ranch in Marin County, Pigeon Point and Pescadero in San Mateo County) and southern (e.g., Abalone Cove in Los Angeles County, Aliso Beach in Orange County, La Jolla in San Diego County) California sites, population size permitting (Sato, 2001). Lastly, this study only examines the bacterial component of the microbiome. Any fungal (Chin et al., 2020; Grice, 2015), archaeal, and viral components of the holobiont that could also be contributing to thermal tolerance or host-microbiome interactions are thus not captured in our current dataset.

At this point it is unknown if the bacterial species we observed in the T. eiseni, T. funebralis, and T. gallina microbiomes are transient (i.e., driven by external factors or stochastic processes) or resident (i.e., explicitly selected for by the host). Resident bacterial species that are a stable, persistent part of the microbiome are more likely to contribute to host thermal tolerance, and thus distinguishing resident vs. transient bacteria in each Tegula host remains an important area of future research. To distinguish these different factions of the Tegula microbiome, further research could follow the approach of Unzueta-Martínez et al. (2022) to differentiate transient vs. resident bacteria in the Eastern oyster Crassostrea virginica. To identify resident bacteria, we could 1) identify any microbiome differences in Tegula species collected directly from multiple different field sites, and 2) determine if these differences persist after several weeks of common garden acclimation in an identical environment. Alternatively, we could determine if certain bacteria persist in each Tegula host in sterile seawater that doesn’t contribute any environmental bacteria (Lokmer et al., 2016, 2016). To identify transient bacteria, we could again collect individuals from multiple different field sites and identify bacterial differences in Tegula individuals held in a common garden environment vs. those sampled directly from the same field site (Unzueta-Martínez et al., 2022).

In addition, to date, we have only characterized the microbiome following a short 5.5-hour heat stress representative of a single low tide period in the field. No assessment has yet been performed of the recovery period following this heat stress; in other words, it is not yet known if the microbial community reverts to the same “baseline” composition after return to non-stressful control conditions. Current marine invertebrate research directly addressing this question is limited, although results in another marine mollusk the green-lipped mussel Perna canaliculus indicate that the microbiome of recovered heat stress samples is more similar to control samples than to sustained heat stress samples (Ericson et al., 2024). Overall, based on the sampling design of this current study, we cannot make any conclusions about the holobiont’s long-term response to heat stress in Tegula species. Thus, more research is needed examining the microbiome at various time points after heat stress to fully understand how the microbiome of Tegula individuals affects thermal tolerance, especially for these intertidal species that regularly experience multiple stressful low tides in a single day.