Decline of a distinct coral reef holobiont community under ocean acidification

Background Microbes play vital roles across coral reefs both in the environment and inside and upon macrobes (holobionts), where they support critical functions such as nutrition and immune system modulation. These roles highlight the potential ecosystem-level importance of microbes, yet most knowledge of microbial functions on reefs is derived from a small set of holobionts such as corals and sponges. Declining seawater pH — an important global coral reef stressor — can cause ecosystem-level change on coral reefs, providing an opportunity to study the role of microbes at this scale. We use an in situ experimental approach to test the hypothesis that under such ocean acidification (OA), known shifts among macrobe trophic and functional groups may drive a general ecosystem-level response extending across macrobes and microbes, leading to reduced distinctness between the benthic holobiont community microbiome and the environmental microbiome. Results We test this hypothesis using genetic and chemical data from benthic coral reef community holobionts sampled across a pH gradient from CO2 seeps in Papua New Guinea. We find support for our hypothesis; under OA, the microbiome and metabolome of the benthic holobiont community become less compositionally distinct from the sediment microbiome and metabolome, suggesting that benthic macrobe communities are colonised by environmental microbes to a higher degree under OA conditions. We also find a simplification and homogenisation of the benthic photosynthetic community, and an increased abundance of fleshy macroalgae, consistent with previously observed reef microbialisation. Conclusions We demonstrate a novel structural shift in coral reefs involving macrobes and microbes: that the microbiome of the benthic holobiont community becomes less distinct from the sediment microbiome under OA. Our findings suggest that microbialisation and the disruption of macrobe trophic networks are interwoven general responses to environmental stress, pointing towards a universal, undesirable, and measurable form of ecosystem change. Video Abstract Supplementary Information The online version contains supplementary material available at 10.1186/s40168-023-01683-y.


