Gene Expression Response to Stony Coral Tissue Loss Disease Transmission in M. cavernosa and O. faveolata From Florida

Since 2014, corals within Florida’s Coral Reef have been dying at an unprecedented rate due to stony coral tissue loss disease (SCTLD). Here we describe the transcriptomic outcomes of three different SCTLD transmission experiments performed at the Smithsonian Marine Station and Mote Marine Laboratory between 2019 and 2020 on the corals Orbicella faveolata and Montastraea cavernosa. Overall, diseased O. faveolata had 2194 differentially expressed genes (DEGs) compared with healthy colonies, whereas diseased M. cavernosa had 582 DEGs compared with healthy colonies. Many significant DEGs were implicated in immunity, extracellular matrix rearrangement, and apoptosis. These included, but not limited to, peroxidases, collagens, Bax-like, fibrinogen-like, protein tyrosine kinase, and transforming growth factor beta. A gene module was identified that was significantly correlated to disease transmission. This module possessed many apoptosis and immune genes with high module membership indicating that a complex apoptosis and immune response is occurring in corals during SCTLD transmission. Overall, we found that O. faveolata and M. cavernosa exhibit an immune, apoptosis, and tissue rearrangement response to SCTLD. We propose that future studies should focus on examining early time points of infection, before the presence of lesions, to understand the activating mechanisms involved in SCTLD.


INTRODUCTION
Florida's Coral Reef (FCR) is a fragile and highly endangered reef system that has been impacted by numerous coral disease outbreaks (Holden, 1996;Richardson et al., 1998;Porter et al., 2001). Since the 1970s, coral disease has reduced coral cover from 30% to 2% (Richardson et al., 1998;Aronson and Precht, 2001;Patterson et al., 2002;Williams, 2005). The decline of coral cover of FCR has been further exacerbated by anthropogenic stressors such as eutrophication (Vega Thurber et al., 2014;Lapointe et al., 2019) frequent hurricanes (Gardner et al., 2005) and intensified bleaching events due to global climate change (Kuffner et al., 2015;Manzello, 2015). In addition to these compounding threats, boulder, and brain coral populations on FCR face significant additional mortality to a newly emerging disease.
At a histological level, SCTLD has been characterized as having multifocal lytic necrosis that starts in the gastrodermis and extends out to the surface epithelia (Landsberg et al., 2020). Disintegration and fragmentation of the tissue, in particular the Symbiodiniaceae-containing cells, as well as the Symbiodiniaceae themselves, have been documented (Landsberg et al., 2020). While it is understood that massive tissue destruction occurs, the mechanisms involved in this response from the coral host remain unknown.
One central question to the cause and spread of SCTLD is whether the coral immune system is affected and thus implicated in the lytic necrosis previously observed. The coral immune system is known to be diverse, and reactive to disease and other environmental perturbations (Palmer andTraylor-Knowles, 2012, 2018). Corals possess a diverse array of pathogen recognition receptors, signaling pathways, and effector responses that have been shown to respond to many different coral diseases (Pinzón et al., 2015;Fuess et al., 2016Fuess et al., , 2018Young et al., 2020). However, the response of the coral immune system to SCTLD is still not known. To address this, we performed several laboratory SCTLD transmission experiments on Orbicella faveolata and Montastraea cavernosa at Mote Marine Laboratory and the Smithsonian Marine Station from 2019 to 2020. We quantified the gene expression using transcriptomic analysis and identified important signatures of immunity and apoptosis. Transcriptomics is invaluable as a diagnostic when the etiology of a disease is not well understood (Fuess et al., 2018;Zhou et al., 2019;Navarro et al., 2020;Young et al., 2020). We found unique differential expression patterns between the two species but found that there were also shared patterns of apoptosis, extracellular matrix rearrangement, and immune response. We hypothesize that apoptosis and immune response are important mechanisms of the host response to SCTLD and propose that future experiments should focus on early time points before the lesions are present, as well as on the interactions of the Symbiodiniaceae within the coral host.

Coral Collection and Transmission Assays
Samples analyzed in this study were collected from three different experiments. Detailed information on the collections and transmission assays are provided below. For all experiments below, "apparently healthy corals" are defined as showing no visual signs of stress or tissue loss. The term "diseased coral" refers to the corals that possessed an active lesion(s) and were used to transmit the disease to experimental corals. "Experimental corals" are defined as apparently healthy ones that are then exposed to a diseased coral. "Control corals" are apparently healthy coral that we exposed to another apparently healthy coral. "Transmission control" refers to controls that acted as a control for coral colony contact, as introducing a different coral species near another can lead to aggression. The transmission controls were not sequenced.

