Inter-domain microbial diversity within the coral holobiont Siderastrea siderea from two depth habitats

Sample collection

Siderastrea siderea colonies were sampled from 27 ± 1 m (n = 8) and 17 ± 1 m (n = 6) from Conch reef surrounding the Aquarius Reef Base (Florida International University) in the Florida reef tract, approximately 8 km offshore from Key Largo. These samples were collected by aquanauts performing a saturation mission in the Aquarius Reef Base habitat during the 20th NASA Extreme Environment Mission Operations (NEEMO) mission in July 2015. Preceding the sampling, the dive team was properly trained and familiarized with the method of collection and the operation of the sampling equipment. Diving from the saturation habitat facilitated greatly extended dives to the deeper collection areas in particular. Fragments of approximately 2 cm² surface area were sampled from healthy colonies from 27 and 17 m using chisel and hammer. Samples were brought to the surface in seawater containing 50 mL tubes within 4 h after collection. Seawater around the sampled corals was immediately removed before fixing the samples in a sterile 20% DMSO, 0.5M EDTA NaCl-saturated solution. It should be noted that collected samples were not actively washed with sterile seawater to remove loosely attached microbes and these microbes may be likely present at background levels in our data set and contribute to the detected patterns or absence of them. Coral tissues were homogenized with high-pressure airbrushing (Rodriguez-Lanetty et al., 2013) and the resulting material was preserved in new 20% DMSO solution at −20 °C until DNA extraction.

DNA extraction, amplification, and sequencing

From each of the eight (27 m) and six (17 m) samples, DNA was isolated using the MP FastDNA spin kit for soil (MP Biomedicals, Santa Ana, CA, USA) with an additional ethanol precipitation step. All DNA samples were divided into three aliquots to allow sequencing of the bacterial 16S rRNA gene, Symbiodinium cp23S rRNA gene, and fungal ITS2 rRNA gene from the same sample.

The cp23S fragment was amplified with primers from Santos, Gutierrez-Rodriguez & Coffroth (2003; 5′-TCAGTACAAATAATATGCTG-3′ and 5′-TTATCGCCCCAATTAAACAGT-3′) and the 16S rDNA V123 region with 27F and 519R primers (5′-AGRGTTTGATCMTGGCTCAG-3′, 5′-GWATTACCGCGGCKGCTG-3′; Turner et al., 1999) using the HotStarTaq Plus Master Mix Kit (Qiagen, Hilden, Germany). The PCR programs for 16S and cp23S consisted of 3 min initial denaturation at 94 °C, 28 cycles of 30 s at 94 °C, 40 s at 53 °C and 1 min at 72 °C, and a final elongation step of 5 min at 72 °C.

To circumvent amplification of DNA from the coral host, the ITS2 region was amplified in a nested PCR using phylum-specific primers (Nikolcheva & Bärlocher, 2004). The first PCR program consisted of 2 min initial denaturation at 95 °C, 35 cycles of 30 s at 95 °C, 30 s at 57 °C and 0:45 min at 72 °C, followed by a final elongation step of 2 min at 72 °C. PCR reactions were conducted in volumes of 20 μL, containing 0.5 μM of each primer, 10 μL 2× PCR Master Mix (Promega, Madison, WI, USA) and 10 ng DNA template. In this first PCR, the universal forward primer ITS3 (5′-GCATCGATGAAGAACGCAG-3′; White et al., 1990) was combined with reverse primers specific for Ascomycota, Basidiomycota, Chitridiomycota, Oomycota, and Zygomycota (Nikolcheva & Bärlocher, 2004). The second PCR program comprised 2 min initial denaturation at 95 °C, 10 cycles of 30 s at 95 °C, 30 s at 60 °C and 30 s at 72 °C, and a final elongation step of 2 min at 72 °C. PCR reactions were identical to the first PCR program with 1 μL of 10⁻¹ diluted amplicon as template. The final amplicons were generated by using the universal ITS3 forward primer and the universal ITS4 reverse primer (5′-TCCTCCGCTTATTGATATGC-3′; White et al., 1990). For all sequencing assays except the first round of the nested ITS2 PCR, primers contained Illumina adapters. Barcodes were located between adapter and forward primer. After amplification, samples were pooled per marker in equal proportions, purified using calibrated Ampure XP beads and DNA libraries were prepared following the Illumina TruSeq DNA library preparation protocol for 300 bp paired-end MiSeq sequencing at MR DNA (Shallowater, TX, USA). De-multiplexed reads were deposited in the SRA database under the accession number PRJNA356144.