Background
Increased atmospheric greenhouse gas concentrations are leading to increased partial pressure of CO 2 (pCO 2 ) and reduced pH in the surface water of the oceans, which have absorbed 25% of all anthropogenic CO 2 to date [1].The impacts of this phenomenon, known as ocean acidification (OA) [2], are predicted to be particularly severe for coral reefs due to declines in net calcification by organisms, which impacts reef structure and diversity [3].Severe impacts on coral reefs are concerning as these ecosystems host vast biodiversity and provide significant ecosystem services to humanity (e.g.coastal protection, food security, and new medicines), and other significant terrestrial and pelagic ocean ecosystems are also directly dependent on their functioning [4].Being able to predict how coral ecosystems change under OA, and what this means for their ability to continue providing ecosystem services to society, is therefore a global environmental and social priority.
Microbes play important roles on coral reefs, from carrying out nutrient cycling, which contributes to ecosystem productivity in nutrient poor waters [5], to their roles in immunity and defence for a wide range of reef invertebrates, including cnidarians, sponges, molluscs, and echinoderms [6][7][8][9].These interactions between macrobes and microbes are known to shift in response to environmental change, such as OA.Such shifts include changing microbial associations with particular macrobes, such as corals [10,11] and sponges [12], and changes in the functional profile of specific microbiomes, such as alterations in nitrogen (N 2 ) fixation by coral-associated bacteria [13] and shifting metabolic activity of free-living bacteria in the water column [14].However, the traditional focus of research on individual macrobes and their microbiomes (holobionts) makes it challenging to scale up our understanding to ecosystemlevel shifts in macrobe-microbe interactions, which, though often overlooked, can play a significant role in driving ecosystem-level change [11,15,16].
Ecosystem-level impacts of OA are expected to result from interaction-mediated changes at the community level, driven by different physiochemical effects of OA at the level of organism metabolism [2,17].Wellestablished organism-level effects of OA include benefits to fleshy algae and other photosynthetic organisms from the resource effect of increased pCO 2 [2,18], which can result in enhanced net dissolved organic carbon (DOC) release [19].In contrast, calcifiers suffer from increased costs of calcification due to reduced pH [2].These differentiated organism-level impacts can shift ecological interactions between taxa [2,17], with cascading effects on energy flows through an ecosystem via altered nutrition and metabolism (i.e.altered trophodynamics [20,21]).Such indirect cascading effects can impact entire trophic networks [22], including macrobe-microbe trophic interactions [23].However, our understanding of the indirect effects of ocean acidification on coral reefs remains limited [24].
Here, we propose that a previously observed ecosystem-level impact of specific stressors on coral reefs -microbialisation -may be generalisable to OA. Microbialisation refers to an increase in microbial biomass resulting from a reallocation of energy from macrobes to microbes [23,25].On coral reefs, the proximal causes of microbialisation have been proposed to be overfishing and eutrophication, which facilitate the enhanced growth of fleshy algae and cause an increased release of dissolved organic carbon (DOC) [26].Elevated DOC has been proposed to increase microbial biomass and disease (the DDAM (DOC, disease, algae, microbes) positive feedback loop) [26].
We build on a commonality between proposed organism-level mechanisms of microbialisation, namely that macrobe communities are expected to be more vulnerable to stressors than microbial communities [26,27], an observation both supported in general [28] and in the case of OA in particular [29].This greater vulnerability of macrobes should lead to declines in some benthic taxa, such as calcified algae, soft and hard corals, sponges [18], and calcified grazers (e.g.gastropods [30]), already documented with OA, and to the disruption of associated macrobe trophic pathways.Examples of these trophic pathways include calcified grazers consuming algal communities (which controls algal proliferation) [31], and sponges removing vast quantities of dissolved organic matter (DOM) from the water column and converting it into food for higher trophic levels (e.g.polychaetes and brittle stars) via the sponge loop [32].With declines in benthic taxa, some of the free energy cycled through such pathways can become available to less impacted taxa, such as environmental microbes [20,33], and is expected to drive the ecological release of such taxa in both density and niche expansion [34].
We hypothesise that we will observe a decline in the compositional and functional distinctness between the benthic holobiont community microbiome and the environmental sediment-dwelling microbiome as a result of microbialisation occurring under OA.Almost all macrobes are holobionts with a symbiotic microbiome [35,36], and therefore, microbialisation has the potential to impact the microbiome of the entire holobiont community (recently referred to as the eco-holobiont [37]).The process of microbialisation should result in decreased compositional and functional distinctness between the benthic holobiont community microbiome and the sediment microbiome through two mechanisms.Firstly, direct impacts on macrobes may alter host metabolism and reduce the resources or habitat available to the holobiont community microbiome, leading to opportunistic environmental microbes displacing holobiont-specialised microbes (e.g.[38]).Secondly, increased abundances of environmental microbes and increased microbe trophic interactions with macrobes (expected due to the loss of macrobe competitors [39]) should lead to increased opportunities for colonisation of the holobiome by environmental microbes [40], including those in the sediment microbiome.
To test whether this decline in the distinctness of the benthic holobiont community microbiome occurs, we generated a unique multiomic dataset from autonomous reef monitoring structures (ARMS) deployed on a natural OA gradient caused by CO 2 seeps.ARMS are threedimensional, artificial settlement structures designed to mimic the structural complexity of coral reef environments, which are increasingly used to monitor coral reefs across the globe (e.g.[41][42][43]).They enable the non-destructive and standardised sampling of a large proportion of reef diversity that is often not studied, including algae and cryptic benthic invertebrates such as sponges, cnidarians, bryozoans, and annelids [44], alongside their associated microbes, and sediment-dwelling environmental microbes.ARMS allow us to investigate OA using a holistic microbial ecosystem approach that integrates across scales from individual microbes and benthic holobionts, to neighbouring holobionts that, in turn, interact with and influence successively larger and more complex communities.We studied the effects of OA using an in situ experimental approach at a location with naturally occurring CO 2 seeps that produce pH and pCO 2 gradients.These seeps have been intensively studied because they can help predict future ocean conditions under OA [45].
Multiomics, in this case metabarcoding and metabolomics, provides a powerful toolkit to investigate the effects of stressors across entire communities.Metabarcoding provides data on genetic community composition and diversity (i.e.compositional metrics) [46], while metabolomics provides comparable data on biochemical composition and diversity of the metabolome (i.e.functional metrics) [47].We first confirm that the expected photosynthetic community shifts take place (as previously documented [18]), consistent with microbialisation.We then analyse the ecosystem-level effects of OA by comparing the distinctness of the sediment microbiome to: (1) the benthic holobiont community microbiome and (2) individual sponge microbiomes.We expected OA to cause a decline in distinctness as a result of ecosystem microbialisation.