Mote Marine Laboratory
The first experiment was conducted from April to June 2019 using apparently healthy and diseased (subacute tissue loss) M. cavernosa (Eaton et al., in prep). For this experiment seven apparently healthy M. cavernosa colonies were collected from the Airport Coral Heads reef in Key West (24.53919 • , −81.77270 • ) which, at the time, was ahead of the known SCTLD invasion zone. Four diseased coral colonies were collected at Looe Key, in the Florida Keys (24.54767 • , −81.45697 • ). An example of the experimental setup is shown in Figure 1.
The second experiment was conducted from July to September 2019 using seven apparently healthy O. faveolata and four diseased M. cavernosa, also at Mote. At the time of collection, only M. cavernosa colonies were found to have active disease and was therefore used as the disease coral in the transmission experiment. Apparently healthy and diseased colonies were collected in the same locations as colonies used for the first experiment.
After collections were completed, all apparently healthy and diseased colonies were immediately transported to the Coral Health and Disease lab at Mote Marine Laboratory in Sarasota, FL for experimentation. Upon arrival, the diseased and apparently healthy colonies were placed in separate bins with ambient seawater overnight, and the transmission experiment was set up the next day. Each bin was equipped with multiple powerheads to maintain circulation, and temperature was regulated by a recirculating temperature-controlled water bath. One 3 × 1.2 m raceway holding twenty 18.9-L aquaria was used for each experiment. Using sterilized gloves, each diseased M. cavernosa colony (n = 4 × 2 experiments) was fragmented on a clean and sterile table. Coral colonies were initially fragmented using an angle grinder (Dewalt DWE402) and trimmed with a diamond blade band saw (Gryphon Corporation, C-40 CR Aquasaw XL).
For both experiments, one of the apparently healthy colonies were randomly chosen to act as a transmission control. The remaining six apparently healthy colonies of each species were fragmented into approximately ten pieces each (∼10 cm long × 3 cm wide × 2 cm high in size) and served as the experimental or control fragments. Apparently healthy colonies were fragmented first to minimize contamination, and all blades and tools were sanitized with a 10% bleach solution before fragmentation. For each study, four of the 20 aquaria were set up as controls with only experimental corals in contact with transmission controls to account for any aggression from direct contact. The other 16 aquaria had experimental coral fragments in contact with diseased coral fragments ( Figure 1A).
Once disease transmission was visible on the experimental samples, they were collected by scraping the surface of the lesioned area with a sterile razor blade. Each experimental sample was paired with a control sample of the same genotype and collected at the same time. The experimental samples had lost approximately 25-95% tissue at the time of collection (Supplementary Table 1). Samples were flash frozen, stored at -80 • C, and sent to the University of Miami, Rosenstiel School of Marine and Atmospheric Sciences on dry ice for 3' RNA-seq library preparation and sequencing.

Smithsonian Marine Station
The third experiment was performed at the Smithsonian Marine Station (SMS) in Ft. Pierce, FL from February to March 2020. For this experiment, five apparently healthy M. cavernosa colonies were collected outside of Dry Tortugas National Park ahead of the SCTLD invasion zone and were brought back to holding tanks at SMS where they were maintained in UVC-sterilized, filtered (0.22 µm) seawater (aquaria system previously described in Ushijima et al., 2020). The apparently healthy colonies were fragmented using a masonry saw (Husqvarna MS 360) fitted with a 35.56 cm diamond, continuous rim circular saw blade, and left to recover for 17 days before the start of the transmission experiment.
Diseased corals were collected off patch reefs near Sand Key in the Florida Keys on February 20, 2020. The corals were transported to SMS on the same day and maintained in filtered seawater overnight, until being fragmented with the masonry saw and used for transmission experiments. Each experimental coral fragment was put into contact with either a diseased fragment (Diploria labyrinthiformis or Pseudodiploria strigosa) or control fragment ( Figure 1B).
Once disease transmission was visible on the experimental fragments it was collected, along with its corresponding control of the same genotype and preserved in RNAlater (Invitrogen) following manufacturer specifications. After incubating for 24 h at room temperature, the samples were stored at -80 • C for later processing. Samples were transported on dry ice to the University of Miami, Rosenstiel School of Marine and Atmospheric Sciences for 3' RNA-seq library preparation and sequencing. Disease transmission ranged from approximately 10-75% tissue loss at the time of collection (Figure 1 and Supplementary Table 1).