Raw sequence reads were joined, de-multiplexed, quality filtered, denoised, purged of chimeras, clustered into operational taxonomic units (OTUs) using the nearest neighbor algorithm and classified with Mothur 1.36.1 software according to the MiSeq SOP (Schloss et al., 2009; Kozich et al., 2013; see Supplemental Information for the executed command steps in Mothur).

To classify the Symbiodinium cp23S sequences we used the sequences from Pochon, Putnam & Gates (2014), who conducted a multigene clade level analysis of the Symbiodinium genus and included the cp23S region using 33 unique sequences representing different types from currently all known clades (clades A–I) with the dinoflagellates Gymnodinium simplex and Polarella glacialis as outgroups. We aligned all sequences with MAFFT (Katoh & Standley, 2013) and formatted the dataset to be compatible with Mothur 1.36.1 (see Supplemental Information for the aligned, unaligned, and taxonomy files). All cp23S contigs were aligned against this dataset for denoising and OTU clustering. The same set of sequences was used to classify the contigs and remove sequences not classified as Symbiodinium. For the bacterial sequences, denoising and classification was conducted using the SILVA reference alignment (release 123; Quast et al., 2013). To effectively eliminate chloroplast, mitochondrial 16S rDNA, and unclassified sequences from the dataset, both the SILVA and GreenGenes (DeSantis et al., 2006) alignments were used. We selected the SILVA reference dataset to align and classify the 16S reads since it is both more recently updated (2015 compared to 2013 for the GreenGenes dataset) and of superior alignment quality (Schloss, 2010). Due to the high variance on the ITS2 region, a high quality reference alignment covering the full fungal kingdom does not exist. Instead, similarity distances were calculated from pairwise alignments. Sequences were initially classified with the User-friendly Nordic ITS Ectomycorrhiza (UNITE) database (Abarenkov et al., 2010), however, due to overall low similarity scores all OTU sequences were additionally compared against the full GenBank record with the basic local alignment search tool (BLAST; Altschul et al., 1990).

For each of the datasets, extremely rare OTUs (represented by 10 or less reads in the total dataset) were excluded from further analysis to prevent diversity inflation with false OTUs resulting from erroneous reads (Bokulich et al., 2013). Classification was conducted with Mothur’s implemented Wang-method and 100 bootstraps per sequence and using a criterion of a minimum of 60 bootstraps per sequence (Wang et al., 2007). 16S and cp23S reads were clustered into OTUs based on the traditional 3% difference criterion. Due to the relative high intragenomic variance on the fungal ITS2 region (Schoch et al., 2012), ITS2 reads were clustered into OTUs using a more conservative criterion of 5% difference.

Statistical analyses

Operational taxonomic unit tables were normalized with cumulative sum scaling and log transformed using the R package “metagenomeSeq” (Paulson et al., 2013). Shannon entropy indices were compared with one-way type-III ANOVAs. To test the hypothesis that microbial communities within the same depths are more similar to each other than between depths, non-parametric analysis of similarities (ANOSIMs) and permutational analysis of variances (PERMANOVAs) were calculated, using Bray–Curtis distances and 9,999 permutations. Non-metric multidimensional scaling (nMDS) plots were produced to visualize similarities in community structures.

