Skip to main content

Functional plasticity in oyster gut microbiomes along a eutrophication gradient in an urbanized estuary



Oysters in coastal environments are subject to fluctuating environmental conditions that may impact the ecosystem services they provide. Oyster-associated microbiomes are responsible for some of these services, particularly nutrient cycling in benthic habitats. The effects of climate change on host-associated microbiome composition are well-known, but functional changes and how they may impact host physiology and ecosystem functioning are poorly characterized. We investigated how environmental parameters affect oyster-associated microbial community structure and function along a trophic gradient in Narragansett Bay, Rhode Island, USA. Adult eastern oyster, Crassostrea virginica, gut and seawater samples were collected at 5 sites along this estuarine nutrient gradient in August 2017. Samples were analyzed by 16S rRNA gene sequencing to characterize bacterial community structures and metatranscriptomes were sequenced to determine oyster gut microbiome responses to local environments.


There were significant differences in bacterial community structure between the eastern oyster gut and water samples, suggesting selection of certain taxa by the oyster host. Increasing salinity, pH, and dissolved oxygen, and decreasing nitrate, nitrite and phosphate concentrations were observed along the North to South gradient. Transcriptionally active bacterial taxa were similar for the different sites, but expression of oyster-associated microbial genes involved in nutrient (nitrogen and phosphorus) cycling varied throughout the Bay, reflecting the local nutrient regimes and prevailing environmental conditions.


The observed shifts in microbial community composition and function inform how estuarine conditions affect host-associated microbiomes and their ecosystem services. As the effects of estuarine acidification are expected to increase due to the combined effects of eutrophication, coastal pollution, and climate change, it is important to determine relationships between host health, microbial community structure, and environmental conditions in benthic communities.


Coastal ecosystems serve as habitat for highly diverse communities that contribute up to 77% of worldwide ecosystem services [1, 2]. Humans directly rely on these environments for activities like tourism and fisheries, and for their ecosystem services [3,4,5]. Environmental conditions in coastal estuarine ecosystems fluctuate rapidly due to changes in nutrient loading, river runoff, and other physical, chemical, and biological factors [6,7,8]. For example, average pH values in coastal waters can vary by as much as one pH unit over both daily and seasonal cycles, reflecting changes in biological outputs like microbial activity, ambient dissolved oxygen (DO), and pCO2 [9, 10]. These frequent changes in estuarine water chemistry are also affected by human activity and coastal geomorphology, and these influences will likely increase over the coming decades [11].

Estuaries such as Narragansett Bay, Rhode Island, USA, provide a natural gradient to study the impacts of eutrophication. The head of the Bay, located in a highly urbanized area, is highly eutrophic while trophic levels at the mouth are more similar to those found over the continental shelf [12, 13]. Previous studies have shown that eutrophication in the headwaters of Narragansett Bay affects many local communities and ecosystem processes at locations downstream, including nitrification rates [14], primary productivity [15], animal physiology, and benthic biodiversity [16, 17]. Over the last 20 years, Narragansett Bay has undergone dramatic changes as a result of targeted efforts to reduce nutrient loading, providing a dynamic model for the study of gradients in estuarine eutrophication [18].

Marine microbial communities play a central role in ecosystem function as the engine for carbon and nutrient cycling. Microbial communities in coastal seawater and sediment exhibit plastic responses to environmental changes or gradients [19,20,21,22]. This may lead to changes in primary productivity, and therefore coastal ecosystem functioning [23]. Studies of bacterial community structures and nitrogen cycling in several coastal lagoons found that physical gradients and nutrients affect sediment microbial interactions and function [19, 24]. In marine sediments exposed to high nutrients, studies reported dramatic changes in ecological function, but no significant differences in microbial community structure [25,26,27].

Host-associated microbiomes are gaining importance as major contributors to ecosystem services and host functioning [28, 29]. Various studies have found that environmental conditions affect microbial community structures in marine hosts, including corals [30], sponges [31], eelgrass [32], seagrass [33], oysters [34] and mussels [35]. Varying pollution levels alter Manila clam, Ruditapes philippinarum, microbiome composition and host susceptibility to chemicals, as well as kelp bacterial community composition [36, 37]. Studies that examine host-associated microbial functional responses to environmental change, however, are limited to survey studies and focus on model organisms (i.e. corals or zebrafish) in lab-based studies [38]. For example, a study of foundational corals found that nutrients did not affect the host fitness or health, but caused shifts in certain microbial taxa that may influence microbial function [30].

In Narragansett Bay, as in other temperate coastal estuaries, the eastern oyster, Crassostrea virginica, is an integral part of the local history, culture, and seafood industries. Moreover, oysters provide many ecosystem functions, including clearing of overlying waters, coastal erosion prevention, and nutrient cycling [39]. Oyster-associated microbiomes significantly contribute to these ecosystem services, as the oyster host retains and provides a habitat for specific bacteria that perform denitrification and assimilate excess phosphorus [40, 41]. Microbes may also aid in maintaining oyster health and homeostasis by controlling infection, performing nutrient removal, or providing metabolites [42,43,44].

Previous studies of microbial ecology in oysters and other hosts have been limited to surveys of microbial community structures in different compartments of the host. There are, however, very few studies investigating the impact of environmental change on the function of host-associated microbiomes, and how those functional changes may affect coastal ecosystem function. The microbiomes of adult oysters, as determined by 16S rRNA gene amplicon sequencing or other genetic markers, vary with location [45, 46], season [34], tissue type [47], disease status [43, 48, 49], and environmental conditions [50]. Some studies have attempted to infer oyster-associated microbial function from 16S rRNA gene amplicon sequencing, but this method relies on phylogenetically conserved function and is largely speculative [51, 52]. In this study, we evaluated the structure and function of eastern oyster, C. virginica, gut microbiomes at five sites along the eutrophication gradient in Narragansett Bay using 16S rRNA gene amplicon sequencing and metatranscriptomics. This survey provides a snapshot of the oyster microbiomes in a relatively small geographic area in a temperate coastal estuary affected by eutrophication, and how these host-associated microbiomes are affected by their local environment.


Sampled sites showed variability in environmental conditions

Five sites were selected along Narragansett Bay: 1.PVD (Providence River: Bold Point Park), 2.GB (Greenwich Bay: Goddard Memorial State Park), 3.BIS (Bissel Cove: Rome Point), 4.NAR (Narrow River), and 5.NIN (Ninigret Pond) (Fig. 1). These sites are representative of a diversity of environmental conditions (i.e. nutrients, dissolved oxygen, pH, salinity) within a coastal estuary and varying levels of anthropogenic inputs (Table 1). A North-South estuarine gradient was detected, especially in nutrient concentrations. Salinity, pH, and DO increased down the Bay from Providence (1.PVD; North) to Ninigret Pond (5.NIN; South), as coastal eutrophication and the influence of river inputs decreased (Table 1, Spearman’s Correlation Coefficients, SCC = -0.8, − 0.8, − 0.9). Nutrient concentrations decreased along the North-South gradient (Table 1, SCC = 0.7, 0.6, 0.9), with 1.PVD showing significantly higher concentrations of nitrite, nitrate, and phosphate than all other sites (p < 0.01). A Principal Component Analysis (PCA) of these measured environmental conditions (each averaged per site), with eigenvalues in components one and two representing 80% of the environmental variation between sites is shown in Fig. 2. Each site was characterized by a subset of environmental factors over the sampling period. The 1.PVD site was characterized by the highest nutrient levels (nitrite, nitrate, phosphate, p < 0.01, compared to all other sites), 2.GB by the highest chlorophyll-a, 3.BIS by the highest ammonium concentrations (p < 0.001), 4.NAR by a higher temperature (NS) and significantly lower salinity than all other sites (p = 0.045), and 5.NIN by significantly higher pH than all other sites (p = 0.023). The average mass, length, and width of oysters at each site decreased down the Bay, with the exception of oysters from 3.BIS, which were significantly heavier and larger than oysters collected at other sites (Table 1, SCC = 0.7,0.7,0.8; p < 0.001).

Fig. 1

Map of study area with 5 sampling locations. A schematic of the samples collected from each site is shown in the bottom right

Table 1 Summary of all measurements collected per site. Environmental values are daily averages ± standard deviation measured at each site 1 day during week of collection. Nutrient values are averages of three-point water samples collected from each site at time of oyster collection. Spearman’s correlation coefficient (SCC, − 1 to 1) was calculated for the association between each parameter and Latitude. The most significant SCC values (| ≥ 0.8|) are shaded green. A value closer to 1 indicates that the parameter decreases from North-South (1.PVD to 5.NIN) and a value closer to − 1 indicates that the parameter increases from North-South. A correlation coefficient of 0 means there is no linear association and that the value does not consistently change along the estuarine gradient. Significant values as compared to the other sites are indicated in bold. *Spearman’s correlation coefficient for Salinity without 4.NAR: -0.8
Fig. 2

Environmental conditions characterizing each site. Principal component analysis of each site is represented by a colored symbol and each environmental factor is represented with an arrow. Orange arrows indicate average environmental values measured in situ during the sampling week (n = 2); light blue arrows represent nutrient concentrations measured from water samples (n = 3)

Differences in microbial community structures were observed between sites and sample types

A total of 2,217,804 quality-controlled, bacterial 16S rRNA gene sequences were analyzed from 50 gut samples and 10 water samples from 5 sites (Table S1). The sequenced mock community and blank control were analyzed to confirm absence of contamination and adequate sequencing proportions (Fig. S1B). Sequence variant analysis and taxonomic classification resulted in the detection of 304 bacterial Orders across 45 Phyla across all samples, which sufficiently covered the estimated diversity in the samples (Fig. S2). The most dominant phyla in the eastern oyster gut samples, averaged for all oysters at all sites, were Cyanobacteria (38 ± 18%) Proteobacteria (21 ± 13%), Tenericutes (6 ± 12%) and Actinobacteria (3 ± 2%) (Fig. S1A). The most dominant phyla in the water column, averaged from all sites, were Proteobacteria (62 ± 10%), Cyanobacteria (15 ± 12%), Bacteroidetes (15 ± 7%), and Actinobacteria (3 ± 2%) (Fig. S1A). Differences in bacterial community structures were observed between the oyster gut and water samples, in addition to between sites for both sample types (gut and water) (Figs. 3 and 4, and S3).