RNA Extraction and TagSeq Library Preparation
Total RNA was extracted using the Qiagen RNeasy Minikit including the recommended 15-min DNase digestion. A Qubit fluorometer v3.0, with the HS-RNA assay kit, was used to assess the quantity of total RNA following the manufacturer's protocol. Total RNA was then converted into 3' complementary DNA (cDNA) libraries using the Lexogen Quantseq FWD Library Preparation kit following manufacturer's protocol. For all samples, library quality and quantity were assessed using an Agilent D1000ScreenTape analysis. A total of 76 samples passed quality control and were sequenced. Samples were sequenced for 100 base pair single-read sequencing using the NOVAseq at the University of Miami, Miller School of Medicine, John P. Hussman Institute for Human Genomics.

Sequence Analysis
All sequence analysis scripts with parameters used are available at https://github.com/cnidimmunitylab/sctld_jamboree. For both species, raw read quality was assessed using FASTQC (Brown et al., 2017). Low quality reads, up to Q10 using the Phred algorithm and adapter sequences, were then trimmed using bbduk.sh in the BBtools suite 1 with the recommended parameters from Lexogen 2 . Trimmed reads from each species were then aligned to their respective publicly available genomes; O. faveolata (NCBI accession GCA_002042975.1) and M. cavernosa (Matz, 2018) using STAR version 2.7.7 (Dobin et al., 2013) and parameters recommended by Lexogen (see text footnote 2). Sequences that did not align to reference genomes were excluded from analysis. Aligned sequence reads were quantified at the gene level using featureCounts in the Subread package (Liao et al., 2014) before being imported into R (v3.6.1) and RStudio (v1.2.1335). For O. faveolata there were many isoforms present for some of the genes. Due to all transcripts for a specific gene having identical annotations the first isoform annotation was selected for downstream analysis. There were no isoforms present in the M. cavernosa genome. Gene ontology (GO) and KEGG pathway annotations were identified for each gene's protein product using the eggNOG-mapper online tool 3 (Huerta-Cepas et al., 2017). These annotations were used in enrichment tests for identified gene lists from the analysis.

Differential Gene Expression Analysis
Pre-filtering of the genes was performed for both species to increase computational speed and reduce memory requirements. For O. faveolata, genes with counts less than or equal to five in four or fewer samples were removed. For M. cavernosa, genes with counts less than or equal to three in two or fewer samples were removed. This differential filtering is due to the varying sample number across species; 26 samples for O. faveolata and 17 samples for M. cavernosa. For initial sample visualization, read counts were transformed using the variance stabilized transformation (VST) from the DESeq2 R package (Love et al., 2014). Transformed data was used in a principal component analysis (PCA) and visualized using ggplot2 (Wickham, 2016). For PCAs, the experimental batch effect was removed using the R package limma (Ritchie et al., 2015). To identify differentially expressed genes between control and diseased fragments, a Wald test with the formula "∼ treatment (control vs. diseased) + experiment (Mote experiment 1 or 2, Smithsonian experiment)" was run in DESeq2 using the unnormalized counts with the default median of ratios normalization applied (Love et al., 2014). For both species, genes with a False Discovery Rate (FDR) adjusted p-value (p adj ) of less than 0.05 were considered significantly differentially expressed.

Orthofinder-Shared Orthologous DEGs Between O. faveolata and M. cavernosa
Orthofinder, using the default settings, was used to identify orthologous genes that were differentially expressed in both O. faveolata and M. cavernosa. To perform this analysis, the predicted protein sequences were collected from the reference genomes of O. faveolata and M. cavernosa and analyzed to identify "orthogroups" (Emms and Kelly, 2019). These orthogroups were then compared between O. faveolata and M. cavernosa and the resulting orthography Venn diagram was made using the R package Venn Diagram.

Weighted Gene Co-expression Network Analysis
For O. faveolata, weighted gene co-expression network analysis (WGCNA) was performed to identify modules of co-expressed transcripts that correlated to that disease treatment (Langfelder and Horvath, 2008). WGCNA was not conducted for M. cavernosa due to the limited sample number. Co-expression analysis identifies genes with similar expression patterns and groups them into modules. Once modules are generated, "eigengenes" that represent the expression patterns of the module are calculated and correlated with additional traits (i.e., control vs. disease). A single signed network was constructed by following the step-by-step network construction method using VST-normalized counts with the bi-midweight correlation and a soft-thresholding power of 20 (Langfelder and Horvath, 2008). Experimental batch effect was also removed using the R package limma prior to network construction (Ritchie et al., 2015). Gene co-expression modules were identified using a cut height of 0.99 on the topological overlap matrix and a minimum module size of 30 genes. Modules with > 85% similar expression profiles were merged. The expression of each network module "eigengene" was correlated against a binary matrix containing categorical treatment information to identify modules with expression patterns linked to a given treatment. Genes with the highest module membership, defined as having highest eigengene correlation values, were identified as the top hub genes for each module.

