Sample collection took place during the boreal spring, between May 26th and June 7th 2018, aboard the R/V Alis, as described in Meunier et al. (2021b). At a hydrothermal vent area at Tutum Bay (Ambitle Island, New Ireland Province, Papua New Guinea), samples were collected from one site exposed to high pCO₂ (4°3′44.82′′S 153°34′56.09′′E) and one reference (control) site with ambient pCO₂ (4°4′0.60′′S 153°34′41.43′′E) located ca. 1.8 km south of the main vent site (Biscéré et al., 2019; Pichler et al., 2019). Collections took place at the same depths (between 3–8 m) for both sampling sites. Coral fragments were cut with a bone cutter while SCUBA diving and transported in individual zip-lock bags to the laboratory within 30 min of collection. We collected a total of 24 samples from 4 coral species at the high pCO₂ site: 7 samples of Acropora secale, 5 of Galaxea fascicularis, 7 of Pocillopora damicornis, and 5 of Porites rus. Similarly, we sampled 27 colonies from 4 coral species at the control site: 10 fragments of A. secale, 5 of G. fascicularis, 9 of P. damicornis, and 3 of P. rus. Thus, a total of 51 samples were collected from both sites. Samples were immersed in DNA/RNA Shield (Zymo, Freiburg, Germany) and shipped on dry ice for DNA extraction, nifH amplification, and sequencing. The P. damicornis samples are the same as those analyzed in Meunier et al. (2021b).

DNA was isolated from 51 coral fragments and 1 negative control (null input sample: to assess for DNA extraction kit contamination) using the Qiagen DNeasy 96 Blood & Tissue kit (Qiagen, Hilden, Germany) with slight modifications to the manufacturer’s instructions. Briefly, coral samples were fragmented with a bone cutter and transferred to Eppendorf tubes with 400 μl of ATL buffer and 20 μl of protease. Samples where then incubated in a thermal mixer (Eppendorf, Hamburg, Germany) for 2 h at 56°C and 600 rpm. Afterward, for each sample, 200 μl of the lysate were transferred into microtubes. To aid the disintegration of coral mucus, tubes were incubated at 56°C for 1 h. DNA extraction then proceeded according to the manufacturer’s instructions. DNA concentrations were quantified using Qubit and the dsDNA High Sensitivity Assay Kit (Invitrogen, Carlsbad, United States). To amplify the nifH gene, the IGK3 (5′-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGGCIWTHTAYGGIAARGGIGGIATHGGIAA-3′) and the DVV (5′-GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGATIGCRAAICCICCRCAIACIACRTC-3′) primers were used (Illumina adaptor overhangs underlined) (Ando et al., 2005). Triplicate PCRs were performed according to Lesser et al. (2018) using the Qiagen Multiplex PCR kit, with a final primer concentration of 1.2 μM in a final reaction volume of 10 μl. Thermal cycler conditions for the nifH PCR amplifications were as follows: 95°C for 15 min, 40 cycles of 95°C for 45 s, 57°C for 45 s, and 72°C for 60 s, followed by a final extension step of 72°C at 10 min. A total of 5 μl of the PCR reactions were run on a 1% agarose gel to confirm successful amplification. Triplicates for each sample were pooled, and samples were purified using Agencourt AMPure beads (Agencourt Bioscience Corporation, Beverly, MA, United States). Samples were then indexed and Illumina sequencing adaptors were added using the Nextera XT Index Kit v2. Successful addition of indexes was confirmed by comparing the length of the initial PCR product to the corresponding indexed sample on a 1% agarose gel. Samples were then cleaned and normalized to 1 ng/μl using the SequalPrep Normalization Plate Kit (Invitrogen, Carlsbad, United States). All samples were pooled in a microtube (6 μl per sample) and concentrated using a CentriVap Benchtop Vacuum Concentrator (Labconco, Kansas City, MO, United States). Next, the quality of the library was assessed using the Agilent High Sensitivity DNA Kit on the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, United States). Quantification was done using Qubit (Qubit dsDNA High Sensitivity Assay Kit, Invitrogen, Carlsbad, United States). Sequencing was performed at 6 pM with 20% phiX on the Illumina MiSeq platform at 2 × 301 bp paired-end V3 chemistry according to the manufacturer’s specifications. Sequencing data determined in this study are available at NCBI BioProject PRJNA668126 at https://www.ncbi.nlm.nih.gov/bioproject/PRJNA668126.

