Metabolomics of Healthy and Stony Coral Tissue Loss Disease Affected Montastraea cavernosa Corals

Stony coral tissue loss disease, first observed in Florida in 2014, has now spread along the entire Florida Reef Tract and on reefs in many Caribbean countries. The disease affects a variety of coral species with differential outcomes, and in many instances results in whole-colony mortality. We employed untargeted metabolomic profiling of Montastraea cavernosa corals affected by stony coral tissue loss disease to identify metabolic markers of disease. Herein, extracts from apparently healthy, diseased, and recovered Montastraea cavernosa collected at a reef site near Ft. Lauderdale, Florida were subjected to liquid-chromatography mass spectrometry-based metabolomics. Unsupervised principal component analysis reveals wide variation in metabolomic profiles of healthy corals of the same species, which differ from diseased corals. Using a combination of supervised and unsupervised data analyses tools, we describe metabolite features that explain variation between the apparently healthy corals, between diseased corals, and between the healthy and the diseased corals. By employing a culture-based approach, we assign sources of a subset of these molecules to the endosymbiotic dinoflagellates, Symbiodiniaceae. Specifically, we identify various endosymbiont- specific lipid classes, such as betaine lipids, glycolipids, and tocopherols, which differentiate samples taken from apparently healthy corals and diseased corals. Given the variation observed in metabolite fingerprints of corals, our data suggests that metabolomics is a viable approach to link metabolite profiles of different coral species with their susceptibility and resilience to numerous coral diseases spreading through reefs worldwide.


INTRODUCTION
Corals are holobionts consisting of the coral animal, photosynthetic dinoflagellate endosymbionts (family Symbiodiniaceae) commonly referred to as zooxanthellae, and a complex microbiome (Rohwer et al., 2002;Maire et al., 2021). Collectively, members of the holobiont contribute to the fitness of the coral animal (Rohwer et al., 2002;Mieog et al., 2009;Gordon and Leggat, 2010;Parkinson et al., 2015). Coral reefs provide a variety of important services, including hosting diverse marine ecosystems (Plaisance et al., 2011;Fisher et al., 2015), supplying pharmacophores for drug development (Bruckner, 2002;Sang et al., 2019), protecting shorelines (Reguero et al., 2021), and contributing to coastal economies with an estimated collective worth of $3.4 billion in the US (Brander and Van Beukering, 2013). Factors including increased ocean temperatures, ocean acidification, pollution, and disease outbreaks present continual challenges to the health of corals (Richardson, 1998;Rosenberg and Ben-Haim, 2002;Rogers and Weil, 2010;Pandolfi et al., 2011;Hoegh-Guldberg et al., 2017;Grottoli et al., 2018;Montilla et al., 2019;Howells et al., 2020). A combination of environmental stressors can increase the vulnerability of corals to disease outbreaks (Maynard et al., 2015;Benkwitt et al., 2020). Coral diseases are notoriously difficult to characterize and the pathogenic agents usually remain unknown (Mera and Bourne, 2018;Vega Thurber et al., 2020). One of the challenges with pathogen-identification is the inherent withinspecies (intra) and between-species (inter) variability to disease susceptibility. Understanding intraspecies and interspecies variability in genotypes, metabolomes, microbiomes, and other physiological factors such as bleaching history of field corals is critical to investigating their contribution to disease resistance and susceptibility, which are essential to effective management and restoration efforts.
Caribbean coral reefs are experiencing an ongoing outbreak of a tissue loss disease named stony coral tissue loss disease (SCTLD), which first appeared along the Florida Reef Tract in 2014 (Miller et al., 2016;Walton et al., 2018;Muller et al., 2020) and continues to spread throughout various countries in the Caribbean (Precht et al., 2016;Alvarez-Filip et al., 2019;Kramer et al., 2020;Muller et al., 2020;Estrada-Saldívar et al., 2021;Heres et al., 2021;Thome et al., 2021). At least 24 different coral species are susceptible to SCTLD, however, intraand interspecific differences have been observed (Precht et al., 2016;Aeby et al., 2019;Ushijima et al., 2020;Walker et al., 2020;Meiling et al., 2021). Reefs affected by SCTLD can have coral community composition altered within a few months of initial observation (Walker, 2018;Estrada-Saldívar et al., 2020), highlighting the imminent threat of this disease. While the etiological agent remains unknown (Landsberg et al., 2020;Neely et al., 2020), host-endosymbiont harmony seems to be disrupted Landsberg et al., 2020), disease progression is suggested to be mediated by bacteria , and treatment with antibiotics is demonstrated to halt lesion progression for a majority of corals in Florida (Neely et al., 2020;Shilling et al., 2021;Walker et al., 2021). Histological studies indicate some characteristic features of SCTLD (Landsberg et al., 2020); however, without knowledge of the etiological agent there is currently no feasible diagnostic test to confirm SCTLD in diseased corals.
The emergence of -omics studies offers a unique opportunity to elucidate the complex relationships within the members of the coral holobiont at the genetic and chemical level (Joyce and Palsson, 2006;Gordon and Leggat, 2010;Ramos-Silva et al., 2013;Ainsworth et al., 2015;Dixon et al., 2015;Kenkel and Matz, 2016;Drake et al., 2018). Through microbiome sequencing, recent studies described shifts in corals' microbial community composition upon SCLTD infection Rosales et al., 2020;Becker et al., 2021). However, the metabolomes, the collection of small molecules called metabolites, of these holobionts have not been evaluated. Metabolomics-based approaches provide insight into an organism's physiological and biochemical state at the time of sample collection (Goodacre, 2007;Viant, 2008;Patti et al., 2012;Wishart, 2019). Thus, metabolomics profiling has been employed to study normal physiology and pathophysiology of diseases. When studying corals, this endeavor is complicated because metabolomic profiles are influenced by the host genotype, the endosymbiont genotype (Symbiodiniaceae), the associated microbiome, the environment, previous history of bleaching, and disease (Muscatine, 1990;Baker, 2003;Roth, 2014;LaJeunesse et al., 2018). Recent application of metabolomics studies on coral holobionts and their environments have demonstrated the applicability of this approach (Sogin et al., 2014;Ochsenkühn et al., 2018;Lohr et al., 2019;Matthews et al., 2020;Roach et al., 2020;Sikorskaya and Imbs, 2020). Direct metabolomics of symbionts have provided insights into their function in the holobiont (Gordon and Leggat, 2010;Hillyer et al., 2017;Sogin et al., 2017;Rosset et al., 2019;Sikorskaya et al., 2021).
Metabolomics studies conducted to date have enabled exploration of changes in metabolite profiles of coral holobionts to specific stressors including thermal bleaching (Quinn et al., 2016b;Sogin et al., 2016;Hillyer et al., 2017;Roach et al., 2021;Williams et al., 2021), but coral metabolomics largely remains an underdeveloped field (Sogin et al., 2014). Two analytical techniques commonly utilized in coral metabolomics include nuclear magnetic resonance (Andersson et al., 2019) and liquid chromatography (LC) coupled to mass spectrometry (MS) (Gordon et al., 2013). Due to the high sensitivity of LC-MS-based approaches and the resulting increased breadth of detected metabolites, LC-MS is a particularly attractive approach toward characterization of metabolome profiles (Lohr et al., 2019). As threats to coral reefs grow, it is important to understand differences in intraspecific chemotypes and genotypes as well as the factors driving these differences. Metabolome studies enhance the understanding of pathogenesis, illuminate the chemical relationships between the coral host, endosymbionts, and other associated microbes, and aid in discovery of biomarkers. There are currently no published studies that compare healthy and SCTLD-infected coral metabolomes or describe the intraspecies variability. Hence, we employed untargeted LC-MS metabolomics to reveal differences in the metabolomic profiles of apparently healthy Montastraea cavernosa corals as well as differences between apparently healthy coral colonies and colonies with disease lesions from tissue samples collected at an endemic site in Broward County, FL. Additionally, we used a culture-based approach to coanalyze the metabolomes of three Symbiodiniaceae species to reveal that at least a third of the top metabolites separating the field coral sample groups belonged to Symbiodiniaceae. We identify a number of Symbiodiniaceae metabolites using a suite of recently developed metabolite feature annotation tools. This study highlights that the culture-based approach, when combined with metabolomics, facilitates an understanding of the contribution of different holobiont components to the metabolomic profiles, which can aid in teasing apart the source of metabolites that dictate the separation in metabolomic profiles. We propose that in addition to the prominent contribution of Symbiodiniaceae to metabolomes of field corals, careful cataloging of host genotype and past environmental insults is important to understand intraspecies differences of field corals.