Gene Ontology Enrichment Analysis
Gene ontology (GO) enrichment analysis was run for both species using the significantly differentially expressed genes (FDR p adj < 0.05) to identify overrepresented GO terms using the plugin BiNGO (Maere et al., 2005) in Cytoscape 4 (version 3.8.2). 2,005 genes were used for O. faveolata, and 582 genes were used for M. cavernosa. The discrepancies in these numbers in comparison to the total number of DEGs is due to gene annotation limitations and duplicate protein identifiers assigned to certain genes. Additionally, GO enrichment analysis was run on two modules of interest (METan and METurquoise) derived from WGCNA for O. faveolata. These modules were chosen because they possessed significant correlations for module membership and gene significance. The reference set of genes, or gene universe, included all 18,066 O. faveolata genes from our data set. The hypergeometric statistical test was used with the Benjamini and Hochberg FDR multiple testing correction. GO terms were deemed significantly enriched if they had an FDR p adj of < 0.05. Biological process, cellular component, and molecular function GO terms were tested for statistical enrichment within DEG lists and WGCNA module lists.

KEGG Enrichment Analysis
For Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis, we looked for enrichment of KEGG terms within lists of differentially expressed genes (2005 for O. faveolata and 582 of M. cavernosa due to limited gene annotations and duplicate protein identifiers assigned to certain genes). Protein identifiers that mapped to unique KEGG functional ontology terms (KO) were used as the reference list (17,287 KO terms for O. faveolata and 14,606 for M. cavernosa) using a hypergeometric test with the "enricher" function in the clusterProfiler package in R Studio (Yu et al., 2012). Significantly enriched KEGG terms were identified as those with FDR p adj < 0.05. These significant KO terms were mapped to KEGG pathways using "Search Pathway" in the KEGG Mapper tool using the best available references, the human genome and coral Acropora digitifera, because KEGG mapper does not have ontologies available for either M. cavernosa or O. faveolata (Kanehisa and Sato, 2020).

M. cavernosa and O. faveolata Sequencing Outcome
Seventeen out of the forty-two M. cavernosa samples (11 experimental and 6 controls) had a minimum read depth greater than two million reads per sample and were then used for downstream analyses. These samples had an average of 7,636,834 reads per sample, and 45.51% uniquely aligned to the M. cavernosa genome (Matz, 2018). Twenty-six out of twentynine O. faveolata samples (13 experimental and 13 controls) sequenced samples had a minimum read depth greater than two million reads per sample and were then used for downstream analyses. These samples had an average read depth of 10,531,729 reads per sample and 42.41% uniquely aligned to the O. faveolata genome (Prada et al., 2016). After filtering out genes with low counts, as detailed in section "materials and methods", the count data included 18,066 genes for O. faveolata and 13,816 genes for M. cavernosa.

Principal Component Analysis and Differential Gene Expression Analysis in Both O. faveolata and M. cavernosa
For all analyses, the batch effect of experiment (Mote experiment 1, Mote experiment 2, or Smithsonian experiment) was removed prior to PCA visualization to distinguish expression patterns due to treatment. Principal component analysis of O. faveolata revealed that experimental samples were different in their gene expression from control samples. The first (PC1) and second (PC2) principal components captured 37% and 13% of variance, respectively, with experimental and control samples clustering apart roughly along PC1 (Figure 2A). Like O. faveolata, principal component analysis of M. cavernosa revealed that experimental and control samples were unique in their expression profiles. 19% of the variance was explained by PC1 and 16% of the variance was explained by PC2, with experimental and control samples clustering apart along PC1 (Figure 2B) Figure 3A). In total, 582 genes were differentially expressed in M. cavernosa, of which 351 DEGs were upregulated and 231 were downregulated ( Figure 3B). The 30 most significantly differentially expressed genes ranked by most significant FDR p adj for O. faveolata included carbonic anhydrase 2-like, tetratricopeptide repeat protein 4-like, collagen alpha-3(VI) chain like, and rhodopsin GQ-coupled like among others ( Figure 4A and Supplementary Table 2). The 30 most significantly differentially expressed genes ranked by most significant FDR p adj for M. cavernosa included collagen alpha chain, galaxin, and DNA-binding protein RFX6 among others ( Figure 4B and Supplementary Table 2).