Quality control of MiSeq nifH amplicon paired reads was done using the TaxADivA (TAXonomy Assignment and DIVersity Assessment) perl wrapper script (Gaby et al., 2018), as detailed in Lesser et al. (2018) with the following settings: -y -k -r 26 -l 29 –keepc4 -g 500. Subsequent to TaxADivA-based read merging, quality control and filtering, samples with <500 sequences were removed from subsequent analyses (Lesser et al., 2018). The remaining 41 samples (including the negative control) were combined into one fasta file and minimum entropy decomposition (MED) was employed at default settings to cluster reads into oligotypes in an unsupervised fashion (Eren et al., 2015). Representative sequences of the resultant 197 final MED nodes were translated to the protein level using FrameBot (Wang et al., 2013). As a reference dataset, the “nifh_prot_ref.fasta” file included with the FrameBot program was used. To separate nifH-like genes belonging to clusters IV and V, from true functional nifH genes belonging to clusters I-III (Raymond et al., 2004), we employed a three-criteria model sensu (Angel et al., 2018). First, after alignment in MAFFT (Katoh et al., 2002), representative sequences for each final MED node were analyzed using a Classification and Regression Tree (CART) model (Frank et al., 2016). Second, we constructed an Evolutionary Placement Algorithm (EPA) Tree (Berger et al., 2011) using the corresponding RaxML implementation (Stamatakis, 2014). Minimum entropy decomposition node representatives were aligned with a reference dataset from the NifMAP pipeline (Angel et al., 2018) that included representatives of all nifH-like gene clusters (“cluster.rep.nifH-chL-bchX.0.9_noDot.fasta”). Query sequences within this combined reference alignment were then placed on the reference tree (“RAxML_bipartitions.MultipleOriginal.tree” also from the NifMAP pipeline). Third, using node representative sequences as queries, a BLASTP (Camacho et al., 2009) search was conducted against the nifH gene database taken from the Zehr lab website (¹ accessed October 2^nd 2020). In addition, we looked for matches in the NCBI non-redundant (nr) protein database (excluding uncultured/environmental sample sequences and automatically generated model annotations) to further confirm the results of the initial BLASTP search. In case of mismatching results but similar percentages of protein identity between the two queries, the result to the curated Zehr nifH database was preferred. In order to determine the cluster designation for the representative sequences of all final MED nodes, a node was considered to be a true nifH-representative (i.e., belonging to clusters I-III) if at least two of the three criteria produced consistent results.

In order to obtain putative species-level annotations and to collapse the MED count matrix for downstream analyses, an OTU (operational taxonomic unit) clustering at the nucleotide level was conducted via CD-HIT-EST (Fu et al., 2012). The clustering cutoff was set to 91.9%, corresponding to the species-level cutoff previously reported (Gaby et al., 2018). Reference sequences of final MED nodes that only differed by their length, i.e., by means of additional or missing nucleotides at the 5′ or 3′ of the sequence, were combined in this clustering step, which led to the removal of 42 redundant MED nodes from an initial total of 197 MED nodes. Lastly, all samples that contained <100 sequences of cluster I-III nifH representatives were removed from subsequent analyses, retaining 25 of the 40 coral samples and 1 negative control. A maximum likelihood phylogenetic tree of all final OTUs was constructed using RAxML from 20 starting trees using the PROTGAMMAAUTO method to determine an optimal substitution model (LG + G in this case). Additionally, 1,000 bootstrap trees were calculated to support the bipartitions of the best ML tree. Trees were visualized using the R package ggtree (Yu et al., 2017).