Site Selection
An experimental site was set up on March 23, 2020 located 475 m from shore (26 • 9.05 N, 80 • 5.75 W), on the Broward-Miami Ecoregion Ridge-Shallow habitat (Walker et al., 2008;Walker, 2012) in about 7 m depth (Figure 1). Nineteen visibly diseased Montastraea cavernosa colonies were tagged, photographed, and mapped across a 100 m × 35 m area of hardbottom using a floating GPS to obtain coordinates for all colonies. Healthylooking colonies were not mapped. Collections were delayed due to the coronavirus pandemic. This site was revisited August 19, 2020 to evaluate disease progression and sample corals. At this time, 6 newly diseased colonies were additionally sampled and three previously diseased colonies no longer showed signs of disease (apparently recovered). Thus, 22 diseased colonies and 3 recovered colonies were sampled in August. On October 14, 2020, three new diseased colonies and ten colonies completely covered in seemingly healthy tissue were added to the site and sampled for metabolomics.

Sampling and Extraction Procedure for Coral Samples
The samples were collected by divers on SCUBA with a 10 mL needleless syringe by agitating 4-5 adjacent polyps with the tip of the syringe. Dislodged tissue fragments and mucus were FIGURE 1 | Site map of sample collection site in Broward County, FL. (A) Shows the southern coast of FL with an arrow indicating the sampling location. (B) Shows the sample site (star) in relation to local benthic habitats and topography as described in Walker et al. (2008). (C) Describes the location of the sampled corals, as well as their health condition at the time of sampling. slowly drawn into the syringe until 10 mL of sample was collected. Diseased colonies were sampled along the disease lesion. Recovered colonies that displayed signs of SCTLD in March 2020, but no longer appeared to be diseased were sampled at the furthest point away from the former disease lesion to sample an area that was least likely to be impacted by SCTLD. For healthy corals, colonies that were completely covered in apparently healthy tissue were chosen to ensure that they were less likely to have had SCTLD or other tissue loss. All samples were then placed in 15 mL conical tubes, frozen and lyophilized.
After weighing, the lyophilized samples were transferred to a 20 mL scintillation vial for extraction. To each conical tube, 6 mL of 2:1 methanol (MeOH): H 2 O was added to dissolve any residual material left in the tube and transferred to the scintillation vial containing the lyophilized sample. To this vial, 4 mL of ethyl acetate (EtOAc) was added, and the samples were extracted (10 ml of 2:2:1 EtOAc: MeOH: H 2 O) for 4 h at 23 • C. The extracts were spun in a Savant SpeedVac for 10 min, and the supernatant was transferred to a 7 mL glass vial (precleaned three times with MeOH). The solvent was removed in vacuo using the Savant SpeedVac set at 35 • C and then lyophilized at −45 • C to remove any remaining water. The concentrate was transferred back into the initial falcon tube in 4 mL of 3:1 MeOH: H 2 O and centrifuged at 4,500 × g for 10 min to remove particulates. For metabolomics, 0.5 mL of the supernatant was transferred to a pre-weighed 1.5 mL microcentrifuge tube, centrifuged, and lyophilized. The remaining supernatant was transferred back to the 7 mL scintillation vial and dried in vacuo. The samples were stored at −20 • C until LC-MS data was acquired.

