Ocean acidification alters the diversity and structure of oyster associated microbial communities

Host‐associated microbial communities are fundamental to host physiology, yet it is unclear how these communities will respond to environmental disturbances. Here, we disentangle the environment‐linked and host‐linked effects of ocean acidification on oyster‐associated microbial communities. We exposed adult oysters (Crassostrea virginica) to CO2‐induced ocean acidification (400 vs. 2800 ppm) for 80 d. We measured the oyster extrapallial fluid pH and sampled the gills for microbial analysis at six time points. We found that different subsets of microbes were linked to acidification (n = 34 amplicon sequence variants [ASVs]) and to host response (n = 20 ASVs) with little overlap (n = 8 ASVs), suggesting that some members of the oyster microbiome were more responsive to environmental conditions while others were more tightly linked to host condition. Our results provide insight into which members of the oyster microbiome may contribute to the health and resistance of their host, and which members are the most vulnerable to changing environmental conditions.

pH (Caldeira and Wickett 2003). Prolonged exposure to OA can have negative effects on calcification, growth, and metabolism of marine calcifying organisms (Anthony et al. 2008;Ries et al. 2009;Talmage and Gobler 2010;Melatunan et al. 2013). Therefore, changes in seawater chemistry will have profound effects on marine life at both macro and micro scales.
Marine microorganisms play key roles in the health of marine invertebrates including cnidarians, sponges, and molluscs (Rosenberg et al. 2007;Roeselers and Newton 2012;Webster and Taylor 2012). Despite the importance of microbes to host physiology, the effects of OA on bivalveassociated microbiomes are inconclusive. Some studies found no effect of elevated pCO 2 exposure on the microbial community structure of scallop tissues (Alma et al. 2020), while others found a significant effect of elevated pCO 2 on species diversity of the oyster hemolymph microbiome (Scanes et al. 2021). To clarify our understanding of how bivalveassociated microorganisms respond to OA, it is essential to disentangle how host response to OA may alter microbial community dynamics.
Microbial community structure often shifts as a result of disturbance. In host-associated microbial communities, the host itself presents an additional set of factors that can drive microbial community dynamics. For example, exposure to OA can caused acidosis in oyster hemolymph (Lanning et al. 2010;Scanes et al. 2018), which in turn can affect the associated microbial community structure (Jones et al. 2015). Additionally, compensatory responses to acidosis may affect the host's energy allocations and direct a greater proportion of energy towards acid-base modulation and away from critical fitness-sustaining processes such as protein synthesis and immune response (Porter et al. 2004;Melzner et al. 2009;Gazeau et al. 2013). Disentangling what portion of the microbiome response is linked to environmental disturbance and what portion is linked to the host's response to that disturbance is essential for predicting the effect of global-change drivers on animal microbiomes.
We used the Eastern oyster (Crassostrea virginica) as model organism to study the influence of environmental stress, host response, and the associated microbiome. C. virginica has been used to study the effects of OA on DNA methylation Venkataraman et al. 2020), transcriptomic (Chapman et al. 2011), proteomic (Tomanek et al. 2011), and physiological responses (Waldbusser et al. 2011). This oyster species has also been the subject of numerous microbiome studies, which have revealed a complex microbial community structure within the oyster digestive system (King et al. 2012), gill (Wegner et al. 2013), and hemolymph (Lokmer and Wegner 2015). Despite the eastern oyster's critical importance to the structure and function of estuarine habitats (Beseres Pollack et al. 2013) and to commercial fisheries (Ekstrom et al. 2015), it is still unknown how the C. virginica microbiome responds to OA.
Bivalve responses to OA are highly variable among individuals (Parker et al. 2011;Thomsen et al. 2017), even within the same species (Parker et al. 2011), which has consequences for bivalve calcification (Miller et al. 2009) and physiology (Gobler and Talmage 2014). The calcifying fluid of bivalves, termed extrapallial fluid (EPF), plays a critical role in shell formation (Weeler and Sikes 1984). Changes in EPF chemistry due to OA can hinder the ability of bivalves to form their shells (Ries et al. 2009;Waldbusser et al. 2011), and can alter oyster metabolic rates and energy budgets by diverting energy towards maintenance of homeostasis rather than towards basic processes, such as immune response, reproduction, and growth (Porter et al. 2004;Melzner et al. 2009;Gazeau et al. 2013). As a result, the ability of bivalves to maintain EPF chemistry can influence bivalves' physiological responses to their external environment (Ries et al. 2009;Cameron et al. 2019;Liu et al. 2020).
To disentangle the environment-linked and host-linked effects of OA on oyster microbiomes, we exposed adult oysters to CO 2 -induced OA and measured their extrapallial fluid pH (pH EPF ). We then linked gill-associated microbial community composition and richness to these environment-linked (OA treatment) and host-linked (pH EPF ) drivers. We tested the following hypotheses: both (1) OA exposure and (2) experimental duration affects composition and structure of microbial communities on oyster gills; and (3) oyster physiological response to OA (measured by ΔpH, a metric previously used by Ries (2011) and Liu et al. (2020), among others, as an indicator of pH regulation at the site of calcification) correlates with taxonomic diversity of microbial communities on oyster gills.