Experimental design and sampling
This study was carried out at CO 2 seeps and adjacent control sites in Milne Bay Province, Papua New Guinea (Fig. 1), located at 9° s latitude in the heart of the Coral Triangle.The two studied seep localities (Upa-Upasina and Dobu) are located along an active tectonic fault where > 99% CO 2 gas has been streaming though the reef substrata at ambient temperature (28.6-29.7 °C) for at least 100 years and probably much longer [18].The reefs surrounding the seeps are under low anthropogenic pressure and have been used to study ocean acidification for the last decade (e.g.[18,48]).The six study sites (two localities, each with three pH levels) exhibit similar geomorphology, temperature and salinity, but contrasting pCO 2 and pH [18].Water temperature, pH, salinity, and pressure at the study sites have been monitored regularly (2010-2016), making this an ideal location to study the isolated impacts of OA.
In April, 2012, eighteen autonomous reef monitoring structures (ARMS; Fig. 2A) were deployed at 3 m depth adjacent to coral reefs at 3 pH levels at each of the 2 localities (mean pH: control 7.99 & 8.01, medium 7.85 at both localities, and low 7.64 & 7.75; n = 3 per pH level [48]).ARMS were collected from the seafloor (Fig. 2B) after 31 months, in November 2014.A 106-μm nitex-lined crate was placed over each ARMS on the seafloor, and they were together returned to the surface, after which each ARMS was placed in an individual holding tank with 45-μm filtered aerated seawater.ARMS were then transported to shore where they were sampled rapidly to minimise molecular and chemical degradation; transportation time of ARMS was < 20 min, and processing of each sample type was completed within 1.5 h.

Sample extraction and multiomics
The standard ARMS processing protocol [44] was modified to test our specific hypothesis.From each ARMS unit, five fractions were collected: the benthic photosynthetic community, the benthic holobiont community, the sediment, Halisarca sp.sponge, and Tethya sp.sponge.To do this, ARMS were removed from their holding tanks, and the 9 plates (17 plate surfaces as there is no accessible bottom surface to the bottom plate) were separated and rinsed to dislodge loosely attached organisms.The water and previously trapped sediment in the holding tank were retained.
First, the benthic photosynthetic community was sampled by randomly subsampling (4 × 1 cm 2 ) the top surface of the top ARMS plate (e.g. Figure 2C), which resembles the algal community found on exposed rocky substrates (e.g.macroalgae, algal turf, calcified algae, and cyanobacterial mats).Second, the two sponge fractions (Halisarca sp. and Tethya sp., e.g. Figure 2D) were generated by sampling morphologically identified sponges from the internal plates.Both are low-microbial abundance sponges [49,50], and Tethya sp. was only found at Upa-Upasina.Thirdly, the benthic holobiont community fraction was generated by scraping, and blending the scrapings from all 17 surfaces, including the remainder of the lightexposed top surface of the top plate.While this fraction will include both photosynthetic and non-photosynthetic organisms, the shaded plate surfaces constitute ~ 16 × the area of the light-exposed top surface of the top plate.The most abundant phyla found in this fraction on other reefs around the world are fairly consistent and include the Porifera, Cnidaria, Bryozoa, Chordata (Ascidiacea), and Annelida [42][43][44].An example of an ARMS plate from which this fraction was collected can be seen in Fig. 2E.This homogenised bulk sample was then subsampled (50 ml).Finally, the sediment fraction was generated by passing the water and sediment from the holding tank through a 500 μm and then a 100-μm sieve and collecting the material which did not pass through the 100-μm sieve.This drained sediment sample was subsampled (10 g).This fraction was therefore primarily composed of sediment, microbes, including free-living sediment dwelling microbes (e.g.bacteria) and microplankton (e.g.single-celled algae such as diatoms and dinoflagellates).From each ARMS unit, each fraction was split in two; one-half was snap frozen for metabolomic analysis by being dropped in liquid nitrogen in a dry shipper, and the other was placed in RNA later for metabarcoding.All samples were returned to the Smithsonian National Museum of Natural History (Washington, DC, USA) in a liquid nitrogen dry shipper.Total DNA was extracted from 10 g of the benthic holobiont community samples and 5 g of the sediment samples using a MO-BIO PowerMax Soil DNA Isolation Kit according to the manufacturer's protocol with the addition of 400 μg/ ml proteinase K and an overnight lysis step at 56 °C and 200 rpm.The benthic photosynthetic community and sponge subsamples were extracted with the DNeasy Blood and Tissue Kit (Qiagen), according to the manufacturer's instructions.All DNA extracts were purified using MO-BIO PowerClean DNA Clean-Up Kits, quantified Qubit dsDNA HS Kit, run on an agarose gel, and DNA quality investigated using ImageJ software.All sample types are known to contain relatively high bacterial biomass, and thus, we do not expect contamination from the lab environment or equipment to be a major issue in DNA libraries.However, extraction and PCR amplification controls were included for all sample types; these were all negative and so were not sequenced.Each sample was analysed with 16S rRNA gene metabarcoding (for the microbe community) and mass spectrometry (for metabolomics).All 16S rRNA gene libraries were prepared for sequencing using the original Earth Microbiome Project protocol using primers 515 F and 806 R, which are designed to amplify prokaryotes (Bacteria and Archaea) [51].To investigate the benthic photosynthetic community, 23S rRNA gene libraries were also prepared using the protocol described by Marcelino and Verbruggen [52] using a two-step PCR procedure that first amplifies the gene fragment followed by ligation of the barcoded Illumina adaptors to the amplicons in a second PCR reaction; this protocol is designed to target both eukaryotic algae and cyanobacteria.
Metabolites were extracted from all fractions in 70% methanol [47].The 70% methanol extraction was chosen to select for slightly polar molecules, encompassing a broad range of the chemosphere [53].Metabolites were separated and identified via liquid chromatography-tandem mass spectrometry using a Bruker Daltonics Maxis qTOF mass spectrometer equipped with a standard electrospray ionisation source according to the methods of Quinn and colleagues [53] .Briefly, the mobile phase was pumped through a Kinetex 2.6 μm C18 (30 × 2.10 mm) ultra-performance liquid chromatography (UPLC) column for a 15-min run.The resulting LC-MS/MS data files were processed through the MZmine2 workflow.The subsequent metabolite feature table was then processed through the GNPS feature-based molecular networking workflow with the default parameters, except that a minimum cosine of 0.65 and a minimum matched peaks of 4 were used for network construction.