Statistical analyses within and between sites and species as well as plotting were done in R using the packages phyloseq, ape, and vegan (R Development Core Team, 2015). Pairwise PERMANOVA tests within and between coral host species for a given sampling site were done using the custom R function “pairwise_adonis.R”² based on a Bray-Curtis dissimilarity matrix. Comparisons of Shannon diversity, inverse Simpson’s Index, and the Chao1 estimator between the two sampling sites over all species were based on Wilcoxon rank sum tests. In addition, we conducted Wilcoxon rank sum tests on those coral host species with samples available from both sites (i.e., P. damicornis and G. fascicularis).

To determine OTUs that were differentially abundant between the control and seep site, linear discriminant analysis (LDA) was conducted in the LEfSe implementation on the Galaxy web application of the Huttenhower lab (³accessed October 2nd 2020). LEfSe statistical tests for differential abundance were carried out with alpha values for factorial Kruskal-Wallis and subsequent pairwise Wilcoxon tests considered significant at ≤0.05. For the LDA analysis, a threshold score of ≥2.0 for OTUs was considered indicative for a significant difference in the relative abundance between sites. This analysis was done across samples from all species, for samples of P. damicornis, and for samples of G. fascicularis.

Patterns of relative abundance based on nifH amplicon sequencing data were confirmed using qPCR (quantitative polymerase chain reaction). For this, OTU_31 was chosen due to its pervasive presence among P. damicornis samples. A specific primer pair for OTU_31 was designed based on its MED representative sequence using Primer-BLAST (Ye et al., 2012). The OTU_31 primer pair (forward: 5′-TTGCGATCCCAAAGCTGACT-3′; reverse: 5′-GCCGGACTCAACACAGAGAA-3′) amplifies a 157 bp long region of the nifH amplicon determined from the MiSeq sequencing. Primer products were checked on a 1.2% agarose gel for single bands after PCR amplification using the following protocol: 90°C for 10 min, followed by 40 cycles of 95°C, 57°C, and 72°C for 30 s each, and a final extension step of 72°C for 10 min. Upon confirmation of a single amplification product using gel electrophoresis, amplicon gel slices from different samples were cut out, and DNA was extracted and purified using the QIAquick Gel Extraction Kit (Qiagen, Germany) following the manufacturer’s protocol. The purified product was Sanger-sequenced to confirm that the amplified sequence was identical to the OTU_31 representative sequence (Supplementary Figure 1). The 16S rRNA gene amplicon primers 784F and 1061R (Andersson et al., 2008; Bayer et al., 2013) were used to produce a standard for relative abundance calculation. To confirm that primer amplification efficiency for the nifH and 16S primer pairs were similar and close to 100%, a 10-fold serial dilution of DNA from a sample with a known presence of OTU_31 was prepared for the determination of qPCR standard curves (Supplementary Figure 2). All qPCR reactions were run on a qTOWER³ 84 using the innuMIX qPCR DSGreen Standard master mix (both Analytik Jena GmbH, Germany) with final primer concentrations of 0.5 μM in a reaction volume of 10 μl with a cycling protocol of 95°C for 10 min, 50 cycles of 95°C for 30 s, 57°C for 30 s, and 72°C for 30 s, followed by 72°C for 10 min, and a subsequent melting curve to assess uniformity of amplification and to confirm the absence of primer dimers. The qPCR samples were run in triplicate in addition to a no-template and an extraction kit negative control for both primer pairs, respectively. Fold change differences were calculated using the ΔΔCT method. For this, the mean ΔCT (CT_nifH – CT_16S) across all samples was subtracted from the ΔCTs of individuals samples and the calculated ΔΔCTs were expressed as relative fold changes of nifH abundance using 2^–ΔΔCT. Pearson’s correlation coefficient was calculated to confirm that qPCR and nifH amplicon sequencing showed similar distribution patterns of OTU_31 in all tested samples.