A Varied Response for Extracellular Matrix Genes in Both Coral Species
In both O. faveolata and M. cavernosa, SCTLD transmission induced the transcription of many genes involved in extracellular matrix building and tissue integrity, including the transcription factor CP-2/Grainyhead, collagens, laminins, integrins, fibronectins, hemicentins, and proteoglycans. In O. faveolata, a total of twenty-seven collagen related genes were identified that had significant L2FC (FDR p adj < 0.05). Nine of these genes had significant upregulation in response to SCTLD and eighteen had significant downregulation in response to SCTLD (Supplementary Table 2). O. faveolata also had eleven hemicentins that were significantly differentially expressed, with seven upregulated and four down regulated (Supplementary Table 2). Additionally, Gene ID: LOC110058607: fibronectin type III domain-containing protein-like (L2FC = 2.692, FDR p adj = 0.00000085); Gene ID: LOC110050720: laminin subunit gamma-1-like isoform X1(L2FC = 2.952, FDR p adj = 0.00004216); Gene ID: LOC110062631, proteoglycan 4-like isoform X1 (L2FC = -1.76, FDR p adj = 0.0003602); Gene ID: LOC110040404: fibronectin type 3 and ankyrin repeat domains protein 1-like isoform X1 (L2FC = 2.164, FDR p adj = 0.0029294); and Gene ID: LOC110050028 integrin alpha-7-like (L2FC = 1.184, FDR p adj = 0.0100922) were all differentially expressed in response  Two distinct clusters separate control fragments (left cluster) from experimental fragments (right cluster), regardless of experiment. The top 30 differentially expressed genes had FDR p adj < 6.5 × 10e -5 . (A) Heatmap of log2 counts per million for O. faveolata top 30 differentially expressed genes. Gene counts for 26 samples (13 control, 13 experimental) are shown. Red represents experimental samples, blue indicates control. Gray represents experiment 1 and black represents experiment 2, of the two separate transmission experiments were run at Mote Marine Laboratory. Three distinct clusters separate gene expression patterns of all control fragments (center cluster), experiment 1 disease exposed fragments (far right cluster), and experiment 2 disease exposed fragments (far left cluster). Top 30 differentially expressed genes had FDR p < 6.5 × 10e -11 . (B) Heatmap of log2 counts per million for M. cavernosa top 30 differentially expressed genes. Gene counts for 17 samples (6 control, 11 experimental) are shown. Red represents experimental samples, blue indicates control. Gray represents experiment run at Mote Marine Laboratory; black represents transmission experiment run at the Smithsonian Marine Station. Table 2). In M. cavernosa, a total of eleven collagens had significant differential expression in response to SCTLD, five upregulated and six downregulated (Supplementary Table 2). M. cavernosa also had four significantly differentially expressed hemicentins, with three upregulated and one downregulated. Additionally, Gene ID: Mcavernosa21960: fibronectin type 3 and ankyrin repeat domains protein 1 (L2FC = 4.076, FDR p adj = 0.03647272); and Gene ID: Mcavernosa03959: Transcription factor CP2 (L2FC = -3.462, FDR p adj = 0.01687218) were all significantly differentially expressed in response to SCTLD transmission (Supplementary Table 2).

Orthology Analysis Reveals Shared Gene Response Between O. faveolata and M. cavernosa
In total seventy Orthofinder assigned orthogroups were found to be differentially expressed in both coral species. Fifty-nine out of the seventy, had congruent differential expression responses, where forty were upregulated and nineteen were downregulated (Figure 5 Table 3). Orthogroups that were up-or

KEGG Enrichment Analysis Found Enrichment for Terms Associated With Hormone Synthesis and Cancer
Of the 2005 DEGs analyzed for O. faveolata, three KO terms, K21626 (FDR p adj = 0.0052), K03068 (FDR p adj = 0.0138), and K07997 (FDR p adj = 0.028), were all significantly enriched. KEGG pathways, p-values, and the number of genes associated with each KO term can be found in Table 1. In general, these terms were associated with KEGG pathways involved in hormone synthesis and cancer. These included pathways for multiple cancers, neurodegenerative diseases, Wnt signaling pathway, and the mTOR signaling pathway. There were no significantly enriched KEGG terms for M. cavernosa.