To identify significant correlations between distributional patterns of microbes we constructed a co-occurrence network (Newman, 2006; Rodriguez-Lanetty et al., 2013). To reduce the complexity of the combined datasets, bacterial OTUs (bOTUs) were merged at the family level and Symbiodinium OTUs (sOTUs) at the type level. OTUs lacking family level classification were grouped by order or higher if all OTUs in that taxonomic level lacked a family classification. From the resulting table, we computed all possible pairwise spearman rank correlation values (ρ) and corresponding p values, corrected for a false discovery rate (Benjamini & Hochberg, 1995). Correlations were considered valid when p < 0.05 and ρ > |0.6|. The network was visualized with the software Gephi (Bastian, Heymann & Jacomy, 2009). The species indicator analysis (Caceres & Legendre, 2009) was computed from the same table as the network (i.e., containing bacteria merged at the family level and Symbiodinium at the type level) with the R package “indicspecies,” using the “multipatt” command, to test for significant associations to one of the depths and assigning bacterial families, Symbiodinium types and fungal OTUs (fOTUs) with p < 0.05 to be depth indicators.

Alignments for additionally investigated OTUs were constructed with MAFFT v. 7. (Katoh & Standley, 2013). Optimal evolutionary models were calculated with MrModelTest v. 2.3 and selected based on the AIC (Nylander, 2004). Maximum-likelihood analyses were conducted on the CIPRES Science Gateway portal (Miller, Pfeiffer & Schwartz, 2012), using RAxML v.8. (Stamatakis, 2014) with a 1,000 bootstrap iterations.

Symbiodinium, bacterial, and fungal community composition

In total, 54 sOTUs were resolved from the 1,455,725 sequences remaining after quality filtering treatment with an average read count of 103,980 ± SD 90,086 per sample. OTUs were classified to types from clades A (7), B (26), C (6), D (1), F (2), G (2), and H (1). The remaining nine sOTUs lacked clade level classification by 60 bootstrap values or more. The overall most abundant, sOTU1 (39.5% average relative abundance) was classified to clade C, sOTU2 (21.7%) to clade A (A3) and sOTU3 (16.4%) and sOTU4 (10.2%) to clade B and type B1, respectively. sOTU5, classified as an H1 type, was the fifth most abundant sOTU (4.8%), and dominated in one of the 14 samples (Fig. 1A). sOTU6 was classified as a G2-1 type and was the sixth most abundant sOTU (3.3%). The one sOTU classified to clade D, type D1 (Symbiodinium trenchii), was rare and occupied the 24th position.

The total amount of bOTUs after quality filtering was 714 from a total of 115,095 sequences with an average of 8,221 ± SD 4,926 reads per sample. The phylum Proteobacteria was the most abundant (49.2%), followed by Chloroflexi (14.7%), Cyanobacteria (7.4%), and Bacteriodetes (6.1%; Fig. 1B). The most abundant bOTU, bOTU1 (10.2%), was classified to the SAR116 clade (Alphaproteobacteria), followed by bOTU2 (5.7%) classified as Synechococcus (Cyanobacteria) and then bOTU3 (4.4%), classified to the clade SAR202 (Chloroflexi). bOTU4 (3.5%) and bOTU5 (3.5%) were classified to Candidatus Branchiomonas (Nitrosomonadaceae) and Candidatus Lariskella (Rickettsiales), respectively.

Despite using fungal-specific primers (Nikolcheva & Bärlocher, 2004), a relative high amount of ITS2 reads from other eukaryotes was retrieved. Only 22.2% of the 790,398 quality filtered ITS2 reads were fungal. Seventeen percent of the reads were from Cnidaria and the remaining reads from other eukaryotes. A group of fungi-like Opisthokonta (Mesomycetozoa) was also represented with 0.15% of the reads. Using the UNITE reference database in the Mothur pipeline, 184 fOTUs were classified to at least the fungal kingdom. However, additional comparison of the complete fOTU sequence set along the GenBank database showed 145 OTUs to yield high similarity scores with non-fungal sequences. These fOTUs were removed from the OTU table, resulting in 39 fOTUs with only similarity hits to fungi or Mesomycetozoa, constituting 141,706 reads in total and by average 10,122 ± SD 9,042 reads per sample. Most fOTUs belonged to Ascomycota (34). Other fOTUs were most similar to Basidiomycota (2), Entomorphthoromycota (1), and Mesomycetozoa (2). ITS2 similarity with sequences from GenBank was notably low (91.0 ± SD 8.2%). The overall most abundant fungus (fOTU1; 48.7% average relative abundance; Fig. 1C) was most similar (76%) to sequences belonging to Lulworthiales (Sordariomycetes). fOTUs 2 (7.2% average relative abundance) and 3 (6.9%) were most similar to sequences from Hyalorhinocladiella (73%; Ophiostomatales) and Physalospora (86%; Xylariales), respectively. A core set of 21 fOTUs, out of 39, was detected in >90% colonies, 15 of which were present in all samples.