We collected 51 samples (plus 1 negative control, see Methods) from four different coral species: 7 and 10 samples of Pocillopora damicornis, 5 and 5 samples of Galaxea fascicularis, 7 and 9 samples of Acropora secale, and 5 and 3 samples of Porites rus from the high and control CO₂ site, respectively. For each species, samples were collected from a high carbon dioxide volcanic seep site with a partial pressure pCO₂ = 1179 μatm (pH_T 7.73) and a nearby control site with an ambient pCO₂ = 436 μatm (pH_T 8.12). After amplifying nifH-like genes using the IKG3 and DVV primer pair and MiSeq sequencing, a total of 7,143,933 read pairs (on average 137,383.3 read pairs/sample) were available for analysis. After quality filtering and read pair merging, 10 samples were excluded due to low sequence counts. With the remaining 41 samples, 6,240,993 sequences were assessed using MED. Post-MED, a total of 4,531,377 sequences remained, of which 637,413 (14.1%) were annotated as functional nitrogenase genes of nifH clusters I-III, while 3,893,964 (85.9%) were identified as belonging to cluster IV and V nitrogenase gene homolog representatives. On average, coral samples were represented by 110,521 sequences after MED. Of the 40 coral host samples that were considered during MED, 15 were discarded due to having <100 sequences annotated to nifH clusters I–III. Consequently, 25 samples remained for downstream analysis: 4 A. secale, 11 P. damicornis, 7 G. fascicularis, 3 P. rus (Table 1). These 25 samples comprised 2,662,462 sequences in total. On average 23.5% of sequences from all samples could be assigned to nifH clusters I-III. Although there were no statistically significant differences in the ratio of functional nitrogenase to nifH homolog annotated sequences between species (ANOVA, p = 0.08), P. damicornis samples showed the highest presence of cluster I–III to cluster IV–V sequences (84% of assesses sequenced) (Figure 1). In addition, samples from the control site had higher relative numbers of cluster I–III sequences (33% of assessed sequences) in comparison to the seeps site (15% of assessed sequences, t test, p = 0.1). The negative control (no input DNA extraction) comprised 6,984 sequences after MED and showed a higher relative number of functional nitrogenase sequences, with 92.7% of retained sequences annotated as true nifH (Table 1 and Supplementary Table 1).

Clustering of the 197 determined MED nodes into OTUs using a 91.9% sequence identity threshold in CD-HIT-EST resulted in 88 OTUs, of which 61 were designated as cluster IV nifH homologs or cluster V bacteriochlorophyll synthesis genes and not considered further (Supplementary Data 1). A maximum likelihood phylogenetic analysis confirmed that all cluster IV/V nifH homolog OTUs were monophyletic and highly conserved (Figure 2). The retained 27 OTUs were assigned to nifH cluster I–III. Of these, 25 were cluster I Mo-Fe nitrogenase representatives assigned as follows: Gammaproteobacteria (9 OTUs), Alphaproteobacteria (8 OTUs), Cyanobacteria (6 OTUs), and Epsilonbacteraeota (2 OTUs) (Waite et al., 2017) (Figure 2 and Supplementary Data 1). The majority of Gammaproteobacteria OTUs were identified as Alteromonadales, some closely matched Pseudomonadales and Vibrionales (Figure 2). All Alphaproteobacteria were annotated as Rhodobacterales according to BLAST in discordance with the Zehr database that annotated all OTUs except OTU_46 as Sphingomonadales. Notably, two OTUs belonging to the order of Rhizobiales were detected in the negative extraction (OTU_70: 65.52% of sequences in negative control, OTU_77: 27.19% of sequences in negative control; Supplementary Data 1), thus considered contaminants (red font annotation in Figure 2), and removed, leaving 23 cluster I OTUs for downstream analyses. Only one OTU (OTU_21) was assigned to nifH cluster II (anaerobe Mo-Fe nitrogenases), a bacterium of the class Opitutae. Likewise, only one OTU (OTU_13) was assigned to Cluster III (Mo-independent nitrogenases), a bacterium from the group Bacteroidetes/Chlorobi.