Overall GO Enrichment Analysis
Fourteen "Biological Process" and five "Cellular Component" gene ontology terms were significantly enriched within the list of 2005 DEG for O. faveolata. Significant "Biological Process" terms included "the regulation of metabolic processes, " "regulation of transcription, " "tissue development, " and "regulation of cell fate specification" ( Table 2). Significant "Cellular Component" terms included: "cellular projection, " "cilium, " "anchoring junction, " and "adherens junction" ( Table 3). Significant GO terms for Biological Process' contained an average of 199.5 protein identifiers per GO term, while significant "Cellular Component" terms contained an average of 73.2 protein identifiers per GO term (Figure 6). Because the relative frequency of a given GO term within our DEG list compared to the entire transcriptional dataset determines significance, there is no relationship between number of protein identifiers and GO term significance. There were no significantly enriched GO terms for M. cavernosa.   Table 5).

GO Enrichment Analysis on METan and METurquoise Modules
GO enrichment analysis was also applied to two significantly correlated modules identified from WGCNA related to disease exposure in the experimental treatment (Supplementary Table 6). METurquiose has 291 significantly enrich terms for "Biological Process, " 57 significantly enriched terms for "Molecular Function, " and 83 significantly enriched terms for "Cellular Components" (Supplementary Table 4). For "Biological Process" the term with the most significant FDR p adj value, was "regulation of cellular component organization" which had 126 gene IDs associated with it (GO ID: 51128, FDR p adj = 1.00E-08). For "Molecular Function" the most significant term was "binding" which had 278 gene IDs associated with it (GO ID: 5488, FDR p adj = 4.93E-09). Lastly, for "Cellular Component" there were two terms that were the most significant: non-membrane-bounded organelle (GO ID: 43228, FDR p adj = 2.74E-14) and intracellular non-membrane-bounded organelle (GO ID: 43232, FDR p adj = 2.74E-14). Both had 175 genes IDs associated with them. METan had no significant GO terms enriched. In this study, we found that there were transcriptional changes in response to SCTLD transmission regardless of the experiment or species. Overall O. faveolata had a much stronger gene expression response than M. cavernosa. While this could be biologically significant, more likely it is due to differences in sequencing success. O. faveolata had an overall higher success with sequencing despite both species having the same preservation methods and similar RNA input quality. We identified many immune, apoptosis, and extracellular matrix genes that were significantly differentially expressed in response to SCTLD and found that that there shared orthologs between O. faveolata and M. cavernosa with a similar gene expression directionality. The exact reason for the quick and deadly spread of SCTLD is unknown. However, based on predictions using hydrodynamic modeling a waterborne pathogen outbreak is suspected (Dobbelaere et al., 2020). An immunocompromised state, defined as a weakened or reduced immune response, has been hypothesized to play a role in the spread of coral diseases (Lesser et al., 2007;Work et al., 2012). Within corals, a weakened immune response would be indicated by attenuated differences in gene expression of immune related genes or proteins, or lower constitutive expression of immune genes or proteins (Palmer et al., 2010;Palmer, 2018). The constitutive immune response of these corals has not been measured, and thus could be reduced indicating immunocompromised state, however, we did find that both O. faveolata and M. cavernosa had significantly differentially expressed immune genes after SCTLD transmission (Supplementary Table 2). Previously, Caribbean corals of modern lineages, such as O. faveolata and M. cavernosa, were identified as being more susceptible to disease in comparison with older lineages due to differences in constitutive expression of antimicrobial activity, melanin, and antioxidants (Pinzón et al., 2014). However, in response to SCTLD, both O. faveolata and M. cavernosa have been found to be less susceptible than many other species . This contradiction may indicate that the types of immune responses that are deployed by different Caribbean corals are not all equally effective against novel diseases and may vary depending on the disease type. Further investigation into the effects of SCTLD on the constitutive immune response will help to clarify our understanding of this difference. This lineagespecific difference in immune response along with variation in species diversity among coral communities may help explain differences in spatial distribution of SCTLD Sharp et al., 2020).

