Methods
Transmission experiment
The ex-situ transmission experiment was conducted in the Experimental Reef Laboratory at the University of Miami’s Cooperative Institute for Marine and Atmospheric Studies over a period of two weeks in March 2020. The disease transmission apparatus consisted of 80, 0.5 L coral vessels with independent flow-through water sources maintained at local ambient temperature (24 ℃) and light conditions (PAR of 250 µ mol m-2 s-1; Studivan et al., 2022). Twenty fragments each of disease-naive O. faveolata and M. cavernosa were sourced from the land-based nursery at Mote Marine Laboratory. Each fragment was split into two equal subfragments for healthy/disease exposure genotype pairs (N = 80 total) and allowed to recover for eight days prior to the disease challenge. A disease donor colony of M. cavernosa with visible SCTLD lesions was collected from Broward County, Florida (26.1479, -80.0939) two days prior to the start of the experiment, and was cut into ~1 x 4 cm fragments using a diamond band saw with each fragment containing a portion of active lesion. Coral fragments in the disease-exposed group were maintained in direct contact with diseased tissue fragments, where disease donor fragments were replaced as needed following total loss of tissue. Genotype pairs in the healthy group were not exposed to any other corals. Corals were observed daily over the course of the experiment to quantify the number of days to onset of disease lesions, as well as to observe the gross signs of lesion formation and tissue necrosis. Following progression of lesions across approximately 50% of the fragment area, coral fragments and their corresponding healthy genotype pair were sampled via a small tissue scrape from near the lesion, which was preserved in TRIzol and flash-frozen in liquid nitrogen. Corals remaining at the end of the experiment were processed in the same manner. Herein, coral fragments are characterized as “healthy” for all samples in the control group, “diseased” for samples in the disease-exposed group exhibiting lesions, and “NAI” (no active infection) for samples in the disease-exposed group that showed no visible signs of disease lesions by the end of the experiment.
Intervention experiment
The in-situ intervention experiment was conducted at long-term monitoring site BC1 (mean depth 5.8 m) in Broward County, Florida (26.1479, -80.0939; Combs et al., 2021; Shilling et al., 2021) over thirteen days in April–May 2020. Fifteen colonies of M. cavernosa with visible SCTLD lesions were tagged and small tissue scrapes were collected 5 cm from the disease lesion and preserved as described previously. Nineteen visually healthy colonies of the same species were also tagged and sampled. Following the initial sampling, diseased colonies were treated with amoxicillin in CoreRx Base 2B as described in Shilling et al. (2021), with the exception that trenching of the coral skeleton in advance of the disease lesion was not conducted. Colonies were then revisited thirteen days later, and tissue sampling was repeated.
Tag-Seq library preparation, sequencing, and bioinformatics
Total RNA was extracted from samples from both experiments using modified versions of protocols for the Zymo Direct-zol RNA Miniprep and Zymo RNA Clean & Concentrator-5 kits (Studivan, 2022d). RNA was quantified on an Agilent Bioanalyzer 2100 and analyzed for purity on a Nanodrop 1000, where samples with RIN scores <4, or poor purity ratios, were re-extracted. All successfully-extracted samples were normalized to 60 ng/µL and shipped to the University of Texas at Austin Genome Sequencing and Analysis Facility for library preparation and sequencing. Tag-Seq library preparation was performed in duplicate using 600 ng of initial RNA for M. cavernosa samples and 100 ng for O. faveolata samples due to initial issues with inhibitory compounds. Following successful amplification of all samples, libraries were pooled and sequenced on a NovaSeq 6000 S2-XP flow cell generating 1 x 100 bp reads with 20% PhiX spike-in, for a target read depth of ~10 million reads per library (~20 million reads per sample).
Raw sequences were processed according to custom Perl scripts in a GitHub repository release (Studivan, 2022d) to deduplicate, trim (FASTX-Toolkit; Gordon & Hannon, 2010), align (Bowtie 2; Langmead & Salzberg, 2012), and quantify gene counts for each sample. Dominant symbiont genera were determined by aligning sample reads to a reference containing Symbiodiniaceae 28S sequences across the four main genera (Symbiodinium spp., Breviolum spp., Cladocopiumspp., and Durusdinium spp.). Orbicella faveolata samples contained a majority of Durusdinium spp. alignments (93.4%), andM. cavernosa samples were dominated by Cladocopium spp. (99.7%); hence, the respective species datasets were aligned to the dominant symbiont references. Cleaned reads were mapped to separate host (M. cavernosa, Kitchen et al., 2015; O. faveolata, Pinzón et al., 2015) and algal symbiont (Cladocopiumspp., Davies et al., 2018; Durusdinium spp., Shoguchi et al., 2021; respectively) reference transcriptomes sequentially to remove cross-contaminated sequences in both host and symbiont transcriptomes. Transcriptome annotations for host and symbiont references were created according to GitHub repositories for O. faveolata(Studivan, 2022b) andM. cavernosa(Studivan, 2022a) using eggNOG-Mapper (Huerta-Cepas et al., 2016, 2018).
Following deduplication, the mean number of reads across all samples was 24.6 ± 0.9 million, with 5.8 ± 0.2 million reads remaining following trimming and cleaning. Alignment rates for O. faveolata samples were 60.3 ± 1.6% for host and 15.6 ± 0.8% for symbionts, resulting in 81,960 and 39,854 isogroups, respectively. For M. cavernosasamples, alignment was 61.2 ± 0.7% for host and 6.5 ± 0.3% for symbionts, with 82,594 and 32,862 isogroups, respectively. The original version of the M. cavernosa transcriptome was chosen over the updated, genome-derived version (https://matzlab.weebly.com/data–code.html) due to low alignment in the latter (28.0 ± 0.3% for host and 7.1 ± 0.4% for symbiont) and poorer representation of isogroups (24,850 and 32,932 isogroups, respectively) for downstream differential expression analysis.
Statistical analysis
All statistical analyses were conducted in the R statistical environment (R Core Team, 2019). Analysis scripts and outputs can be found in a GitHub repository release associated with this manuscript (Studivan, 2022c). Transmission metrics from the transmission experiment (time to onset of disease lesions) were analyzed for an effect of species using an ANOVA of Box-Cox transformed data, and with a fit proportional hazards regression model in the packages survival(Therneau, 2021) andsurvminer(Kassambara et al., 2021) to determine the relative risk of developing disease lesions between species.
Gene counts were first separated by experiment, species, and host/symbiont datasets, then pre-filtered to remove isogroups with a cumulative count <10 across all samples. Data were subsequently transformed using a variance stabilizing transformation in the DESeq2 package (Love et al., 2014), and examined for outlier samples using the packagearrayQualityMetrics(Kauffmann et al., 2009). Samples violating the distance between sample arrays criterion (S a) were removed from the respective experiment and host/symbiont datasets (Table 1; Table 2). Significance testing for each experiment dataset was performed with PERMANOVAs using 1e 6 permutations, and sample dispersion was visualized using PCoAs of Manhattan distance in the package vegan(Oksanen et al., 2015). Differential expression analyses were conducted usingDESeq2 models for the respective experiments. For O. faveolata in the transmission experiment, the model genotype + exposure group (control, sctld) was used, whereas for M. cavernosa , the model genotype + disease status (healthy, diseased, NAI) was used since not all fragments exhibited disease signs. For the intervention experiment, a single-factor concatenated treatment (control, SCTLD) and time (0,1) model was used. Contrast statements to conduct pairwise comparisons for each of the experiment datasets as described in the GitHub repository release (Studivan, 2022c). Visualization of differentially expressed genes (DEGs) was done using the pheatmap package (Kolde, 2019).
Functional enrichment analyses with Mann-Whitney U tests were conducted for each experiment’s pairwise comparisons using the packagesGO-MWU (Wright et al., 2015) with an alpha cutoff of 0.05 and KOGMWU(G. B. Dixon et al., 2015) for full gene ontology (GO) and eukaryotic orthologous group (KOG) annotations, respectively. To identify correlational relationships between gene modules (genes with similar expression patterns) and experimental variables, weighted gene co-expression network analyses (WGCNA; Langfelder & Horvath, 2008) were conducted as described in the associated GitHub repository release (Studivan, 2022c). Module-associated genes were then exported for GO enrichment analysis to examine any relationships between gene expression patterns and additional sample traits beyond diseased versus healthy status.
Finally, DEGs unique to and shared between species and experiments were determined using the VennDiagram package (Chen & Boutros, 2011). For the transmission experiment, common DEGs (orthogroups) and associated GO/KOG annotations for diseased/healthy comparisons inO. faveolata and M. cavernosa were identified by blast annotations using OrthoFinder (Emms & Kelly, 2015, 2019), and visualized by heatmaps and KOG classes. For comparison of DEGs in M. cavernosa samples from the transmission and intervention experiments (e.g., to compare DEGs in disease-exposed transmission samples to antibiotic-treated samples in the intervention experiment), matching gene and KOG annotations were examined in the same fashion as the cross-species comparisons.