Exposure experiment
We collected 84 adult C. virginica (80-100 mm shell length) from three different locations within Plum Island Sound, Massachusetts, United States (42 45 0 05.9 00 N, 70 50 0 13.3 00 W; 42 43 0 30.7 00 N, 70 51 0 18.1 00 W; and 42 40 0 54.4 00 N, 70 48 0 48.6 00 W). Once in the lab, we placed the oysters in separate 50-liter flowthrough tanks for 8 d. We installed a luer-lock port on the right valve of each oyster , returned them to their respective tanks, and waited four days before randomly assigning them to one of two nominal pCO 2 treatments, 400 or 2800 ppm. These treatment levels represent present day atmospheric concentrations (400 ppm) and a pCO 2 level (2800 ppm) that favors the dissolution of inorganic calcite (Ω < 1). Although this value is higher than near future predictions of OA for the open ocean, coastal waters, including the site from which the oysters were collected (Downey-Wall et al. 2020), experience short-term acidification events that exceed even the highest pCO 2 conditions employed in this experiment, resulting from daily tidal cycles, input of meteoric waters, and seasonal oxidation of organic matter. Detailed exposure and seawater Table 1. Summary of the seawater chemistry parameters and extrapallial fluid pH at each time point. Average measured seawater parameters: salinity (Sal), temperature (Temp), pH on the total scale (pH (total scale)), total alkalinity (TA), and dissolved inorganic carbon (DIC). Average calculated seawater parameters: Seawater pH on the seawater scale (pH (SW scale)), partial pressure of CO 2 (pCO 2 ), carbonate ion concentration (CO 3 2À ), bicarbonate ion concentration (HCO 3 À ), dissolved CO 2 (CO 2 ), calcite saturation state (Ω calcite ), and pCO 2 of the mixed gases in equilibrium with experimental seawaters (pCO 2 (ppm-v)).
A randomized block design was used to ensure that equal numbers of oysters from each population were distributed across all treatments. Each pCO 2 treatment was replicated over six 42-liter tanks, with four oysters from each population per tank. The oysters were acclimated at present-day pCO 2 conditions for 30 d, the high-pCO 2 treatment was then ramped up to target pCO 2 over a 12-h period and maintained for 80 d. To produce the target pCO 2 levels, we mixed compressed CO 2 with compressed CO 2 -free air (control treatment) or with compressed air (high OA treatment) using solenoid-valve-controlled mass flow controllers at flowrates proportional to the target pCO 2 conditions. The mixed gases were then bubbled into the experimental seawater treatments. We measured the pH EPF of the oysters during tissue sampling through the installed luer-lock ports (Downey-Wall et al. 2020), as described in greater detail in the Supporting Information (Methods: EPF pH measurements).