Based on the retained 25 OTUs, the highest number of distinct diazotrophs was found among samples of P. damicornis (14 OTUs), followed by G. fascicularis (12 OTUs), P. rus (8 OTUs), and A. secale (5 OTUs). Generally, there was little overlap in diazotroph community composition between coral species (Figure 3A). For each coral species, at least 25% of the OTUs were not observed in any other species, and none of the 25 OTUs was present in all four host species. Apart from the differences between coral species, diazotroph assemblage between colonies of the same species was also largely variant (Figure 3B). We could not detect any OTU that was consistently present across all samples of a given coral species, site, or host species-site combination with one exception. OTU_05 (Gammaproteobacteria) was found in all but one control site sample of G. fascicularis (Figure 3C). Further, we found three OTUs that were shared between three coral host species: OTU_31 (Cyanobacteria) was present in P. damicornis, G. fascicularis, and A. secale, while OTU_18 and OTU_23 (both Gammaproteobacteria) were present in P. damicornis, G. fascicularis, and P. rus. Thus, we found little similarity in diazotroph abundance between coral species as well as between colonies of a given coral species. A principal coordinate analysis of the Bray-Curtis dissimilarity matrix of relative OTU abundances between samples did not reveal a pattern of distinction between the seep and control site (Figure 3B; PERMANOVA: F = 0.88, p = 0.53). Rather, coral host species explained more of the variation between samples than site (coral species R² = 0.22; site R² = 0.03), with G. fascicularis and P. rus samples clustering closer together. On the other hand, samples of P. damicornis and A. secale appeared more spread out on the ordination. PERMANOVA analysis further supported this observation (F = 1.93, p = 0.0041) and pairwise comparisons indicated that the difference stems from the difference between G. fascicularis and P. damicornis samples (Supplementary Table 2; F = 2.84, FDR-adjusted p = 0.03) as well as between G. fascicularis and A. secale (F = 2.38, FDR-adjusted p = 0.03). There was no difference between sets of samples when site and coral host species were considered as factors (F = 0.75, p = 0.70).

Coral samples were associated with an average of only 3.08 ± 0.26 OTUs. In general, each sample was dominated by one or two highly abundant OTUs (Figure 3C). The OTU with the highest relative abundance in every sample made up on average 82.22 ± 4.02% of nifH sequences. Alpha diversity was generally very low among all samples (Table 1). There were no differences between samples or species from the control and seep site when looking at Shannon diversity (Wilcoxon = 45.5, p = 0.081), Inverse Simpson’s index (Wilcoxon = 44.5, p = 0.073), and Chao1 estimator (Wilcoxon = 54.5, p = 0.195) for all host species (Supplementary Figure 3), or when only samples of P. damicornis (Wilcoxon, all p > 0.05) or G. fascicularis (Wilcoxon, all p > 0.05), the two coral species with samples from both sites, were considered (Supplementary Table 3).

Testing for differential relative abundance of OTUs (LefSe analysis) between the high pCO₂ and ambient site identified no OTUs that were indicative for either the seep or control site across all four coral host species. When only P. damicornis samples were considered, however, there were clear differences in the distribution of diazotroph taxa between samples from both sites. At the class level, Gammaproteobacteria were more prevalent at the ambient site, while Alphaproteobacteria, specifically OTU_46 (closest matching to family Rhodobacteraceae), were more prevalent at the high pCO₂ site (LDA effect size = 5.41, p = 0.013). For control samples, OTU_03 (family Vibrionaceae) was identified as differentially abundant (LDA effect size = 5.13, p = 0.036), with Gammaproteobacteria as a whole being discriminant (LDA effect size = 5.40, p = 0.015). By contrast, in G. fascicularis, none of the OTUs exhibited a significant differential abundance between the high pCO₂ and control site.

