Core Community Persistence Despite Dynamic Spatiotemporal Responses in the Associated Bacterial Communities of Farmed Pacific Oysters

A breakdown in host-bacteria relationships has been associated with the progression of a number of marine diseases and subsequent mortality events. For the Pacific oyster, Crassostrea gigas, summer mortality syndrome (SMS) is one of the biggest constraints to the growth of the sector and is set to expand into temperate systems as ocean temperatures rise. Currently, a lack of understanding of natural spatiotemporal dynamics of the host-bacteria relationship limits our ability to develop microbially based monitoring approaches. Here, we characterised the associated bacterial community of C. gigas, at two Irish oyster farms, unaffected by SMS, over the course of a year. We found C. gigas harboured spatiotemporally variable bacterial communities that were distinct from bacterioplankton in surrounding seawater. Whilst the majority of bacteria-oyster associations were transient and highly variable, we observed clear patterns of stability in the form of a small core consisting of six persistent amplicon sequence variants (ASVs). This core made up a disproportionately large contribution to sample abundance (34 ± 0.14%), despite representing only 0.034% of species richness across the study, and has been associated with healthy oysters in other systems. Overall, our study demonstrates the consistent features of oyster bacterial communities across spatial and temporal scales and provides an ecologically meaningful baseline to track environmental change.


