Three-Dimensional Molecular Cartography of the Caribbean Reef-Building Coral Orbicella faveolata

All organisms host a diversity of associated viruses, bacteria, and protists, collectively defined as the holobiont. While scientific advancements have enhanced the understanding of the functional roles played by various components of the holobiont, there is a growing need to integrate multiple types of molecular data into spatially and temporally resolved frameworks. To that end, we mapped 16S and 18S rDNA metabarcoding, metatranscriptomics, and metabolomic data onto three-dimensional reconstructions of coral colonies to examine microbial diversity, microbial gene expression, and biochemistry on two colonies of the ecologically important, reef-building coral, Orbicella faveolata and their competitors (i.e., adjacent organisms interacting with the corals: fleshy algae, turf algae, hydrozoans, and other corals). Overall, no statistically significant spatial patterns were observed among the samples for any of the data types; instead, strong signatures of the macroorganismal hosts (e.g., coral, algae, hydrozoa) were detected, in the microbiome, the transcriptome, and the metabolome. The 16S rDNA analysis demonstrated higher abundance of Firmicutes in the coral microbiome than in its competitors. A single bacterial amplicon sequence variant from the genus Clostridium was found exclusively in all O. faveolata samples. In contrast to microbial taxa, a portion of the functionally annotated bacterial RNA transcripts (6.86%) and metabolites (1.95%) were ubiquitous in all coral and competitor samples. Machine learning analysis of microbial transcripts revealed elevated T7-like cyanophage-encoded photosystem II transcripts in O. faveolata samples, while sequences involved in bacterial cell division were elevated in turf algal and interface samples. Similar analysis of metabolites revealed that bacterial-produced antimicrobial and antifungal compounds were highly enriched in coral samples. This study provides insight into the spatial and biological patterning of the coral microbiome, transcriptome, and metabolome.