We validated the nifH amplicon sequencing data using qPCR to assess and rule out potential biases stemming from differences in sequence amplification (see Figure 1) and specificity of degenerate primers (see Table 1) (Gaby and Buckley, 2012; Angel et al., 2018; Lesser et al., 2018; Tilstra et al., 2019). We chose OTU_31 as a representative because it was broadly and abundantly present among samples of P. damicornis (Figure 3C). Based on a qPCR standard curve, the specific primer pair for OTU_31 showed high efficiency (104.84%), comparable to the primer pair used for 16S rRNA gene amplification that was used as a reference (100.92%). Sanger sequencing of the single PCR product obtained confirmed that the amplicon is identical to the representative sequence of OTU_31 (Supplementary Figure 1). We assayed all samples of P. damicornis using the OTU_31 specific primer pair and the 16S rRNA gene primer pair (see Methods). The comparison of fold change values relative to 16S and the relative abundance in the nifH sequencing dataset showed high correlation between qPCR and amplicon sequences (Figure 4) with a Pearson correlation coefficient of r = 0.87 (p < 0.001).

To analyze nifH amplicon sequencing data, we employed a pipeline based on (Lesser et al., 2018) and modified using methods from Angel et al. (2018) to confidently assign sequences to nifH clusters. Amplification of the nifH gene is commonly achieved through the use of highly degenerate primers at high concentrations using a high number of PCR cycles. The consequential amplification and sequencing of homologous sequences that do not represent a functional nitrogenase subunit is inevitable. Thus, it is important for downstream analyses to identify and remove such false positives. Employment of three different methods to identify such sequences of nifH cluster IV and V (nifH homologs without known function and bacteriochlorophyll synthesis genes) (Raymond et al., 2004), we found that only 14.1% of sequences could be considered a functional nitrogenase gene in the present dataset. This is in line with what (Lesser et al., 2018) and (Angel et al., 2018) observed using the same primers. Notably, we found that the amount of functional nifH sequences was highly variable across individual samples, which varied widely in their overall read depths as well as the read depths of cluster I-III nifH sequences.

After removal of all cluster IV/V homologs of nifH from the initial 88 OTUs, 27 of them were identified as belonging to true, functional nitrogenase genes. Most of these nitrogenase sequences belonged to cluster I Mo-Fe nitrogenases, while cluster II and III were only represented by one OTU each. These cluster I-III nifH OTUs covered a wide range of phylogenetic backgrounds, with the commonly found groups of diazotroph associates of corals being represented by Cyanobacteria as well as Alpha- and Gamma proteobacteria. Notably, samples in this study contained two cluster I nifH OTUs belonging to the Rhizobiales order, closely related to Bradyrhizobium sp. and Rhodopseudomonas oryzae, respectively. Rhizobiales are commonly reported in diazotroph assemblages of corals (Santos et al., 2014; Lesser et al., 2018; Bednarz et al., 2019), sometimes comprising a large relative abundance of up to 70% (Lema et al., 2012). In this study, two Rhizobiales OTUs of the family Bradyrhizobiaceae were identified across samples, but they also amplified in the negative DNA extraction control (no sample input) where they were the most abundant. This finding cautions the high abundance and importance of Rhizobiales in coral diazotroph communities proposed previously and highlights the importance of implementing negative controls, as noted previously (Laurence et al., 2014; Salter et al., 2014; Weyrich et al., 2019). Following the above, we excluded these two OTUs (OTU_70, OTU_77) from consideration in downstream analysis.

