Multiple stressors interact primarily through antagonism to drive changes in the coral microbiome

Perturbations in natural systems generally are the combination of multiple interactions among individual stressors. However, methods to interpret the effects of interacting stressors remain challenging and are biased to identifying synergies which are prioritized in conservation. Therefore we conducted a multiple stressor experiment (no stress, single, double, triple) on the coral Pocillopora meandrina to evaluate how its microbiome changes compositionally with increasing levels of perturbation. We found that effects of nutrient enrichment, simulated predation, and increased temperature are antagonistic, rather than synergistic or additive, for a variety of microbial community diversity measures. Importantly, high temperature and scarring alone had the greatest effect on changing microbial community composition and diversity. Using differential abundance analysis, we found that the main effects of stressors increased the abundance of opportunistic taxa, and two-way interactions among stressors acted antagonistically on this increase, while three-way interactions acted synergistically. These data suggest that: (1) multiple statistical analyses should be conducted for a complete assessment of microbial community dynamics, (2) for some statistical metrics multiple stressors do not necessarily increase the disruption of microbiomes over single stressors in this coral species, and (3) the observed stressor-induced community dysbiosis is characterized by a proliferation of opportunists rather than a depletion of a proposed coral symbiont of the genus Endozoicomonas.

The two therapies likely disrupt the gut microbiome through different mechanisms and are independently associated with dysbiosis. In one case, the stressors produced opposite responses in the abundance of a single bacterial genus, Alistipes. Yet, an antagonistic interaction was not tested for but easily could be with a crossed design with patients receiving both therapies. In an environmental example, warm-or cold-stressed oysters crossed with bacterial infection by vibrios showed evidence of synergy as warm-stressed oysters experienced the highest mortality following infection 13 . When evaluating the oyster hemolymph microbiome, an interaction term of stress × infection in the univariate analysis of alpha diversity and multivariate analysis of beta diversity was not included, but if included in statistical methods, would clarify the type of interaction between the two stressors.
Using robust statistical methods and interaction models benchmarked in the microbiome field 14 , we investigated how a global stressor, thermal stress, interacts with local stressors, nutrient pollution and predation, to alter the coral microbiome. Corals, currently experiencing major threats of climate change and nutrient pollution, can function as environmental sentinels and are thereby prime candidates for multiple stressor experiments. Previous studies of the coral microbiome have shown that stress tends to increase species richness 8,9,[15][16][17] and cause shifts in community composition from potentially beneficial symbiotic bacteria that dominate healthy corals to potentially opportunistic or pathogenic bacteria that dominate stressed corals 8,9,15,[17][18][19][20] . Beta diversity, or species turnover between samples, has also been reported to increase with stress 7,9,21,22 , and stressed corals have microbial communities distinct from control corals [23][24][25][26] .
Therefore, we designed a fully-crossed experiment to investigate biological responses including alpha and beta diversity indices and differential abundance modeling of individual taxa with stress. For the purposes of this study, we define a stressor to be any external disturbance from the host's environment that causes a quantifiable change in microbial community structure. We utilized univariate and multivariate statistical techniques to parse the main effects and interactions among stressors. The coral Pocillopora meandrina was exposed to increased seawater temperatures, pulse nitrate and ammonium enrichment, and simulated predation in a factorial mesocosm tank experiment with all possible combinations of these stressors. We hypothesized that local stressors like nutrient pollution and predation would interact synergistically with thermal stress to reduce the host's ability to regulate its microbial community which would be manifested by: (1) an increase in the compositional heterogeneity and variability (beta diversity) among stressed corals compared to the controls, and (2) an increase in community evenness in stressed corals as a result of (3) shifts from few dominant symbiotic bacterial taxa to a myriad of potentially opportunistic bacterial taxa that bloom and become overrepresented in stressed corals.