Shared Immune Activity in Both O. faveolata and M. cavernosa
In this study, we found many shared immune-related genes that were differentially expressed in both O. faveolata and M. cavernosa. In both species, peroxidases (ex. peroxidasins and peroxiredoxins) were identified as significant DEGs (LOC110065690 and Mcavernosa31629, Supplementary Table 2), as shared orthogroups (OG0014907: Peroxiredoxin activity, OG0002872: peroxidase activity) (Supplementary Table 3), and as having module membership in METan, the module most significantly correlated to disease transmission (Supplementary Table 5). Overall, peroxidases are important enzymes involved in innate immune defense. Within scleractinians and gorgonians peroxidases have previously been identified as strong immune responders to heat stress and disease (Mydlarz and Harvell, 2007;Palmer, 2018). Peroxidasin is linked to both extracellular matrix function and innate immunity; and peroxiredoxins are a family of peroxidases that are involved in scavenging reactive oxygen species (ROS) such as peroxides (Nelson et al., 1994;Péterfi et al., 2009;Abbas et al., 2019). Previous studies on the sea anemone model, Aiptasia, found that peroxidasins expression is affected by symbiotic state and heat stress (Lehnert et al., 2014;Cziesielski et al., 2018). The exact role, however, could not be elucidated because there was a variation in response, and indications that there may be important genotypic differences in expression (Lehnert et al., 2014;Cziesielski et al., 2018). In response to a pathogen exposure, including bacteria, viruses and fungi, the transcription of peroxiredoxins are activated and a host immune response is initiated (Abbas et al., 2019). Within corals, ROS production has been implicated in coral bleaching and disease (Venn et al., 2008;Weis, 2008;Hansel and Diaz, 2021). Here the presence of significantly upregulated peroxiredoxins in both coral species, faveolata. Significant enrichment was tested for a list of 2005 uniquely annotated differentially expressed genes between control and disease exposed O. faveolata fragments. Color indicates FDR p adj from hypergeometric overrepresentation test (yellow = most significant). Size of dot indicates number of uniquely annotated genes that were identified and mapped to gene ontology terms by protein identifier, in each significant GO category. may indicated that an antioxidant response is occurring to either an environmental response, immune response, or a combination of the two. This shared response should be further investigated to understand the core mechanisms and the source of ROS production. Overall, given the multiple areas where shared responses were identified, the role of peroxidases in SCTLD response should be further investigated to understand this potential protective and/or destructive mechanisms.
Other annotated orthogroups linked to immunity included: OG0001112: transmembrane receptor protein tyrosine kinase activity, and OG0003082: Transforming growth factor-beta (TGF-beta) family (Supplementary Table 5). In both species the TGF-beta family genes associated with this orthogroup were downregulated (Supplementary Table 5). Previously, a study in O. faveolata found that exposure to exogenous TGF-beta had an immune suppressive effect, while inhibiting TGF-beta maintained the baseline immune response (Fuess et al., 2020). The downregulation of TGF-beta in both corals we hypothesize indicates that an immune response to SCTLD transmission is occurring. Further investigation into the modulation of this pathway in response to SCTLD will be informative due to its potential to be immune suppressive. Protein tyrosine kinase (PTKs) act as promoters of proinflammatory cytokine production (Nag and Chaudhary, 2009). Within O. faveolata, PTK has previously been identified as significantly expressed in response to both exposure to lipopolysaccharides and Vibrio alginolyticus (Fuess et al., 2016). We hypothesize presence of this significantly upregulated orthogroup in our dataset indicates that an immune response is occurring through PTK signaling. Further research on the role of PTK signaling, including the use of proteomics and enzymatic assays during SCTLD infection will help in our understanding of the mechanisms of this signaling pathway.
Unique Immune Responses of O. faveolata and M. cavernosa to SCTLD Within O. faveolata, we found that the most highly expressed immune gene was fibrinogen-like protein 1. This gene was also the second most highly expressed DEG overall for O. faveolata (Supplementary Table 2). Vertebrate fibrinogen related proteins are important modulators of coagulation, and in invertebrates are hypothesized to play a role in pathogen defense, as a pathogen recognition receptor (PRR) (Zhang et al., 2010;Hanington and Zhang, 2011;Adema and Loker, 2015). In a distantly related cnidarian, Hydractinia symbiolongicarpus, fibrinogens are differentially expressed in response to both Gram-positive and Gram-negative bacteria (Zárate-Potes et al., 2019). In corals, the function of this gene and related gene-family members is still not understood, however its highly significant gene expression in response to SCTLD may indicate that it is important PRR for detecting SCTLD-associated pathogens. Within M. cavernosa, the mostly highly expressed gene overall was an immune related gene WD repeat-containing protein 36 (Supplementary Table 2). WD repeat-containing proteins are a large family of beta sheet containing proteins that have a varied function including immunity (Li and Roberts, 2001). In humans, these proteins are associated with many human diseases and can play a role in both immune defense and apoptosis (Li and Roberts, 2001). Again, the functional mechanisms of this gene within coral disease response are not understood but given its high differential expression within M. cavernosa it could be an important immune mechanism for response to SCTLD.