Fig. 3

Effect of site and sample type on diversity indices in present bacterial community structures. a Simpson’s Index of Diversity calculated using ASV-level 16S rRNA gene amplicons for gut samples (left, n = 10) and water samples (right, n = 2). Global p-values were calculated using the Kruskal-Wallis rank-sum test, and pairwise p-values were calculated with the Wilcox rank-sum test. b NMDS plot visualization of Bray-Curtis beta-diversity (k = 2) at the ASV level for gut samples by Site (left) and all samples by Type (right). The ellipse lines show the 95% confidence interval (standard deviation). p-values indicate significance of grouping with adonis2 Permutational Multivariate Analysis of Variance Using Distance Matrices test. c Within site Bray-Curtis dissimilarity index values for gut samples (n = 90; 9 comparisons for each of 10 samples per site). Global p-value was calculated using the Kruskal-Wallis rank-sum test, and pairwise p-values were calculated with the Wilcox rank-sum test

Fig. 4

Effect of site and sample type on present and active bacterial community structures. a Number of bacterial Orders shared between the water 16S rRNA gene amplicons, oyster gut 16S rRNA gene amplicons, and oyster gut metatranscriptomes (vertical bars). The total number of Orders found in each group is shown in the horizonal bar graph on the left. b Relative percent abundances (square root rescaled) of top 30 bacterial Orders associated with seawater samples (n = 2) or oyster gut tissue (n = 10 or 5), per site. The most abundant bacteria (16S rRNA gene amplicons, middle, n = 10) and the most transcriptionally active bacteria (metatranscriptomes, right, n = 5) in the oyster gut are shown

Effect of sample type on microbial community structures

The structure of the gut microbiome was distinct from the water microbial community, regardless of the sampling site (Fig. 3b, stress = 0.19, df = 4, PERMANOVA R2 = 0.46, p = 0.001). Of the 304 Orders detected in the 16S amplicon data, the water and gut samples had 135 (45%) Orders in common, while 8 (2%) were exclusively found in the water and 161 (53%) were found only in the oyster gut, suggesting selection by the host (Fig. 4a). LEfSe analysis revealed 22 Orders significantly more abundant in the water than the gut samples, including Flavobacteriales, Rhodobacterales, Rhodospirillales, and Oceanospirillales (Figs. 4b and S4A, LDA > 2, p < 0.05). Conversely, Corynebacteriales, Propionibacteriales, Desulfobacterales, and Mycoplasmatales were some of the 14 Orders more abundant in gut samples (Figs. 4b and S4A, LDA > 2, p < 0.05). Significantly more unknown bacterial Orders were detected in the oyster gut samples, compared to the water (Fig. 4b, t-test p < 0.001).

Effect of site on microbial community structures

The oyster gut bacterial communities from each site were significantly different at the ASV level (Fig. 3b, Table S2, df = 4, stress = 0.19, PERMANOVA R2 = 0.15, p = 0.001). Samples from 1.PVD and 3.BIS showed significantly lower alpha-diversity (Fig. 3a, Simpson’s Index; p < 0.01) than samples at other sites. This correlated with specific microbial signatures found at each site. For example, Corynebacteriales and Synechococcales were significantly more abundant in the gut samples from 4.NAR and 5.NIN than at other sites, and Verrumicrobia were significantly more abundant at 4.NAR than in others sites (Figs. 4b and S4B, LDA > 2, p < 0.05). The within-site variability of oyster gut microbiomes was significantly higher at the Northern sites (1.PVD, 2.GB, 3.BIS) than at the Southern sites (4.NAR, 5.NIN) (Fig. 3c, Bray-Curtis; p = 0.001).

Oyster gut samples at 2.GB showed significantly higher percentages of chloroplast-associated 16S rRNA gene amplicons (50 ± 27%), which is consistent with high chlorophyll-a levels measured at this site (Figs. 2 and S4B, LDA > 2, p < 0.05). These sequences were removed prior to calculating dissimilarity metrics. Oyster gut samples from all sites shared 105 Orders (34% of 304 total), while 9–31 (3–10%) Orders were distinct to gut samples at certain sites (Fig. S3).

Transcriptionally active microbial community structures differ from ASV community structures

A total of 409 million metatranscriptomic 150 bp-long, quality-controlled paired-end reads were obtained from 25 gut samples (n = 5 per site; Table S1), which sufficiently covered approximate 97.6 ± 0.4% of the diversity in the samples (Fig. S2). Direct taxonomic annotation of these merged paired-end reads classified 2.35 ± 0.02% microbial reads (Table S1). This level of annotation is comparable with other studies in host-associated systems and most probably due to incomplete taxonomic coverage in reference databases and high levels of host nucleic acids [24, 33, 53]. Gene classification using marker ribosomal genes was not possible due to the rRNA depletion performed during library prep and subsequent biased removal of these common taxonomy markers [54]. Other commonly used methods based on house-keeping genes are not useful with environmental samples as the ones in this study, most probably due to the incomplete levels of annotation of these marker genes in these sample types [55]. Sixty-eight bacterial Orders across 29 Phyla were detected in the taxonomically annotated reads, of which 36 (53%) were also detected in the gut 16S amplicon data. The most active annotated phyla in the gut samples (all oysters) were Proteobacteria (46 ± 5%) and Firmicutes (16 ± 8%) (Fig. 4b). The most active taxa (Bacillales, Pseudomonadales, and Enterobacterales; as detected in the metatranscriptomes) were not the most abundant taxa (Cyanobaceria, Mycoplasmatales, and Unknown Proteobacteria; as detected by 16S rRNA gene amplicon analysis) (Fig. 4b). There were 103 (out of 296, 35%) Orders detected in the oyster gut 16S amplicons that were not detected in the metranscriptomes, most likely due to methodological biases (Figs. S2 and Fig. 4a).

While microbiome structures of oyster gut samples (detected by 16S rRNA gene amplicon analysis) showed differences between sites (Fig. 3b), microbiome structures of the active taxa (determined by taxonomic annotation of metatranscriptomic reads) in the gut samples were not different between sites (Fig. S5, species level). This may be due to the higher number of species detected using the metatranscriptomic approach as compared to the more conserved 16S rRNA gene amplicon sequencing, as seen by comparing the “core microbiome” detected with each method (taxa occurring in > 80% of samples, Fig. S6 and Table S3). Metatranscriptomic analysis detected a higher number of conserved species, as compared to the 16S rRNA gene amplicon analysis (72 conserved taxa in metatranscriptomes compared to 15 16S rRNA gene ASVs). Of those, only 2 taxa were identified by both methods, including a Propionibacteriaceae and a Synechococcus sp. The largest portion of the conserved sequences in both datasets corresponded to the Gammaproteobacteria, with 21 and 3 annotated taxa in the metatranscriptomic and 16S rRNA gene amplicon datasets respectively, followed by Actinobacteria (7 and 2 respectively). Gammaproteobacters conserved in oysters from Rhode Island waters included several Aeromonads, Enterobacters, Pseudomonads, and Vibrionales as identified by the metatranscriptomic dataset. Other conserved taxa uniquely detected by the metatranscriptomic dataset included annotations to Acidobacteria (1), Bacteroidetes (2), Chlamydia (1), Firmicutes (12 bacilllales and 5 clostridiales), Fusobacteria (1), Spirochaetaceae (1), and Verrumicrobiales (1). On the other hand, the 16S rRNA gene analysis detected one conserved ASV in the mollicutes (Fig. S6 and Table S3).

Transcriptional responses in the oyster gut microbial community reflect the estuarine gradient in Narragansett Bay

Although no significant differences were detected between sites on the taxonomy of the transcriptionally active microbial taxa in the gut tissue (Fig. S5), their transcriptional responses varied based on the environmental conditions at each site (Fig. 5). In order to increase statistical power in the characterization of environment (i.e. eutrophication) on gene expression, the microbial transcriptional response at the more eutrophic/urbanized northern sites (1–3) was compared to that of the southern sites (4–5; considered as the less urbanized “control group”). This resulted in 11 SEED Level-1 pathways that showed significantly differential gene expression between northern and southern sites. A significant upregulation of bacterial stress responses and general metabolic activities (carbohydrates, respiration, amino acids, fatty acids, lipids etc.) was seen at northern sites, as well as a downregulation of photosynthesis, metabolic transport, and motility and chemotaxis (Benjamini-Hochberg adjusted p < 0.05; Fig. 5a). Select pathways were then further analyzed by comparing gene expression at each site to the mean of all sites.

Fig. 5

Pathways showing significant differential gene expression patterns between sites. a Differential expression of all significant (Benjamini-Hochberg padj < 0.05) Level 1 pathways at the Northern sites (1–3, n = 15), compared to Southern sites (4–5, n = 10). A red bar (fold change> 0) indicates upregulation in the North and a blue bar (fold change< 0) indicates downregulation in the North. b Differential expression of Level 2 Stress Response pathways at each site, relative to the mean of all other sites (n = 5, Benjamini-Hochberg *padj < 0.05, **padj < 0.01)

Stress response

In order to further examine the effects of anthropogenic factors (e.g. eutrophication, urbanization) on oyster microbial community and function, a more in-depth analysis of differences in the expression of genes involved in stress responses and nutrient cycling was performed (SEED level 2 annotation). Differential expression of genes in stress response and nutrient pathways at each of the sites was compared to the mean level of expression at all other sites (Figs. 5 and 6b). A significant upregulation in the expression of microbial genes involved in dealing with osmotic stress was detected in samples from 2.GB (as compared to the mean of all sites), as well as a significant downregulation in genes involved in periplasmic stress (p < 0.05). Conversely, a significant upregulation in genes involved in periplasmic stress (e.g. rseA, degS, deQ) and downregulation in genes involved in osmotic stress (e.g. genes coding for betaine aldehyde dehydrogenase and choline dehydrogenase) and oxidative stress (e.g. genes coding for NAD G3P dehydrogenase) was detected in oyster gut samples from 5.NIN (p < 0.05, Fig. 5b, S7, and S8). Expression of microbial genes involved in other acute stress responses, including acid stress, cold shock, and heat shock, were not significantly different between sites.