INTRODUCTION
The holobiont is the sum of an organism and all of its associated microbes and viruses, which interact through a complex suit of biochemicals. Technologies recently adopted for marine molecular ecology have widened the lens with which ecologists can investigate a holobiont, including standard marker amplicon sequencing, metabarcoding, metagenomics, metatranscriptomics, metaviromics, metaproteomics, metabolomics, and others. Each of these molecular approaches on their own offer a unique view into facets of macroorganismal, microbial, viral, or biochemical dynamics (Hasin et al., 2017;Pinu et al., 2019;Roach et al., 2020Roach et al., , 2021a. For example, 16S and 18S rDNA amplicon sequencing have allowed for targeted taxonomic profiling of prokaryotic and eukaryotic microbiomes in many systems from corals to humans (e.g., Rohwer et al., 2002;Turnbaugh et al., 2007). Metagenomics, or the shotgun sequencing of DNA, has allowed for the identification and quantification of microbial and viral taxa, and functional genes (Handelsman et al., 1998;Breitbart et al., 2002;Edwards and Rohwer, 2005;Franzosa et al., 2018). Metagenomics also provides the ability to assemble genomes de novo from community-level sequencing, and identify and describe new gene families (Parks et al., 2017;Ovchinnikov et al., 2017). Furthermore, purification-based metagenomic methods, such as metaviromics, utilize physical and chemical properties to isolate and sequence specific gene pools (e.g., viral pool, archeal pool) that would otherwise comprise small percentages of total DNA present in the holobiont (Thurber et al., 2009;Poulos et al., 2018). Together, these emerging approaches have improved the capacity for understanding the complex nature of holobiont biology and ecology, by providing the means to investigate all of the holobiont constituents in situ.
Metatranscriptomics, or the shotgun style sequencing of RNA, allows for the quantification of total gene expression in a community (Wang et al., 2009). In addition, amplicon sequencing of conserved coding regions of the genome (of both macroorganisms and microbes) has provided species-level and higher taxonomic identifications of operational taxonomic units (OTUs) or amplicon sequence variants (ASVs) to single nucleotide differences (Callahan et al., 2016;Thompson et al., 2017). Metabolomics provides insight into the physiological and metabolic activity of a holobiont by identifying and quantifying small molecules present in the holobiont (Quinn et al., 2016;Sogin et al., 2016;Hartmann et al., 2017;Matthews et al., 2020). Many of these methods complement each other, for example, open-reading frame protein predictions from metagenomic data can provide information for new proteins or molecules in metabolomic datasets (Wang et al., 2014). High throughput sequencing of amplified DNA or environmental DNA (eDNA) using standard gene markers, known as metabarcoding, has allowed for a massive number of censuses of single and multicellular marine eukaryotes such as algae, invertebrates, and vertebrates (Blaxter et al., 2005;Taberlet et al., 2012;Leray and Knowlton, 2015). Integrating these tools is essential for understanding the roles of microbes, viruses, and biochemicals in large-scale ecological processes. Combining these approaches can improve the capacity for scientists to observe and link states and processes from all aspects of the central dogma of biology in an ecosystem.
In addition to molecular approaches for ecology, developments in the literal 'lens'-i.e., photographic imaging technology-have also increased the amounts and types of information that can be obtained from natural environments (Burns et al., 2015;Roach et al., 2021b). Photogrammetric threedimensional modeling provides a more realistic visualization and representation of environments than traditional twodimensional image collections and has allowed for new types of spatial analysis (Burns et al., 2016;George et al., 2021). For example, analyses of 3D reconstructions of corals and their competitors have demonstrated that the geometric properties of corals 'in part' determine the outcome of competitive coralalgal interactions in several coral genera . Similar techniques for 3D reconstruction have also been used in terrestrial landscapes (Schneider et al., 2014) and freshwater ecosystems (Marazuela et al., 2018).
The combination of molecular data with spatially explicit 3D cartography tools are now allowing for unprecedented insight into pattern and process in many systems. For example, metabolite mapping on the human body demonstrated that the composition of molecules on the human skin varied depending on people's daily routine (Bouslimani et al., 2015). These technologies have yet to be fully utilized in marine environments due to the technological difficulties of in situ sampling and data collection. Here, we combine 16S rDNA amplicon sequencing, 18S rDNA amplicon sequencing, metatranscriptomics, and metabolomics with three-dimensional molecular cartography tools to provide 3D molecular maps of several components of coral and algae holobionts, allowing us to investigate the spatial distribution and ecological drivers of microbial taxa, transcripts, and small molecules.
The goal of this study was to use cutting edge -omics and imaging techniques to explore the processes occuring where different benthic holobionts meet (i.e. interfaces), providing insight into the benthic interaction ecology of the ecologically important, reef-building coral Orbicella faveolata. O. faveolata (previously Montastraea faveolata), or the mountainous star coral, is native to the Caribbean Sea and Gulf of Mexico. O. faveolata is a massive reef-building coral currently listed as endangered by the International Union for Conservation of Nature (IUCN) (Aronson et al., 2008), with an estimated 50% abundance loss in the last 30 years. On healthy reefs, O. faveolata often dominate the reefs at depths of 10-20 m, and these corals face several threats including coral bleaching, algal overgrowth, and diseases such as yellow band disease, black band disease, and white plague (Kimes et al., 2013).
Here, several different molecular techniques were used to investigate specific components of O. faveolata holobionts. 16S and 18S rDNA sequencing data served as a metric to assess total bacterial and protistan taxonomic diversity, and to provide a general overview of the three-dimensional composition of O. faveolata's microbiome. A total RNA sequencing approach was used to establish holobiont gene expression profiles, and metabolomic profiling allowed for the detection of small molecules, yielding a snapshot of the holobiont's biochemistry. We combined these approaches to test the hypothesis that holobiont taxa and function may be driven by spatial patterns (e.g., physical distance), but found no evidence supporting this. Rather, biological variables (e.g., the macroorganismal hosts) proved to be more important for taxonomic, functional, and molecular composition. Integration of multiple data types in a spatial framework provides a better understanding of the distinct components of the holobiont and how they mediate ecological interactions between macroorganisms. Future work should focus on scaling these methods both in terms of higher sampling resolution and sampling across larger spatial scales to further explore unique molecular profiles for specific organisms or physiological states (e.g., stress, disease, competition) and to examine if there are breakpoints where spatial patterns do exist.

Sample Collection
Two individual colonies of Orbicella faveolata were selected as "corals of interest" at two seperate dive sites (Water Factory, 12 • 06 31.7 N 68 • 57 13.4 W and East Point, 12 • 02 35.1 N 68 • 45 43.5 W) on the island of Curaçao in the southern Caribbean (Figure 1). Samples from the O. faveolata colonies and their interacting organisms (turf algae, Millepora complanata, Orbicella annularis, and one conspecific O. faveolata) were collected along transects using the cardinal directions (Figure 1). For each North, East, South and West transect, two coral samples (I and II), one interaction/interface (III) and two competing organism (IV and V) samples were collected. The interaction/interface sample was taken from the border where the two organisms touch, and includes tissue from both of the organisms. Three tissue biopsies (1 cm diameter × 1 cm depth) from each location (I, II, III, IV, and V) were sampled for metabolomics, metatranscriptomics and microbial diversity (i.e., 16S and 18S rDNA sequencing). The tissue biopsies were collected by divers on SCUBA using a diamond encrusted core drill bit (Lasco Diamond Products). Samples were stored in ziploc bags until their return to a moored vessel with indoor laboratory space where samples for metatranscriptomics (RNA) were fixed with RNAlater (cat AM7020) in cryovials and flash frozen in liquid nitrogen. Samples for mass spectrometry (metabolomics) were stored in 70% LC/MS grade methanol/water in glass vials and stored at −20 • C. Samples for metabarcoding (DNA) were also flash frozen in liquid nitrogen. All molecular samples were later moved to −80 • C prior to processing and analysis.

Photogrammetry, 3D Model Construction and Molecular Cartography
Over 200 photographs were taken of each O. faveolata colony and their interacting organisms at the two separate sites (Water Factory and East Point) using a Canon Rebel T4i with a 35-mm lens and two Keldan 800 lumen video lights to illuminate the corals uniformly . A ruler, placed along the interface of the coral colonies, was photographed to set the scale for the resulting 3D digital models. Pre-sampling images were used for the 3D coral reconstruction, and post-sampling images were taken to map the sample location on the 3D models. Photos from each coral colony were color corrected and subsequently imported into Autodesk R ReCapTM Photo 21.0 to create spatially accurate 3D coral reconstructions (Burns et al., 2015;Leon et al., 2015). The models were scaled using the in-reef ruler incorporated into the 3D-model and the ReCap R Photo 'set scale and units' tool. The models were exported as a STL file as well as an OBJ file that included corresponding MTL and JPG files.
The STL models were then imported into MeshLab 2020.03, and the spatial coordinates for each sampling point in the 3D models were determined using methods from Protsyuk et al. (2018). The location of each sample was selected on the 3D model surface using the 'reference scene' tool in MeshLab, and the x, y, z coordinates for each sample location were recorded. The corresponding sample data was then added to the spatial coordinates in a CSV file. Three-dimensional molecular heatmaps were generated using the online software platform ili (Protsyuk et al., 2018). These maps were constructed with the OBJ files with corresponding MTL and JPG files as the input for the 3D surfaces and CSV files for the -omics data and spatial coordinates. Colored mesh was overlaid using the JPG file, and a linear 'hot' color scale was used for all 3D molecular heatmaps. The data was displayed using a spot border capacity and radius of 1.

DNA/RNA Extraction
DNA was extracted from tissue biopsies using DNeasy PowerBiofilm kit (Qiagen Cat No./ID: 2400-50) and quantified using Qubit DNA HS assay kit. Total RNA was extracted from the samples using the AllPrep DNA/RNA kit (Qiagen Cat No./ID: 80204). DNA was used for amplicon sequencing and the RNA portion was dedicated for metatranscriptomes. RNA was quantified using Qubit RNA HS assay kit and checked on an FIGURE 1 | Images (A,C) and 3D models (B,D) of the Orbicella faveolata corals with sampling points overlaid for visualization. Inset in panels (A,C) show the geographic location of the sampling efforts. Compass insets indicate the orientation of the coral colonies for cardinal direction. Three samples were taken from each sampling point (I-V), one for 16S and 18S rDNA amplicon sequencing, one for metatranscriptomic sequencing, and one for LC MS/MS metabolomic analysis for a total of 60 samples (24 coral, 12 interface, and 24 competitor) taken at each site (Water factory and East Point). Models were created using Autodesk ReCap Photo and 3D molecular mapping was performed using ili (scale bar = 20 cm). White arrows indicate the interface between O. faveolata and the interacting competitors (turf algae, Orbicella annularis, Millepora complanata). White arrows indicate the interface between O. faveolata and the interacting competitors (turf algae, Orbicella annularis, Millepora complanata).

Amplicon Sequencing
To determine the bacterial diversity of O. faveolata and its competitors, extracted DNA was sent to the Integrated Microbiome Resource (IMR) at Dalhousie University where the V6-V8 region of the 16S rDNA was amplified, and the amplicon samples were sequenced with Illumina MiSeq using 300 + 300 bp paired end V3 chemistry, resulting in 656,698 reads. To study the protist diversity, PCR amplification of 18S rDNA was performed using non-metazoan 18S rDNA primers (UNonMet). The cleaned and purified PCR products were also sent to the IMR at Dalhousie University where the V4 region of the 18S rDNA was amplified, and the amplicon samples were sequenced with Illumina MiSeq using 300 + 300 bp paired end V3 chemistry, resulting in 2,972,176 reads.
The paired-end reads were demultiplexed and denoised using DADA2 as part of the Qiime 2 TM pipeline (Bolyen et al., 2019), and sequences were trimmed using the parameters "-p-trimleft-f 19 -p-trim-left-r 20 -p-trunc-len-f 290 -p-trunc-len-r 290" for 16S amplicon sequences and "-p-trim-left-f 15 -p-trim-leftr 19 -p-trunc-len-f 290 -p-trunc-len-r 290" for 18S amplicon sequences. Taxonomic classification was performed using the q2-feature-classifier (Bokulich et al., 2018) trained against the SILVA database version 132, and phylogenetic trees were built with FastTree as part of the Qiime 2 TM pipeline. Long branching amplicon sequence variants (ASVs) were removed from 16S and 18S rDNA datasets, along with ASVs classified as plastids and mitochondria in the 16S rDNA dataset. For 16S rDNA, ASVs were collapsed into a Qiime 2 taxonomy table at the phylum level, and analyses were conducted on both bacterial ASVs and phyla.
Bacterial transcripts were identified in our metatranscriptomic sequence dataset using the mapping algorithm SUPER-FOCUS against the RAST-SEED database for microbial functions with default parameters (Overbeek et al., 2014;Silva et al., 2016). Annotations for specific eukaryote coding regions were conducted using FRAP, which utilizes SMALT mapping, with a 96% identity over 100% of the length of the qualityfiltered and trimmed sequence read (Ponstingl and Ning, 2010). Sequences mapping to coral and Symbiodinium reference genomes were removed prior to the microbial subsystem annotation (Supplementary Table 2). A compilation of publicly available coral genomes and Symbiodinium genomes from reefgenomics.org were used to identify coral and Symbiodiniumassociated reads (Liew et al., 2016).

Ultra-Performance Liquid Chromatography -Tandem Mass Spectrometry
Methanolic extracts were analyzed as described in Roach et al. (2020). Analysis was done using Ultra Performance Liquid Chromatography (UPLC) with a Kinetex C18 reverse phase column (Phenomenex Inc.) connected to a Maxis Q-TOF Mass Spectrometer (Bruker Daltonics). The raw datafiles from the MS machine were converted into .mzXML files with the Bruker Data Analysis software version 4.1. The .mzXML files are available on the MassIVE database (massive.ucsd.edu) under numbers MSV000080597 and MSV000080632 (same dataset). The .mzXML files were imported into MZmine2 (Pluskal et al., 2010) of which a beta version was used (2.37.1.corr17.7, see Supplementary Methods). Thresholds used were the same as in Roach et al. (2020) only the ADAP chromatogram builder function was used with the following settings: a minimum height intensity was 3.00E + 03, the group intensity threshold 3.00E + 03, minimum group size in number of scans was set to 3, together with a mass tolerance of 25 ppm or 0.05 m/z.
In addition to the MZmine workflow, the MZmine metacorrelate function was used which performed a Pearson correlation analysis with the following settings: Retention time tolerance: 0.1 min, the minimum height: 4.00E + 03 and noise level: 3.00E + 03. The minimum number of samples in all was set to an absolute value of 2, and there was no minimal number of samples per group. 60% of the intensity overlap was set as a minimum and gap filled features were excluded. Minimum number of data points was set to 5, and two points on the edge. The shape correlation was set to 85% minimum. Ion Identity networking function searched for all modifications and all adducts from the positive ion mode within a 0.05 -25 ppm, with a minimum height of 4.00E + 03. The maximum charge was set to 4, and there was a maximum for 3 molecules per cluster. For more details on the MZmine workflow and specific settings see Supplementary Materials.

Metabolomics: Molecular Networking and Spectral Library Searching
The Feature-Based Molecular Networking (FBMN) workflow (Nothias et al., 2020) in the GNPS environment (Wang et al., 2016) was used to create the Ion Identity Network (IIN) (job ID: 415d51c0fd25492187b2d822eaa6f217). All fragment ions within ± 17 Da of the precursor m/z were removed. MS/MS spectra were filtered by selecting the top 6 fragment ions within a ± 50 Da window throughout the spectrum. Mass tolerances were set for both the precursor ion and the MS/MS fragment ion to 0.05 Da. A molecular network was created with edges having a cosine score above 0.7 and more than 4 matched peaks. Nodes had to appear in each other's top 20 of most similar nodes, and the maximum size of a molecular family was set to 500. If a molecular family was larger than 500, the lowest scoring edges were removed.
Library spectra were filtered in the same way as the input data. For a positive library match the spectra needed a minimum of 4 matching peaks and a cosine score above 0.7. Analog search mode was used by searching against MS/MS spectra within a difference of 100.0 Da of the precursor ion. The edges from the IIN were added separately by the output file of MZmine.
GNPS annotated 53 features based on MS2 spectra as known compounds by matching MS2 spectra to library spectra. Another 238 features were highly similar to MS2 spectra and marked as analog hits. Information on feature structures from the different workflows were integrated with MolNetEnhancer and resulted in 564 Classyfire annotations on the level of Kingdom and Superclass, 562 on class level, and 420 on subclass level. A total of 265 runs were used as input for MZmine, which detected 3355 features. Analysis on MS1 level was done on 31 samples and 24 blanks. Contaminant features were flagged and removed if the maximum peak area in one of the blanks was larger than half of the mean peak area of all samples. Peak area background noise level was set on 4929 as was the smallest peak area before the gap filling step in MZmine. Features needed to pass the peak area background noise level threshold in at least 2 samples. In addition, feature peak areas smaller than ± two-fold background noise (peak area of 10000) were set to 0, but features were kept. This resulted in 2315 features.

Statistical Analyses
For the 16S and 18S rDNA amplicon sequencing, the abundance tables were rarefied by subsampling the reads to account for library size and imported using QIIME 2R 1 . The transcriptome abundance data tables were taken using read hits counts normalized by the total number of hits per library, the metabolomic data table consists of areas under the curve. Maximum likelihood trees of 18S rDNA were built in IQ-TREE v1.5.4. (Nguyen et al., 2015).
The spatial analysis was conducted using physical distances measured over the 3D coral's surface. The shortest 3D topological distance between two sample points was calculated using the ReCap R Photo 'measure distance' tool, and a pairwise physical distance matrix was generated for all samples within the coral colonies. The spatial patterns of the different data types were then analyzed using linear regression analyses in R version 3.4.4 (R Core Team, 2017) with ggplotRegression (Wickham and Chang, 2016), and power analyses were conducted using the R pwr package 1.3-0 (Champely et al., 2017).
All statistical analyses were conducted in R version 3.4.4 (R Core Team, 2017). Abundance data tables of each data type were manipulated with dplyr and tidyverse and used to create Bray-Curtis dissimilarity matrices with the Vegan package 2.5-2 (Oksanen et al., 2016). Additionally, the Vegan and pairwiseAdonis packages were used to generate PERMANOVA test statistics on the Bray-Curtis dissimilarity matrices and provided pairwise comparisons between variables 2 . The Principle Coordinates and eigenvectors were generated by the Ape package 5.4-1 with a Cailliez correction for negative eigenvalues (Paradis and Schliep, 2019) and the Non-metric Multi-dimensional Scaling (NMDS) analysis was 1 https://rdrr.io/github/jbisanz/qiime2R/ 2 https://github.com/pmartinezarbizu/pairwiseAdonis conducted using the metaMDS function of the Vegan package 2.5-2 (Oksanen et al., 2016). Principal Coordinate Analysis coordinates with eigenvectors and NMDS were plotted with ggplot2 (Wickham and Chang, 2016), and heatmaps of relative abundances were plotted using the heatmap function in R. Supervised and unsupervised random forests analyses were run using rfPermute 2.1.81 (Archer, 2016), and 3,000 trees were built using a set seed of 25. Shannon entropy and rarefaction species richness was calculated using the Vegan package 2.5-2. F-Tests (two-sample for variances) and t-Tests (two-sample assuming equal variance or assuming unequal variance) were used for diversity metrics.

Co-occurrence Analysis
To integrate the different datasets, the mmvec tool 3 was used to calculate conditional probabilities of the occurrence of transcriptomes or metabolites with the occurrence of specific bacteria on a phylum level (Morton et al., 2019). The standalone version of mmvec was used and co-occurrence probabilities were extracted by applying a softmax transformation. The conditional probabilities of the top predictors were visualized in a two way heatmap using JMP (JMP, Version 14. SAS Institute Inc., Cary, NC, United States, 1989-2021).

Ubiquity Calculations
The methods from Hester et al., 2016 were used to calculate the ubiquity and relative abundance of the 16S and 18S rDNA ASVs, bacterial transcripts, and metabolites. Briefly, the ubiquity of an individual ASV/transcript/metabolite was determined by the number of samples the ASV/transcript/metabolite was present in divided by the total number of samples. For example, a ubiquity of 1 resulted from the ASV/transcript/metabolite appearing in 100% of the samples. The relative abundance was calculated as a proportion of the entire community (e.g., proportion a specific ASV/transcript/metabolite to all ASVs/transcripts/metabolites from all samples). Ubiquity-abundance was plotted using ggplot2 (Wickham and Chang, 2016).

Three-Dimensional Photomosaics and Spatial Analysis
Analysis of the 3D photomosaic models allowed for the calculation of geometric properties of the corals and revealed that the East Point O. faveolata colony had a maximum height of 1.80 m and a maximum diameter of 1.70 m with a total surface area of 196,388 cm 2 , whereas the Water Factory colony was smaller with a maximum height of 1.06 m and a maximum diameter of 0.70 m with a total surface area of 8,973 cm 2 . The East Point O. faveolata colony had an interface length of 969 cm (colony perimeter length), and the main competing organism interacting with the O. faveolata colony included turf algae (46.81% of the interface), along with a closely related species of coral, Orbicella annularis (17.80% of the interface). Coral overhangs where competitors were absent made up 24.12% of the interface, and interactions with crustose coralline algae (CCA) were limited to 11.27%. The Water Factory O. faveolata colony had an interface length of 545 cm, and the main interactions involved turf algae (27.36%), along with a direct interaction between the coral colony and Millepora complanata (fire coral/hydrozoan, 32.24%). The remainder of the interface consisted of coral overhangs (32.36%) and interactions with CCA (8.04%). The coral colonies and interacting organisms were within a 2.5 m × 2.5 m (6.25 m 2 ) reef plot at East Point and a 1 m × 1 m (1 m 2 ) reef plot at Water Factory.
Molecular cartography and subsequent linear regression analysis of spatial patterning demonstrated no evidence for spatial autocorrelation or isolation by distance, as there was no significant correlation between pairwise physical distance and pairwise sample similarity for any data type within the coral colonies (16S rDNA, 18S rDNA, metatranscriptomes, or metabolomes) (all R 2 s < 0.04 and all p values > 0.38, Figure 2). Additional sampling may have improved these correlations; however, a power analysis revealed that 28 samples would be required to see a significant trend (p < 0.05) using a strong effect size (0.5) and a power of 80% (Supplementary Table 1). Since 21-29 samples for each data type were used in this study, additional sampling would not likely improve the weak correlations observed between physical distance and sample similarity (Supplementary Table 2).

Diversity of Bacterial Communities: 16S rDNA
A total of 656,698 16S rDNA sequences were generated with an average of 24,322 sequences per sample (n = 27). A total of 3,034 amplicon sequence variants (ASVs) were identified, and ASVs identified as chloroplasts, mitochondria, and long branching chimeras were removed. Overall, 16S rDNA sequencing demonstrated that samples were not spatially structured (p = 0.42, Figure 2A) but were rather found to cluster by community type (Figures 3A,B; Supplementary Figure 1). At the ASV level, there was significant clustering of coral samples separate from all other samples (Figures 3A,B PERMANOVA p < 0.003) as well as significant differences between corals and interface samples (Figures 3A,B, PERMANOVA p = 0.003). This was largely driven by the differences in the O. faveolata microbiome compared to the microbiome of turf algae and interface samples ( Figure 3B, PERMANOVA p = 0.02 and 0.02, respectively). Similar trends were observed at the phylum level, though PERMANOVA p-values were generally less significant than at the ASV level (Supplementary Figure 2). Furthermore, when comparing between colony versus within colony variance, there was a significant difference in 16S rDNA annotations between colonies at the ASV level ( Figure 3B PERMANOVA p = 0.04), but not at the phylum level (Supplementary Figure 2 PERMANOVA p = 0.73). Interestingly, the bacterial microbiome of the two hydrozoan competitor samples did not cluster separately as an outgroup, but were rather more similar to turf algae than to O. faveolata (Figures 3A,B).
At the ASV level, there was an overall trend for the coral microbiome to be less diverse than interfaces or competitor (turf algae and M. complanata) microbiomes ( Figure 3D), with coral samples having the lowest Shannon entropy (4.43 ± 0.92), followed by the interface (5.11 ± 2.34), and finally the competitor samples (5.31 ± 1.55). It is also notable that turf algal samples had the highest entropy (6.24 ± 0.65) of all sample types. An even stronger trend of diversity was observed for rarefaction species richness (referred to as richness), with coral samples being significantly less rich than interface (p = 0.002) or competitor samples (p = 0.01). This pattern was also observed at the phylum level (Supplementary Figure 2).
Random forests analysis demonstrated that the bacteria in the phyla Firmicutes (specifically Clostridium), Bacteroidetes (specifically Flavobacterium) and Proteobacteria (specifically Alphaproteobacteria) were among the most important bacterial ASVs in predicting the sampling area (i.e., coral, interface, competitor). This was further supported by 9 of the top 10 most important random forests predictors at the ASV level coming from these same three bacterial phyla (i.e., Firmicutes, Bacteroidetes, and Proteobacteria) (Supplementary Figure 3). Interestingly, many of the top BLAST hits were to previously sampled O. faveolata (Supplementary Figure 3).
The ubiquity-abundance analysis revealed that a single ASV, in the genus Clostridium was uniquely ubiquitous in O. faveolata (i.e., a Clostridium sp. was found in all O. faveolata samples and only in O. faveolata samples, Figure 3C). The relative abundance of the Clostridium ASV (d6d6d57d15072e10d222a16f3e6d97e9) in O. faveolata samples ranged from 1.31% to 0.03% with an average of 0.47% ± 0.45%. The Clostridium ASV was also one of the top ten predictors in the random forest analysis, indicating the potential importance for this ASV in the core microbiome of O. faveolata (Supplementary Figure 2). BLAST results showed that the Clostridium ASV shared 100% sequence identity with an uncultured bacterium clone collected from O. faveolata in Panamá in 2008 (Sunagawa et al., 2010). The Clostridium ASV was highly abundant and ubiquitously distributed across all O. faveolata samples, but absent from any and all interface or competitor samples (Video S1). The phylum Firmicutes, which includes Clostridium, also showed a similar distribution across the coral colonies and their competitors, however Firmicutes was highly abundant in all corals including O. annularis (Figure 4). Other phyla, such as Cyanobacteria, displayed the opposite pattern, with low relative abundances in O. faveolata and higher relative abundances in turf algal competitors and interfaces ( Figure 4C).
Several rare bacterial taxa were also present at both sites. Four unique ASVs belonging to Endozoicomonas, a symbiont of marine animals, were found at low relative abundances (0.02-0.99%) in the O. faveolata colonies from both sites, and the bacterial symbiont was present in 60% of the coral samples. Five different Endozoicomonas ASVs were present in the hydrozoan competitor, M. complanata, at higher relative abundances (0.65-7.49%), suggesting O. faveolata and M. complanata specific symbionts. ASVs belonging to the understudied TM6 phylum were also identified at low relative abundances (0.01-0.26%), specifically from two interface samples and one turf

Diversity of Microbial Eukaryotes: Non-metazoan 18S rDNA
A total of 2,972,176 18S rDNA sequences were generated with an average of 102,488 sequences per sample (n = 29). A total of 541 ASVs were identified and ASVs identified as bacteria and long branching chimeras were removed. Symbiodiniaceae, specifically Cladocopium spp., was the most abundant taxon in coral samples, whereas red (Florideophyceae) and green (Ulvophyceae) algae were the dominant taxa in turf algae samples (Supplementary Figure 4). ASVs belonging to these taxa were also the main drivers of the clustering observed in the PCoA and NMDS (Supplementary Figures 1, 4). The majority of the ubiquitous ASVs unique to O. faveolata were also Symbiodiniaceae (Supplementary Figure 5), and a phylogenetic analysis showed no clustering of ASVs based on coral colony or coral species (Supplementary Figure 6). However, Symbiodiniaceae ASVs from M. complanata samples formed a well-supported clade that was distinct from the coral Symbiodiniaceae clade, suggesting hydrozoan-specific and coral-specific eukaryotic symbionts (Supplementary Figure 6). BLASTn results showed that the coral Symbiodiniaceae ASVs shared a 99.01% identity with several Clade C symbionts (Cladocopium spp.) whereas M. complanata ASVs had a 98.77% identity to Clade A symbionts (Symbiodinium spp.).
Corallicolid, an intracellular apicomplexan symbiont of coral, was present at low abundances in three coral samples including the competing O. faveolata colony at East Point (northern transect, sample V, relative abundance = 0.075%), the competing O. annularis colony at East Point (southern transect, sample V, relative abundance = 0.29%) and the interface between O. faveolata and M. complanata at Water Factory (western transect, sample III, relative abundance = 0.36%). This symbiont has been reported in O. faveolata, O. annularis, and several other cnidarians (Kwong et al., 2019), but its presence in some coral colonies and absence from nearby conspecific colonies suggests a sporadic association with specific coral species.

Function of Bacterial Communities: Metatranscriptomes
Similar to 16S and 18S rDNA, the functional annotation of bacterial metatranscriptomes also demonstrated significant clustering of samples by host organism, with O. faveolata samples being significantly different than interface (PERMANOVA p < 0.01) and competitor samples (PERMANOVA p < 0.001) (Figures 5A,B and Supplementary Figure 1). In stark contrast to 16S rDNA data where the hydrozoan outgroup (pink triangles in Figure 5A) grouped with turf algae, the bacterial transcriptomes group hydrozoan samples with coral, suggesting a potential taxonomic, functional decoupling in these microbiomes. Of the 5,206 functionally annotated bacterial RNA transcripts from all host organisms, 357 (6.86%) were ubiquitous to all samples, and no functions were unique to O. faveolata; a strong contrast to the 16S rDNA results (Figure 3C). Of the ubiquitous transcripts, the NAD(P)H-quinone oxidoreductase function was the most abundant, with a mean relative abundance of 12.4% in all samples; furthermore, NAD(P)H-quinone oxidoreductase was also shown to be the most important predictor for determining sample type in the random forests classification analysis (Supplementary Figure 7). FIGURE 4 | Heatmap and 3D cartography of bacterial phyla from Orbicella faveolata colonies, competitors and interfaces. (A) Heatmap of phyla selected from the supervised random forests analysis. Relative abundances are log-transformed and colors are normalized to columns with red representing the lowest relative abundance and white representing the highest. Asterisks highlight phyla shown in 3D models. 3D coral models mapped with the relative abundances of (B) Firmicutes and (C) Cyanobacteria from both sampling sites. The relative abundance of Firmicutes ranges from 0.21-93.5% at East Point and 0.0-89.7% at Water Factory. For Cyanobacteria, the relative abundance ranges from 0.0-14.4% at East Point and 0.0-0.49% at Water Factory.
Random forests analysis revealed that microbial transcripts were enriched in phage-encoded photosystems in coral samples and enriched for transcripts involved in bacterial cell division, motility, and chemotaxis in turf algal and interface samples (Figure 6). At all levels of the SEED hierarchy (i.e., Levels 1,2,3, and function) phage related genes were found to be significantly enriched in coral samples and a driver of the clustering seen in the data (Figure 6). Specifically, we found the photosystem II protein D1 (PsbA) transcript to be driving the T7-like cyanophage pattern at higher levels (Supplementary Figures 8, 9). This has been visualized in a two way dendrogram in Figure 6A and via a 3D molecular heatmap in Figure 6B. Other random forests top predictors included bacterial respiration along with cell regulation and signaling, both of which were more abundant in coral than the interface or competitors. The opposite pattern was observed for functions involving membrane transport, cell wall and capsule, iron acquisition and metabolism, virulence, dormancy and sporulation, and phages (prophages, transposable elements, plasmids) ( Figure 6A).

Associated Biochemicals: Metabolomes
The metabolomic data set consisted of 2,315 unique features of which 53 (2.28%) had a spectral match to a known compound in the GNPS database with another 230 analog hits (i.e., contains a known mass shift and has a similar spectrum). Similar to 16S and 18S rDNA, and transcriptomic data, metabolites   (Quinn et al., 2016). The C18 Lyso-PAF metabolite (m/z 510.39) was also the second top predictor of coral, turf algae, and interface groups in the random forests analysis on all metabolites (Figure 8 and Supplementary Figure 10).
The ubiquity-abundance analysis revealed that out of 2,315 metabolites, 45 were found in all samples (1.94%) while 219 and 39 metabolites were unique to coral (O. faveolata and O. annularis) and turf algae, respectively. Of the coral specific metabolites, 136 were only present in O. faveolata. (Figure 7C). Of those compounds found only in O. faveolata, the FIGURE 6 | Molecular heatmap of bacterial metatranscriptomes from Orbicella faveolata colonies, competitors and interfaces. (A) Heatmap of SEED Level 1 categories selected from the supervised random forests analysis. Phage Group 1 includes phages, prophages and transposable elements while Phage Group 2 also includes plasmids. Colors are normalized to columns with red representing the lowest relative abundance and white representing the highest. Asterisks highlight Level 1 categories shown in 3D models. 3D coral models mapped with the relative abundance of bacterial transcripts belonging to (B) Phage Group 1 and (C) Cell Cycle and Cell Division at both sampling sites. The relative abundance of Phage Group 1 transcripts range from 0.06-0.74% at East Point and 0.002-0.42% at Water Factory. For Cell Cycle and Cell Division, the relative abundances range from 0.03-0.53% at East Point and 0.02-0.28% at Water Factory.
cyanobacterial-produced antifungal compound lobocyclamide (dereplicator hit) and the Streptomyces-produced antimicrobial coleimycin (analog hit) were among the most ubiquitous compounds (Supplementary Figure 10). Five unknown metabolites with mass to charge ratios of 467.31, 240.10, 432.34, 1149.74, and 397.21 were also uniquely ubiquitous in O. faveolata (ubiquity = 0.92). Unlike 16S, 18S, and transcriptome samples, the random forests analysis of metabolites did not find the ubiquitous molecules among the top compounds for predicting sample type. Instead, the random forests analysis identified two fatty acid esters (ClassyFire annotation) connected in the network, a tyrosine-proline dipeptide (analog hit), and a ceramide (analog hit) to be the top predictors of sample type (Figure 8 and Supplementary Figure 9).

Relationship Between Diversity of Molecular Data Types
There were no significant correlations between the entropy (H') (all R 2 s ≤ 0.07 and p-values ≥ 0.33, Supplementary Figure 11) or richness (all R 2 s ≤ 0.02 and all p-values ≥ 0.52, Supplementary  Figure 11) among different data types, indicating potential decouplings in taxonomic, functional, and biochemical profiles. An analysis of all data types sampled within the 2.5 × 2.5 m reef site at East Point (6.25 m 2 ) and a 1 m × 1 m reef area at Water Factory (1 m 2 ) also showed differences in entropy between sites. The Water Factory site had the most diverse 18S rDNA (entropy = 3.71) and bacterial function (transcriptomes, entropy = 6.22), whereas East Point had the highest diversity of 16S rDNA (entropy = 4.33) and metabolites (entropy = 6.75) (Supplementary Figure 12). Furthermore, when all samples were considered, metabolomes were seen to be the most diverse of all data types (entropy = 6.76), while 18S rDNA was the least diverse (entropy = 3.83) (Supplementary Figure 12).
The conditional probabilities of several microbial transcripts and metabolites varied between Firmicutes and Bacteroidetes (Figure 9), two of the main phyla identified in the previous analyses. Dormancy and sporulation along with cell division and cell cycle had a greater conditional probability with Bacteroidetes (0.270 and 0.039) relative to Firmicutes (0.010 and 0.007), whereas the conditional probability of phosphorus metabolism, RNA metabolism, and virulence was higher with Firmicutes (0.134, 0.060 and 0.097) than with Bacteroidetes (0.027, 0.004 and 0.002). The conditional probability of Phage group 1 (phages, prophages and transposable elements) with Firmicutes (0.040) was also higher than with Bacteroidetes (0.009). For metabolites, glycerophosphoethanolamine and cis-9hexadecenoic acid (8S-HETrE) had a higher probability with Firmicutes (0.019 and 0.001) than with Bacteroidetes (0.002 and 7.04 × 10 −4 ). Firmicutes and these metabolites were also more abundant in corals than in other organisms (Figures 4,  8), highlighting a potential association between the coral microbiome and metabolome.

DISCUSSION
Coral and algal competition has been studied extensively utilizing molecular techniques, such as amplicon sequencing, metagenomics, and metabolomics (Barott et al., 2012;Haas et al., 2016;Quinn et al., 2016;Pratte et al., 2018;Silveira et al., 2019;Clements et al., 2020;Roach et al., 2020). Here, we integrate bacterial and protistan metabarcoding, metatranscriptomics, and metabolomics across 3D models of the endangered reef-building coral Orbicella faveolata and its competitors. This approach allowed us to integrate and visualize molecular data in a spatial framework to investigate spatial and biological patterns within and between O. faveolata colonies and their competitors (Figure 1).
All molecular profiles (DNA, RNA, and metabolites) revealed ecologically relevant differences between the O. faveolata colonies, the interface, and the corresponding benthic competitors (Figures 2-8). Holobionts clustered by microbial composition (16S and 18S), bacterial gene function (transcriptomes), and metabolites (Figures 3, 5, 7). Many uniquely ubiquitous microbes and metabolites, and sporadic microbial symbionts were found to be associated with O. faveolata. No spatial patterns were observed, suggesting that a single coral colony is a diverse and dynamic consortium of host, microbes, and associated functions/metabolites (Figure 2). However, larger sampling efforts may be warranted to investigate this further and the scale of sampling should be taken into consideration. We also utilized machine-learning algorithms to compare thousands of variables from different data types (Figures 4, 6, 8) and found potential associations between FIGURE 10 | Summary of the Orbicella faveolata holobiont and its main competitor, turf algae. Core microbial members of O. faveolata are shown along with their potential functions and interactions with the coral host (red circle). Question marks represent functions that are unknown. The change in abundance of specific bacterial groups, bacterial and phage functions, and metabolites across coral and turf algae competitors are illustrated by orange and green triangles. molecular functions and various members of the coral holobiont (Figures 9, 10).

Uniquely Ubiquitous and Sporadic Coral Symbionts
A putative Clostridium sp. ASV was uniquely ubiquitous to the O. faveolata colonies, pointing to it as a core member of the O. faveolata holobiont (Figure 3C and Supplementary Video 1). Our findings are corroborated by a previous study that reported this same Clostridium sp. ASV to be present in O. faveolata from Panamá in 2008 (Sunagawa et al., 2010), which suggests stability of the ASV as an O. faveolata symbiont in the Caribbean Sea. Clostridium are commonly found in corals, likely inhabiting hypoxic portions of the O. faveolata holobiont. They have been found in a wide variety of coral species across large geographic spatial scales, where they may occupy low-oxygen mucosal microniches generated by the breakdown of complex carbon (McKew et al., 2012). Our findings highlight the need for future studies on compartmentalization of the coral holobiont, for example, using multi-omics with respect to different anatomical regions of the coral (Sweet et al., 2011).
Other microbial symbionts showed differential associations with the coral holobionts and their respective competitors. Endozoicomonas, a bacterial symbiont of marine animals (Neave et al., 2016), was found at low abundances in the majority of O. faveolata samples at both sites. Conversely, Endozoicomonas were highly present in the hydrozoan competitor, M. complanata. Endozoicomonas ASVs were unique to each cnidarian host, suggesting at least two species of Endozoicomonas, each specific to O. faveolata and M. complanata. Endozoicomonas is likely involved in the cnidarian sulfur cycle (Tandon et al., 2020), and may play a relevant functional role in the holobionts. Corallicolid, an apicomplexan coral symbiont, was also present in some coral colonies and absent from nearby conspecific colonies, showing a potential sporadic association with specific coral species. At East Point, corallicolids were found in the O. faveolata and O. annularis colonies interacting with the coral of interest, while it was present in only one O. faveolata interaction sample at the Water Factory site. Corallicolids are known to infect a diverse range of cnidarian hosts, however, their ecological role remains unknown (Kwong et al., 2019), and additional sampling is required to determine the apicomplexan's association with various coral species.
Several bacterial patterns at the phylum level were observed to be specific to the coral holobionts. Broadly, abundances of the phylum Firmicutes, which contains the genus Clostridium, were higher in the coral samples, and abundances of Bacteroidetes were higher in turf algae samples resulting in significantly higher Bacteroidetes-to-Firmicutes ratios at the interface between coral and turf algae (Figure 4). Interestingly, Bacteroidetes-to-Firmicutes ratios have been linked to coral health; in specific, this ratio has been shown to be elevated in coral holobionts that are losing against turf algal competitors and also in disease states (Closek et al., 2014;Roach et al., 2017Roach et al., , 2020. The differential abundances of these phyla have been studied more in depth with regards to human health, where the gut microbiome profile of obese individuals is linked to a decrease in Bacteroidetes-to-Firmicutes ratios (Turnbaugh et al., 2006). This suggests that the correlation between the relative abundances of Bacteroidetes and Firmicutes, and host physiology is potentially specific to the animal host, and likely has to do with the delicate balance of anaerobic and aerobic respiratory processes.

Relevant Transcripts and Metabolites
Other bacterial phyla, such as Cyanobacteria, displayed relatively low abundances in coral samples compared to the interface and competitor samples. Conversely, transcripts from a T7-like cyanophage (family Podoviridae, Cyanopodovirus), specifically, photosystem II psbA transcripts, were found in high abundance in coral samples compared to interface and competitor samples (Figures 4, 5). Other coral studies have found abundant cyanobacterial photosystem genes, including psbA and psbD, that were potentially acquired through phage infections (Weynburg et al., 2017). Viral-mediated horizontal gene-transfer of photosystems genes, and other genes, from and to the bacterial host has been observed (Chénard and Suttle, 2008;Thompson et al., 2011;Little et al., 2020); furthermore, a bacteriophageadherence-to-mucus model suggests that phage interact with mucosal surfaces, such as the surface of a coral, through outer capsid protein domains (Hoc) previously identified in Cyanophage (Sullivan et al., 2005;Barr et al., 2013;Silveira and Rohwer, 2016). Our observation of the expression of T7like cyanophage photosystem genes may suggest these viruses are active in coral mucus and possibly a widespread feature of coral holobionts.
Metabolomic analysis highlighted 136 known molecules as being unique to the O. faveolata holobiont (Figure 7C). Among the metabolites unique to O. faveolata, the bacterial produced metabolites, lobocyclamide, a known antifungal compound produced by Cyanobacteria (MacMillan et al., 2002), was one of the most ubiquitous (Supplementary Figure 10). Bacterial-produced antifungals have been found previously in Gorgonian octocorals, where Pseudoalteromonas sp. have been shown to produce them in a light-dependent manner (Moree et al., 2014). Another uniquely ubiquitous O. faveolata metabolite was the Streptomyces-produced antimicrobial coleimycin (Supplementary Figure 10). The significant enrichment of these bacterial-produced antimicrobials highlights the potential for bacterial modulation of the host microbiome.
Other compounds found to be elevated in the O. faveolata holobiont included several putatively bioactive lipids such as lyso-platelet activating factor and a ceramide analog (Figures 7,  8, 9). These lipids can serve as modulators of the coral immune system, and have been shown to play a role in corals' response to various stressors including turf algae overgrowth and temperature stress (Quinn et al., 2016;Galtier d'Auriac et al., 2018;Roach et al., 2020Roach et al., , 2021a. Our findings highlight the importance of these and other bioactive lipids and bring into question the role of the coral immune system in ecological interactions (Quistad et al., 2014).

Diversity and Spatial Patterns
Overall, corals had lower diversity of microbial composition and function, but displayed higher metabolite diversity compared to interfaces and turf algae. This is consistent with previous studies on coral and algae interactions where the bacterial consortium within corals was found to be less taxonomically rich and diverse than algal holobionts ( Figure 3D; Barott et al., 2011). In terms of spatial patterns of the taxonomic, functional, and metabolite data, we observed no evidence for significant isolation by distance or spatial autocorrelation (Figure 2). Several other environmental factors may impact spatial patterns within a coral colony, such as water movement over the coral surface (e.g., boundary layers) or incident light; molecular factors such as gene regulation and post-translational modifications could also affect this spatial patterning.

Co-occurrence Networks Linking Omics Datasets
Understanding the distribution of various components of the holobiont in context to each other remains a largely under investigated aspect of molecular ecology in host-associated systems. To this end, an artificial neural network approach (mmvec) was used to calculate the conditional probabilities of the multi-omics annotations produced in this study (Morton et al., 2019). The network analysis revealed strong cooccurrences between Cyanobacteria and photosynthesis along with Planctomycetes and nitrogen metabolism transcripts (Figure 9). Planctomycetes is one of the only known bacterial groups to perform anaerobic ammonium oxidation (anammox), a process that significantly contributes to the nitrogen cycle (Jetten et al., 2009). Interestingly, increased levels of bacterial nitrogen fixation in the coral holobiont have been hypothesized to contribute to coral disease states and disrupture in the stability between corals and their endosymbiotic partners (Shashar et al., 1994;Santos et al., 2014;Rädecker et al., 2015). Because most Planctomycetes have strong associations with macroalgae, lack peptidoglycan walls (e.g., exhibit resistance to antimicrobials produced by other bacteria), and possess holdfasts for attachment, it is probable that members of this phylum contribute to processes associated with algal competition and coral bleaching. Additionally, the highest co-occurrence found was between transcripts for dormancy and sporulation and the bacterial phylum, Chlamydiae (Figure 9). All members of Chlamydiae are obligate intracellular symbionts of eukaryotes and infect a wide range of hosts including protists and animals (Horn et al., 2000;Roulis et al., 2013). The relationship between these symbionts and dormancy and sporulation functions provides insight into the lifecycle of these intracellular bacteria. The human pathogen, Chlamydia, can exist in a metabolically reduced state known as an elementary body that allows for longterm survival in nutrient-poor conditions (Hogan et al., 2004). This 'spore-like state' would likely be common on coral reefs where host conditions vary due to the fluctuating environment. However, the potential hosts of these specific Chlamydiae symbionts remain unknown, as is the case for the majority of Chlamydiae (Dharamshi et al., 2020).
The neural network analysis also linked the presence of specific transcripts and metabolites with Firmicutes and Bacteroidetes. As previously mentioned, Firmicutes were highly abundant in coral samples (Figure 4), and this group co-occurred with phage, prophage, and transposable element transcripts (Figure 9) which were also highly abundant in coral samples (Figure 6). On the other hand, cell division and cell cycle transcripts were more abundant in interface and algal samples, and the transcripts co-occurred with Bacteroidetes (Figure 9). These patterns where Bacteroidetes and cell division/growth associate with the interface and algae, whereas Firmicutes and phage related functions associated with coral provide further evidence of the functional roles of these taxonomic groups. The co-occurrence of glycerophosphoethanolamine and cis-9-hexadecenoic acid (8S-HETrE), two highly abundant coral metabolites, with Firmicutes (Figure 9) also suggests a potential association between the coral microbiome and metabolome (Figure 10).

Current Limitations and Future Directions
Here, we highlight the novelty and robust power of incorporating multi-omics into natural history and ecological frameworks (Figure 10). However, among the limitations of our findings is the use of database-dependent analyses. In the future, tools such as artificial neural networks (Mendez et al., 2019;Cantu et al., 2020) can be used in bioinformatic and chemoinformatic approaches to shed light on sequence and mass-spectrometry data that is absent in repositories. Furthermore, less invasive sampling techniques could be used to capture molecular dynamics over relevant temporal scales (e.g., before and after stress events). Our data reveals important in situ molecular ecology of the coral holobiont, and follow-up studies should extend these sampling schemes to entire reefs to investigate holobionts over multiple spatial and biogeographic scales.
This study focused on investigating the spatial and ecological drivers of two coral colonies in unprecedented detail using four different omics techniques combined with photogrammetric spatial reconstruction. The lack of spatial patterns observed is potentially a result of other uninvestigated variables such as environmental factors or boundary layer processes. The lack of spatial structuring at this scale may also be due to the resolution and scale of sampling. For example, samples could be collected more densely or distributed over larger areas to investigate these trends at other scales such as host-associated microenvironments or across whole reef systems. Because these methodologies are becoming more affordable (i.e., the cost of DNA/RNA sequencing and mass spectrometry), scaling these approaches to address holobiont spatial aspects at the whole reeflevel or even the level of island chains or ocean basins is now feasible. In addition, 3D image reconstruction has been applied at the reef-scale and is relatively straightforward to implement (e.g., Edwards et al., 2017;Roach et al., 2021b). Future molecularmapping methods combined with environmental and temporal data will elucidate additional functions performed by the diverse members of the coral holobiont that, in turn, influence molecular to reef-scale processes.

CONCLUSION
The data herein revealed that microbial composition, function, and molecular profiles are strong drivers of differences between coral and algae holobionts, and at the spatial scales examined, no evidence of spatial patterns within the O. faveolata colonies were observed for the various data types. This dataset revealed uniquely ubiquitous microbes and metabolites, along with potentially sporadic microbial symbionts associated with this endangered reef-building coral. The 16S rDNA data from O. faveolata demonstrated that a single ubiquitous Clostridium ASV was present in all colony samples, and abundant T7like cyanophage transcriptional functions were also found in O. faveolata, adding new insight into the viral fraction of the coral holobiont. Metabolomics profiling displayed high levels of bioactive lipids in the O. faveolata samples, highlighting their role in coral immune responses. In addition, the co-occurrence of bacterial taxa, microbial transcripts, and small molecules revealed specific bacterial phyla with strong associations to various metabolic functions, such as Planctomycetes and nitrogen metabolism along with Chlamydiae and dormancy functions. These methodologies provide a new framework to examine holobiont structuring and the molecular ecology of the coral holobiont.

DATA AVAILABILITY STATEMENT
For metabolomes, the .mzXML files are available on the MassIVE database (massive.ucsd.edu) under numbers MSV000080597 and MSV000080632 (same dataset). 16S rDNA amplicon, 18S rDNA amplicons, and total RNA transcriptome sequence data are deposited and publicly available on the Sequence Read Archive (SRA -https://www.ncbi.nlm.nih.gov/sra) under BioProject: PRJNA692629.

AUTHOR CONTRIBUTIONS
ML, EG, MA, FR, and TR designed research. FR, TR, EG, SB, and BM conducted fieldwork and collected samples. ML, EG, MA, JS, SB, JH, ZQ, VB, BW, DP, TR, and FR performed research. ML, EG, MA, RQ, PK, PD, FR, and TR contributed new reagents and analytic tools. ML, EG, MA, JS, ZQ, and TR analyzed data. ML, EG, and TR wrote the manuscript. MA, JS, SB, JH, ZQ, VB, MR, AG, DP, CS, AH, LK, MV, RQ, PK, PD, and FR provided edited versions of the manuscript. All authors contributed to the article and approved the submitted version.

FUNDING
We would like to thank the following organizations for funding this research: the National Science Foundation for Grants G00009988 (to TR), the Gordon and Betty Moore Foundation for Grants #3781 (to FR), #9702 (to FR), and #9201 (to PK), and the University of British Columbia (International Doctoral Fellowship to EG).

ACKNOWLEDGMENTS
We would like to thank the staff of CARMABI research station.

SUPPLEMENTARY MATERIAL
The Supplementary Figure 11 | Regression analyses of metabolomic, 16S rDNA amplicons, and microbial functional Shannon Entropy (H') and Richness (s). 16S amplicon relative abundances for these analyses were conducted from ASV relative abundances. Microbial functions for these analyses were from SEED subsystems functional level relative abundances. No significant correlations were found across the three data types. For metabolomes, relative abundance of all hits were used. Blue lines represent the fit line and blue areas represent the fit confidence (95%).
Supplementary Figure 12 | The range of Shannon entropy of all sample types from East Point, Water Factory and both sites combined. The size of the circles represent the reef area sampled at East Point (6.5 m 2 ), Water Factory (1 m 2 ), and total area (7.5 m 2 ).
Supplementary Figure 13 | Heatmaps from the neural network (mmvec) analysis. The conditional probability (dark red for low and yellow/white for high) that the presence of a metabolite co-occurs with a bacterial phylum or transcript is shown. (A) All metabolites related to all bacterial phyla. (B) All metabolites related to all SEED level 1 transcripts.
Supplementary Table 1 | Power analysis of the sampling size (n) required to be statistically significant (p ≤ 0.5) with variable effect sizes (r) and a power of 80%.
Supplementary Table 2 | Power analysis using sample sizes (n) from this study to determine the power at variable effect sizes (r) that are statistically significant (p ≤ 0.05).
Supplementary Table 3 | Sequence mapping identification results including the percentage of reads assigned to microbial subsystems, coral reference genomes, and Symbiodinium reference genomes.