Microbiome sample collection
We destructively sampled the oyster microbiomes at seven time points, one during acclimation and six during exposure. The exposure time points were at 24 and 36 h, and 9, 22, 50, and 80 d. At each time point, we sampled six oysters from each treatment. We chose these time points because we wanted to capture both short-and long-term effects. We shucked and dissected the oysters to aseptically remove gill tissues, placed them individually in 2 mL cryovials, immediately froze them in liquid nitrogen, and stored them at À80 C until we extracted the DNA. We chose to examine gill tissues because they are constantly in contact with the surrounding seawater through filtering and are thought to provide an entry point to the oyster for parasites and other pathogens (King et al. 2019).

Nucleic acid preparation and sequencing
We extracted nucleic acids from approximately 0.25 g of each oyster gill using the DNeasy PowerLyzer PowerSoil kit (Qiagen) following the manufacturer's protocol. To amplify the V4 region of the 16S rRNA gene in triplicate 25 μL PCR reactions, we used the primers 515FY: 5 0 TATGGTAATTGTGT GYCAGCMGCCGCGGTAA 3 0 (Parada et al. 2016) and 806RB: 3 0 AGTCAGTCAGCCGGACTACNVGGGTWTCTAAT 5 0 (Apprill et al. 2015). For each sample, we combined and purified the triplicate product using the Agencourt AmPure Magnetic Beads system (Beckman Coulter), and then ligated Illumina paired-end adapters with unique Nextera XT V2 indexes. We quantified our libraries using the Quant-iT PicoGreen dsDNA Assay Kit (Invitrogen) and pooled them at equimolar concentrations. We sequenced our libraries in two sequencing runs using paired end, 2 Â 250 bp kits on an Illumina MiSeq with V2 chemistry at Northeastern University.  Library preparation and quantification details are described in the Supporting Information (Methods: Nucleic acid preparation and sequencing).

Sequence and statistical analysis
For each sequencing run, we used the DADA2 (v1.7.0) workflow with default parameters (Callahan et al. 2016), implemented in R Studio (v4.0.0), to quality-filter, merge pairedend reads, remove chimeric sequences, group the sequences into amplicon sequence variants (ASVs), and assign taxonomy against the Silva database (version 132; Quast et al. 2013). We then merged the ASV tables from both sequencing runs for downstream analyses (Supporting Information Methods: Sequence analysis).
We investigated the environment-linked effects of OA on the oyster gill microbiome by analyzing the effects of: (1) OA exposure and (2) experimental duration. We also considered the influence of population (Supporting Information Methods: Statistical analysis; Results: Oyster population and microbial community composition; Table S1). We ran PER-MANOVA with the "adonis2" function in Vegan (Dixon 2003) and used DESeq2 to identify significantly differentially abundant taxa between the high and control treatment levels at each time point (p < 0.05 and Benjamini and Hochberg 1995 adjusted p < 0.01). To investigate the relationship between the microbes and the host physiological response, we tested whether the difference between pH EPF and pH seawater (ΔpH = pH EPF -pH seawater ), as well as pH EPF alone, correlated with taxonomic richness of microorganisms using linear regression analyses. The ΔpH roughly indicates how much the oysters regulated their pH EPF relative to their surrounding pH seawater in support of calcification. A positive ΔpH indicates that the oyster is maintaining pH EPF relative to pH seawater (most favorable for calcification); a ΔpH of zero indicates that the oyster is maintaining pH EPF equivalent to pH seawater ; and a negative ΔpH indicates that the oyster is maintaining pH EPF below that of pH seawater (least favorable for calcification). We modeled the association between ASV abundance and ΔpH using the "randomForest" R package (Liaw and Wiener 2002) and confirmed model performance by examining the out-ofbag error rate (a method that uses bootstrap aggregation to assess performance without a training set). We also performed leave-one-out cross-validation with 999 permutations in the "caret" R package (Kuhn 2008). Finally, we used set theory to identify ASVs that were either shared or were exclusive to the ASVs that were linked to treatment (differentially abundant ASVs between treatment levels) and ASVs that were linked to host response (Supporting Information Methods: Statistical analysis).