Bioinformatics
A bioinformatic pipeline was implemented in R for the 16S and 23S rRNA gene libraries.Amplicon sequence variants (ASVs) were generated from raw sequencing data using the Divisive Amplicon Denoising Algorithm (DADA2 v1.24.0 [54]).Reads were quality filtered to maintain Q30 scores while maintaining at least 50 base pair overlap and removing any base pair below Q2 [55].Default maxEE (2) and truncQ (2) parameters were used, 16S rRNA gene sequences were truncated at a length of 150 base pairs on both strands, 23S rRNA gene sequences had the first 20 base pairs (nonbiological primers) trimmed from both strands and were truncated at a length of 249 base pairs on the forward strand and 212 base pairs on the reverse strand (see Table S1 for numbers of sequences passing denoising steps).Taxonomy was assigned using the DECIPHER v2.24.0 R package [56] -which has been shown to have higher accuracy than popular classifiers including BLAST and the RDP classifier [56] -and the GTDB (16S rRNA gene [57]) and microgreen (23S r RNA gene [58]) databases.16S rRNA gene ASVs identified as plastids were subsequently removed.Samples were not rarefied as part of bioinformatic processing [59] but only when required for specific statistical analyses (see ASV-level Shannon diversity below).

Statistical analyses
PERMANOVAs were performed for each fraction separately, with the 16S rRNA gene sequencing and metabolomic data subdivided into community fractions (benthic photosynthetic community, benthic holobiont community and sediment) and organism fractions (Halisarca sp. and Tethya sp.sponges), resulting in 10 PERMANOVAs (Table S2).Prior to fitting PERMANOVAs, a multivariate analogue of Levene's test for homogeneity of variances (betadisper) was applied to ensure PERMANOVA tests could be applied.PERMANOVAs fit locality and ordinal pH as explanatory variables, and locality was treated as a blocking factor, except in the case of the Tethya sp.sponge which was only found at one locality and so was fit with pH as the only explanatory variable.An eleventh PERMANOVA was run on the 23S rRNA gene data, following the same approach.Bonferroni corrections were applied to all p-values obtained from the PERMANOVAs to account for multiple testing.All PERMANOVAs and supporting NMDS visualisations were based on Morisita dissimilarities between samples, as they have been shown to be most reliable in the case of under sampling [60].
Total ASV richness and phylum level Shannon diversity, each accounting for unobserved ASVs, were estimated for each metabarcoding sample using a breakaway model [61] and a DivNet model, treating all samples as independent observations, respectively [62].The estimated richness and estimated Shannon diversity of all metabarcoding samples were then modelled using a single betta hierarchical mixed model for each metric (including all fractions).This modelling approach was chosen to account for explanatory variables, richness variance, and richness estimation error [61].Fraction, ordinal pH, and the interaction of fraction and pH were included as fixed effects and locality as a random effect.Compound richness and Shannon diversity were calculated from untransformed data for metabolomic samples from all fractions.Each metric was modelled with a single linear mixed model.Fraction, pH, and the interaction of fraction and ordinal pH were included as fixed effects and locality as a random effect.In addition, ASV level Shannon diversity was calculated (no statistical estimation procedure was applied) for all samples rarefied to even depth (n = 50,000) to test the sensitivity of the results found for phylum level Shannon diversity.
The change in abundance of phyla and metabolites with OA was analysed using DeSeq2 differential abundance analysis with a negative binomial distribution [63] .Differential abundance and dispersions were calculated for each community fraction (benthic photosynthetic community, benthic holobiont community, and sediment) and multiomic analysis type separately using a DESeq2 design formula with variables of locality and pH.This enabled change within each community fraction to be examined.However, abundance and dispersions were calculated for both sponge fractions together using a design formula with variables of species, locality, and pH.This enabled shared change occurring across sponges under OA to be examined.Wald significance tests were conducted for changes in differential abundance under OA, with a parametric fit of dispersions [63].
Microbiome/metabolome distinctness was calculated for each ARMS as the proportion of unique sequences found within the benthic holobiont community microbiome/ metabolome which were not also found in the sediment microbiome/metabolome.Microbiome/metabolome distinctness was modelled using a linear mixed model with ordinal pH as a fixed effect and locality as a random effect.A likelihood ratio test was used to infer the significance of ordinal pH as a fixed effect.Note that this analysis was not conducted for the benthic photosynthetic community as we are testing whether distinctness is reduced for the general community of macro-organisms, and the benthic holobiont community is a more general sample including both photosynthetic and non-photosynthetic organisms.The same approach was taken to calculate and model microbiome/ metabolome distinctness for individual holobionts with the additional random effects of (i) ARMS identity (nested within locality), as multiple individual holobionts were collected from the same ARMS, and (ii) sponge species.
Benthic holobiont community microbiome distinctness from the sediment microbiome was also modelled with a modified mixture Sloan neutral community model (MSNCM [64]).This additional modelling approach captures the contribution of each of the sediment and benthic holobiont community microbiome metacommunities to the composition of benthic holobiont community microbiomes from individual ARMS, thus providing an alternative abundance-based test of whether benthic holobiont community microbiomes become more distinct from the sediment microbiome under OA.The original Sloan neutral community model describes the frequency of occurrence of ASVs in a community as a function of their abundance in the metacommunity, with a single free parameter (m: migration) which can be interpreted as the probability of neutral dispersal or alternatively inverse dispersal limitation.The MSNCM used here models ASV frequency in sampled benthic holobiont community microbiomes from each pH regime as a function of its abundance in two metacommunities: (1) all benthic holobiont community microbiomes from the same pH regime as the sample and (2) all sediment microbiomes from the same pH regime as the sample.Each metacommunity is fit with its own migration/ inverse dispersal limitation parameter (mholo, menv), and a mixture parameter (mix) is fit describing the contribution of each metacommunity.The model is fit to samples from each pH regime separately using non-linear least-squares fitting as detailed in Burns et al. [65] .
Fifteen 23S rRNA gene metabarcoding libraries were generated across the pH gradient, to confirm the expected effect of OA on the benthic photosynthetic communities.The composition of the benthic photosynthetic community differed significantly by pH (F = 5.5, p < 0.05), with significant declines in phylum Shannon diversity (95% CI [-0.76, -0.55], p < 0.05) and ASV Shannon diversity (95% CI [-2.27, -0.33], p < 0.05) at lower pH.Lower pH was associated with significantly increased differential abundance of the dominant phylum Ochrophyta (of which 99.7% of reads were from the class Phaeophyceae, and 71.4% were from the genus Sargassum; Fig. 3; Figure S1).
Forty-seven 16S rRNA gene metabarcoding libraries were generated from the microbiomes of the two sponge species: 17 from Halisarca sp.individuals and 30 from Tethya sp.individuals.Microbiome composition was only significantly different across the pH gradient for Tethya sp.sponges (F = 5.9, p < 0.05; Table S2; Figure S2).While there was no significant effect of pH on ASV richness or ASV level Shannon diversity, reduced pH was associated with significantly lower phylum level Shannon diversity in Tethya sp.sponge microbiomes (95% CI [− 0.44, − 0.12], p < 0.05), but no such effect was found in Halisarca sp.sponge microbiomes (Table S4; Figure S3; see Table S6 for a summary of all significant 16S rRNA patterns).At the phyla level, lower pH was associated with significantly higher abundances of Marinisomatota, SAR324, and Cyanobacteria and significantly lower abundances of Firmicutes, Desulfobacterota, and Bacteroidota (Fig. 3) in both sponge species microbiomes.