Assemblage patterns in Symbiodinium, bacterial, and fungal communities

Differences in Symbiodinium community structure were detected between the two depths as visual in the nMDS plot (Fig. 2A) and statistically supported by ANOSIM and PERMANOVA analyses (p = 0.001 and 0.002, respectively). This was opposite for bacteria (p^ANOSIM = 0.450, p^PERMANOVA = 0.232) and fungi (p^ANOSIM = 0.261, p^PERMANOVA = 0.118) where differences between depths were not significant. While based on Shannon entropy indices no significant differences were detected in α-diversity between depths in any of the three microbial domains, the mean diversity tended to be higher in the deep habitat (27 m) across the communities (Fig. 3).

Co-occurrence patterns and depth indicators

The resolved network, conducted to explore co-occurrence patterns, contained 30 nodes (fOTUs, Symbiodinium types or bacterial families) and 21 edges (significant spearman correlation values above 0.6 between nodes), three of which were negative (Fig. 4). The average degree (the average number of edges per node) was 1.4 and the average weighted degree, weighted based on spearman correlation values, was 0.94. The modularity of the network or degree to which a network tends to unfold into modules or sub-communities (Blondel et al., 2008) was 1.04.

Within the network, one large module of co-occurring nodes was identified, containing six nodes. The module was connected with three negative edges to a node classified as S-70 (Planctomycetes). The node in the module with the most edges was an uncultured bacterial group that is designated in the SILVA database as “PAUC34f” and had three positive correlations and one negative correlation with the Planctomycete group. Other nodes were classified to Proteobacteria (Desulfurellaceae, Nitrospinaceae, and Nitrosococcus), Chloroflexi (SAR202), and Acidobacteria (Subgroup 11). In the overall network, the majority of nodes were bacterial. Symbiodinium type A13 and an unclassified clade A type co-occurred in a pair. Fungal nodes in the network were fOTUs 23 and 33, most similar to Infundibulomyces (Ascomycota; incertae sedis) and Conicomyces (Chaetosphaeriales), respectively. fOTU 23 co-occurred with the Bacteriodetes family Flammeovirgaceae and fOTU33 with Rhodobiaceae (Alphaproteobacteria; Rhizobiales).

The species indicator analysis identified Symbiodinium types C1 and C90 as indicators of the deep habitat (27 m), (Table 1). Among the bacterial groups, RSB-22 (Verrucomicrobia) was indicative for 27 m and an unclassified alphaproteobacterium and Desulvovibrio (Deltaproteobacteria) were for 17 m. fOTU 12 (with 85% ITS2 similarity to a sequence from Basidiobolus, phylum: Entomorphthoromycota) was significant for 17 m. None of the depth indicators appeared in the network.

Siderastrea siderea associated Symbiodinium

Along the depth cline, light intensity decreases exponentially and is a major variable structuring reef diversity composition (Van den Hoek et al., 1978). In addition to overall intensity, longer wavelengths attenuate faster and as a result, the light spectrum narrows with depth. It is therefore not surprising that changes in the dominantly associated Symbiodinium type have been shown to occur over relative short depth ranges (Frade et al., 2008; Bongaerts et al., 2015). Despite the narrow range covered in this study, an exceptional diversity of Symbiodinium was detected in association with S. siderea, though it is not the first work to report high diversity of Symbiodinium associated with S. siderea (Finney et al., 2010; Bongaerts et al., 2015). Contrastingly, other work has indicated limited diversity of Symbiodinium hosted by S. siderea, though the dominant symbiont may differ (Baker, 2001; Warner et al., 2006). Siderastrea congeners have been previously shown to follow a similar reef site driven pattern of low diversity within site, but with differences between sites (Thornhill, Fitt & Schmidt, 2006; Monteiro et al., 2013).