Microbial community composition
Oyster gill microbial communities were highly variable among individuals throughout the entire experiment in both the control and high pCO 2 treatments (Fig. 1). There were noticeable changes over time in control samples, with bacterial orders such as Vibrionales and Alteromonadales having greater relative abundances at the beginning of the experiment than at the end (Fig. 1a). Despite the large variability within treatments and time points, ecologically meaningful treatment effects were observed at individual time points and these effects varied over the duration of the experiment. For example, 24 h after initial exposure, Oceanospirillales had a considerably higher relative abundance in the high-treatment  Table S2 in the Supporting Information for more detailed taxonomic and abundance information.
samples compared to the controls. In contrast, 9 d after initial exposure, Oceanospirillales appears to be less abundant in the high-treatment samples relative to the controls (Fig. 1).

OA treatment effects on microbial community composition
We found that exposure to OA changed significantly the composition and structure of oyster gill-associated microbial communities. When examining Bray-Curtis dissimilarity indices with an NMDS (stress = 0.21), there was a clear separation between treatment centroids (Fig. 2a,b). This observation was further supported by a PERMANOVA that showed significant differences in microbial community β diversity by treatment (Table S1: F = 1.72, df = 1, p = 0.045), time point (Table S1: F = 1.79, df = 5, p = 0.001), and their interaction (Table S1: F = 1.43, df = 5, p = 0.030).

Differential abundance of ASVs by OA treatment
Differential abundance analyses revealed 34 ASVs that had significantly different abundances (DeSeq2, p value < 0.05, Microbe abundance linked to host physiological response (using ΔpH = pH EPF À pH seawater ). Oysters with a ΔpH greater than zero were actively maintaining a higher pH EPF than their surrounding seawater, the majority of these oysters were found in the in the high OA treatment and were associated with increased ASV richness. (a) Scatter plot showing the relationship between ASV richness and ΔpH within each treatment. The dark blue line shows the best fit linear regression (R 2 = 0.2, p = 0.00972) for oysters in the high-pCO 2 treatment. There was no significant relationship between ASV richness and ΔpH of oysters in the control treatment (R 2 = 0.065, p = 0.165). (b) Heat map of the 20 ASVs that are the most important contributors to a randomForest regression model that was trained to predict ΔpH from microbial community composition in oysters exposed to OA. The heat map shows the relative abundance of the ASVs in samples spanning the ΔpH range (left). Bars on the right denote the importance of each ASV in improving model performance and the color of the bars represents the order of the ASV. See Table S3 in the Supporting Information for more detailed taxonomic information.
adjusted p value < 0.01) between control and high OA treatments at any given time point (Fig. 2c; Table S2). At 24 h after initial exposure, there were a number of Vibrionales and Oceanospirillales ASVs that were less abundant in the hightreatment samples relative to the controls. However, at later times, some ASVs were more abundant under OA conditions. For example, at 9 d after exposure began, ASVs from Alteromonadales, Cellvibrionales, and Campylobacteriales, among other orders, had a higher abundance in high pCO 2 conditions relative to the controls (Fig. 2c). Similarly, at 80 d after exposure began, there was a higher abundance of Vibrionales, Alteromonadales, and Campylobacteriales ASVs in samples exposed to OA relative to samples in the control treatment.

Correlation between ΔpH and taxonomic diversity
We used ΔpH (i.e., pH EPF À pH seawater ) as a proxy for host physiological response to OA. We found that ASV richness in oyster gills was positively related to ΔpH of oysters exposed to the high pCO 2 treatment ( Fig. 3a; p < 0.01, R 2 = 0.2), with ASV number increasing as ΔpH increases. In contrast, we found no relationship between ASV richness in oyster gills and ΔpH of oysters exposed to the control treatment ( Fig. 3a; p = 0.165). Similarly, we found that ASV richness in oyster gills was positively related to pH EPF of oysters exposed to the high pCO 2 treatment ( Fig. S1; p < 0.05, R 2 = 0.15), with ASV number increasing as pH EPF increases. Furthermore, we found no relationship between ASV richness in oyster gills and pH EPF of oysters exposed to the control treatment ( Fig. S1; p = 0.18).