Sampling and Extraction Procedure for Cultured Endosymbionts
Symbiodiniaceae isolates, Breviolum, Symbiodinium, and Durusdinium, from the University of Buffalo Undersea Reef Research (BURR) collection were received from Richard Karp and Andrew Baker, University of Miami. These endosymbionts were isolated by Mary Alice Coffroth from Orbicella faveolata corals collected in the Florida Keys between 2002 and 2005. All isolates were cultured in F/2 media, incubated at 27 • C, with 20 µE of light on a 14:10 diurnal cycle. To prepare extracts, 500 µl of each culture was transferred to a 1.5 mL microcentrifuge tube. Cells were pelletized at 500 × g for 5 min to remove media and then washed twice with seawater filtered through a 0.22 µm -pore membrane filter (FSW). A sample of 0.2 µm-filtered seawater (500 µL) was extracted in tandem as a control. The washed pellet was resuspended in 100 µL of HPLC grade water for transfer into a pre-cleaned 7 mL scintillation vial and then frozen and lyophilized. Samples were extracted overnight at −20 • C in 4 mL of 2:2:1 EtOAc: MeOH:H 2 O. Solvents were removed on a SpeedVac at 35 • C for 4 h. The dried samples were resuspended in 3:1 MeOH:H 2 O and transferred into a 1.5 mL microcentrifuge tube. The samples were centrifuged at 4,500 × g for 10 min, and the supernatant was transferred to the pre-weighed microcentrifuge tube. This supernatant was removed via SpeedVac for 3 h after which the concentrated extract was frozen, then lyophilized. The weight of lyophilized extract was recorded.

