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.