Results stressors drive symbiont decreases and opportunist increases in relative abundance.
To test the individual and interactive effects of increased temperature, nutrient enrichment, and predation on coral microbiomes, replicates of the coral Pocillopora meandrina were exposed to each combination of individual, double, and triple stressors in a fully crossed tank experiment ( Supplementary Fig. S1). Patterns in the relative abundance of different microbial taxa were assessed with generalized linear mixed-effects models (GLMM) and clearly revealed a dominant member of the coral community (Fig. 1) of the genus Endozoicomonas are proposed bacterial symbionts with coral, and are typically underrepresented in stressed corals 8,[27][28][29] . All control corals were dominated by OTU-Endo with an average relative abundance of 94.88 ± 0.86% (Fig. 1). The main effects of high temperature, scarring, and NO 3 − all significantly decreased the relative abundance of OTU-Endo in the GLMM (p < 0.001, p < 0.001, p < 0.05; Supplementary Table S1, Fig. 1). Compared to control corals, the relative abundance of the dominant OTU-Endo decreased by more than 50% in the high temperature treatment (41.19 ± 16.61%) and decreased nearly 50% in the scarred treatment (48.21 ± 16.61%). High temperature and scarring interact to reduce the sum of the independent effects on the dominant taxon (p < 0.01; Supplementary Table S1), decreasing OTU-Endo by only ~25% (70.67 ± 5.44) compared to controls. In fact, all two-way interactions were antagonistic and significantly reduced the magnitude of the response predicted by main effects (Supplementary Table S1). Alternatively, both three-way interactions significantly increased the response of the individual stressors on OTU-Endo after accounting for all main effects and two-way interactions (Supplementary Table S1).
The decreases observed in OTU-Endo closely mirror increases in the second most abundant OTU of the family Desulfovibrionaceae (OTU-Desulfo). For instance, all four main effects significantly increased the relative abundance of OTU-Desulfo (Supplementary Table S2) with high temperature causing the greatest change in this OTU from control corals with 0.09 ± 0.09% to 25.56 ± 13.44%, followed by NO 3 − (10.62 ± 10.18%) and scarring (9.92 ± 6.30%). Alternatively, the interaction of high temperature and scarring only increased the relative abundance of OTU-Desulfo to 1.51 ± 0.78%. Two-way and three-way interactions showed similar significant, but opposite directional changes as those observed for OTU-Endo (Supplementary Table S2). For the third most abundant OTU of the family Enterobacteriaceae, the main effect of scarring caused significant increases in abundance and significant interactions between scarring and nutrient treatments (Supplementary Table S2). The OTU from family Amoebophilaceae significantly increased in scarring and high temperature and the three-way interaction with NH 4 + was significant (Supplementary Table S2). And lastly, the OTU from family Moraxellaceae significantly increased in the NH 4 + treatment and two-way interactions with scarring and nutrient treatments were significant (Supplementary Table S2). For all OTUs, significant two-way interactions were antagonistic, or less than the sum of the main effects, and three-way interactions were synergistic when accounting for main effects and two-way interactions. When evaluating changes in relative abundance, the microbial community under stress is marked both by a decrease in the dominant symbiont and increases in lower abundant opportunistic taxa.
High temperature and scarring drive changes in community diversity. Alpha diversity metrics were assessed with linear mixed effects models (LMM) and showed the absence of synergisms. On average, Chao1 estimates found that control corals had a mean of 40.42 ± 3.48 unique OTUs and had the lowest standard error across all treatments ( Supplementary Fig. S2). Neither main effects nor interactions were significant predictors in the LMM with Chao1 index as the response variable (Supplementary Table S3). Faith's phylogenetic diversity ranged from 4.09 ± 0.55 in control corals to 8.99 ± 1.81 for corals in high temperature. The main effects of high temperature and scarring significantly increased Faith's phylogenetic diversity (p < 0.05, p < 0.05; Supplementary  Table S3). Additionally, the interaction between high temperature and NH 4 + showed a significantly antagonistic response on Faith's phylogenetic diversity (p < 0.05, Supplementary Table S3). Although evidence suggests that stressors generally increase microbial species richness in coral microbiomes 8 , we found no significant differences in species richness, but rather changes in phylogenetic diversity with stressors.
All treatments tended to increase Simpson's index or community diversity (3-to 7-fold) compared to the controls, with an index of 0.10 ± 0.02 (Fig. 2a). Patterns in community diversity closely match those of the dominant OTU-Endo and reflect changes in community 'evenness' which is accounted for when calculating the Simpson's index. Both high temperature and scarring significantly increased Simpson's Index compared to the controls (p < 0.01 and p < 0.01; Supplementary Table S3). High temperature and scarred treatments produced the highest mean community diversity (0.62 ± 0.16 and 0.67 ± 0.11, respectively), greater than 6 times that of the controls. Interestingly, scarring and high temperature interact to reduce the independent effects in the regression analysis (p < 0.01; Supplementary Table S3). This interaction can be observed with interaction plots and described in two ways: (1) high temperature reverses the effect of increased community diversity from control to scarred corals (Fig. 2b), or (2) scarring decreases the difference in community diversity between corals in ambient to high temperature seawater (Fig. 2c). Enrichment with NH 4 + or NO 3 − also significantly interacts to reverse the main effect of scarring, thereby decreasing community diversity in scarred corals compared to controls (p < 0.01, p < 0.05; Supplementary Table S3; Supplementary Fig. S3). The two three-way interactions between scarring, high temperature, and nutrient enrichment also produced a significant result (p < 0.05, p < 0.05; Supplementary Table S3), suggesting that the interaction between any two stressors depends on the level of the third stressor. In the linear model framework, interaction type is not directly interpretable with the less-than-or greater-than-additive definition since three-way interactions in the model account for all main effects and two-way interactions. Therefore, the three-way interaction increased community diversity after accounting for all antagonistic two-way interactions, although the triple stressor treatments are less than most single stressor treatments (Fig. 2a). Despite the significant main effects and interactions in the LMM, no pairwise treatment comparisons were significant with a Tukey post hoc test.
Beta diversity measures show less-than-additive effects during microbiome exposure to multiple stressors. PERMANOVA results, with treatment as the predicting factor, showed the presence of distinct microbial communities (Fig. 3a p-value, Supplementary Table S4). To visualize differences in community location, mean separation from communities in other treatments (between-group-distance) was calculated (Fig. 3a). Control corals had low average between-group distances and were significantly different from single stressors of NH 4 + , high temperature, and scarring ( Fig. 3a, Supplementary Table S4). However, unlike our hypothesis that with additional amounts of stress microbiomes would become increasingly distinct, single stressors of scarring www.nature.com/scientificreports www.nature.com/scientificreports/ and high temperature produced the greatest mean between-group distances rather than interacting stressors (Fig. 3a, yellow versus purple and red boxes). In fact, high temperature and scarring combined was significantly different from controls but not from either high temperature or scarring alone (Fig. 3a, Supplementary  Table S4). Additionally, triple stressors (red boxes) were not significantly different from single stressors (Fig. 3a, Supplementary Table S4). Therefore, the combination of multiple stressors had less-than-additive effects on the change in composition of the microbial communities. And, while a true synergism is not possible for relative distance measures with a maximum value of 1.0, combined stressors are generally less than single stressors. The linear model PERMANOVA results showed that in corals experiencing the main effect of high temperature, microbial communities were significantly distinct from controls (p < 0.05; Supplementary Table S5). The three-way interaction between temperature, nutrients, and scarring was also significant (p < 0.05; Supplementary www.nature.com/scientificreports www.nature.com/scientificreports/ Table S5), suggesting that the biological response to a single stressor is influenced by the other two. Like double stressor treatments, triple stressors ( Fig. 3a red boxes) interact antagonistically to produce less community distance or distinctness than scarring or high temperature alone. Our results suggest that the high temperature as a main effect and high temperature and scarring as independent treatments have the most influential effect on shifting microbial communities.
PERMDISP results also showed significant differences among treatments in dispersion magnitude (Fig. 3b  p-value, Supplementary Table S6). Distance to group centroid measures showed that control corals harbored microbial communities that were more homogenous, and therefore less dispersed, than those belonging to stressed corals (Fig. 3b). The addition of stressors (with the exception of NO 3 and scarring combined) increased distance-to-centroid and increased variance suggesting that stress causes dispersion or destabilization of the microbiome. These stochastic changes may show evidence of the Anna Karenina Principle (AKP) in which dysbiotic animal microbiomes vary more in community composition than healthy microbiomes 21 . When considering the dominance of OTU-Endo within the community (Fig. 1), however, evidence suggests that the single taxon is driving the dispersion effect. In fact, the relative abundance of OTU-Endo was significantly negatively correlated with sample distance-to-centroid measurements (estimate: −0.35, p < 0.001). The ability to detect statistically meaningful variance in other taxa is marginal due to the dominance of OTU-Endo. Instead, dispersion may be artificially increased in stressed corals due to an increase in the relativized number of rare taxa and the reduction in OTU-Endo dominance.