Mass Spectrometry Data Acquisition and Analyses
All dried extracts were resuspended in 100% MeOH (LC-MS grade) containing 1 µM of sulfadimethoxine as an internal standard. The samples were analyzed with an Agilent 1290 Infinity II UHPLC system (Agilent Technologies) using a Kinetex 1.7 µm C18 reversed phase UHPLC column (50 × 2.1 mm) coupled to an Impact II ultrahigh resolution Qq-TOF mass spectrometer (Bruker Daltonics, GmbH, Bremen, Germany) equipped with an ESI source for MS/MS analysis. MS spectra were acquired in positive ionization mode, m/z 50-2,000 Da. An active exclusion of two spectra was used, implying that an MS 1 ion would not be selected for fragmentation after two consecutive MS 2 spectra had been recorded for it in a 0.5 min time window. The exclusion was reconsidered and additional MS 2 spectra was acquired if five-fold enhancement in intensity was observed. The eight most intense ions per MS 1 spectra were selected for further acquisition of MS 2 data. The chromatography solvent A: H 2 O + 0.1% v/v formic acid and solvent B: MeCN + 0.1% v/v formic acid were employed for separation. Flow rate was held constant at 0.5 mL/min throughout the run. The gradient applied for chromatographic separation was 5% solvent B and 95% solvent A for 3 min, a linear gradient of 5% B-95% B over 17 min, held at 95% B for 3 min, 95% B-5% B in 1 min, and held at 5% B for 1 min, 5-95%B in 1 min, held at 95% B for 2 min, 95-5%B in 1 min, then held at 5%B for 2.5 min at a flow rate of 0.5 mL/min throughout. A mixture of 6-compounds (amitriptyline, sulfamethazine, sulfamethizole, sulfachloropyridazine, sulfadimethoxine, coumarin-314) was run as quality control every 8 samples to ensure consistent instrument and column performance.
The raw data was converted to mzXML format using vendor software and preprocessed on MZmine 2.53 using mass detection, chromatogram building, chromatogram deconvolution, isotopic grouping, retention time alignment, duplicate removal, and missing peak filling (Pluskal et al., 2010). This processed data was submitted to the feature-based molecular networking workflow on the Global Natural Products Social Molecular Networking (GNPS) platform. Herein, the output of MZmine includes information about LC-MS features across all samples containing the m/z value of each feature, retention time of each feature, the area under the peak for the corresponding chromatogram of each feature, and a unique identifier. The MS 2 spectral summary contains a list of MS 2 spectra, with one representative MS 2 spectrum per feature. The mapping information between the feature quantification table and the MS 2 spectral summary was stored in the output using the unique feature identifier and scan number, respectively. This information was used to relate LC-MS feature information to the molecular network nodes. The quantification table and the linked MS 2 spectra were exported using the GNPS export module (Pluskal et al., 2010;Nothias et al., 2020) and the SIRIUS 4.0 export module (Pluskal et al., 2010;Dührkop et al., 2019). Feature-based molecular networking was performed using the MS 2 spectra (mgf file) and the quantification table (csv file). Three quantification tables were generated; one with LC-MS feature list of data acquired on extracts of all coral samples analyzed in the study, the second with LC-MS feature list of data acquired on the extracts of all coral samples and extracts of cultured Symbiodiniaceae, and the third with LC-MS feature list of data acquired on the extracts of coral samples collected in October and extracts of cultured Symbiodiniaceae.
Two molecular networks were generated using this workflow. To generate each network, the data was filtered by removing all MS/MS fragment ions within ± 17 Da of the precursor m/z. MS/MS spectra were filtered by choosing only the top 6 fragment ions in the ± 50 Da window throughout the spectrum. The precursor ion mass tolerance was set to 0.02 Da and the MS/MS fragment ion tolerance to 0.02 Da. A molecular network was then created where edges were filtered to have a cosine score above 0.7 and more than 4 matched peaks. Further, edges between two nodes were kept in the network if and only if each of the nodes appeared in each other's respective top 10 most similar nodes. Finally, the maximum size of a molecular family was set to 100, and the lowest scoring edges were removed from molecular families until the molecular family size was below this threshold. The spectra in the network were then searched against GNPS spectral libraries (Wang et al., 2016). The library spectra were filtered in the same manner as the input data. All matches kept between network spectra and library spectra were required to have a score above 0.7 and at least 4 matched peaks. The molecular networks were visualized using Cytoscape version 3.7.2 (Shannon et al., 2003). The compound annotations follow the "level 2" annotation standard based upon spectral similarity with public spectral libraries or spectra published in the literature as proposed by the Metabolomics Society Standard Initiative (Sumner et al., 2007). Annotations were confirmed with commercial standards when available. The MS2LDA analysis was performed as previously described with default parameters (van der Hooft et al., 2017;McAvoy et al., 2020). The output of this analysis is available at https://gnps.ucsd.edu/ProteoSAFe/status.jsp?task=a7dfdc1606c1 48c1bbda9c0fbc803d2e.
Prior to statistical analysis, blank subtraction was performed on the quantification tables using a Jupyter notebook, available at GitHub -Garg-Lab/Jupyter-Notebook-Blank-Subtraction-Broward-County-Corals. The mean of each feature in the samples is compared separately to the mean of the feature in the solvent controls and LC-MS blanks. A feature is retained if its sample mean is greater than 0.25 × the mean of the sample in the blanks. The entire quantification table is exported, with each feature marked as "true" (feature is retained) or "false" (feature is not retained). Visualization of the molecular network, metadata mapping, and feature filtering was performed using Cytoscape (v 3.7.2) (Shannon et al., 2003). Within Cytoscape, node filtering was performed using metadata and quantification tables and the nodes corresponding to features present in blank were removed. All Euler diagrams were created using EulerAPE (v 3.0.0) (Micallef and Rodgers, 2014). The filtered quantification table were submitted to MetaboAnalyst 5.0 for all multivariate analysis (Chong et al., 2019). Pareto scaling was used for normalization, then principal component analysis and hierarchical clustering using the Euclidean distance measure and Ward algorithm for clustering were performed. Heat maps of discriminating features ranked using non-parametric ANOVA were generated (Supplementary File 1). The software PRIMER version 7 was used to test similarities in metabolite species between sampled corals with differing health conditions. Cluster analyses with SIMPROF tests and non-metric multidimensional scaling (MDS) plots were constructed from a Bray-Curtis similarity matrix of log transformed metabolite data. The metabolite similarities between corals using the coral condition status at the time of sampling were tested using analyses of similarities (ANOISM). The similarity percentages (SIMPER) were calculated for each cluster to identify the metabolites most responsible for the differences between those clusters with 94% similarity.
Additional compound annotations were performed using the XCMS Online, and in silico tools such as MolDiscovery (v.1.0.0) and SIRIUS 4.5 with CSI:FingerID and CANOPUS, and through literature searches for compounds known to be produced by corals and their endosymbionts for untargeted analysis. MolDiscovery, available through the GNPS platform, compares in silico generated mass spectra of small molecules from a variety of databases with experimental MS 2 spectra and includes a similarity score for each reported match (Cao et al., 2020). The LC-MS feature list from the extracts of coral samples collected in October (mgf file) was submitted to the MolDiscovery workflow. The mzXML files of data acquired on extracts of corals collected in October were submitted to the XCMS Online platform and processed using the parameters for UPLC Bruker Q-TOF pos (6675) (Gowda et al., 2014).
CSI:FingerID was employed using default settings. The fragmentation trees were computed for submitted LC-MS features, and predicted chemical formulas and structures that best matched the experimental spectra were ranked using machine learning (Dührkop et al., 2019). The compound classes were assigned using CANOPUS, which uses MS/MS spectra and the chemical structures predicted by CSI:FingerID to propose chemical classes for a feature (Djoumbou et al., 2016;Dührkop et al., 2021).