Depth zonation in S. siderea associated Symbiodinium was reported by Bongaerts et al. (2015) where, Symbiodinium D1 was dominant at 2 m depth. For the site sampled in this study, such shallow depths are unattainable and were not included. However, sOTU24, which was classified as D1, was detected at background levels (0.012% average relative abundance) and could potentially have been more abundant at shallower depths. From 2 to 10 m in depth Bongaerts et al. (2015) detected clade B types to be dominant in a minority of sampled colonies and from 5 to 50 m the type C3 was the predominantly detected primary symbiont. Likewise, in this study we identified an unclassified clade B (sOTU3) and clade C (sOTU1) type to be the most abundant Symbiodinium types at 17 and 27 m, respectively. Although the abundance of sOTU1 was by average more than two-fold higher at 27 m compared to 17 m on average, sOTU1 (being the only unclassified clade C type) was not detected as an indicator of 27 m. However, the other clade C OTUs, C90 (sOTUs 34, 45, and 52) and C1 (sOTUs 38 and 53) were indicators of the 27 m depth.

It should be noted that we did not test for a genotype effect and that the detected Symbiodinium zonation depth patterns might reflect the existence of S. siderea genotypes specific associations. Nevertheless, to our knowledge there is no published data indicating that genotypes within coral species living under the same conditions host Symbiodinium symbionts from differing clades while on the other hand depth has been shown to influence coral-associated microbes (hence Bongaerts et al. (2015) but see also Klaus et al. (2007), Frade et al. (2008), Lee et al. (2012), and Glasl et al. (2017)).

In addition to hosting a diverse background community, varying Symbiodinium types from different clades can be the dominant colonizer. Therefore, we posit that Siderastrea have the physiological capacity for broad host flexibility, but suggest that realized flexibility is winnowed by environmental conditions. Stable conditions could select for a community strongly dominated by one type. The observation of high Symbiodinium diversity may in part result from our study being the first to employ high-throughput sequencing of S. siderea associated Symbiodinium, which has the power to detect the presence of rare symbionts. However, many of the S. siderea symbionts are present in sufficient abundance for detection through the previously employed methodologies such as cloning and DGGE. This then suggests that the historically limited sampling effort on this host species and/or a novelty specific to our study site explain the surprising diversity we report.

Symbiodinium type A3, identified as the second most abundant sOTU, is in the Caribbean known for predominance in shallow specialist hosts, including Acropora and branching Porites from 1 to 10 m (LaJeunesse, Loh & Trench, 2009). Clade A types have physiological adaptations to facilitate light tolerance (Banaszak, LaJeunesse & Trench, 2000; Reynolds et al., 2008). Despite this, Finney et al. (2010) also reported an increasing abundance of A3 with increasing depth within Siderastrea radians and Stephanocoenia intersepta hosts. While Siderastrea spp. are recognized to be more efficient at scattering light than Acropora (Marcelino et al., 2013), it seems unlikely that S. siderea from the deeper habitat are able to mimic the light environment of corals more than 10 m shallower. Thus, the A3 identified in deep Siderastrea and Stephanocoenia may represent a cryptic type or phenotypically divergent genotype as elsewhere in the A3 lineage (Parkinson et al., 2015; Pinzon et al., 2011).

To our knowledge, this is only the second study to report high abundance of Symbiodinium clade H from a scleractinian (after Bongaerts et al., 2015) and one of few to report clade G (van Oppen, Mahiny & Done, 2005; Kimes et al., 2013; Thomas et al., 2014). These clades are known to associate with foraminifera and porifera, respectively (Coffroth & Santos, 2005; Stat, Carter & Hoegh-Guldberg, 2006; Pochon, Putnam & Gates, 2014). Both H1 and G2 Symbiodinium types were present in all 14 samples and the H1 OTU was the dominant sOTUs in one of the deep habitat samples. Similarly, Bongaerts et al. (2015), detected H1 as the dominant Symbiodinium type in one Porites astreoides colony. However, further experimental study is required to characterize the nature of these associations.

Siderastrea siderea associated bacteria