Differential abundance analysis shows stress is marked by an increase in opportunists.
From differential abundance analysis with DESeq2, a total of 56 unique OTUs were differentially abundant in one or more of the treatments or interactions (Fig. 4, Supplementary Table S7). On average, main effects resulted in a 2.80 log2 fold-change in the differentially abundant OTUs (yellow-colored line; Fig. 4). This increase was driven primarily by high temperature and scarred treatments which have an average log2 fold-change of 6.35 and 7.23, respectively (yellow-dashed lines, Fig. 4). In contrast, NO 3 − and NH 4 + treatments on average significantly decreased taxa compared to controls (−5.92 and −10.06, respectively, yellow-dashed line, Fig. 4). The NH 4 + treatment resulted in the fewest differentially abundant OTUs (Fig. 4), suggesting that the microbial communities under NH 4 + enrichment are most similar to the control corals. Ammonium (NH 4 + ) is a fish-derived form of nitrogen that can be beneficial to corals 30 , which likely would reflect a healthy microbial community. A single OTU from the proposed symbiont family of Endozoicimonaceae increased in abundance in high temperature, however, this taxon was not the same dominant OTU-Endo from Fig. 1. In fact, OTU-Endo was not identified as having significantly changed in any treatments. When differential abundance analysis was repeated with an OTU table summarized by Family and Genus, higher-level significant changes generally agree with those of individual OTUs (Supplementary Table S8). Therefore, the average increase in abundance of OTUs from families such as Rhodobacteraceae, Sphingomonadaceae, and Desulfovibrionaceae (Fig. 4), suggests that the changes in relative abundance and dysbiosis resulting from stressors are characterized by an enrichment of pathogenic or opportunistic bacteria rather than a depletion of symbionts.  Table S4). (b) Weighted Unifrac mean distance-to-centroid by treatment. P-values denote the significance of treatment group in a PERMDISP using betadisper. Pairwise betadisper was used to assign significance codes for group distances (Supplementary Table S6). Groups sharing a letter are not significantly different from each other. Box colors represent the type of stressor combination applied to the corals: none = teal, single = yellow, double = purple, triple = red. Distances are ordered by increasing median, and note that red (triple stressor) treatments are not clustered on the far left.
www.nature.com/scientificreports www.nature.com/scientificreports/ Interaction type depends on the type and number of stressors. The response of bacterial taxa in two-way interactions was variable depending on family of OTU, OTU, and type of interacting stressors. Generally, two-way interactions between stressors were antagonistic in nature aside from the interaction between NO 3 − and high temperature which produced a synergistic average increase in taxa (purple-colored box, Fig. 4). In the linear framework model, zero log2-fold change between two main effects (grey lines in purple-colored box, Fig. 4) signifies no interaction or the additive model (sum of the two main effects). Instead, the model shows that the average Figure 4. Differential abundance of OTUs with significant response to stressors using DESeq2 and the generalized linear model framework. Only OTUs present in greater than 15 samples were included in the analysis. Each point represents a single OTU that increased or decreased significantly (FDR corrected p < 0.05) with the stressor or stressor interaction. Each row and dot color corresponds to a microbial family (i.e. multiple OTUs from a single family are increased in multiple treatments). Family name in [] represents a recommended taxonomic annotation by GreenGenes. Box colors represent the stressor type: main effect = yellow, two-way interaction = purple, three-way interaction = red. The colored line within each box represents the mean log2FoldChange for OTUs with significant changes in that stressor type. The dashed colored line within comparisons represents the mean log2FoldChange for OTUs with significant changes within the individual stressor effects and interactions. The gray colored line at 0 log2FoldChange denotes no effect or no interaction. www.nature.com/scientificreports www.nature.com/scientificreports/ log2-fold change in bacterial taxa for two-way interactions was less than the sum of the individual main effects by −4.95 (purple-colored line, Fig. 4). This suggests that when combined, two individual stressors act antagonistically to dampen the main effects. The family Rhodobacteraceae had the most differentially abundant OTUs. A total of 10 and 14 OTUs from this family significantly increased in abundance in the high temperature and scarred treatments, respectively, whereas, 11 of these shared OTUs decreased in the two-way interaction. Likewise, three OTUs of the family Desulfovibrionaceae all increased in the high temperature and scarred treatments. However, in the two-way interaction, two of these OTUs decreased compared to the expected additive model.
The GLM model, however, also shows that three-way interactions between stressors are generally synergistic when considering changes in taxa abundances (red-colored box, Fig. 4). The three-way interaction takes into account the main effects and each of the two-way interactions. The null model predicts that the addition of a third stressor has no effect on the interaction of the other two stressors. Compared to this null model, the three-way interactions of high temperature, scarring, and nutrients resulted in an average 19.55 log2 fold increase in taxa (red-colored line, Fig. 4). Numerous differently abundant OTUs in the two-way and three-way interactions have a 30 log-fold change. These results may, however, be a caveat of the DESeq2 method to calculate a change in abundance associated with the presence of a taxon that was formerly absent in the treatment contrast. Notably, the magnitude of the average log2 fold-change increases with increasing number of stressors (yellow-, purple-, and red-colored lines, Fig. 4), which likely results from reduced power from the consecutive addition of interaction terms, thereby, requiring a larger change for a significant statistical result. Regardless of these differences in the magnitude of the response, the evaluation of factorial interactions with a GLM agree with results from the dominant taxa ( Supplementary Tables 1 and 2) and Simpson's Index (Supplementary Table 3) in characterizing double stressors as antagonistic and triple stressors as synergistic.