Metabolomic Profiles of All Samples
The apparently healthy and diseased M. cavernosa samples were collected in two batches. The first batch was collected in August of 2020 (samples labeled as Diseased_Aug, 22 samples) and the second batch was collected in October of 2020 (samples labeled as Diseased_Oct, 3 samples and Healthy_Oct, 10 samples) (Figure 1 and Supplementary Figure 1). Additionally, three samples that were previously observed to have signs of SCTLD when the site was set up in March 2020, but no longer displayed signs of disease were also collected from the same site (samples labeled as Recovered_Aug). The extracted metabolites were stored frozen at −20 • C in a dried form and then the metabolomics data acquisition was performed in a single run to avoid batch-tobatch variation. The area under the curve for each detected MS 1 feature was obtained from extracted ion chromatograms using MZmine2. The resulting feature table was analyzed with a suite of data visualization and annotation tools, which are described below and detailed in the methods (Figure 2). After removing metabolite features detected in solvent controls, 3,914 metabolite features were detected in the coral extracts. The corresponding data on these features was first subjected to principal component analysis (PCA) to visualize metabolomic variation between samples. When PCA was performed on all samples, 44.3% of the total variation was explained by the PCA components 1 and 2 ( Figure 3A). This PCA indicated relatively separate clustering of the healthy (n = 10, green dots) from the diseased samples (n = 25, pink and orange dots), however, the recovered samples (n = 3, cyan dots) did not all cluster with healthy or diseased samples and showed large withingroup variation similar to the diseased samples. The first PCA component explained the variation within each sample class, whereas the second component explained separation between the healthy and diseased corals. A permutational multivariate analysis of variance (PERMANOVA) was performed with 999 permutations to investigate if the metabolome of diseased and healthy stony corals differed significantly. The PERMANOVA F-test with F (p) value of 4.834 (0.001) supported the observation that the metabolomes differ between healthy and diseased coral samples. Next, we employed partial least square discriminant analysis (PLS-DA), which as opposed to PCA, is a supervised method and maximizes the separation between pre-defined groups. The PLS-DA model revealed clear separation between healthy and diseased samples, while still highlighting withingroup variation with lower Q 2 value of 0.64 (Supplementary  Figure 2). Such within-group variation was previously shown for white syndrome disease in Acropora and Platygyra (Ochsenkühn et al., 2018). Thus, within-group variation is a hallmark of metabolomic profiles of coral holobionts collected in the wild and should be carefully considered in study design employing -omics methods. The metabolomics variation among healthy corals in a nursery has also been investigated for at least one species of coral and was attributed partly to the host genotype, while other factors were proposed for with-in group variation (Lohr et al., 2019). The metabolome profile variation in visually healthy corals in the ocean may be attributed to multiple sources, including host genotype, associated endosymbionts (Symbiodiniaceae), the microbiomes, the surrounding abiotic environment, past bleaching history, and whether these corals are under stress before visual signs of disease appear. The variation in the diseased samples may additionally arise from severity of disease, duration of infection, or variation in groups of colonizing microorganisms.

Metabolomic Profiles of Samples Collected in October
Assessing the metabolomes of apparently healthy and diseased corals collected in October alone revealed tighter clustering of the diseased samples with relatively separate clustering from healthy samples suggesting that the disease progression may lead to convergence of metabolomes (pink dots, Figure 3B). The clustering of the diseased samples collected in October was tighter than the August diseased samples (Figure 3A). In this regard, the three principal components for October samples explained 77% of the total variation with component 1 accounting for the majority of the variation (56.3%) as compared to 23.3% of the variation explained by PC1 when all samples  were considered (Figure 3 and Supplementary Figure 3A). Furthermore, the healthy samples 3-6 displayed greater spread across all components suggesting large within-group variation in these apparently healthy samples. Among these, the healthy 3 showed the largest variation along all three components in PCA and displayed gray pigmentation of the entire colony (Supplementary Figure 1). Thus, the phenotypic variation between the coral colonies of the same species in the field is captured at the metabolome level suggesting metabolomics approaches could be productive in the future for linking specific phenotypes to chemotypes as databases for compounds detected in coral holobionts are populated. The M. cavernosa corals have been suggested to maintain high levels of genetic diversity compared to other coral species (Budd et al., 2012). The relative intraspecies variation in metabolomic profiles of different species of corals in the future will enable insights into whether genetic diversity can be a driver of observed metabolomic variation.

Geographic Proximity and Metabolomic Variation
To visualize the relation of spatial proximity of coral colonies to observed clustering in PCA, we categorized the corals by the 94% similarity MDS clusters in a map and circled colonies with close relationships in the hierarchical cluster analysis (HCA) (Figure 4), wherein samples are clustered in a dendrogram in which the linkages and their relative proximity equates to the similarity of metabolomic profile. The closer the branches, the more similar are the given object's variables, herein metabolite abundances. Two major clusters were observed in HCA (branches 1 and 2, Figure 4A and Supplementary Figure 3B). Each cluster was further subdivided in two subclusters ( Figure 4A, branches 3, 4 and branches 5, 6). In each of these sub-clusters, all apparently healthy samples clustered together whereas diseased samples separated into additional subclusters. The Recovered_Aug samples did not follow a defined pattern; Recovered_1 and Recovered_3 samples clustered with healthy samples whereas the Recovered_2 sample clustered with the diseased samples. To query where this clustering is dictated by the proximity of coral colonies, the individual branches in HCA were colored based on their location shown in the site map ( Figure 4A). Indeed, a number of neighboring branches revealed geographic proximity, however, this observation was not consistent across all colonies within the HCA groups. The yellow branches in the HCA and corresponding yellow circles in the map highlight connected colony pairs that were close neighbors. The purple and orange highlighted branches and circles were colonies in the same HCA sub-cluster, but were not directly connected. Some colonies in these branches were also in close spatial proximity. Based on the two clusters observed in HCA (Figure 4A branches 4 and 6), the healthy samples were colored in green and blue in the PCA plot in Figure 3B. The clustering patterns were reproducible for diseased corals when analyses was performed on data generated on diseased samples only (Supplementary Figure 3C).
The non-metric multi-dimensional scaling (MDS) plots divided the samples into eight distinct clusters ( Figure 4B). Categorizing the corals by the 94% similarity clusters in the map revealed no obvious large-scale spatial clustering. Many colonies with similar metabolite communities were spread throughout the site and spatially comingled with dissimilar colonies. However, some neighboring colony pairs were close in similarity indicating potential within-MDS-cluster spatial patterns. For example, MDS cluster A indicated colonies 1 and 2 were most similar to each other and furthest away from the other colonies in the similarity group, which was also the case in the map. So was the case with cluster E colonies 7 and 8. However, there were many more cases where colonies with dissimilar metabolite communities were in close spatial proximity. Analyses of Similarities (ANOSIM) were used to test for similarities in metabolites between corals using the coral condition status at the time of sampling. Since three distinct clusters were evident at 94.6% similarity in the Diseased_Aug samples, similarity percentages (SIMPER) were calculated for each cluster to identify the metabolites most responsible for the differences between those clusters. A number of features describing the differences between the healthy and diseased samples include metabolites produced by Symbiodiniaceae such as glycolipids, betaine lipids, and tocopherols are described below.
The variation in coral metabolomes is a hallmark of complex interactions between the holobiont members that differ in genotypes, chemotypes, as well the varying levels of disease progression observed in field samples. The observations from a smaller number of field samples in October in our study require validation in future studies that employ time series design post the onset of SCTLD. Metabolomics paired with description of host and endosymbiont genotype as well as microbiome structure and bleaching history is essential to explain underlying factors resulting in observed clustering. When coupled with metaproteomics and metatranscriptomics, the underlying mechanisms of such variation can be further explained at the biochemical level.