The composition of the bacterial community associated to S. siderea was not structured differently between the depths covered in this study, nor did we detect differences in total Shannon diversity. The stability in alpha diversity is in line with previous environmental sequencing studies that sampled coral microbiomes from different depths (Klaus et al., 2007; Lee et al., 2012; Glasl et al., 2017).

The 16S sequence of the most abundant bacterial phylotype (bOTU1) was classified to the SAR116 clade which despite its classification under the Rickettsiales in the SILVA record, is part of the Rhodospirillaceae (Alphaproteobacteria; Rhodospirillales; Grote et al., 2011). Further investigation of this bOTU showed it was 99% similar to 16S sequences from coral holobionts (Lopez et al., 2008; Sunagawa, Woodley & Medina, 2010; Yang et al., 2013), and identical to that of a bacterium which was recently shown to be a calcifying endosymbiont of a sponge, where it lacked a cell wall and was vertically transmitted (Uriz et al., 2012; Garate et al., 2017). Garate et al. (2017) opted the presence of a “calcibacterium” clade within SAR116 mainly including sequences from sponge and coral holobionts. Based on phylogenetic analysis including sequences from Oh et al. (2010) and Garate et al. (2017) as well as all SAR116 sequences in the SILVA record (Fig. S1), bOTU1 was resolved within a highly supported clade within SAR116 that was indeed dominated by sequences from sponges and corals. As discussed in Garate et al. (2017), the presence of bacterial calcite spherules has not yet been investigated in cnidarians but the high abundance of a calcibacterium related bOTU within the S. siderea microbiome and repeated recovery of related sequences from other corals (Lopez et al., 2008; Sunagawa, Woodley & Medina, 2010; Yang et al., 2013), set the stage for further research to explore whether bacterial calcification might occur in corals as well.

The second most abundant bOTU was classified to Synechococcus, a group of cyanobacteria found in other studies as an abundant coral symbiont that fixes nitrogen utilized by Symbiodinium (Lesser et al., 2004, 2007; Carlos, Torres & Ottoboni, 2013; Morrow et al., 2015). Due to the use of different 16S rDNA regions in this study and Lesser et al. (2004), sequences cannot be directly compared. Nevertheless, the high abundance of Synechococcus in S. siderea might suggest a potential contribution of nitrogen fixation within this holobiontic association, similar to that of Montastraea cavernosa and Synechococcus (Lesser et al., 2004). Another high abundant bacterium (bOTU3) was classified to Chloroflexi SAR202. Chloroflexi are not commonly reported in coral microbiomes, however have been detected on a few occasions (Kimes et al., 2013; Lee et al., 2012), and in particular from S. siderea (Cardénas et al., 2012; Kellogg et al., 2014). Finally, the abundant bOTU5 was assigned to the Rickettsiales, and further phylogenetic investigation showed this sequence to cluster with sequences acquired from other Cnidaria (Fraune & Bosch, 2007; Sunagawa, Woodley & Medina, 2010; Kimes et al., 2013) under a fully supported node (see Fig. S2), sister to the vertically transmitted intracellular Candidatus Lariskella (Matsuura et al., 2012) and Candidatus Midichloria, a candidate genus of ecologically wide spread intracellular symbionts originally isolated from Acanthamoeba (Montagna et al., 2013).

Siderastrea siderea associated fungi

Although fungi have been studied primarily in terrestrial ecosystems, it is known that their presence in the coral holobiont is ubiquitous (Kendrick et al., 1982; Le Campion-Alsumard, Golubic & Priess, 1995; reviewed in Yarden (2014) and Raghukumar (2012)). The ecology of the coral-associated fungi, however, is poorly understood (Golubic, Radtke & Le Campion-Alsumard, 2005). While fungi within the skeleton (or endolithic fungi) may invade coral tissue during times of stress (Yarden, 2014), fungi have also been found to parasitize on endolithic algae (Le Campion-Alsumard, Golubic & Priess, 1995).