Introduction
Bacterial communities associated with animals and plants play an important role in mediating processes at the individual host and wider community and ecosystem scales [1][2][3][4][5]. In marine systems, the surfaces of host organisms are in direct contact with seawater, providing a constantly changing bacterial community [6]. As such, marine hosts often have thousands of bacteria associations that are extremely responsive across various spatial and temporal scales, which makes them challenging to understand [7]. Recent attempts to simplify this complexity have focussed around identifying signatures of stability through space and time [8][9][10][11]. Here, by identifying core stable taxa, it may be possible to define components of the community that are important for host function, whilst also facilitating comparisons across individuals, populations and ecological contexts [12,13]. Such considerations are important for understanding and interpreting future disturbances and environmental changes [14,15], and there is increasing interest in developing microbial indicators to assess coastal ecosystem health [16].
A range of stressors can disrupt host-bacterial relationships [10,[17][18][19], with wide-ranging implications for host performance, including enhanced susceptibility to pathogens [15]. Bivalve molluscs seem particularly vulnerable to stress-induced disease outbreaks, which often result in mass mortality events that can have serious and far-reaching ecological and economic ramifications [20,21]. For example, the summer mortality syndrome (SMS) in the Pacific oyster, Crassostrea gigas, has become more common in recent decades and is an ever-increasing constraint to the expansion of the aquaculture sector [22,23]. Whilst the underlying mechanisms causing SMS are multifaceted and complex [24], the breakdown of the host-bacteria relationship has been frequently observed [23,[25][26][27]. However, microbial communities associated with oyster hosts are most often sampled when they are already in a state of dysbiosis, usually during a disease event [28]. This limits current ability to develop microbially based monitoring approaches, as our understanding of natural seasonal dynamics (without periods of disease) are generally lacking (but see [29]).
In the face of recent and projected ocean warming trends, a better understanding of host-bacteria relationships is even more pressing, given the increasing likelihood of disruption and detrimental host impacts [19,27,30,31]. For C. gigas, warming can push previously safe areas into prevailing climates associated with SMS. Even in areas already affected, the additional stress burden imposed by more intense summer temperatures now elicits SMS with increasing frequency and intensity and without an obvious aetiological agent [27,28,32]. In the Northeast (NE) Atlantic, the Irish and Celtic Seas represent a transition zone for SMS (~ 19 °C summer SST). Mortality events here are infrequent, but risk increases dramatically with decreasing latitude and increasing temperatures [33]. Therefore, warming trends and associated mortality events may threaten the expansion of the industry over the coming decades [34]. With this in mind, we characterised the host-bacterial relationship at two C. gigas farms in the Irish/Celtic Seas over the course of a year. In doing so, we aimed to (i) establish important baselines against which dysbiosis can be measured and (ii) examine patterns of spatiotemporal variability in bacterial communities to elucidate possible drivers of structure and richness.

Sampling Approach
Sampling took place at two commercial oyster farms in Ireland ( Fig. 1) (Bannow Bay and Cromane) in spring (March/ May), summer (June/July/August) and winter (December) 2018. Increased sampling frequency was conducted during summer to capture shifts in bacterial communities during periods of thermal stress and potential dysbiosis. During each sampling event, eight oysters were randomly taken from trestles at low shore height and three 1-l replicates of seawater were collected in sterile Nalgene bottles. Oysters and seawater were frozen at − 20 °C until DNA extraction.
After defrosting, each oyster was shucked into a 50-ml Falcon tube, mixed with an equal volume (v:v) of sterile artificial seawater and blended using a tissue homogeniser. Defrosted seawater was concentrated by filtering through a 0.22-μm nitrocellulose filter. DNA was extracted from 1 ml of oyster homogenate and whole seawater filters using Qiagen DNeasy PowerSoil extraction kits, following the manufacturer's instructions. DNA was then weighed using a qubit fluorimeter and re-suspended to 2 ng/μl. Library preparation and sequencing of the V4 region of the 16S rDNA gene using primers (515f GTG CCA GCMGCC GCG GTAA + 806r GGA CTA CHVGGG TWT CTAAT) was conducted by using StarSEQ (StarSEQ GmbH, Mainz, Germany) following an optimised protocol of [35]. At least one negative PCR control was run on each plate and demonstrated runs were free from contamination.

Sequence Processing
All processing and analysis was conducted in the R statistical environment. Paired-end reads were processed according to the BIOCONDUCTOR workflow for microbiome data analysis [36]. Sequences were trimmed and truncated using the "filterAndTrim" function in DADA2 with the following parameters: truncLen, f = 240, r = 160; truncQ = 2; and trimLeft, f = 20, r = 19, to remove primers and low-quality reads. Amplicon sequence variants (ASVs) were resolved using DADA2 [36]. Chimeras (0.97% of sequences) were removed using the "removeBimeraDenovo" function in DADA2. Sequence taxonomy was assigned using the RDP naïve Bayesian classifier against the SILVA release 132 database [37] using the "assignTaxonomy" function in DADA2. Sequence read counts, taxonomic assignments and metadata were assembled as an object in the R package ""PHYLOSEQ" and was used in downstream analysis [38]. Samples containing < 10,000 reads, taxa contributing < 0.01% of the reads in the dataset and ASVs identified as mitochondria, chloroplasts or Archaea were then removed from the PHYLOSEQ object. Sequence counts were then expressed as relative abundance (in proportion to the total sample count). Sequences are accessible through the EMBL database (accession no. PRJEB52444). The ASV table and metadata are available at https:// figsh are. com/s/ b36ed 8e187 2f496 d437a.

Statistical Analysis
After sequence processing, three seawater samples had to be discarded from the dataset. Due to this, replication was too low for some of the interaction terms between sample type, month and season. As the focus of this study was to track the shifting oyster bacterial community through time, we made initial comparisons between oyster and seawater samples as a single dataset to determine overall differences between the two sample types. We then based subsequent analyses on differences between sites and months solely on the oyster samples. To account for differences in sequence depth between samples in alpha diversity estimates, the dataset was rarefied to the minimum sample depth (12,430 reads), using the "rar-efy_even_depth" function in PHYLOSEQ. Alpha diversity for each sample was estimated through the Chao1 index [39] implemented through the "estimate_richness" function in PHYLOSEQ. The Chao1 index estimates ASV richness, and the standard error surrounding this estimate, based on the observed number of ASVs, the observed number of ASVs occurring only once and the observed number of ASVs occurring only twice [39]. Alpha diversity was compared using a two-way analysis of variance (ANOVA). Model factors consisted of site (fixed factor; two levels, Bannow, Cromane) and month (fixed factor; six levels, March, May, June, July, August, December). Differences in community structure were determined using PERMANOVA (Anderson, 2001) based on Bray-Curtis dissimilarity and implemented through the "Adonis" function in the package "VEGAN" [40]. This analysis was repeated across all taxonomic levels. Post hoc pairwise comparisons were performed to determine where the differences in community structure lay (at p < 0.05). Model design was the same as that for alpha diversity. Differences in multivariate dispersion between communities were examined using the "betadisper" function in "VEGAN". A similarity of percentage (SIMPER) procedure was conducted to determine which taxa contributed the most to any observed dissimilarities.
We defined core taxa at the ASV level and used a compositional dataset. There is no consistent definition of a "core" in the literature with authors setting prevalence thresholds from 50 to 100%. Here, we used a prevalence threshold of 80% (of the total dataset with all months and sites included).

Alpha Diversity
Overall, alpha diversity (Chao1 index) ranged from 377 ± 15.4 in oysters to 413 ± 9.5 in seawater but there was no significant difference between the two sample types (F (1, 123) = 2.6, p = 0.14). When oysters were analysed separately, alpha diversity differed by month but there was no significant effect of site or the interaction term (Table 1). Post hoc analysis showed that March (561 ± 27.9) and December (519 ± 24) were significantly greater than all other months (May 258 ± 31, June 311 ± 18, July 309 ± 13 and August 312 ± 23) but were similar to one another (Fig. 3). Alpha diversity also differed significantly by month in seawater samples (F (5,27) = 3.27, p = 0.02). Post hoc analysis showed that the only significant difference was observed between March and May (Fig. S1).

Shared ASVs
Out of the 17,533 ASVs recorded across this study, 2014 ASVs (11.5%) were shared between seawater and oysters. The majority of ASVs were found in oysters, which hosted 12,945 ASVs compared to 2574 ASVs recorded solely in seawater. Across the oyster samples, 2189 ASVs (15%) were shared between sites, and these were similar irrespective of the sampling month. Here, shared ASVs between sites ranged from 226 ASVs (7.6%) in August to 314 ASVs

Community Structure
Initial comparisons between oyster and seawater samples showed bacterial communities to be clearly differentiated (Fig. S2), and further analysis focussed solely on oysterassociated communities. PERMDISP showed no significant differences in within-factor multivariate dispersion for either month (F (5, 86) = 2.24, p = 0.06) or site (F (1, 90) = 0.05, p = 0.81). Bacterial community structure exhibited a month × site interaction (Table 1) suggesting the magnitude of difference between sites was not consistent between months or vice versa. This pattern was evident when the dataset was aggregated to coarser taxonomic resolutions (Table S1). Post hoc analysis showed all pairwise comparisons within this interaction term to be significant. Nonmetric multidimensional scaling (nMDS) ordination showed a clear division between sites and differentiation between March and December and all other months (Fig. 4). SIMPER analysis revealed that this was largely due to a shift in the most dominant taxa. In warmer months (May, June, July and August), ASVs belonging to the genus Mycoplasma were far more abundant, whilst AS1 from the family Spirochaetaceae was consistently present in far lower abundance (Table S2). Together, these ASVs accounted for up to 41% of observed dissimilarity between monthly comparisons. nMDS ordination also showed communities at Bannow in May to be more structurally dissimilar than those sampled in other warmer months (June, July and August), a pattern not observed at Cromane. SIMPER analysis revealed that this was largely driven by a dominance (~ 45%) of ASV-7 from the genus Pseudomonas at Bannow during May.

Discussion
The host-bacteria relationship is important for the healthy functioning of benthic organisms, but understanding the complex and dynamic nature of microbial communities remains challenging. For the Pacific oyster, Crassostrea gigas, such understanding is crucial given this species' vulnerability to disease outbreaks and subsequent dysbiosis of this relationship. Here, we characterised the associated bacterial communities at two Irish C. gigas farms, over the course of a year. We found C. gigas harboured spatiotemporally variable bacterial communities that were distinct from bacterioplankton in surrounding seawater. However, despite high variation, we observed clear patterns of stability in the form of a small core component that was persistent across space and time.

Spatial Structuring
We found clear structuring between sites that was evident in every sampling period. This is consistent with the spatial structuring observed in C. gigas [29] and Sydney rock oyster, Saccostrea glomerata [41], in Australian farms and the eastern oyster, Crassostrea virginica, across the east coast of the USA [42]. Whilst we do not have the necessary environmental data to examine potential underlying mechanisms, a range of processes operating over the geographic scale covered here (i.e. ~ 400 km) may be important drivers of variability. These include deterministic factors such as temperature [43], pH [19], nutrient loads [44] and sediment characteristics [45] or neutral processes (e.g. isolation) that facilitate ecological drift [46]. The clear structuring between sites is in contrast to the high within-site variability and small (albeit significant) site-level differences observed between "wild" C. gigas populations in the north and south Wadden Sea [47]. Greater similarity among open coast "wild" Wadden Sea populations, compared to that for our estuarine sites, may be related to greater connectivity, homogenising any selection pressures or drift at that spatial scale. Future studies coupling high-resolution in situ environmental data with microbial community structure are needed to better understand mechanistic drivers of pattern.

Temporal Structuring
The structure of bacterial communities associated with C. gigas varied markedly across sampling periods. In particular, communities sampled in the cooler months (March and December) were distinct and harboured greater diversity compared with those in the warmer months (May, June, July and August). Patterns in alpha diversity mirrored the prevailing patterns of the bacterioplankton (seawater), where greater mixing of the water column and elevated coastal runoff may elevate diversity in winter, as has been recorded in many temperate seas [48][49][50][51]. The clear structuring between warm and cool months is similar to that reported in the oysters, C. virginica [52] and S. glomerata [41]. A large proportion of variation between months was consistently driven by ASV4 (genus Mycoplasma) and ASV1 (family Spirochaetaceae) that were also members of the core component (see below). These taxa displayed clear seasonal dynamics with ASV4 overrepresented in warmer months and ASV1 in cooler months. Differences between months were also driven by a transient component of bacterial communities that was unique to any given sampling month. It is likely that temporal structuring is driven by both stochastic effects associated with passively acquiring ASVs from a dynamic and shifting planktonic community, and deterministic processes imposed by the host on resident/core taxa.

Core Component
Despite clear differences between site and sampling month, there was a small temporally and spatially stable "core", which contributed disproportionately to overall sample abundance. Many of these core taxa have also been observed in C. gigas from other systems [29,47,53]. For example, taxa from the Spirochaetaceae family were found across six estuaries spanning ~ 500 km along Australia's east coast and the genera Mycoplasma and Synechococcus were consistently found across 12 sites in Port Stephens estuary [28,29]. Moreover, they have also been observed in other oyster species around the world [41,44,52,54,55] and may, therefore,

Fig. 5
Relative abundance of the core bacteria community (defined at taxa present in > 80% of samples at a relative abundance > 0.1%) associated with the Pacific oyster, Crassostrea gigas. Abundance is expressed as proportion of entire sample represent part of a consistent core component for oysters generally. Core taxa are hypothesised to be associated with a "healthy" microbiome and may be critically important to the host [9,11]. Whilst the role they play for function of the host remains largely unknown, they are consistently associated with healthy (when compared to diseased) individuals. Lasa et al. (2019) [26] compared healthy C. gigas to those symptomatic with SMS during mortality events in populations across Europe. They found taxa within our persistent core (Mycoplasma, Synechococcus and Spirochaetaceae) to dominate healthy (non-infected) individuals, suggesting a role in oyster health and fitness. Similarly, a decrease in Mycoplasma has been associated with infection of the protozoan Martelia sydneyi in the Sydney rock oyster, S. glomerata [56], and a reduction in the wider class Mollicutes in eastern oysters, C. virginica, infected with the alveolate Perkinsus marinus [55]. However, further studies incorporating other hosts and a greater understanding of the functional profiles of these core taxa is required before the ubiquity and utility of this core can be determined.

Conclusion
In summary, we identified stable and variable features of the host-bacteria relationship of Pacific oysters, which, with their extensive introduced distribution (> 50 countries) and commercial dominance in many regions, are perhaps the world's most globalised bivalve. Microbial communities are increasingly recognised for their role in mediating host resilience to environmental perturbations, and there is increasing interest in developing microbial indicators to assess ecosystem health. Importantly, no mortalities associated with SMS were reported by farms during the study, which means these communities represent "healthy" and "normal" baselines. This represents a crucial first step towards identification of microbial indicators to assess the health of oyster farms. Future studies may build upon this and document how the breakdown of this relationship may impact host condition. This may lead to robust microbial indicators in response to a range of climatic and local stressors.