Description of Metabolite Features Underlying Variations in Metabolomic Profiling
We used feature-based molecular networking, and various other in silico approaches listed in Figure 2 and described in the methods sections, to annotate metabolite features showing within group and between group variations. The distribution of metabolite features across healthy, diseased, and recovered samples was visualized as UpSet plot (Figure 5). Of 3,914 metabolite features, 2,231 features were observed in all coral samples. Among the metabolite features unique to each group, the largest number of unique features were detected in diseased corals collected in August (366 features) followed by apparently healthy samples (234 features), recovered samples (29 features) and the least number of unique features were detected in diseased coral samples collected in October (11 features). An important component to understanding the response of corals to disease is the ability to assign the source of the metabolites; that is whether the metabolites originate from the host, endosymbionts, the associated microbiome, or the environment. We hypothesized that the large proportion of the metabolome unique to healthy samples may be attributed to the endosymbiotic dinoflagellates, Symbiodiniaceae, as these symbionts display signs of stress during disease and some diseased coral colonies, including all three collected in October, did display bleaching along the tissue margin of the lesion while others did not. The M. cavernosa has been reported to display bleaching along the margin of the lesion in response to SCTLD . To examine metabolite sources, we cultured three genera of Symbiodiniaceae, extracted the metabolites from the cultures, and acquired untargeted metabolomics data as described for coral samples. The cultured Symbiodiniaceae include Symbiodinium sp., Durusdinium sp., and Breviolum sp. The MS/MS spectra from the coral tissue extracts and the cultured Symbiodiniaceae were preprocessed together using MZmine. Of the 3,992 features retained after blank subtraction, 20.3% of these features were shared between the coral holobiont and the cultured symbionts, with 74.3% of the features unique to the coral holobiont, and 5.4% of features unique to cultured Symbiodiniaceae ( Figure 6A). The metabolites produced in culture by three species of cultured Symbiodiniaceae were also analyzed and significant overlap in their metabolomes was observed, as well as unique metabolites detected in culture extracts of each species (Figure 6B). To visualize features driving the differences between various sample types in PCA, a heat map was generated on the top-50 metabolite features ranked using non-parametric ANOVA (Figure 7). The MS 2 spectra for these top 50 features that showed differential abundance between healthy, recovered, and diseased samples from the output of ANOVA were searched using GNPS, NIST 2017, MolDiscovery, XCMS Online, and SIRIUS with CSI:FingerID. Where annotations were not possible, the chemical class was assigned using CANOPUS. These database matches, proposed compound class by CANOPUS, and putative molecular formulas assigned by SIRIUS with CSI:FingerID are reported in Supplementary Table 1. Many of the features driving the separation among individuals within the healthy coral category as well as between the healthy and the diseased coral samples were detected in both the cultured Symbiodiniaceae and field corals (Figure 7, indicated by±). Even though these Symbiodiniaceae were not isolated from the M. cavernosa coral colonies analyzed in the present study (see section "Materials and Methods") and the cultured Symbiodiniaceae can differ phenotypically from when they are present in the coral host (Maruyama and Weis, 2021), the significant overlap observed between field corals and cultures suggests that a culturebased approach is a viable step to tease apart at least a portion of Symbiodiniaceae metabolomic contributions to the coral holobiont. Indeed by manipulating culturing conditions, various different metabolites were shown to be produced by Symbiodiniaceae (Nakamura et al., 1998). Similar approaches have been used to identify microbiome-specific metabolites in infectious diseases in humans (Quinn et al., 2016a;Garg et al., 2017;Melnik et al., 2019). In addition, spatial patterning of microbes and metabolites has also been used to assign sources of detected molecules (Little et al., 2021). The variability in the Symbiodiniaceae metabolites within healthy samples was apparent from the heat map and may arise from different host genotypes or differences in associated microorganisms. It is now well known that some coral colonies can display resilience to disease and the role of Symbiodiniaceae in this resilience has been proposed (Correa et al., 2009;Littman et al., 2010;Rouzé et al., 2016;Rosset et al., 2019). Thus, in the absence of mass spectral databases for coral symbionts and microbiomes, a culture-based approach can be employed to further assess their role in metabolome variation in field corals. An alternative approach to culturing allowing direct comparison is to acquire metabolomics data on isolated in hospite pure Symbiodiniaceae fractions from the coral tissue (Sikorskaya et al., 2021).
To gain chemical insights into the overlapping features detected in extracts of coral and cultured Symbiodiniaceae, we attempted annotation of the metabolite features scored in the MDS analysis in Figure 4B, PLS_DA in Supplementary  Figure 2 and ANOVA selected features in Figure 7 using FBMN and in silico approaches. The features in this category belonged to lipids, including betaine lipids such as diacylglycerylcarboxyhydroxymethylcholine (e.g., DGCC, m/z 490.374), monogalactosyl-and digalactosyldiacylglycerols (MGDG and DGDG), sphingolipids, lyso-PC, steroid-like molecules and tocomonoenol (Supplementary Table 1). Among these, betaine lipids are characterized by the presence of trimethylated hydroxyamino acid connected to the diacylglycerol moiety through an ether bond. To generate a comprehensive list of betaine lipids detected, feature-based molecular networking ( Figure 8A) was paired with discovery of mass structural motif specific to betaine lipids using MS2LDA analysis ( Figure 8B). The annotation of mass structural motif in Figure 8B allowed annotation of additional betaine lipids not present in the cluster in Figure 8A (and Supplementary Table 2). A comprehensive network of betaine lipids including the ones produced in the cultured Symbiodiniaceae is shown in Supplementary Figure 4. Of note, some ether linked higher chain length acyl tails were observed only in culture extracts of Symbiodiniaceae observed via presence of nodes with higher m/z (Supplementary Figure 4). The distribution of a subset of annotated betaine lipids is shown as heat maps in Figure 8C and distribution of all annotated betaine lipids is shown in Supplementary Figure 5. Notably, the healthy sample 3 displayed the largest variation in metabolomics data from PCA plot and also displayed the largest diversity of betaine lipids (Supplementary Figure 5). The betaine lipids have been previously reported in endosymbionts from various coral species (Rosset et al., 2019;Sikorskaya, 2020;Sikorskaya et al., 2021) and are reported as biomarkers of bleaching history (Roach et al., 2021) and thermal tolerance (Leblond et al., 2015;Rosset et al., 2019). It has been previously suggested that betaine lipids being the component of cell membrane are influenced by the host genotype (Sikorskaya et al., 2021). The temporal changes in profiles of these lipids as well as differences between different species undergoing SCTLD in the future will provide important insights into their functions and potential role in resilience. Glycolipids are present in thylakoid membranes of Symbiodiniaceae and their importance in heat tolerance has been previously reported (Tchernov et al., 2004;Leblond et al., 2015;Rosset et al., 2019). The glycolipid cluster in FBMN ( Figure 9A) was annotated using MolDiscovery, which compares user uploaded MS/MS spectra with in silico generated MS/MS spectra of small molecules, and the annotations were propagated using FBMN and MS2LDA and compared  to mass spectrometry data in the literature. The feature with m/z 797.519 was annotated as MGDG, potentially heterosigma-glycolipid II, a diacylglycerol glycolipid that contains omega-3 polyunsaturated fatty acid residues (Kobayashi et al., 1992). The detection of specific fragments with m/z 333.243 and 259.206 in the MS 2 spectra supports the presence of 18 carbons with 4 desaturations as one of the acyl tails ( Figure 9B). Therefore, feature m/z 797.519 was annotated as a monogalactosyldiacylglycerol, MGDG (20:5/18:4). Additional MGDG glycolipids were annotated using propagation of FBMN (Supplementary Table 3). The fragments with m/z 333.243 and 259.206 resulted in a MS2LDA substructure motif specific to galactosyldiacylglycerols containing an 18:4 acyl tail. Using this chemical information, the features m/z 771.504 and 769.488 were annotated as MGDG (18:4/18:4) and MGDG (18:5/18:4), respectively. These Symbiodiniaceae produced glycolipids displayed different patterns than Symbiodiniaceae produced betaine lipids in healthy as well as diseased samples ( Figure 9C and Supplementary Figure 6). Indeed, glycolipids and not betaine lipids were shown to be variable when Cladocopium C3 was exposed to thermal stress (Rosset et al., 2019). These differences suggest that in addition to the variability in the abundance of the Symbiodiniaceae, their chemotypes are also variable across diseased samples. In recent work by Landsberg and colleagues, histopathological examination of several species of corals with tissue loss attributed to SCTLD revealed disruption in physiology of host and Symbiodiniaceae (Landsberg et al., 2020). The differences in Symbiodiniaceae-specific metabolites further support this finding. In algae and dinoflagellates, MGDGs are one of the main classes of lipid that comprise thylakoid membranes (Li-Beisson et al., 2019;Sikorskaya et al., 2021), on which light-dependent photosynthesis occurs, and their functions as antitumor, anti-inflammatory, anticancer, and appetite-suppressing agents have been described from plant sources. Their variability among corals undergoing different diseases should be compared and correlated to Symbiodiniaceae genotypes to link their functions to holobiont health. Since changes in glycolipids in thylakoid membranes have been linked to thermotolerance in corals and in other microalgae (Huang et al., 2017) and increased frequency of mass bleaching events are occurring, deeper understanding of these lipid distributions in different species of corals either resilient or susceptible to disease as well as in culture extracts of isolated symbionts is necessary. Previously, GC-MS based metabolomics did not reveal correlation between symbiont composition and the host metabolite profiles (Matthews et al., 2020), but LC-MS based metabolomics of the entire holobiont in our case, clearly reveals metabolite differences that are attributed to the symbionts via a culture-based strategy.
In addition to various lipids such as LysoPC, sphinogolipids, and steroid-like molecules detected at higher levels in healthy corals (Supplementary Table 1), an α-tocopherol analog, tocomonoenol also detected in the Symbiodiniaceae cultures and was among the top ranked features (Figure 10). Detection of tocomonoenol in culture extracts suggests that the origin of this molecule is endosymbiotic Symbiodiniaceae, again supporting that culture-based strategy could facilitate description of individual contributions of each partner of the holobiont. Indeed, the chloroplast-associated antioxidant alpha-tocopherol was observed to be accumulated in intracellular metabolite pools of symbionts during thermal bleaching (Hillyer et al., 2017) and the antioxidant activity of tocomonoenol and other tocopherols have been previously explored (Yamamoto et al., 1999;Beppu et al., 2020). The annotation of the tocomonoenol was enabled by the in silico tool MolDiscovery and supported by the MS 2 spectra for α-tocopherol available in the NIST database ( Figure 10B). Thus, in the absence of mass spectral databases, such in silico tools serve as a valuable strategy enabling a glimpse into structural scaffolds detected in the coral holobiont, and culture-based strategies allow assignment of the biological partner contributing the metabolite (Figure 10). FBMN also allowed visualization of various analogs of α-tocopherol ( Figure 10A). The annotation of α-tocopherol was confirmed using commercial standard (Supplementary Figure 7). The mirror plots in Supplementary Figure 8A support that metabolites with m/z 424.332 (orange node, Figure 10A) and m/z 410.317 (orange and cyan-colored node, Figure 10A) are tocotrienols. While tocopherol with a fully saturated polyprenyl side chain and tocomonoenol with a single unsaturation of the side chain were observed at higher abundances in healthy samples, tocotrienol with three unsaturations of the side chain (m/z 410.318) accumulated in diseased corals (Supplementary Figure 8B).
Lastly, the metabolites observed at higher abundances in the diseased samples were queried based on a heat map generated on features ranked via non-parametric ANOVA for log transformed data (Supplementary Figure 9). Several metabolites with similar MS 2 spectra in this class belonged to the chemical class of unsaturated fatty acids such as docosahexaenoic acid (Supplementary Figure 10 and Supplementary Table 4). Accumulation of unsaturated fatty acids upon stress induced by exposure to high concentrations of pollutant, ethylhexyl salicylate, in the coral Pocillopora damicornis has been reported and their potential as biomarker for stress in corals has been suggested (Stien et al., 2020). These observations warrant targeted analysis of polyunsaturated fatty acids and their downstream effectors such as eicosanoids in corals undergoing disease.