Relationship between ΔpH and individual microbial taxa
The randomForest model that determined which members of the microbial community were most strongly associated with ΔpH, indicated 51.86% of the variation in ΔpH explained by ASVs, with mean squared residuals of 0.006. Leave-one-out cross-validation confirmed model performance  with a moderate coefficient of determination (R 2 = 0.46). The accuracy of the model was measured using the mean absolute error (MAE), calculated as the deviation of the predicted from observed values (MAE = 0.05). We identified the 20 most important ASVs in improving model performance and found that the orders Campylobacterales and Oceanospirillales were among the taxa with the highest relative importance in explaining the relationship between ΔpH and ASV richness (Fig. 3b). ASVs from Campylobacteriales, Alteromonadales, and Rhodobacterales, were present in high abundances in oysters that had a ΔpH > 1 (Fig. 3b).
Set theory to compare ASVs that were linked to OA treatment and linked to host response Comparison of ASVs linked to treatment (differentially abundant ASVs between treatment levels) and ASVs linked to host response (ASVs associated with ΔpH) found that only eight of the ASVs overlapped between the two sets (Fig. 4). The overlapping ASVs belonged to the orders Campylobacterales, Cellvibrionales, Vibrionales, Alteramonadales, Rickettsiales, Oceanospirillales, and Pseudomonadales. Twenty-six ASVs, within the orders Holosporales, Flavobacteriales, Parvibaculales, Sphingobacteriales, and Clostridiales were exclusive to the list of ASVs that were linked to OA treatment (i.e., were differentially abundant between treatments). In contrast, 12 ASVs within the orders Bacteroidales, Entorobacteriales, Betaprotobacteriales, and Desulfovibrionales were exclusive to the list of ASVs that were linked to oyster physiological response (e.g., were associated with ΔpH in oysters exposed to OA) (Fig. 4).

Discussion
The oceans will continue to acidify as we continue to increase the partial pressure of atmospheric CO 2 . It is unclear, however, how host-associated microbes will respond to this continued acidification of seawater. Our results show that some members of oyster gill microbial communities are linked to OA treatment, while others are linked to the oyster's physiological response to OA (i.e., ΔpH).