Fig. 6

Differential gene expression in nitrogen and phosphorus metabolism at each site. a Differential expression of Nitrogen pathways at all sites, relative to the mean of the others (n = 5, Significance: Benjamini-Hochberg *padj < 0.05, **padj < 0.01). (top) Total differential expression of overall nitrogen metabolism, indicated with the blue-green colors. (bottom) Relative log fold change in nitrogen metabolism pathways, indicated with yellow-red colors. b Differential expression of Phosphorus pathways at all sites, relative to the mean. (top) Total differential expression of overall phosphorus metabolism, indicated with the blue-green colors. (bottom) Relative log fold change in phosphorus metabolism pathways, indicated with purple colors

Nitrogen metabolism

Nutrient cycling is central to ecosystem services provided by oysters. Nitrogen and phosphorus are especially important, since they are the major components of eutrophication and often limiting factors to primary production [11, 56]. No significant changes in expression of genes involved in overall nitrogen metabolism (SEED level 2 annotation) was observed in the gut oyster microbiome from the different sites (Fig. 6a top), despite the significant differences in levels of nitrate, nitrite, and ammonium levels detected between sites (Table 1). High levels of variability in the expression of genes involved in the different pathways involved in nitrogen metabolism were observed between oysters within sites. Significant differences between sites were observed, however, in the patterns of expression of genes from specific pathways involved in nitrogen metabolism (SEED level 1 annotation; Fig. 6a bottom), reflecting differences in the responses of oyster gut microbes to the environmental conditions at each site. At the northernmost site (1.PVD), there was a significant downregulation of denitrification and NO detoxification-related genes (e.g. nosF and cytochrome c-dependent nitric oxide reductase (cNor)) compared to the mean of all sites, while at the southernmost site (5.NIN), a significant upregulation of genes involved in ammonia pathways (e.g. genes coding for NR(I), GNLE, and nitrate reductase) and a downregulation of nitrilase genes was observed (p < 0.05, Fig. S9).

Phosphorus metabolism

Expression of genes in the oyster gut microbiome involved in phosphorus metabolism decreased down the Bay, with microbial communities in the guts of oysters from the most southern (5.NIN) and northern (1.PVD) sites respectively showing significantly lower and higher levels of expression of genes involved in phosphorus metabolism than the mean of the sites (p < 0.01; Fig. 6b top). An upregulation of genes involved in the phosphate pathway (e.g. alkaline phosphatase) was observed in the gut microbiome of oysters from the northernmost site (1.PVD), as well as an upregulation of genes in the phosphonate pathway in oysters from 2.GB (e.g. phosphonoacetaldehyde hydrolase) (Fig. S6b bottom and S10). These two sites also showed the highest ambient phosphate concentrations (Table 1). Conversely, there was a significant downregulation of phosphate and phosphonate pathways at the southernmost site compared to the mean of all sites (5.NIN, p < 0.01, Fig. S6b bottom).


A better understanding of the effect of environmental conditions on both the structure and function of oyster-associated microbes is important for the management of oyster populations and optimization of the ecosystem services they provide. This study characterized the composition and function of eastern oyster-associated microbiomes at sites within a temperate, urbanized estuary. Oyster gut microbiomes during the summer were diverse in composition and differed between sites. Differences between the structure of microbiomes between water and oyster gut were consistent with selection and amplification of taxa from the water environment by the oyster host, as described in previous work [50]. Significant differences in expression of several gene pathways (stress response, nutrient utilization) were observed between sites, reflecting the environment at each of the sites. In particular, the gut microbial community of oysters collected at the northern sites, experiencing higher levels of nutrients and anoxia, showed upregulation of genes associated with stress response and phosphorus metabolism, as compared to the less eutrophic southern sites. Microbes in the gut of oysters from the less eutrophic sites showed a relative upregulation of genes associated with nitrogen metabolism. These responses varied according to the eutrophication gradient, indicating that the responses of oyster gut-associated microbiomes reflect the local environment, despite the fact that they are located within the host (i.e. the oxygen and nutrient status of the water is pervasive in the oyster gut). This is the first study combining 16S rRNA gene amplicon and metatranscriptomic data to look at both the composition and function of microbes associated with bivalves in an estuarine eutrophication gradient, with findings that are relevant to the development of restoration projects geared to maximize ecosystem services provided by oysters.

Surprisingly, expression of certain pathways involved in nitrogen metabolism in oyster-associated microbiomes was significantly higher at sites with the lowest levels of nutrients (NO2, NO3, NH4+) in the water at the southern range of the estuarine gradient than at the more eutrophic northern sites. These results are consistent with previous findings showing that oxygen conditions control nitrogen and phosphorus cycling in the sediments by limiting nutrient availability [57], with high dissolved oxygen concentrations in water promoting nitrogen removal [58]. This interaction between oxygen concentration (or redox state) and nitrogen metabolism has been well-documented in marine sediments: higher DO and low NO3 concentrations stimulate denitrification, while the opposite occurs with high NO3 concentrations [59, 60]. Our findings indicate that the environment in the oyster gut, as it relates to N and P cycling, reflects the overall environmental conditions at the site, consistent with expectations from sediment and water column observations. In the more eutrophic, anoxic, and acidic waters of the Providence River, where there is no available ammonia, oyster-associated microbiomes would upregulate pathways using the phosphate available from the sediment [61, 62]. Oysters in these northern sites may have been enriched in microbes adapted to prefer phosphate over nitrogen, due to selection driven by prevailing environmental conditions. Alternatively, higher levels of expression of genes involved in ammonia and nitrate utilization at the southern sites, in which less eutrophic conditions were found for nitrite, could reflect a physiological compensation driven by a need to scavenge the lower amounts of nitrate at these sites [63]. This is, however, not observed for phosphate. The potential similarity in function of microbes involved in N and P cycling between oyster gut and seawater/sediment may be explained by unique aspects of the physiology and anatomy of bivalves like oysters. Bivalves are ectotherms with an open circulatory system, and perform high levels of filtration when actively feeding [64, 65]. Therefore, we hypothesize that the particular nutrients (N, P compounds) and conditions (DO, pH) driving N and P cycling, must be, overall, similar in the oyster gut to the ambient water environment.

Microbiomes in the gut of eastern oysters collected at each of the sites also reflected potential stressors at each of the sites. For example, the increase in microbial oxidative stress observed at 1.PVD and 2.GB has been widely observed in microbial communities in response to anoxia, pollution, and toxins [66, 67]. The upregulation of periplasmic stress response (due to stressors within the inner bacterial membrane) observed in samples collected at site 5.NIN is likely coupled with increased nitrogen metabolism and transport [68,69,70]. In general, as eutrophic conditions worsen, bacteria will expend more energy on stress response and metabolic activities, a trend that has also been shown in marine sediment microbiomes [71, 72]. Other stressors, including pathogens, toxins, or chemical pollutants may have contributed to the differential expression in stress response pathway between sites and require further study.

Comparison of the microbial composition between water and oyster samples suggest that oysters select and amplify certain bacterial species through feeding selection, niche colonization, and/or evasion of immunity, as shown in oysters and other invertebrate host species [73, 74]. The fact that bacterial composition in gut samples does not completely reflect the community in the water samples may indicate that oysters amplify rare members in the water community and/or retain bacteria previously acquired through time horizontally from the water or vertically from parents. Consistent with the hypothesis of amplification, certain bacterial taxa that are known intracellular anaerobes were relatively more abundant in gut than water samples. These include members of the Mycoplasmatales, Mollicutes, and Clostridiales. Mycoplasmatales have been identified as common invertebrate symbionts and are avid biofilm-formers, allowing them to survive and replicate in the host [75, 76]. Besides these intracellular bacteria, Proteobacteria formed the most abundant and active phylum in the overall community, consistent with published literature in oysters [34, 48, 77]. High variability in the relative abundances of certain taxa (i.e. Mycoplasmatales or Caulobacterales) among oysters within sites suggests that host acquisition of bacteria from surrounding waters is shaped not only through exposure during feeding, but also factors like host health, characteristics, origin, and/or genetics [42, 78]. Further examination of within-site variability and its relationship with other host parameters (e.g. health and physiological status, genetics) may reveal how certain taxa are promoted in each oyster [79]. Decreased microbial diversity has been associated with health-compromised hosts, which may limit their plasticity and ability to respond to environmental change [48, 80]. Conversely, microbiomes of unhealthy hosts may be associated with an increase in diversity in tissues other than the usually more microbially diverse gut, as opportunistic and pathogenic bacteria proliferate in these tissues when compromised by disease [81, 82]. The interplay between the environment, oyster-associated microbiomes, and host health will be the focus of further study.