Previous studies suggest a mutualistic relationship between corals and their associated nitrogen-fixing bacteria that could provide adaptive potential to the holobiont under a changing environment (Lesser et al., 2007; Cardini et al., 2014; Benavides et al., 2016; Voolstra and Ziegler, 2020). The functional importance of diazotrophs could thus affect the coral holobiont in two opposing ways: (i) diazotroph community assemblages are flexible, thereby allowing to dynamically change in response to changing environmental conditions (Roder et al., 2015; Ziegler et al., 2017; Voolstra and Ziegler, 2020), or (ii) the dependence of the holobiont on fixed nitrogen results in a near obligate association with a low level of flexibility under changing environmental conditions.

We observed that most coral samples harbored only two to five distinct diazotrophs, and each sample was dominated by one or two taxa. The overall validity of the patterns observed in the amplicon sequencing results was confirmed by the high correlation between qPCR-based cycle threshold (C_T) fold-change difference and relative abundance in the amplicon sequencing data of OTU_31. Thus, considering that we identified on average only one or two OTUs of high relative abundance in most samples, coral-associated diazotroph communities may be even less diverse than reported in previous studies (Lema et al., 2012, 2014; Lesser et al., 2018; Bednarz et al., 2019), at least for the coral species in this study.

Considering all samples analyzed, coral host species were predominantly associated with Alpha- and Gammaproteobacteria, as well as Cyanobacteria (Figure 3C). This is largely in line with previous studies of nifH-based analyses in corals that found pre-dominance of Alpha-, Delta-, Gammaproteobacteria (Zhang et al., 2016; Lesser et al., 2018), Alpha-, Delta-, Gammaproteobacteria as well as Cyanobacteria and Bacteroidetes (Bednarz et al., 2019), or Alphaproteobacteria (majoritively Rhizobiales, but see above) (Lema et al., 2014). Notably, Deltaproteobacteria, although found in other studies (Zhang et al., 2016; Bednarz et al., 2019), were not found in samples from this study.

In P. damicornis, diazotroph communities exhibited prevalent association with Gammaproteobacteria at the control site in comparison to samples from the high pCO₂ site that were predominantly associated with an Alphaproteobacterium of the Rhodobacteraceae family (OTU_46, Figure 3C). This resembles findings by Yang et al. (2015) that observed a shifting dominance of Alpha- to Gammaproteobacteria in P. damicornis in a coast to open sea gradient. Bacteria of the family Rhodobacteraceae are commonly associated with stress, injury, or disease in corals (Sekar et al., 2008; Sunagawa et al., 2009; Roder et al., 2014a,b).

The absence of diazotroph diversity or conserved coral species-specific assemblages are striking and warrant further investigation regarding the underlying reasons and the broad validity of these observations. From the current data, it seems that diazotrophs may be opportunistically acquired (i.e., diazotroph assemblage is haphazard) from the surrounding environment, based on their functional contribution. As such, it may be unlikely to find a pattern of phylosymbiosis and cophylogeny as demonstrated previously for coral-associated bacteria (Pollock et al., 2018). This may provide one explanation for the de facto absence of conserved diazotroph community patterns. Considering the notion that corals fix on average less nitrogen than other benthic reef organisms (Cardini et al., 2014), diazotrophy of coral-associated microbes may not be as essential for holobiont functioning as previously thought.

Although volcanic CO₂ seeps offer a unique opportunity to examine the consequences of OA in a natural setting, they offer only a steady-state view of microbiome composition under high pCO₂ conditions and do not allow teasing apart the adaptive potential underlying microbiome assemblage differences. Further studies will be required to complement nifH sequencing results with physiological measurements to disentangle the contribution and implications of diazotroph microbiome composition for holobiont functioning under ambient and high pCO₂ conditions. Ideally, isolating and culturing approaches of differentially abundant OTUs combined with genome sequencing will clarify a potential functional transition from heterotrophic Gammaproteobacteria to autotrophic Alphaproteobacteria under OA, as observed in P. damicornis. In addition, in situ transplantation experiments of coral fragments between high and ambient pCO₂ sites offer a promising approach to study diazotroph assemblage dynamics.