The overall low ITS2 similarity of fOTU sequences to GenBank (91.0 ± SD 8.2% average similarity) corroborates the presence of a high novel diversity associated with S. siderea’s microbiome. This agrees with Amend, Barshis & Oliver (2012), where Acropora hyacinthus associated fungi had an average similarity of 97% to sequences from GenBank, using the more conserved large ribosomal subunit (LSU) rDNA marker. A notable pattern within the fungal community was the dominance of the Lulworthiaceae related fOTU1 as the most abundant in nine of the 14 samples. In each of the remaining samples, different fungi (fOTUs 2–6) occupied this position.

Amend, Barshis & Oliver (2012) did not detect differences between fungal communities of A. hyacinthus from sites on the same reef that differed by 1.5 °C. Over this range the dominant Symbiodinium type shifted. Similarly, in this study where samples were collected from the same reef but with difference in depth instead of temperature, the dominant Symbiodinium type shifted (from the unclassified C type to a combination of B types) while such a pattern was not detected for fungal communities.

In Amend, Barshis & Oliver (2012) a core set of 11 OTUs was identified, present in >90% of the samples. Notably, 21 fOTUs in our study were detected above this threshold. Due to the use of different markers, sequences of these core-members cannot be compared between studies. However, Hortaea werneckii and Lindra sp., highlighted in their study, are related to core fOTUs from S. siderea. Among the well-classified fOTUs, fOTU35 is highly similar to H. werneckii (99%), present in 13 of the 14 samples. Lindra belongs to the Lulworthiales, of which 7 fOTUs (4 in more than 90% of the samples) yielded closest similarity hits including the most abundant fOTU1 and fOTU2. The Lulworthiales is an order of marine fungi, predominantly isolated from algae or detritus (Kohlmeyer, Spatafora & Volkmann-Kohlmeyer, 2000). Given the low similarity, the sequences acquired from S. siderea likely reflect one or more novel groups within or related to this order.

The Internal transcribed spacer region has become the standard DNA marker in fungal environmental sequencing assays (Abarenkov et al., 2010). Despite this, our data emphasize that coral-associated fungal ITS sequences are underrepresented in sequence databases, resulting in poor classification of detected OTUs. This problem was also encountered by Wainwright et al. (2017) and corroborates that marine fungi are poorly studied with respect to their terrestrial relatives (Richards et al., 2012; Peay, Kennedy & Talbot, 2016). For future attempts to characterize fungal diversity on coral reefs more conserved markers such as the SSU or LSU rRNA might be recommendable. However, the overall high degree of detected fungal novelty warrants a need for fungal isolation efforts from corals to fill the gaps in the sequence databases. Finally, the fungal-specific primers from Nikolcheva & Bärlocher (2004) used in this study failed to prevent amplification of non-fungal sequences, a problem also encountered by Amend, Barshis & Oliver (2012) on the LSU marker.

Inferences from network analysis

The network of the combined Symbiodinium, bacterial, and fOTU tables resolved predominantly co-occurrence patterns involving bacteria. However, two co-occurrences were observed between bacteria and fungi. Co-occurrence of Symbiodinium A3 and an unclassified type A likely indicates either intragenomic variability within the A3 type or niche-overlap between closely related types. The network contained one sub-community, or module, consisting of six bacterial groups and connected by three negative correlations with a group of Planctomycetes of the class Phycisphaerae. The central node, having three positive correlations and one negative included bOTUs classified in the SILVA database to the unrecognized phylum PAUC34f. This enigmatic group of bacteria has not yet been cultured and sequences in the SILVA database are solely from metagenomic studies (Quast et al., 2013).

A noteworthy member of this node is the Chloroflexi SAR202 clade, as this includes the third most abundant bOTU. The node classified as Subgroup 11 of the Acidobacteria was connected to the PAUC34f node. Acidobacteria are specialized aerobic oligotrophs (Kielak et al., 2016) obtaining energy from denitrification and known to engage in syntrophic associations with different proteobacteria (Spring et al., 2000; Meisinger et al., 2007). Proteobacteria in the module were most similar to the genus Nitrosococcus and the family Nitrospinaceae, both of which are groups of nitrifying bacteria (Kowalchuk & Stephen, 2001; Lücker & Daims, 2014). Although little is known about the majority of organisms from this module, its members have distinctive metabolisms and it seems unlikely that the co-occurrence patterns are solely a result of niche-overlap.