The differences in microbial composition observed in oyster gut samples between the 16S rRNA gene and the metatranscriptomics data were not unexpected, considering that metatranscriptomes reflect actively transcribing taxa, while the 16S rRNA gene amplicon method detects both live and recently dead microbes [83]. Although a small percentage of the metatranscriptomic reads were annotated, this method was able to capture the annotated microbial diversity at the species level. There are other potential technical biases that could have affected differences in microbial composition as determined by these two different methods, including PCR bias introduced in the process of amplification of the 16S rRNA gene portion, RNA degradation during storage, and issues related to uneven annotation between taxa [84, 85]. Although we have tried to minimize the impacts of these biases through the inclusion of several quality controls commonly used in microbiome research (e.g. mock communities, the use of technical and biological replicates and complementary analyses methods) [86], these biases should also be considered in the interpretation of these results. Despite technical challenges, comparisons between the 16S rRNA gene amplicon data with metatranscriptomic analysis of the oyster-associated microbial community, however, may provide some initial insights into identification of which of the core microbes show a symbiotic relationship with the oyster host, versus those that are transient food in the gut (i.e. accumulate in the oyster gut through association with food selectively ingested by oysters through filter feeding) [50, 87]. In particular, of the selected taxa shown to be relatively more abundant in the gut samples as compared to the water, a subset was detected to be transcriptionally active, particularly Bacillales, and Vibrionales, suggesting that these bacteria are not immediately digested in the gut or eliminated by the immune system. These results are consistent with the fact that these taxa are commonly cultured from oyster samples [48, 50, 88, 89]. Conversely, despite the high relative abundance of Synechococcales and other Cyanobacteria detected both in oyster gut samples and water from some of the sites through 16S rRNA gene amplicon sequencing, the gut metatranscriptomes did not show a relative enrichment in levels of expression of genes involved in recent photosynthesis, suggesting this higher abundance was transient and a reflection of recent feeding activity.


This study has implications for quantification of ecosystem services provided by eastern oyster restoration and aquaculture. In Narragansett Bay, oyster fisheries were a dominant industry in the late 1880s, but a combination of pollution, overfishing, and dredging lead to the collapse of oyster populations in the 1940s [90]. In recent years, numerous efforts have been made to renew oyster reefs and restore their ecosystem services in Narragansett Bay. A common goal of oyster restoration projects is improvement of water quality by stimulation of environmental denitrification [39, 40]. Our findings support that removal of bioavailable nitrogen by denitrification, an important ecosystem service provided by oysters, declines in low oxygen, nutrient rich environments [62, 91, 92]. Enhanced denitrification would occur at high dissolved oxygen and nutrient rich environments, such as the conditions observed at 4.NAR during the summer. This implies that if the environmental microbial community does not have the genes necessary for the nitrogen pathway and/or the environmental conditions do not favor the process, then the addition of oysters to the site will not promote the ecosystem service. The prevailing environmental conditions and function of the resident environmental microbial community should be considered when selecting sites for oyster farming and restoration. In this study, 4.NAR and 5.NIN would provide the greatest return on investment for a restoration project, if only the benefits of denitrification are considered.

In summary, the estuarine gradient affected eastern oyster-gut associated microbial communities through changes in community composition, microbial stress responses, and microbial nutrient utilization. Combined, these results have implications for environmentally-driven changes in oyster microbial acclimation and potential ecosystem services. As the effects of estuarine acidification are expected to increase due to the combined effects of eutrophication, coastal pollution, and climate change, it is important to determine relationships between host health, microbial community structure, and environmental conditions. The results presented here form a baseline for future studies exploring how human-driven estuarine acidification affect overall oyster health and their associated ecosysem services.


Sample collection

Five sites were selected along Narragansett Bay, Rhode Island, USA and wild eastern oysters, C. virginica, were collected from the northern 4 sites, and farmed oysters were collected from 5.NIN (no wild oysters were found). Environmental data for temperature, pH, DO, salinity, and chlorophyll-a were collected using a YSI 6 Series Multiparameter Water Quality Sonde (Model 6920VS) every 30 s for 15 min at each site during the morning and afternoon hours on 1 day of the week of sampling, and averaged to account for diel variation in these parameters [93]. Sample collections were completed at low tide at each site from August 17–25, 2017 and consisted of oyster and water samples at each of the 5 sites with scientific collector’s permit #212 granted by the RI Department of Environmental Management. A total of 150 oysters (30 per site) were randomly collected within a 10 m2 transect at each of the 5 sites. On the day of collection, whole oysters (with shell) were weighed, shell width and length was measured, and samples of gut tissues (around 300 mg) were dissected and immediately preserved in RNAlater (Invitrogen) for RNA/DNA extractions. Preserved tissue samples were stored at − 80 °C until nucleic acid extractions. Up to 1 L of seawater was collected from the location and depth at which the oysters were collected in each site. Duplicate seawater samples were filtered using a peristaltic pump onto a 0.22 μm Sterivex filter (Millipore Sigma), filled with 2 mL of RNAlater, and then stored at − 80 °C until DNA extraction. An additional sample of seawater (30 mL) was filtered through a 0.22 μm syringe-top Polyether-sulfone (PES) filter and frozen at − 80 °C for nutrient analyses. Nutrient concentrations (nitrite, nitrate, ammonium, and phosphate) were analyzed in triplicate using a Lachat QuickChem QC8500 automated ion analyzer operated by the University of Rhode Island Marine Sciences Research Facility.

Gut DNA and RNA extraction

Total nucleic acids were extracted from 150 to 200 mg of gut tissue (n = 10 random oysters per site; 50 total) using the Qiagen Allprep PowerViral DNA/RNA extraction kit with modifications as follows. The gut tissue sample was added directly to a 0.1 mm glass bead tube (Qiagen), along with 600 μL of Solution PV1 and 6 μL of sterile β-mercaptoethanol to minimize RNA degradation. The samples were subjected to bead beating for 5 min, followed by proteinase K digestion at 55 °C for 1 h in a shaker at 300 rpm. The supernatant was transferred to a new microcentrifuge tube and the protocol continued according to the manufacturer’s recommendations. Following nucleic acid extraction, the concentration was quantified using a Nanodrop 2000 instrument (ThermoFisher).

RNA purification from a 5 μL total nucleic acid aliquot was performed using the DNase Max I kit (Qiagen) according to the manufacturer’s protocol in a 50 μL reaction volume for 5 oyster gut samples per site. DNA purification of a 30 μL total nucleic acid aliquot was performed using an adapted version of the DNeasy PowerLyzer PowerSoil Kit. In brief, the total nucleic acids were transferred to a new 2 mL microcentrifuge tube, 1200 μL of Solution C4 was added, then vortexed. Next, 4 μL of RNase A solution was added to the sample and incubated for 2 min at room temperature. The treated DNA was loaded to a spin column, washed with Solution C5, and eluted in 50 μL of Solution C6. DNA and RNA concentrations were quantified with both a Nanodrop 2000 instrument (ThermoFisher) and Qubit Fluorometer High-Sensitivity reagents (Invitrogen).

Seawater DNA extraction

Total DNA from water samples was extracted from the duplicate Sterivex filters using the Qiagen Allprep PowerViral DNA/RNA and DNeasy PowerLyzer PowerSoil kits with modifications as follows. The RNAlater was flushed out of the filters using a sterile syringe, and filters were rinsed with 2 mL of 1X sterile nuclease-free Phosphate Buffer Saline (PBS, pH 7.4, Invitrogen). Solution PV1 (1800 μL) and β-mercaptoethanol (18 μL) were added to the filter cartridge and incubated at 37 °C for 30 min. Next, 20 μL of proteinase K was added to the filter and digested at 55 °C for 1 h. The supernatant was flushed from the filter, and the protocol continued according to the manufacturer’s recommendations. DNA was purified from the entire total nucleic acid product and quantified using the methods described above.

Nucleic acid amplification and sequencing

In order to obtain a comprehensive representation of the gut microbial community and their activities, 2 types of sequencing were performed: 16S rRNA gene amplicon of the V6 region (DNA, a measure of overall composition) and whole shotgun metatranscriptomes (RNA, a snapshot of functional activity at the time of collection) [94]. Amplicons of the V6 region of the 16S rRNA gene in the 50 gut DNA samples (10 per site), 5 water samples, a sample of a mock community (Zymo Research DNA standard I), and blank PCR control were prepared using 967F/1064R primers. A two-step PCR reaction using 300 ng of gut DNA or 10 ng of water DNA was performed in triplicate 33 μL reactions as previously described [95, 96]. The PCR products were analyzed with 75 bp paired-end sequencing to obtain overlapping reads on an Illumina MiSeq at the Genomics and Sequencing Center at University of Rhode Island.

The metatranscriptomic libraries were prepared from 2 μg of gut RNA (n = 5 randomly selected per site), fragmented at 500 nt using Covaris ultrasonification, and treated with the Illumina Ribo-Zero Gold rRNA Removal Epidemiology Kit prior to library prep to remove both host and bacterial rRNA. Illumina TruSeq PCR-free library kits were used to prepare the libraries, and then verified using both KAPA library quantification kits and Agilent Bioanalyzer. The resulting metatranscriptomic libraries were sequenced on an Illumina NovaSeq S4 to obtain 2 × 150 bp paired-end reads at the Yale Center for Genome Analysis.

Processing and analysis of sequencing data