Oyster gill microbes linked to OA treatment
The influence of OA on bivalve-associated microorganisms is inconsistent in the literature. Some studies have found no effect of OA on bivalve microbiomes (Alma et al. 2020), while others have found significant shifts in bivalve microbiome composition and structure as a result of OA (Scanes et al. 2021). Our findings agreed with the latter, showing that C. virginica microbiomes are significantly altered by OA exposure over time.
The inherent variability of animal microbiomes (Mengoni et al. 2013;Utter et al. 2016) was also evident in our study and was exacerbated by the inter-individual variability of bivalve responses to OA (Parker et al. 2011;Thomsen et al. 2017).
We identified taxa that responded to OA exposure over both short and long intervals. At 24 h, most differentially abundant ASVs were more abundant in the control than in the high OA treatment, indicating taxa that were immediately susceptible to OA disturbance. Among the susceptible ASVs were members of the genus Vibrio, which contains multiple species known to decrease their attachment and replication capabilities when exposed to acidic ocean conditions (Huq et al. 1984;Larsen 1984), which could explain their low abundance on oyster gills under acidified conditions. Oceanospirillales were also initially susceptible to OA in this experiment. Members of this order also had significantly lower abundances in marine biofilms under acidification conditions (Ng and Chiu 2020) and on two different coral species near CO 2 seeps (Morrow et al. 2015). Cultured species of Oceanospirilalles exhibit optimal growth under alkaline pH conditions (Wang et al. 2015), which is consistent with their low abundance after the initiation of the high OA treatment.
At intermediate and later time points, there were new and returning ASVs that were abundant in the high OA treatment that were not present or were present in low abundances in the control samples. Among the ASVs that were more abundant in the high-treatment samples at later time points were members of the order Campylobacterales. Most species of this order are capnophiles that thrive in high-CO 2 environments (Al-Haideri et al. 2016;Waite et al. 2017). ASVs in the genus Vibrio initially had low abundances in the high-CO 2 treatment, but at later time points they were significantly more abundant in the high-CO 2 treatment than in the control. When exposed to environmental disturbance, bacteria can adapt and persist by re-adjusting cellular functions that are better suited for the new environmental conditions (Hottes et al. 2013;Mocali et al. 2017), which could explain the resurgence of some ASVs after initially being depleted in the OA treatment.
Oyster gill microbes linked to oyster physiological response We used the difference between pH EPF and pH seawater (i.e., ΔpH) to investigate the relationship between oyster physiological response and microbial communities on oyster gills. Although this relationship is significant, the variance explained by this model is low and should be interpreted with caution (Fig. 3A). However, considering the high interindividual variability in oyster response to OA (Parker et al. 2011), compounded by the high microbiome variability that results from complex microbe-microbe, microbe-host, and microbe-environment interactions in each one of the oysters, it is remarkable that there is any significant amount of microbial community diversity variance that could be explained by a single variable (ΔpH). Additionally, there was no relationship between ΔpH and ASV richness in the control treatment; however, the oysters that maintained a ΔpH above zero were almost entirely found in the high-OA treatment and had higher ASV richness than in the control treatment, suggesting that the OA treatment may affect microbial richness in a way that is mediated by the host treatment response (ΔpH).
Previous studies have shown that marine bivalves can regulate their calcifying fluids in response to OA (Cameron et al. 2019;Liu et al. 2020), this response likely impacts the energy budgets of bivalves as they expend additional energy in trying to maintaining homeostasis (Porter et al. 2004;Melzner et al. 2009;Lanning et al. 2010;Gazeau et al. 2013). We observed that oysters that maintained their pH EPF above pH seawater (i.e., elevated ΔpH), harbored a higher microbial species richness on their gills. This may be because oyster hosts were allocating more energy towards regulating their EPF and not enough energy towards their immune response which can maintain a healthy microbiome. This is consistent with previous literature showing that exposure to elevated pCO 2 can increase the richness of oyster-associated microbiomes (Scanes et al. 2021) and demonstrated that reduced ability of animals to regulate their microbial communities can result in highly variable and speciesrich microbiomes (Foster et al. 2017;Zaneveld et al. 2017).
We identified microbes that were associated with a positive ΔpH, which suggests their abundances may result from the host's physiological response to OA treatment. For example, one of the ASVs associated with a high ΔpH is from the order Campylobacterales genus Arcobacter, a genus commonly found on stressed oysters in postharvest storage (Fernandez-Piquer et al. 2012;Madigan et al. 2014). Other ASVs associated with elevated ΔpH were in the order Vibrionales genus Vibrio. Species in this genus are known oyster pathogens (Richards et al. 2015). Much more work is needed to understand the occurrence of microbes on hosts exposed to OA and how they are modified via host response, to better understand whether OA will increase host susceptibility to disease.
We compared the 34 ASVs that were linked to treatment (differentially abundant ASVs between treatments) and 20 ASVs that were linked to host response (ASVs associated with elevated ΔpH) and found that only eight ASVs overlapped between the two sets (Fig. 4). This suggests that some members of the oyster microbiome are more responsive to the environment while others are more tightly linked to the host. Understanding how changing environments affect ecologically and economically important species like the Eastern oyster and their associated microbiome is essential for predicting how these organisms will respond to anthropogenic disturbances in this era of global change.