Apoptosis Signals in Response to SCTLD Transmission
Apoptosis is the cellular process of programmed cell death through activation of caspases and proteases (D'Arcy, 2019). It is often initiated in response to stressors such as heat and hypoxia. Here, we found that both O. faveolata and M. cavernosa had significantly differential expression of apoptosis-related genes in response to SCTLD transmission, both up-and downregulated, as well as high module membership within METan, a module significantly correlated with disease transmission (Supplementary Table 5). The hub gene, with the highest module membership in METan did not have an annotation, but notably the apoptosis regulator, Bax-like was found to have high membership in this module (Supplementary Table 5). Bax with another protein called Bak, are pro-apoptotic proteins that, when activated, form holes in the membranes of mitochondria causing a proapoptotic cascade (Westphal et al., 2014). The high module membership of this gene within this module, along with other pro-apoptotic genes present within this module may imply that a coordinated apoptosis response is occurring in response to SCTLD. Previous studies on coral immunity have found that apoptosis is a critical response to disease, and propose that it may be a critical mechanism for removing foreign pathogens for the coral host and is activated in more disease susceptible species (Fuess et al., 2017;Roesel and Vollmer, 2019;Zhang and Lieberman, 2020). Previous histopathological analysis of corals with SCTLD found that lytic necrosis (LN), which starts from the gastrodermis, was the primary attribute of this disease (Landsberg et al., 2020). We believe that the discrepancy between the previous histopathology observations of LN and our current observations of apoptosis may be linked to timing of sampling. If a cell doesn't have enough ATP for completing apoptosis, then secondary necrosis can occur through the swelling and lysis of the cells (D'Arcy, 2019). We hypothesize that this differential gene expression and complex varied apoptotic response is due to the late term effects of SCTLD infection, which could lead to LN, much like what was observed previously in histological studies.

Extracellular Matrix Genes as Primary Responders to Disease
Tissue rearrangement as indicated through the differential expression of extracellular matrix genes such as collagens and other wound healing factors has previously been documented as an important component of the heat stress and disease response in corals (Traylor-Knowles, 2019; Young et al., 2020). For example, in a previous study on the effects of disease on the gene expression of Acropora palmata found that fragments with no disease transmission had an enrichment of genes including "collagen alpha chains, " and "protocadherin-like, " indicating that extracellular matrix genes may play an important role in preventing disease progression (Young et al., 2020). In the current study, we observed that the critical extracellular matrix genes: collagen and hemicentin were differentially expressed in response to SCTLD in both species (Supplementary Table 2). It should also be noted that two orthogroups associated with collagens were also identified as shared between O. faveolata and M. cavernosa (OG0010810, OG0010929, Supplementary  Table 3). We hypothesize that the differential expression of collagens and other wound healing related genes is being activated in response to the tissue degradation occurring during apoptosis and subsequent LN. In the future, investigation into the role of tissue rearrangement and integrity as an important mechanism for protection against SCTLD may be useful. If SCTLD is a waterborne transmissible disease, as has been hypothesized, then it could presumably be transmitted if the epithelial of the coral is compromised. In the future, testing if resistant genotypes have an enrichment of cell adhesion and wound healing factors that act as protection against pathogen invasion will be informative for our understanding of how this disease is transmitted (Young et al., 2020).

SUMMARY
In this study we analyzed O. faveolata and M. cavernosa samples that were exposed to SCTLD through three laboratory transmission assays. We found that treatment was a greater driver of variance than time or location of the experiment. In both species, we identified many different DEGs and orthogroups in response to SCTLD. We observed that genes from apoptosis, immunity, and extracellular matrix gene pathways were differentially expressed, and found that apoptosis genes had high module membership within METan; the module most correlated with disease transmission. Additionally, we propose that immune reactivity related to peroxidases, PTK, and TGF-Beta pathways be further investigated in the context of SCTLD. Together these data suggest that SCTLD is causing both an immune and apoptotic response, as well as tissue rearrangement in corals. In the future, studies should focus on testing earlier time points of infection before the physical manifestation of the lesion, as this may give more information into the activating mechanisms involved in SCTLD.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found on the National Center for Biotechnology Information, NCBI Project Number: BioProject ID: PRJNA737001. All codes used for the analysis of this work can be found on GitHub: https://github. com/cnidimmunitylab/sctld_transcriptomics_2021.

AUTHOR CONTRIBUTIONS
NT-K, KE, EM, KR, VP, and BU designed the experiments and facilitated the experiment. KR, BU, and KE conducted the experiments. NT-K, MC, BY, AD, MD, AG, NK, GS, and CM processed, analyzed, and interpreted the data. All authors participated in the writing and editing of the manuscript.

FUNDING
The collection, processing and analysis of this work was funded by the Florida Department of Environmental Protection Office of Resilience and Coastal Protection and Coral Reef Conservation Program (DEP CRCP) and by the National Science Foundation Grant number: 1951826.