CONCLUSION
In this study, we employed metabolomic profiling to visualize intraspecies variation of M. cavernosa corals with no visual sign of disease, corals that apparently recovered, and corals presumably with SCTLD. The unsupervised PCA, HCA, and MDS revealed metabolomic variation within healthy corals, diseased corals, and between healthy and diseased corals, which was attributed to various metabolites including lipids and tocopherols produced by Symbiodiniaceae. In many instances, the healthy samples that clustered closer in PCA and HCA plots were also immediate neighbors in the site map. Thus, the proximity resulted in similar metabolomes in most cases with a few exceptions. This observation was not true for most diseased colonies, which showed wider spread in metabolomic profiles in the PCA plot. In the future, time series experimental design for disease progression, along with characterization of host genotype, symbiont genotype, and microbiome structure is important to explain the source of this variability in the metabolomic profiles.
The emerging approaches in metabolomics technology allow analyses of hundreds to thousands of metabolites. The major bottleneck in data analyses in metabolomics stems from our inability to annotate the metabolites detected and to assign the biological source of the metabolite to an individual holobiont partner. In the absence of comprehensive mass spectral databases for coral metabolites, annotations were made using in silico methods and literatures searches. These compound annotations were propagated to additional unknown metabolite features through a feature-based molecular networking approach and by discovery of sub-structural motifs using MS2LDA. Thus, we demonstrate that use of in silico tools serve as a practical approach in the absence of databases to aid and improve metabolite identifications for profiling of coral holobionts. Better understanding of holobiont response to disease can be gleaned if metabolomics variation can be further assigned to the respective partner in the holobiont. In this work, we highlight a feature-based molecular networking approach to metabolomics data acquired on cultured Symbiodiniaceae and the holobiont to assign sources of detected metabolites for the first time. Using this strategy, we annotated various top features describing metabolomics variation to Symbiodiniaceae demonstrating the usefulness of this approach. A caveat of this culture-based approach is that not all metabolites relevant in the holobiont may be produced in culture. A tedious approach to overcome this limitation is to culture the symbionts under different growth conditions such as temperature variation, nutrient limitation, etc., to induce production of metabolites otherwise observed in a holobiont. In addition, direct extraction of metabolites from the relatively pure fraction of in-hospite Symbiodiniaceae isolated from the coral should be applied for direct comparison. Indeed, this approach has been used previously to describe lipid composition of different species of Symbiodiniaceae under temperature stress. Once the source of a metabolite is assigned to the holobiont partner, genomics can be employed to decipher the biochemical pathways and their regulatory mechanisms. As additional metabolomic datasets become available on cultured coral holobiont partners including cultured bacteria, such analyses will continue to provide biochemical insights into the workings of a coral holobiont. Thus, linking metabolites to their biological producer through comparative metabolomics in a coral holobiont is an important next step in understanding symbiosis in corals and its disruption during disease and environmental stress.

DATA AVAILABILITY STATEMENT
The LC-MS/MS datasets generated in this study are available in the public repository, Mass spectrometry Interactive Virtual Environment (MassIVE) at https://massive.ucsd.edu/ with the identifier MSV000087471.

AUTHOR CONTRIBUTIONS
JD conducted metabolomics data acquisition, analyses, data interpretation, and assisted with manuscript editing. OJ assisted with metabolomics data analyses. KP and JH conducted the experiments. NG helped with project conception, supervised the project, designed the experimental approach, and assisted with various methodological, data and statistical analyses, and wrote the manuscript. VP helped design the sampling scheme, supervised sample handling and extractions, and contributed to interpretation of data and manuscript editing. BU helped with project conception and editing of manuscript. BW helped with project conception, experimental design, data analyses, and editing of manuscript. All authors contributed to the article and approved the submitted version.

FUNDING
This study was supported by the Environmental Protection Agency to NG, VP, and BU (award number X7-01D00120 -0), the Florida Department of Environmental Protection Office of Resilience and Coastal Protection-Southeast Region to VP (award number B7C0F5), BU (award number B7F150) and BW (award number B558F2), and NSF CAREER award to NG (award number 2047235).