Biochemical diversity and composition
One-hundred and seven metabolome libraries were analysed from the same 18 ARMS.Each fraction (benthic photosynthetic community metabolome; benthic holobiont community metabolome, sediment metabolome, and sponge metabolomes) had between 18 and 34 samples across the pH gradient (Table 1).These samples produced 1211 compounds, of which 4.62% were identified, representing 6% of all molecules.
Fifty-five metabolome libraries were generated from the sediment metabolome, the benthic holobiont community, and the benthic photosynthetic community metabolomes (n = 18, 18, 19, respectively).There was no significant compositional difference in these community metabolomes at lower pH (Table S2; Figure S2).There was also no significant effect of pH on compound richness for any community metabolome (Table S7; Figure S3).However, Shannon diversity was significantly lower at lower pH in the benthic holobiont community metabolome (95% CI [− 0.07, − 1.18], p < 0.05) and significantly higher at lower pH in the sediment metabolome (95% CI [0.02, 0.80], p < 0.05; Table S8; Figure S3).Decreased pH was associated with significant differences in the abundance of several identified compounds (note that more than 95% of compounds were not identifiable).In the sediment metabolome, glycerophospholipids (lysophosphatidylcholines (LPCs) and phosphocholines) and pheophorbide A (a chlorophyll-derived compound) were less abundant, and beta-carotene was more abundant.In the benthic holobiont community benzene derivatives, chondramide B and mesoporphyrin IX (the latter two both anticarcinogens) were less abundant; a range of glycerophospholipids (again LPCs and phosphocholines), sucrose, and beta-carotene were more abundant.In the benthic photosynthetic community, pheophorbide A was more abundant, and various glycerophospholipids had significantly different abundances, with LPCs mostly having higher abundances and phosphocholines lower abundances (Fig. 4).Fifty-two metabolome libraries were generated from the two sponge species, with 18 from Halisarca sp. and 34 from Tethya sp.There was no significant compositional difference in the metabolomes of either sponge at different pH (Table S2; Figure S2).There was also no significant effect of pH on compound richness or Shannon diversity for the sponge metabolomes (Figure S3).Among identified compounds (note that more than 95% of compounds were not identifiable), decreased pH was associated with significantly different abundances in a variety of glycerophospholipids and a significantly increased abundance of a benzene derivate across both sponge species (Fig. 4).
The dispersal probability from the sediment microbiome into the benthic holobiont community microbiome, as measured by the MSNCM (menv weighted by the mixing parameter), was also higher at lower pH.The frequency of occurrence of ASVs across the benthic holobiont community microbiome samples was well described by their abundance in the benthic holobiont community and sediment metacommunities through the MSNCM.However, the fit of the model was poorer at medium and low pH (control: 0.64, medium: 0.29, low: 0.38).The ratio between mholo and menv (both weighted by the mixing parameter) was also lower at lower pH, reflecting an overall increased contribution of sediment microbial abundance to determining the microbiome composition of benthic holobiont communities under OA (Table 2).
Change in microbial and chemical distinctness of individual sponge holobionts (Tethya sp. and Halisarca sp.) with ocean acidification was also analysed to examine whether community-level patterns were observed in individual holobionts.Individual sponge holobiont microbiome distinctness from sediment was lower at lower pH (95% CI [− 7.87, − 1.75], p < 0.05; Fig. 5).This individual holobiont microbiome distinctness model explained 61% of variation (R 2 ) compared to 9% of variation (R 2 ) explained by the equivalent community-level microbiome distinctness model described above.There was no significant effect of pH on sponge metabolome distinctness.