Discussion
Contrary to our hypothesis, our overall results suggest the global and local stressors tested in this tank experiment generally do not act synergistically to induce dysbiosis in the coral microbiome of Pocillopora meandrina. In fact, we find that the biological response in the microbial community to stress does not scale positively with increasing number of stressors. We predicted that when sequentially adding stressors to the system, we would see a concurrent increase in deterministic changes to the microbiome (Fig. 5a,c,e). For beta diversity, deterministic changes would produce clusters with increasingly distant locations from the control community (Fig. 5a). Stochastic changes would likewise produce communities that were more dispersed or variable (Fig. 5c). Instead, we found that the greatest deterministic changes in the microbial community resulted from single stressors, while interactions produced an intermediate level of change resulting in antagonisms that decreased the individual effects (Fig. 5b). For stochastic changes, any environmental stressor was sufficient to induce dispersion around the centroid of healthy corals (Fig. 5d), although this dispersion was likely a result of the single dominant taxon (Fig. 1). The changes in alpha diversity, however, did not scale positively with the number of stressors, and single stressors appear to increase community diversity more than two or three stressors combined (Fig. 5e,f). We also found that stress induced a myriad of opportunists to invade the community, shifting species dominance away from coral symbionts. The dynamics observed in species' abundance profiles of the microbial community following a perturbation may be explained by each particular microbes' nutrient preference and competitive ability 31 .
In contrast to more heterogeneous communities that may be more robust to changes in community evenness, the control corals were dominated by a single taxon initially and thus exhibited low evenness. We would predict then that any perturbation to the system would only increase evenness reflected in higher Simpson's diversity (e.g., Fig. 2). As such, when stressors were applied to the coral host, bacterial community evenness increased when the dominance shifted from the OTU-Endo symbiont to other taxa such as Desulfovibrionaceae 8,19,32 , Enterobacteraceae 20,23 , Amoebophilaceae 33,34 , Moraxellaceae 35,36 , and Rhodobacteraceae 8,16,37,38 (Fig. 1). Contrary to previous work 8 , we did not see an increase in species richness with stress. This may be a result of the mesocosm tanks restricting natural presence/absence dynamics on the reef. Instead, the results suggest a reshuffling of microbial members rather than an increase of new species. Many microbiome studies seek to understand whether dysbiosis, or an imbalance in microbiota, is marked by an invasion or proliferation of pathogens or a depletion of beneficial bacteria 39,40 . Yet taxon relative abundance measures alone do not provide enough information to answer this question.
Individual responses of taxa to stress and their contribution to microbiome dysbiosis can be assessed with differential abundance analysis 37,41,42 . Unlike diversity measures which are driven by the changes of dominant taxa in the community, differential abundance testing can identify changes in minor players in the community. DESeq2 can be used to model the abundance of each taxon independently, while accounting for the discrete positive nature of count data and the compositionality of the community using a generalized linear model (GLM) 43 . Using the linear model framework, we expected main effects to increase opportunistic bacteria, and interactions to produce synergistic effects as the community becomes increasingly compromised (Fig. 5g). Instead, we found that two-way interactions produced antagonistic responses among opportunistic taxa. This apparent antagonism may be a result of a dominance effect, in which one stressor accounts for most or all of the biological response, changing susceptible taxa such that the second stressor has no additional effect 5 . Stressors may provide some benefit or resource that normally limits the abundance of opportunistic taxa. For instance, high temperature may increase bacterial reproduction and metabolic rate, while scarring may increase free nutrients in the form of amino acids or open niches. These results suggest that opportunists such as Rhodobacter or Desulfovibrio spp. are not co-limited by the resources provided by high temperature and scarring. For instance, opportunistic taxa may be proliferating at such a high rate due to increased temperature and increased reproductive rates, that additional free nutrients from scarring do not compound the effect. In contrast, three-way interactions resulted in synergies as invading taxa continued to increase in abundance (Fig. 5h). This suggests that opportunistic taxa that had maximized their biological response under two stressors, were in fact co-limited by some resource provided by www.nature.com/scientificreports www.nature.com/scientificreports/ a third stressor. For instance, the addition of nitrogen may have allowed some opportunistic taxa to surpass the maximum threshold of reproductive or metabolic potential set by high temperature and scarring. Alternatively, the difference in interaction type between two-way and three-way interactions may be a result from the coral host's compromised immune system [44][45][46] . The coral host effectively regulates it's associated microbial community under two stressors with a heightened immune response. However, with the addition of a third stressor, innate immunity could be overwhelmed, and the host could no longer regulate its community, thereby allowing a synergistic proliferation of opportunists.
Despite the current bias in interaction literature toward identifying synergies 5 , our study highlights multivariate and univariate statistical tools that can be combined with standard methods in microbial ecology to more accurately characterize interaction types to host-microbiome systems. Community diversity measures are standardly conducted in microbiome research 14 , however, they have rarely been used to explicitly test for antagonisms or synergisms between environmental stressors using a microbiome dataset 7,9,12,13,37 . Although there is  Location (a,b) and dispersion (c,d) effects represent measures of beta diversity. Community evenness (e,f) represents patterns in Simpson's index. Patterns in taxa differential abundance in log2FoldChange using the generalized linear model framework are displayed in g and h with the gray line denoting no effect or no interaction. Colors represent the type of stressor combination applied to the corals: none = teal, single = yellow, double = purple, triple = red.
www.nature.com/scientificreports www.nature.com/scientificreports/ no evidence that these measures respond linearly to stress, these analyses revealed unexpected patterns of community response to increasing amounts of stress. This study presents an initial evaluation of the utility of these community diversity measures in characterizing interactions between different combinations of stressors that are known to damage the coral host and produce compositional changes in its microbiome.