16S rRNA gene amplicon sequences were demultiplexed and quality filtered using DADA2 (v1.6.0) implemented in QIIME2 (v2018.4.0) with additional parameters --p-trunc-len-r 65 --p-trunc-len-f 76 --p-trim-left-r 19 --p-trim-left-f 19 to determine analysis sequence variants (ASVs) [97, 98]. All ASVs were summarized with the QIIME2 pipeline (v2018.4.0) and classified directly using the SILVA database (99% similarity, release #132) [99, 100]. Processed ASV and associated taxonomy data was exported as a count matrix for analysis in R (v3.4.1). Significant differences in Orders between sample types or gut samples at each site were calculated using linear discriminant analysis effect site (LEfSe) implemented on the Galaxy server [101, 102]. Non-bacterial and chloroplast sequences were then removed, and the data was normalized by percentage to the total ASVs in each sample for further dissimilarity metric analysis.

All descriptive and statistical analyses were performed in the R statistical computing environment with the vegan v2.5.5 and phyloseq v1.28.0 packages [103, 104]. Rarefaction curves and sequencing coverage estimates were generated using the rarecurve() and rareslope() commands with sample = [number of reads in smallest sample] in vegan v2.5.5 [105]. Simpson’s diversity values were calculated for each sample at the ASV level using the vegan package and analyzed using the non-parametric Kruskal–Wallis rank sum test in R. Non-metric dimensional analysis (NMDS) was used to determine the influence of sample type or field site on the ASV-level composition, implemented using vegan. The Bray-Curtis dissimilarity metric was calculated with k = 2 for max 50 iterations and 95% confidence intervals (standard deviation) were plotted. Statistical testing of beta-diversity was done using the PERMANOVA adonis2 test implemented in vegan (method = “bray”, k = 2) [106, 107]. Within-site variability was calculated using the command vegdist (method = “bray”, k = 2) and the matrix was simplified to include samples compared within each site. Additional visualizations were computed using the ComplexHeatmap v3.9 and UpSetR v1.4.0 packages [108, 109].

Raw reads from the microbial community metatranscriptomes were first quality controlled with Trimmomatic software v0.36 with parameters PE -phred33 SLIDINGWINDOW:4:15 MINLEN:70 [110]. Metatranscriptomic analysis was performed using scripts from the SAMSA2 pipeline [111], with the following modifications. The quality-controlled paired-end reads were combined using PEAR v0.9.10 and then rogue rRNA reads were removed from the merged reads using SortMeRNA v2.1 [112, 113]. Taxonomic and functional annotation of the data were performed against RefSeq and SEED Subsystem databases, respectively, using DIAMOND v0.9.23 [114]. Core microbiomes were calculated using a custom script in R to determine which taxa (16S rRNA amplicons: ASVs; metatranscriptomes: species) were present in > 80% of the samples in each dataset.

Custom scripts using DESeq2 v1.14.1 were used to calculate differential expression between sites using the command DESeqDataSetFromMatrix(), with values transformed to account for differences in read abundance between samples [115]. The resulting changes in expression at all pathway levels were exported to R for analysis and visualization using ggplot2 v3.2.1 and cowplot v1.0.0 [116,117,118]. Significant differences in gene expression are reported using Benjamini-Hochberg adjusted p-values. NMDS analysis was used to describe differences in gene expression between sites, and was calculated at the species and order level as described above. All processed sequencing files, bash scripts, QIIME2 artifacts, and Rscripts to reproduce the figures in the manuscript are available on Zenodo [119].

Environmental statistical analysis

All statistical analyses of environmental and sequencing data were performed in R (v3.4.1 R Development Core Team, 2011) as follows. The environmental principal component analysis (PCA) was calculated using the prcomp (scale = TRUE) command implemented in base stats v3.6.1, and then plotted using autoplot() enabled by ggfortify v0.4.7 [120]. Significant differences in environmental parameters between sites were determined using all raw data subset by site and parameter. The data per site was compared using the compare_means() command from the ggpubr v0.2.2 package [121]. The method for each comparison was defined as “anova” for initial testing, then “t.test” for pairwise comparisons. Adjusted p-values were calculated using the Benjamini-Hochberg method by adding “p.adjust.method = BH” to the command [115].

Availability of data and materials

The raw sequences generated for this study can be found in the NCBI Sequence Read Archive under BioProject no. PRJNA598635. Table S4 lists the corresponding BioSample and accession numbers for each individual sample. All processed sequencing files, bash scripts, QIIME2 artifacts, and Rscripts to reproduce the figures in the manuscript are available in the Zenodo repository, [119].



Providence River, Bold Point Park


Greenwich Bay, Goddard Memorial State Park


Bissel cove, Rome Point


Narrow River


Ninigret Pond


Amplicon sequence variant




Dissolved oxygen


Linear discriminant analysis


Linear Discriminate Analysis Effect Size


Phosphate buffer saline


Principal component analysis


Polymerase chain reaction


Quantitative Insights Into Microbial Ecology


Reference Sequence Database


Spearman’s correlation coefficient


  1. 1.

    Costanza R, D’Arge R, De Groot R, Farber S, Grasso M, Hannon B, et al. The value of the world’s ecosystem services and natural capital. Nature. 1997;387:253–60.

    CAS  Article  Google Scholar 

  2. 2.

    Martínez ML, Intralawan A, Vázquez G, Pérez-Maqueo O, Sutton P, Landgrave R. The coasts of our world: ecological, economic and social importance. Ecol Econ. 2007;63:254–72.

    Article  Google Scholar 

  3. 3.

    Bulleri F, Chapman MG. The introduction of coastal infrastructure as a driver of change in marine environments. J Appl Ecol. 2010;47:26–35.

    Article  Google Scholar 

  4. 4.

    Liquete C, Piroddi C, Drakou EG, Gurney L, Katsanevakis S, Charef A, et al. Current status and future prospects for the assessment of marine and coastal ecosystem services: a systematic review. PLoS One. 2013;8:e67737.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  5. 5.

    Firth LB, Knights AM, Bridger D, Evans AJ, Mieszkowska N, Moore PJ, et al. Ocean sprawl: Challenges and opportunities for biodiversity management in a changing world. In: Oceanography and Marine Biology: An Annual Review; 2016. p. 193–269.

    Google Scholar 

  6. 6.

    Nixon SW. Coastal marine eutrophication: a definition, social causes, and future concerns. Ophelia. 1995;41:199–219.

    Article  Google Scholar 

  7. 7.

    Sunda WG, Cai WJ. Eutrophication induced CO2-acidification of subsurface coastal waters: interactive effects of temperature, salinity, and atmospheric P CO2. Environ Sci Technol. 2012;46:10651–9.

    CAS  Article  Google Scholar 

  8. 8.

    Waldbusser GG, Salisbury JE. Ocean acidification in the coastal zone from an Organism’s perspective: multiple system parameters, frequency domains, and habitats. Annu Rev Mar Sci. 2014.

  9. 9.

    Alexandre A, Silva J, Buapet P, Björk M, Santos R. Effects of CO2 enrichment on photosynthesis, growth, and nitrogen metabolism of the seagrass zostera noltii. Ecol Evol. 2012;2:2625–35.

    Article  PubMed  PubMed Central  Google Scholar 

  10. 10.

    Baumann H, Wallace RB, Tagliaferri T, Gobler CJ. Large natural pH, CO2 and O2 fluctuations in a temperate tidal salt marsh on Diel, seasonal, and interannual time scales. Estuar Coasts. 2014;38:220–31.

    Article  Google Scholar 

  11. 11.

    Wallace RB, Baumann H, Grear JS, Aller RC, Gobler CJ. Coastal Ocean acidification: the other eutrophication problem. Estuar Coast Shelf Sci. 2014;148:1–13.

    CAS  Article  Google Scholar 

  12. 12.

    Oczkowski A, Nixon S, Henry K, DiMilla P, Pilson M, Granger S, et al. Distribution and trophic importance of anthropogenic nitrogen in Narragansett Bay: an assessment using stable isotopes. Estuar Coasts. 2008;31:53–69.

    CAS  Article  Google Scholar 

  13. 13.

    Calabretta CJ, Oviatt CA. The response of benthic macrofauna to anthropogenic stress in Narragansett Bay, Rhode Island: a review of human stressors and assessment of community conditions. Mar Pollut Bull. 2008;56:1680–95.

    CAS  Article  PubMed  Google Scholar 

  14. 14.

    Berounsky VM, Nixon SW. Rates of nitrification along an estuarine gradient in Narragansett Bay. Estuaries. 2006;16:718.

    Article  Google Scholar 

  15. 15.

    Oviatt CA. Impacts of nutrients on Narragansett Bay productivity: a gradient approach. In: Science for ecosystem-based management. New York: Springer New York; 2008. p. 523–43.

    Google Scholar 

  16. 16.

    Hale SS, Cicchetti G, Deacutis CF. Eutrophication and hypoxia diminish ecosystem functions of benthic communities in a New England estuary. Front Mar Sci. 2016;3:249.

    Article  Google Scholar 

  17. 17.

    Pelletier M, Ho K, Cantwell M, Perron M, Rocha K, Burgess RM, et al. Diagnosis of potential stressors adversely affecting benthic invertebrate communities in Greenwich Bay, Rhode Island, USA. Environ Toxicol Chem. 2017;36:449–62.

    CAS  Article  PubMed  Google Scholar 

  18. 18.

    Oviatt C, Smith L, Krumholz J, Coupland C, Stoffel H, Keller A, et al. Managed nutrient reduction impacts on nutrient concentrations, water clarity, primary production, and hypoxia in a north temperate estuary. Estuar Coast Shelf Sci. 2017;199:25–34.

    CAS  Article  Google Scholar 

  19. 19.

    Highton MP, Roosa S, Crawshaw J, Schallenberg M, Morales SE. Physical factors correlate to microbial community structure and nitrogen cycling gene abundance in a nitrate fed eutrophic lagoon. Front Microbiol. 2016;1691.

  20. 20.

    Meyer J, Riebesell U. Reviews and syntheses: responses of coccolithophores to ocean acidification: a meta-analysis. Biogeosciences. 2015;12:1671–82.

    Article  Google Scholar 

  21. 21.

    Paerl HW, Dyble J, Twomey L, Pinckney JL, Nelson J, Kerkhof L. Characterizing man-made and natural modifications of microbial diversity and activity in coastal ecosystems. Antonie van Leeuwenhoek. Int J Gen Mol Microbiol. 2002;81:487–507.

    CAS  Article  Google Scholar 

  22. 22.

    Nogales B, Lanfranconi MP, Piña-Villalonga JM, Bosch R. Anthropogenic perturbations in marine microbial communities. FEMS Microbiol Rev. 2011;35:275–98.

    CAS  Article  PubMed  Google Scholar 

  23. 23.

    Paerl HW, Dyble J, Moisander PH, Noble RT, Piehler MF, Pinckney JL, et al. Microbial indicators of aquatic ecosystem change: current applications to eutrophication studies. In: FEMS Microbiology Ecology; 2003. p. 233–46.

    Google Scholar 

  24. 24.

    Kieft B, Li Z, Bryson S, Crump BC, Hettich R, Pan C, et al. Microbial community structure-function relationships in Yaquina Bay estuary reveal spatially distinct carbon and nitrogen cycling capacities. Front Microbiol. 2018;9:1282.

    Article  PubMed  PubMed Central  Google Scholar 

  25. 25.

    Bowen JL, Ward BB, Morrison HG, Hobbie JE, Valiela I, Deegan LA, et al. Microbial community composition in sediments resists perturbation by nutrient enrichment. ISME J. 2011;5:1540–8.

    Article  PubMed  PubMed Central  Google Scholar 

  26. 26.

    Bulseco-McKim AN, Giblin AE, Tucker J, Murphy AE, Sanderman J, Hiller-Bittrolff K, et al. Nitrate addition stimulates microbial decomposition of organic matter in salt marsh sediments. Glob Chang Biol. 2019;25:3224–41.

    Article  Google Scholar 

  27. 27.

    Chen J, McIlroy SE, Archana A, Baker DM, Panagiotou G. A pollution gradient contributes to the taxonomic, functional, and resistome diversity of microbial communities in marine sediments. Microbiome. 2019;7:104.

    Article  PubMed  PubMed Central  Google Scholar 

  28. 28.

    Beinart RA. The Significance of Microbial Symbionts in Ecosystem Processes. mSystems. 2019;4:e00127–19.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  29. 29.

    Nyholm SV, Graf J. Knowing your friends: invertebrate innate immunity fosters beneficial bacterial symbioses. Nat Rev Microbiol. 2012;10:815–27.

    CAS  Article  PubMed  Google Scholar 

  30. 30.

    Shaver EC, Shantz AA, McMinds R, Burkepile DE, Thurber RLV, Silliman BR. Effects of predation and nutrient enrichment on the success and microbiome of a foundational coral. Ecology. 2017;98:830–9.

    Article  PubMed  Google Scholar 

  31. 31.

    Cleary DFR, Swierts T, Coelho FJRC, Polónia ARM, Huang YM, Ferreira MRS, et al. The sponge microbiome within the greater coral reef microbial metacommunity. Nat Commun. 2019;10:1644.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  32. 32.

    Lin HJ, Nixon SW, Taylor DI, Granger SL, Buckley BA. Responses of epiphytes on eelgrass, Zostera marina L., to separate and combined nitrogen and phosphorus enrichment. Aquat Bot. 1996;52:243–58.

    Article  Google Scholar 

  33. 33.

    Crump BC, Wojahn JM, Tomas F, Mueller RS. Metatranscriptomics and amplicon sequencing reveal mutualisms in seagrass microbiomes. Front Microbiol. 2018;388.

  34. 34.

    Pierce ML, Ward JE, Holohan BA, Zhao X, Hicks RE. The influence of site and season on the gut and pallial fluid microbial communities of the eastern oyster, Crassostrea virginica (Bivalvia, Ostreidae): community-level physiological profiling and genetic structure. Hydrobiologia. 2016;765:97–113.

    CAS  Article  Google Scholar 

  35. 35.

    Li Y-F, Xu J-K, Chen Y-W, Ding W-Y, Shao A-Q, Liang X, et al. Characterization of gut microbiome in the mussel Mytilus galloprovincialis in response to thermal stress. Front Physiol. 2019;10:1086.

    Article  PubMed  PubMed Central  Google Scholar 

  36. 36.

    Milan M, Carraro L, Fariselli P, Martino ME, Cavalieri D, Vitali F, et al. Microbiota and environmental stress: how pollution affects microbial communities in Manila clams. Aquat Toxicol. 2018;194:195–207.

    CAS  Article  PubMed  Google Scholar 

  37. 37.

    Marzinelli EM, Qiu Z, Dafforn KA, Johnston EL, Steinberg PD, Mayer-Pinto M. Coastal urbanisation affects microbial communities on a dominant marine holobiont npj Biofilms. Microbiomes. 2018;4:1.

    Article  Google Scholar 

  38. 38.

    Rocca JD, Simonin M, Blaszczak JR, Ernakovich JG, Gibbons SM, Midani FS, et al. The microbiome stress project: towards a global meta-analysis of environmental stressors and their effects on microbial communities. Front Microbiol. 2018;9:3272.

    Article  PubMed  Google Scholar 

  39. 39.

    Grabowski JH, Brumbaugh RD, Conrad RF, Keeler AG, Opaluch JJ, Peterson CH, et al. Economic valuation of ecosystem services provided by oyster reefs. Bioscience. 2012;62:900–9.

    Article  Google Scholar 

  40. 40.

    Kellogg ML, Smyth AR, Luckenbach MW, Carmichael RH, Brown BL, Cornwell JC, et al. Use of oysters to mitigate eutrophication in coastal waters. Estuar Coast Shelf Sci. 2014;151:156–68.

    CAS  Article  Google Scholar 

  41. 41.

    Caffrey JM, Hollibaugh JT, Mortazavi B. Living oysters and their shells as sites of nitrification and denitrification. Mar Pollut Bull. 2016;112:86–90.

    CAS  Article  PubMed  Google Scholar 

  42. 42.

    Wegner K, Volkenborn N, Peter H, Eiler A. Disturbance induced decoupling between host genetics and composition of the associated microbiome. BMC Microbiol. 2013;13:252.

    Article  PubMed  PubMed Central  Google Scholar 

  43. 43.

    Lokmer A, Wegner KM. Hemolymph microbiome of Pacific oysters in response to temperature, temperature stress and infection. ISME J. 2015;9:670–82.

    CAS  Article  PubMed  Google Scholar 

  44. 44.

    Trabal Fernández N, Mazón-Suástegui JM, Vázquez-Juárez R, Ascencio-Valle F, Romero J. Changes in the composition and diversity of the bacterial microbiota associated with oysters (Crassostrea corteziensis, Crassostrea gigas and Crassostrea sikamea) during commercial production. FEMS Microbiol Ecol. 2014;88:69–83.

    Article  Google Scholar 

  45. 45.

    King WL, Siboni N, Kahlke T, Dove M, O’Connor W, Mahbub KR, et al. Regional and oyster microenvironmental scale heterogeneity in the Pacific oyster bacterial community. FEMS Microbiol Ecol. 2020;96:1–12.

    Article  Google Scholar 

  46. 46.

    Lokmer A, Goedknegt MA, Thieltges DW, Fiorentino D, Kuenzel S, Baines JF, et al. Spatial and temporal dynamics of pacific oyster hemolymph microbiota across multiple scales. Front Microbiol. 2016;1367.

  47. 47.

    Lokmer A, Kuenzel S, Baines JF, Wegner KM. The role of tissue-specific microbiota in initial establishment success of Pacific oysters. Environ Microbiol. 2016;18:970–87.

    CAS  Article  PubMed  Google Scholar 

  48. 48.

    King WL, Jenkins C, Seymour JR, Labbate M. Oyster disease in a changing environment: decrypting the link between pathogen, microbiome and environment. Mar Environ Res. 2019;143:124–40.

    CAS  Article  PubMed  Google Scholar 

  49. 49.

    King WL, Jenkins C, Go J, Siboni N, Seymour JR, Labbate M. Characterisation of the Pacific Oyster Microbiome During a Summer Mortality Event. Microb Ecol. 2018:1–11.

  50. 50.

    Pierce ML, Ward JE. Microbial ecology of the Bivalvia, with an emphasis on the family Ostreidae. J Shellfish Res. 2018;37:793–806.

    Article  Google Scholar 

  51. 51.

    Arfken A, Song B, Bowman JS, Piehler M. Denitrification potential of the eastern oyster microbiome using a 16S rRNA gene based metabolic inference approach. PLoS One. 2017;12:e0185071.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  52. 52.

    Murphy AE, Kolkmeyer R, Song B, Anderson IC, Bowen J. Bioreactivity and microbiome of biodeposits from filter-feeding bivalves. Microb Ecol. 2019;77:343–57.

    CAS  Article  PubMed  Google Scholar 

  53. 53.

    Antczak M, Michaelis M, Wass MN. Environmental conditions shape the nature of a minimal bacterial genome. Nat Commun. 2019;10:3100.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  54. 54.

    Petrova OE, Garcia-Alcalde F, Zampaloni C, Sauer K. Comparative evaluation of rRNA depletion procedures for the improved analysis of bacterial biofilm and mixed pathogen culture transcriptomes. Sci Rep. 2017;7:41114.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  55. 55.

    Herbert ZT, Kershner JP, Butty VL, Thimmapuram J, Choudhari S, Alekseyev YO, et al. Cross-site comparison of ribosomal depletion kits for Illumina RNAseq library construction. BMC Genomics. 2018;19:199.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  56. 56.

    Howarth RW. Nutrient limitation of net primary production in marine ecosystems. Annu Rev Ecol Syst. 1988;19:89–110.

    Article  Google Scholar 

  57. 57.

    Testa JM, Brady DC, Di Toro DM, Boynton WR, Cornwell JC, Kemp WM. Sediment flux modeling: simulating nitrogen, phosphorus, and silica cycles. Estuar Coast Shelf Sci. 2013;131:245–63.

    CAS  Article  Google Scholar 

  58. 58.

    Alzate Marin JC, Caravelli AH, Zaritzky NE. Nitrification and aerobic denitrification in anoxic-aerobic sequencing batch reactor. Bioresour Technol. 2016;200:380–7.

    CAS  Article  PubMed  Google Scholar 

  59. 59.

    Rysgaard S, Risgaard-Petersen N, Niels Peter S, Kim J, Lars PN. Oxygen regulation of nitrification and denitrification in sediments. Limnol Oceanogr. 1994;39:1643–52.

    CAS  Article  Google Scholar 

  60. 60.

    Smith MS, Tiedje JM. Phases of denitrification following oxygen depletion in soil. Soil Biol Biochem. 1979;11:261–7.

    CAS  Article  Google Scholar 

  61. 61.

    Gomez E, Durillon C, Rofes G, Picot B. Phosphate adsorption and release from sediments of brackish lagoons: pH, O2 and loading influence. Water Res. 1999;33:2437–47.

    CAS  Article  Google Scholar 

  62. 62.

    Lam P, Kuypers MMM. Microbial nitrogen cycling processes in oxygen minimum zones. Annu Rev Mar Sci. 2010;3:317–45.

    Article  Google Scholar 

  63. 63.

    Sowell SM, Wilhelm LJ, Norbeck AD, Lipton MS, Nicora CD, Barofsky DF, et al. Transport functions dominate the SAR11 metaproteome at low-nutrient extremes in the Sargasso Sea. ISME J. 2009;3:93–105.

    CAS  Article  PubMed  Google Scholar 

  64. 64.

    Robledo JAF, Yadavalli R, Allam B, Pales-Espinosa E, Gerdol M, Greco S, et al. From the raw bar to the bench: Bivalves as models for human health. Dev Comp Immunol. 2018;92:260–82.

    Article  Google Scholar 

  65. 65.

    Ward JE, Shumway SE. Separating the grain from the chaff: particle selection in suspension- and deposit-feeding bivalves. J Exp Mar Biol Ecol. 2004;300:83–130.

    Article  Google Scholar 

  66. 66.

    Lesser MP. Oxidative stress in marine environments: biochemistry and physiological ecology. Annu Rev Physiol. 2006;68:253–78.

    CAS  Article  PubMed  Google Scholar 

  67. 67.

    Alves de Almeida E, Celso Dias Bainy A, Paula de Melo Loureiro A, Regina Martinez G, Miyamoto S, Onuki J, et al. Oxidative stress in Perna perna and other bivalves as indicators of environmental stress in the Brazilian marine environment: Antioxidants, lipid peroxidation and DNA damage. Comp Biochem Physiol A Mol IntegPhysiol. 2007;146:588–600.

    CAS  Article  Google Scholar 

  68. 68.

    Spiro S. Nitrous oxide production and consumption: regulation of gene expression by gassensitive transcription factors. Philos Trans R Soc B Biol Sci. 2012;367:1213–25.

    CAS  Article  Google Scholar 

  69. 69.

    Reyes C, Schneider D, Lipka M, Thürmer A, Böttcher ME, Friedrich MW. Nitrogen metabolism genes from temperate marine sediments. Mar Biotechnol. 2017;19:175–90.

    CAS  Article  PubMed Central  Google Scholar 

  70. 70.

    Raivio TL, Silhavy TJ. Periplasmic stress and ECF sigma factors. Annu Rev Microbiol. 2002;55:591–624.

    Article  Google Scholar 

  71. 71.

    Zhang Y, Wang L, Hu Y, Xi X, Tang Y, Chen J, et al. Water organic pollution and eutrophication influence soil microbial processes, increasing soil respiration of estuarine wetlands: site study in Jiuduansha wetland. PLoS One. 2015;10:e0126951.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  72. 72.

    Meyer-Reil LA, Köster M. Eutrophication of marine waters: effects on benthic microbial communities. Mar Pollut Bull. 2000;41:255–63.

    CAS  Article  Google Scholar 

  73. 73.

    Parfrey LW, Groussin M, Mazel F, Loudon A, Kwong WK, Davis KM. Is Host Filtering the Main Driver of Phylosymbiosis across the Tree of Life? mSystems. 2018:97–115.

  74. 74.

    Apprill A. The Role of Symbioses in the Adaptation and Stress Responses of Marine Organisms. Annu Rev Mar Sci. 2020;12:annurev-marine-010419-010641.

  75. 75.

    McAuliffe L, Ellis RJ, Miles K, Ayling RD, Nicholas RAJ. Biofilm formation by mycoplasma species and its role in environmental persistence and survival. Microbiology. 2006;152:913–22.

    CAS  Article  PubMed  Google Scholar 

  76. 76.

    Fraune S, Zimmer M. Host-specificity of environmentally transmitted mycoplasma-like isopod symbionts. Environ Microbiol. 2008;10:2497–504.

    CAS  Article  PubMed  Google Scholar 

  77. 77.

    Stevick RJ, Sohn S, Modak TH, Nelson DR, Rowley DC, Tammi K, et al. Bacterial community dynamics in an oyster hatchery in response to probiotic treatment. Front Microbiol. 2019;10:1060.

    Article  PubMed  PubMed Central  Google Scholar 

  78. 78.

    King WL, Siboni N, Williams NLR, Kahlke T, Nguyen KV, Jenkins C, et al. Variability in the composition of pacific oyster microbiomes across oyster families exhibiting different levels of susceptibility to OsHV-1 μvar disease. Front Microbiol. 2019;473.

  79. 79.

    Dupont S, Lokmer A, Corre E, Auguet J-C, Petton B, Toulza E, et al. Oyster hemolymph is a complex and dynamic ecosystem hosting bacteria, protists and viruses. Anim Microbiome. 2020;2:12.

    Article  Google Scholar 

  80. 80.

    Kinross JM, Darzi AW, Nicholson JK. Gut microbiome-host interactions in health and disease. Genome Med. 2011;3:14.

    Article  PubMed  PubMed Central  Google Scholar 

  81. 81.

    de Lorgeril J, Escoubas JM, Loubiere V, Pernet F, Le Gall P, Vergnes A, et al. Inefficient immune response is associated with microbial permissiveness in juvenile oysters affected by mass mortalities on field. Fish Shellfish Immunol. 2018;77:156–63.

    CAS  Article  PubMed  Google Scholar 

  82. 82.

    Clerissi C, de Lorgeril J, Petton B, Lucasson A, Escoubas J-M, Gueguen Y, et al. Microbiota composition and evenness predict survival rate of oysters confronted to Pacific oyster mortality syndrome. Front Microbiol. 2020;11:311.

    Article  PubMed  PubMed Central  Google Scholar 

  83. 83.

    Li R, Tun HM, Jahan M, Zhang Z, Kumar A, Fernando D, et al. Comparison of DNA-, PMA-, and RNA-based 16S rRNA Illumina sequencing for detection of live bacteria in water. Sci Rep. 2017;7(1):1–11.

  84. 84.

    Pollock J, Glendinning L, Wisedchanwet T, Watson M. The madness of microbiome: attempting to find consensus “best practice” for 16S microbiome studies. Appl Environ Microbiol. 2018;84:e02627–17.

    Article  PubMed  PubMed Central  Google Scholar 

  85. 85.

    Shakya M, Quince C, Campbell JH, Yang ZK, Schadt CW, Podar M. Comparative metagenomic and rRNA microbial diversity characterization using archaeal and bacterial synthetic communities. Environ Microbiol. 2013;15:1882–99.

    CAS  Article  Google Scholar 

  86. 86.

    Hornung BVH, Zwittink RD, Kuijper EJ. Issues and current standards of controls in microbiome research. FEMS Microbiol Ecol. 2019;95:45.

    CAS  Article  Google Scholar 

  87. 87.

    Newell RI, Jordan SJ. Preferential ingestion of organic material by the American oyster Crassostrea virginica. Mar Ecol Prog Ser. 1983;13:47–53.

    Article  Google Scholar 

  88. 88.

    Rampadarath S, Bandhoa K, Puchooa D, Jeewon R, Bal S. Early bacterial biofilm colonizers in the coastal waters of Mauritius. Electron J Biotechnol. 2017;29:13–21.

    CAS  Article  Google Scholar 

  89. 89.

    Riiser ES, Haverkamp THA, Borgan Ø, Jakobsen KS, Jentoft S, Star B. A single vibrionales 16S rRNA oligotype dominates the intestinal microbiome in two geographically separated Atlantic cod populations. Front Microbiol. 2018;1561.

  90. 90.

    Schumann S. Rhode Island’s shellfish heritage: an ecological history. Narragansett: University of Rhode Island; 2015.

    Google Scholar 

  91. 91.

    Zehr JP, Church MJ, Moisander PH. Diversity, distribution and biogeochemical significance of nitrogen-fixing microorganisms in anoxic and Suboxic Ocean environments. In: Past and present water column anoxia. Dordrecht: Kluwer Academic Publishers; 2006. p. 337–69.

    Google Scholar 

  92. 92.

    Howarth R, Chan F, Conley DJ, Garnier J, Doney SC, Marino R, et al. Coupled biogeochemical cycles: Eutrophication and hypoxia in temperate estuaries and coastal marine ecosystems. Coupled biogeochemical cycles: eutrophication and hypoxia in temperate estuaries and coastal marine ecosystems. Frontiers in Ecology and the Environment. 2011;9(1):18–26.

  93. 93.

    Codiga DL, Stoffel HE, Deacutis CF, Kiernan S, Oviatt CA. Narragansett Bay hypoxic event characteristics based on fixed-site monitoring network time series: intermittency, geographic distribution, spatial synchronicity, and interannual variability. Estuar Coasts. 2009;32:621–41.

    CAS  Article  Google Scholar 

  94. 94.

    Graham EB, Knelman JE, Schindlbacher A, Siciliano S, Breulmann M, Yannarell A, et al. Microbes as engines of ecosystem function: When does community structure enhance predictions of ecosystem processes? Front Microbiol. 2016;214.

  95. 95.

    Eren AM, Vineis JH, Morrison HG, Sogin ML. A filtering method to generate high quality short reads using Illumina paired-end technology. PLoS One. 2013;8:e66643.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  96. 96.

    Huse SM, Mark Welch DB, Voorhis A, Shipunova A, Morrison HG, Eren AM, et al. VAMPS: a website for visualization and analysis of microbial population structures. BMC Bioinformatics. 2014;15:41.

    Article  PubMed  PubMed Central  Google Scholar 

  97. 97.

    Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJA, Holmes SP. DADA2: high-resolution sample inference from Illumina amplicon data. Nat Methods. 2016;13:581–3.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  98. 98.

    Caporaso JG, Kuczynski J, Stombaugh J, Bittinger K, Bushman FD, Costello EK, et al. QIIME allows analysis of high-throughput community sequencing data. Nat Methods. 2010;7:335–6.

    CAS  Article  Google Scholar 

  99. 99.

    Bokulich NA, Kaehler BD, Rideout JR, Dillon M, Bolyen E, Knight R, et al. Optimizing taxonomic classification of marker-gene amplicon sequences with QIIME 2’s q2-feature-classifier plugin. Microbiome. 2018;6:90.

    Article  PubMed  PubMed Central  Google Scholar 

  100. 100.

    Bolyen E, Rideout JR, Dillon MR, Bokulich NA, Abnet CC, Al-Ghalith GA, et al. Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nat Biotechnol. 2019;37:852–7.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  101. 101.

    Segata N, Izard J, Waldron L, Gevers D, Miropolsky L, Garrett WS, et al. Metagenomic biomarker discovery and explanation. Genome Biol. 2011;12:R60.

    Article  Google Scholar 

  102. 102.

    Afgan E, Baker D, Batut B, Van Den Beek M, Bouvier D, Ech M, et al. The galaxy platform for accessible, reproducible and collaborative biomedical analyses: 2018 update. Nucleic Acids Res. 2018;46(W1):W537–W544.

  103. 103.

    Dixon P. VEGAN, a package of R functions for community ecology. J Veg Sci. 2003;14:927–30.

    Article  Google Scholar 

  104. 104.

    McMurdie PJ, Holmes S. Phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data. PLoS One. 2013;8:e61217.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  105. 105.

    Chao A, Jost L. Coverage-based rarefaction and extrapolation: standardizing samples by completeness rather than size. Ecology. 2012;93:2533–47.

    Article  PubMed  Google Scholar 

  106. 106.

    Mcardle BH, Anderson MJ. Fitting Multivariate Models to Community Data : A Comment on Distance-Based Redundancy Analysis Published by : Ecological Society of America Stable URL : Ecology. 2010;82:290–7. doi:

  107. 107.

    Warton DI, Wright ST, Wang Y. Distance-based multivariate analyses confound location and dispersion effects. Methods Ecol Evol. 2012;3:89–101.

    Article  Google Scholar 

  108. 108.

    Gu Z, Eils R, Schlesner M. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics. 2016;32:2847–9.

    CAS  Article  PubMed  Google Scholar 

  109. 109.

    Conway JR, Lex A, Gehlenborg N. UpSetR: an R package for the visualization of intersecting sets and their properties. Bioinformatics. 2017;33:2938–40.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  110. 110.

    Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.

    CAS  Article  Google Scholar 

  111. 111.

    Westreich ST, Treiber ML, Mills DA, Korf I, Lemay DG. SAMSA2: A standalone metatranscriptome analysis pipeline. BMC Bioinformatics. 2018;19:175.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  112. 112.

    Zhang J, Kobert K, Flouri T, Stamatakis A. PEAR: a fast and accurate Illumina paired-end reAd mergeR. Bioinformatics. 2014;30:614–20.

    CAS  Article  PubMed  Google Scholar 

  113. 113.

    Kopylova E, Noé L, Touzet H. SortMeRNA: fast and accurate filtering of ribosomal RNAs in metatranscriptomic data. Bioinformatics. 2012;28:3211–7.

    CAS  Article  PubMed  Google Scholar 

  114. 114.

    Buchfink B, Xie C, Huson DH. Fast and sensitive protein alignment using DIAMOND. Nat Methods. 2014;12:59–60.

    CAS  Article  PubMed  Google Scholar 

  115. 115.

    Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  116. 116.

    Wickham H. ggplot2: elegant graphics for data analysis. New York: Springer-Verlag; 2009. ISBN 978-3-319-24277-4.

  117. 117.

    R Development Core Team R. R: A Language and Environment for Statistical Computing. 2011.

    Google Scholar 

  118. 118.

    Wilke CO. cowplot: Streamlined Plot Theme and Plot Annotations for “ggplot2”. 2019.

    Google Scholar 

  119. 119.

    Stevick R. Analyses for functional plasticity in oyster gut microbiomes along a eutrophication gradient in an urbanized estuary. Zenodo. 2020.

  120. 120.

    Tang Y, Horikoshi M, Li W. ggfortify: Unified Interface to Visualize Statistical Result of Popular {R} Packages. {R} J. 2016;8:478–89.

  121. 121.

    Kassambara A. ggpubr: “ggplot2” Based Publication Ready Plots. 2019.

Download references


We are grateful to the Rhode Island Department of Environmental Management for assistance with site permits for oyster collections, and all the students and friends who helped with field collections. We would like to thank the Yale Center for Genome Analysis and Janet Atoyan at the URI Genomics and Sequencing Center for sequencing support and analytics, and Rodrigue Spinette for assistance with the nutrient analysis. 


This work was supported by URI Coastal Institute Grants in Aid, The Nature Conservancy of RI - TNC Global Marine Team, NSF Graduate Research Fellowship 1244657 to RJS, National Science Foundation EPSCoR Cooperative Agreement #OIA-1655221, and the Blount Family Shellfish Restoration Foundation. This material is based upon work conducted at a Rhode Island NSF EPSCoR research facility, the Genomics and Sequencing Center, supported in part by the National Science Foundation EPSCoR Cooperative Agreement #OIA-1655221.

Author information




RJS, AFP, and MG-C contributed conception and design of the study; RJS collected and prepared the samples for sequencing, performed the sequence analysis, and wrote the manuscript; all authors contributed to manuscript revision, read and approved the submitted version.

Corresponding author

Correspondence to Marta Gómez-Chiarri.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Figure S1.

(A) Percent abundances of the 10 most abundant phyla for 16S rRNA gene amplicon sequencing data by sample type and site. All other taxa are grouped into “Others.” Mock community control samples are shown at the left. (B) Percent abundance of the top 10 most abundant ASVs in the mock community samples. All other taxa are grouped into “Others.”

Additional file 2: Figure S2.

Sequencing coverage analysis for 16S rRNA gene amplicon (left) and Metatranscriptomic samples (right). (A) Rarefaction curves for each sample colored by site for gut samples or water samples. The minimum sample size (raremax) is shown and indicated with a dashed line on each plot. (B) The slope calculated at the raremax for each sample is shown. (C) The estimated coverage (100–100*slope) for each sample is shown. Mean coverage and standard deviation for each method is shown in the bottom right.

Additional file 3: Figure S3.

Number of bacterial Orders shared between the oyster gut and seawater 16S rRNA gene amplicons at each site (vertical bars). The total number of Orders found in each group is shown in the horizontal bar graph on the right. Intersections in gray denote comparisons that include the water samples.

Additional file 4: Figure S4.

(A) Linear discriminant analysis Effect Size (LEfSe) analysis of bacterial Orders in seawater (n = 10), compared to gut samples (n = 50) in the 16S rRNA gene amplicons (one-against-all). (B) LEfSe analysis of bacterial Orders in gut samples at each site (n = 10) in the 16S rRNA gene amplicons (all-against-all). Only significantly increased taxa are shown. Significance was determined by LDA score > 2.0, alpha value = 0.05 for factorial Kruskal-Wallis test, and alpha value = 0.05 for pairwise Wilxocon test.

Additional file 5: Figure S5

. NMDS plot visualizations of Bray-Curtis beta-diversity (k = 2) at the (A) Species and (B) Order levels for gut metatranscriptomic samples by Site. The ellipse lines show the 95% confidence interval (standard deviation). p-values indicate significance of grouping with adonis2 Permutational Multivariate Analysis of Variance Using Distance Matrices test.

Additional file 6: Figure S6.

Heatmap of taxa identified as the core bacterial community in the (A) 16S rRNA gene amplicon data and the (B) metatranscriptomic data. Green boxes indicate the presence of the core taxa in each sample per site. Core taxa was defined as occurring in > 80% of the samples per sequencing type.

Additional file 7: Figure S7.

Differential expression (log fold change) of SEED Level 4 gene annotation of Oxidative Stress response groups at each site, relative to the mean of the others. All significantly regulated genes are outlined in red and annotated with an asterisk (n = 5, Benjamini-Hochberg *padj < 0.05, **padj < 0.01).

Additional file 8: Figure S8.

Differential expression (log fold change) of SEED Level 4 gene annotation of Osmotic and Periplasmic stress response groups at each site, relative to the mean of the others. All significantly regulated genes are outlined in red and annotated with an asterisk (n = 5, Benjamini-Hochberg *padj < 0.05, **padj < 0.01).

Additional file 9: Figure S9.

Differential expression (log fold change) of SEED level 4 gene annotation of nitrogen metabolism pathways at each site, relative to the mean of the others. All significantly regulated genes are outlined in red and annotated with an asterisk (n = 5, Benjamini-Hochberg *padj < 0.05, **padj < 0.01).

Additional file 10: Figure S10.

Differential expression (log fold change) of SEED level 4 gene annotation of phosphorus metabolism pathways at each site, relative to the mean of the others. All significantly regulated genes are outlined in red and annotated with an asterisk (n = 5, Benjamini-Hochberg *padj < 0.05, **padj < 0.01).

Additional file 11: Table S1.

Sequencing summary statistics, including the number of reads that passed quality control (QC) in each 16S rRNA gene amplicon and metatranscriptomic sample. No metatranscriptomes were sequenced for water samples.

Additional file 12: Table S2.

Bray-Curtis beta-diversity summary statistics calculated using adonis2 Permutational Multivariate Analysis of Variance Using Distance Matrices test.

Additional file 13: Table S3.

Taxa present in > 80% of all samples in the 16S rRNA gene amplicon and metatranscriptomic datasets. This table corresponds to Fig. S6.

Additional file 14: Table S4.

Raw sequencing data accession information. This spreadsheet includes BioSample information for each oyster with corresponding Accession numbers for 16S rRNA gene amplicon or metatranscriptome files.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Stevick, R.J., Post, A.F. & Gómez-Chiarri, M. Functional plasticity in oyster gut microbiomes along a eutrophication gradient in an urbanized estuary. anim microbiome 3, 5 (2021).

Download citation


  • Microbiome
  • 16S rRNA gene sequencing
  • Metatranscriptomics
  • Coastal eutrophication
  • Oysters
  • Crassostrea virginica