Discussion
Our results are consistent with the simplification of the algal community seen under OA elsewhere [73].Our results further demonstrate that this simplification effect extends across the entire benthic photosynthetic community holobiome as algal-associated microbes also decline in Shannon diversity with OA.While we do not observe an overall compositional shift in the benthic photosynthetic community metabolome, implying that photosynthetic function is largely conserved through these changes, we do see shifts in specific metabolites.Increasing concentration of chlorophyll-derived pheophorbide A suggests increased chlorophyll turnover [74], consistent with the doubling of macroalgal benthic cover previously observed under OA at these [18].We also see significant increases in sucrose in the benthic holobiont community under OA, potentially related to sugar-enriched dissolved organic carbon (DOC) released by these algae.
Both findings are consistent with the initial steps in the DDAM mechanism of microbialisation [26], whereby increased DOC is released from fleshy algae, in this case Sargassum [19].We did not, however, test for the increased microbial biomass previously shown to result from increased DOC stimulating the microbial loop and causing a shift in ecosystem trophic structure [26,75].The main microbial phyla that we observed to increase in abundance in our community samples (i.e.not sponge samples) with low pH (Desulfobacterota and WOR-3) differ from those observed by Haas [26] (Alphaproteobacteria, Gammaproteobacteria, and Bacteroidetes).However, this difference in specific taxonomic patterns in response to different stressors, across different sites, and when sampling different substrates (water vs. sediment) is expected and is a key motivation for seeking to identify more general ecosystem-level effects of stressors, such as microbiome distinctness.
The data presented here provide the first evidence of declining holobiont community distinctness in microbes and metabolites under OA.Our results build on the evidence that OA changes holobiont microbiomes [76] by demonstrating a systematic decline of a distinct benthic holobiont community microbiome, such that it becomes more compositionally similar to the sediment microbiome.The composition of the sediment microbiome in our study is representative of surface sediments from the region [77] and global analyses of marine sediments [78].Specifically, the results of the MSNCM suggest that this effect may result from the benthic holobiont community being increasingly colonised by sediment microbes as pH decreases.A similar effect has previously been observed for specific sponge [79] and coral [80,81] holobionts in response to different environmental stressors, but has not previously been observed at the community level.While organism-level studies provide information about the response of key species (e.g.habitat builders) to environmental change, a holistic approach is needed to accurately evaluate and predict impacts on coral reefs.
The synthesis of knowledge across scales, from individual microbes and holobionts to ecosystem-wide communities and processes, has recently been called for by multiple authors [29,75,82].Autonomous reef monitoring structures (ARMS) provide a novel tool for taking this "nested ecosystem approach" and conducting in situ experiments.
We explain the decline in the distinctness of the benthic holobiont community from the sediment microbiome as being caused by increased opportunities for colonisation of benthic holobiont communities by environmental microbes due to microbialisation.However, individual macrobes under stress may also become less able to regulate their microbiomes, while colonisation opportunities remain constant.In extreme cases, this inability to regulate the microbiome can result in traumatic dysbiosis [80], a more heterogeneous microbiome (Anna Karenina principle, [83], and host death.However, it is unlikely that this process is the major contributor to the observed community-level effect seen here because macrobe community compositions have already shifted under OA conditions at these vent sites, with an increased dominance of taxa that are less impacted by the stressor [18,84].Therefore, the notion that the majority of the macrobe community is experiencing dysbiosis associated with acute organism-level stress seems unlikely.For example, some sponges are known to thrive under OA, and do not exhibit evidence of organismal stress [85][86][87].In this study, the two sponge holobionts individually analysed (Tethya sp. and Halisarca sp.) showed reduced distinctness of their microbiomes from the sediment microbiome under OA.However, neither showed evidence of increased compositional heterogeneity of their microbiomes as expected under dysbiosis by the Anna Karenina principle, which predicts organism-level stress to reduce the ability of macrobes to regulate their microbiomes [83].In addition, metabolomes for these two sponges do not become significantly less distinct from the sediment under OA, which is consistent with the hypothesis that these sponges were not under stress.
Colonisation of holobionts by environmental microbes may support the resilience of macrobe communities, as it has been shown to allow some hosts to acclimatise to new environmental conditions, for example by allowing the host to make use of changing energy sources, and facilitate greater adaptation than can be afforded by host phenotypic plasticity [88,89].Different degrees of microbial restructuring observed among different sponge species indicate that horizontal transmission differs between species, and these variations affect the ability of sponges to persist under OA conditions [12,90].Here, we find significant increases in Cyanobacteria associated with sponges under OA conditions, which can contribute > 50% of a sponge's carbon demand [91] and likely provide at least some sponge species with enhanced scope for growth in these seep environments [12].We also find significant increases in Desulfobacterota in the benthic photosynthetic community, holobiont community, and sediment; this phylum includes many organisms capable of reducing sulphur compounds [92,93].Only one of the two ocean vents (Dobu) is known to release hydrogen sulphide [18], so this does not explain the increase in Desulfobacterota observed across sites.We note that while we do not identify changes in the differential abundance of any metabolites with known roles in sulphur cycling, this is not evidence of their absence due to the low level (less than 5%) of identification of metabolites.While their role here is unknown, coral reefs are important hotspots of marine sulphur and increased sulphate reduction rates of marine microbial communities have been found to occur between a pH of 6 and 7 [94], suggesting that rates may increase with OA.As the marine sulphur cycle is a quintessential example of algal-bacterial interactions [95], it will be important for future studies to investigate the impact of algal-derived microbialisation on the marine sulphur cycle, especially as new components and pathways in the sulphur cycle are still being identified [96].
One would expect the dynamics of microbes and holobionts to be universal to all ecosystems [36,97,98], though they may emerge from different organism-level interactions.Therefore, microbialisation, and the observable property of declining holobiont community distinctness under environmental change, could represent a universal ecosystem stress response.Identifying such a general, undesirable response (microbialised ecosystems typically have lower intrinsic and use values [99]) and a metric of ecosystem change has clear benefits to policy and evaluation.For example, ecosystem change and the associated risk of ecosystem collapse are the underpinning concept leveraged for the IUCN Red List of Ecosystems [100], but defining collapse for each ecosystem individually is a time-consuming and contentiously value-laden task [101,102].Furthermore, as microbial communities respond rapidly to environmental change, microbial bioindicators could provide signatures of change with the speed and resolution to allow real-time responses by ecosystem managers.Generating predictions of ecosystem change based on a mechanistic understanding of all organism-level effects of stressors remains unrealistic [17,103].Therefore, identifying general ecosystem-level changes under stress presents a promising route towards a more efficient predictive ecosystem science, responding to the urgent needs of the biodiversity and climate crisis.

Fig. 1 Fig. 2
Fig. 1 Location of the study localities, Dobu and Upa-Upasina, in Milne Bay Province, Papua New Guinea (A & B)

Fig. 3
Fig.3Heat maps of significant differential abundance of phyla with decreasing pH.Significant change in differential abundance of algal phyla (23S rRNA gene) with decreasing pH are seen in (A); algal phyla are shown on the left, and families are shown on the right.Algal taxonomy was assigned using the microgreen database[58].Significant change in differential abundance of microbial phyla (16S rRNA gene) with decreasing pH are seen in (B); bacterial phyla are shown on the left.Microbial taxonomy was assigned using the GTDB taxonomy[57]

Fig. 5
Fig. 5 Microbiome and metabolome distinctness of benthic holobiont communities as a function of pH.Microbiome (A) and metabolome (B) of the benthic holobiont community versus the sediment microbiome and microbiome (C) and metabolome (D) of the two sponge microbiomes versus the sediment microbiome.Distinctness was calculated as percentage of ASVs/metabolites not shared between microbiomes/ metabolomes.Horizontal dotted line indicates 50% distinct.**p-value < 0.01

Table 1
Summary table of all 16S rRNA gene and metabolomic data generated

Table 2
Mixture Sloan neutral community models.Summary of fitted parameters at each pH regime and model fit statisticsMixture weighted