Methods experimental design and sampling of coral tissue for microbial analysis.
To test the individual and interactive effects of increased temperature, nutrient enrichment, and predation on coral microbiomes, a fully crossed tank experiment (Supplementary Fig. S1) was conducted at the University of California Gump Research Station (17°29′26.04″S, 149°49′35.10″W) in Mo' orea, French Polynesia. The experiment was conducted in September of 2016 using twelve independent 150 L flow-through, temperature controlled mesocosms with natural locally-sourced sea water at a flow rate of 448.1 ± 24.1 mL per minute and run under a 12:12 light:dark cycle with ~700 μmol quanta per m² per second. Each of the twelve tanks was independently regulated for temperature (chilling loop and heater) and for light. Ten Pocillopora meandrina colonies were collected from 3-4 m on the Mo' orea north shore fore reef and transported by boat to the Gump Research Station to be immediately fragmented. Each colony was fragmented into twelve coral nubbins (~5 cm in length), epoxied with Z-Spar to vertically stand on plastic mesh, and distributed to each of mesocosm tanks at 26 °C ± 1 °C for a total of ten nubbins per tank.
After a 24-hour acclimation period, five (half) of the nubbins in each tank were randomly selected and mechanically scarred on the branch tip with 8 mm snub nose pliers. Therefore, each tank contained two treatments: control or scarred with a specific temperature and nutrient regime. (Supplementary Fig. S1). Pliers were sterilized between each scarring to prevent cross-contamination. The resulting scars were ~2 mm deep and removed the tissue layer and portions of the skeletal structure. These scars are representative of corallivory by parrotfishes and pufferfishes, two of the most common scraping corallivores on these reefs. Half of the tanks were then randomly selected and heated to 29 °C (ca. < 1 °C change per hour) to mimic temperatures associated with thermal stress on Mo' orean reefs 47 . Next, each tank was assigned one of three nutrient treatments: pulse of 4 µM of nitrate (NO 3 − ), pulse of 4 µM of ammonium (NH 4 + ), or no enrichment controls. For 21 days, each enriched tank was spiked with either nitrate or ammonium twice a day every 12 hr. Flow in all mesocosms was ceased for an hour immediately following the nutrient additions. This resulted in a twice daily hour-long nutrient pulse, followed by 5 hours of steady dilution and 6 hours of ambient concentration.
For microbial analysis, a subset of the nubbins (n = 6 per treatment) were randomly selected while controlling for source tank and colony to minimize samples while ensuring sufficient replication after 21 days of exposure to stressors (Supplementary Fig. S1, Supplementary Table S9). For microbial analysis, a healthy branch tip was clipped off of each unscarred nubbin and the branch tip around the scar was collected for scarred samples. These samples were frozen and shipped to Oregon State University.
16S library preparation and sequencing, quality control, and initial data processing. DNA was extracted using the MoBio Powersoil ® DNA Isolation Kit according to the manufacturer's recommendations.
Polymerase Chain Reaction (PCR) was performed on the V4 region of the 16S rRNA gene using the primer pair 515F (5′-GTG CCA GCM GCC GCG GTA A-3′) and 806Rb (5′-GGA CTA CHV GGG TWT CTA AT-3′) that targets bacterial and archaeal communities 48,49 . Amplicons were barcoded with barcoding primers with Nextera adapters, pooled in equal volumes for sequencing 50 , and purified with AMPure XP beads. Amplicon pools were paired-end sequenced on the Illumina MiSeq sequencing platform, 2 × 300 bp end version 3 chemistry according to the manufacturer's specifications at the Oregon State University's Center for Genome Research and Biocomputing (CGRB) Core Laboratories. QIIME (v1.9) 51 was used for quality control and selection of operational taxonomic units (OTUs). Demultiplexed raw reads were trimmed of adapter and primer sequences and pair-end sequences merged. Sequences with a total expected error or total sum of the error probabilities >1 for all bases were discarded. Chimeras were removed and 97%-similarity OTUs were picked using USEARCH 6.1 52 , the 97% GreenGenes 13_8 reference database 53 , and QIIME's subsampled open-reference OTU-picking protocol 54 . In this protocol, sequences that failed to hit the reference collection are randomly subsampled for de novo clustering. Taxonomy was assigned using UCLUST, and reads were aligned against the GreenGenes database using PyNAST 55 . The GreenGenes reference database is commonly used in microbiome analysis 14 and has been validated on animal and coral microbiomes in numerous studies 33,[56][57][58] . FastTreeMP 59 was used to create a bacterial phylogeny with constraints defined by the GreenGenes reference phylogeny.
Both a rarefied and unrarefied OTU table were created for downstream analyses. First, OTUs were filtered out of the starting table if their representative sequences failed to align with PyNAST to the GreenGenes database or if they were annotated as mitochondrial or chloroplast sequences. After this step, the number of reads per sample ranged from 1 to 87262 with a median of 9742 per sample and 3383 unique reads. OTUs with less than 100 counts were then removed from the OTU table resulting in a total of 430 unique reads, and a histogram of reads by samples was plotted ( Supplementary Fig. S4). We did not find that any one OTU was associated with one particular sample. Samples with less than 1000 reads were considered undersampled and discarded from the table resulting in an unbalanced experimental design (Supplementary Table S9 www.nature.com/scientificreports www.nature.com/scientificreports/ statistical analyses investigating alpha and beta diversity. To determine how dominant taxa within the community respond to stressors, the rarefied OTU table was first transformed to proportions to identify the five OTUs with the highest relative abundances. These OTUs were then evaluated for changes with treatment using generalized linear mixed-effects models (GLMM) which allow for inclusion of random effects. Raw OTU counts from the unrarefied table were then included as the response variable in the GLMM with an offset by the log of total sequencing depth to reflect relative abundances. Offset variables are commonly used in count models to adjust for differential exposure times and, for these purposes, to adjust for different sequencing depths. Each OTU count distribution was evaluated for normality with quantile-quantile plots and the Shapiro-Wilk test for normality 66 . A logistic regression with unrarefied count data with an offset was selected over a linear regression with proportion data from the rarefied table since the data were not normally distributed and evidence suggests logistic regressions perform better than arcsine-transformed data in linear regressions 67 . For all models, temperature, nutrients, and scarring were assessed as fixed effects and factorial interaction terms and tank and colony as separate random effects. For the most abundant OTU in the dataset, OTU-Endo, a poisson distribution was used since it resulted in normal residuals using lme4 68 (v1.1.15) and lmerTest (v2.0.36) to obtain P-values. For the remaining four OTUs a zero-inflated negative binomial regression was used to account for the excess of zeros in the data using glmmTMB (v0.2.2.0). The negative binomial distribution was also evaluated based on normality of residuals. A single zero-inflation parameter was modeled for all observations.
To improve normality of alpha diversity metrics, Chao1 and Faith's phylogenetic distance were log-transformed, while Simpson's index was arcsine-transformed. Stressor effects on each alpha diversity metric were assessed with linear mixed effect models (LMM) using lme4 with temperature, nutrients, and scarring as fixed effects and factorial interaction terms and tank and colony as separate random effects. P-values were approximated with the lmerTest (v2.0.36).
For beta diversity analyses, the Weighted Unifrac dissimilarity matrix generated in phyloseq (v1.23.1) was used to test for location and dispersion effects acting on the microbial community. First, differences between treatments were assessed with permutational multivariate analysis of variance (PERMANOVA) using the adonis function in the package vegan (v2.4.6) both with treatment as the single factor and with pairwise comparisons between treatments 69 . Next, permutation tests for homogeneity in multivariate dispersion (PERMDISP) were performed using the betadisper function in the package vegan both with treatment as the single factor and with pairwise comparisons between treatments 70 . For visualizing differences in beta diversity, mean distance-to-centroid values by treatment were extracted from the betadisper results, and mean separation from communities in other treatments or between-group distances were averaged by treatment from the pairwise dissimilarity matrix generated in phyloseq. Together, PERMANOVA and PERMDISP provide a comprehensive analysis of deterministic and stochastic changes to the microbiome by evaluating between-group differences and within-group dispersions, respectively. For comparability, the two tests were run with treatment as the single factor, since PERMDISP does not allow for a specified model formula. Pairwise analysis of variance for the PERMANOVA test were conducted using a modified version of pairwiseAdonis 71 (v0.0.1), and for the PERMDISP test with the permutest command in vegan. Distance-to-centroid was also regressed against the arcsine-transformed relative abundance of OTU-Endo to determine the correlation between community dispersion and the dominant taxon. Lastly, stressor main effects and interactions were evaluated with PERMANOVA by fitting a linear model to the distance matrices using factorial interaction terms. All diversity analyses were conducted in R (v3.4.0). Differential OTU abundance analysis. To analyze differences in abundance of bacterial taxa across stressors and stressor interactions, a negative binomial generalized linear model (GLM) was fitted with the R package DESeq2 (v1.16.1) using the unrarefied OTU table 43 . DESeq2 does not support random effects and therefore a mixed model was not used. All rare taxa that were present in fewer than 15 samples were excluded from this table to create a pre-filtered, unrarefied OTU table. In the DESeq2 method, raw counts are modeled with a negative binomial distribution which is commonly used for overdispersed count data 72 . Additionally, "size factors" or normalization factors are calculated with a median-of-ratios method to normalize differences in sequencing depth between samples 43 . The GLM design specified nutrient regime, temperature, scarring, and their interactions as factors with beta priors set to false. Wald post-hoc tests were used to identify factors in the model that significantly affected each taxon compared to the control level by building a results table from a specified treatment contrast. To control the rate of false positives due to multiple comparisons, differentially abundant taxa within each treatment contrast were identified as significant with Benjamini-Hochberg FDR p-values less than 0.05 73 . The pre-filtered, unrarefied OTU table was summarized to Family and Genus level using the tax_glom command in phyloseq. The DESeq2 method was then applied to these two tables to demonstrate accordance between significant changes at the OTU level with significant changes at higher taxonomic levels.