Sample collection

Coral colonies showing active Stony Coral Tissue Loss Disease (SCTLD) as well as nearby completely healthy colonies were targeted for sampling on February 11 and 13, 2020 on Buck Island (“Outbreak”, 18.27883°, -64.89833°), and Black Point (“Existing”, 18.3445°, -64.98595°) reefs, respectively, in St. Thomas, USVI (Fig. 1). Buck Island was considered a recent outbreak site, and will be referred to as “Outbreak.” Coral species sampled included Montastraea cavernosa, Colpophyllia natans, Meandrina meandrites, and Orbicella franksi (Table 1). Black Point had been experiencing SCTLD for approximately one year, and will be referred to as “Existing.” Coral species sampled at the Existing site were Montastraea cavernosa and Colpophyllia natans (Table 1). SCTLD was identified by single or multiple lesions of bleached or necrotic tissue with epiphytic algae colonizing the recently dead and exposed skeleton (Fig. 2). At both reefs, some paling of colonies was apparent, especially on Orbicella spp., as a result of a recent bleaching event in October 2019. Due to this, on Orbicella spp. it was challenging to distinguish SCTLD from White Plague-type diseases, which generally occur following bleaching events (Miller et al., 2009). As a result, we avoided sampling Orbicella spp., except when it was clear the colony had regained full coloration and the disease lesion was consistent with SCTLD infection.

To investigate if putative pathogens were recoverable from seawater surrounding diseased colonies, near-coral seawater was sampled from 2-5 cm away from each coral colony prior to tissue sampling. This occurred via negative pressure with a 60 ml Luer-lock syringe (BD, Franklin Lakes, NJ, USA). Two seawater samples were collected over each colony displaying SCTLD lesions: one sample was taken directly above healthy tissue approximately 10 cm away from the lesion, when possible, and a second sample over diseased tissue. Syringes were placed in a dive collection bag for the duration of the dive, and then the seawater was filtered through a 0.22 mm filter (25 mm, Supor, Pall, Port Washington, NY, USA) while onboard the boat and the filter with holder was placed in a Whirl-pak bag and kept on ice until returning to the shore. While onshore, filters were placed in sterile 2 ml cryovials (Simport, Beloeil, QC, Canada) and frozen in a liquid nitrogen dry shipper.

After near-coral seawater sampling, a single tissue sample was taken from the healthy colonies and two samples were taken from each diseased colony. Diseased tissue samples were collected at the interface between healthy and newly bleached tissue (Fig. 2), and then visually healthy tissue was targeted approximately 10 cm away from the disease interface (when possible, sometimes diseased colonies had very little healthy tissue remaining). For some colonies, limited healthy tissue remained and tissue approximately 3 to 5 cm away was targeted. Coral tissue and mucus “slurries” were collected by using 10 ml non-Luer lock syringes (BD). The syringe tip was used to agitate and remove tissue from the colony, while suspended tissue, mucus, and seawater were simultaneously aspirated into the syringe. To control for the significant amount of seawater and seawater-associated microbiota introduced into the tissue sample, a total of nine 10 ml syringes were used to capture only ambient reef seawater greater than 1 m off the reef benthos. This would allow us to assess how much seawater-based bacteria were in the slurry samples, and are hereafter referred to as “Syringe Method Control” samples. Immediately after collection, the syringes were placed in a Whirl-pak bag underwater to prevent the loss of tissue or mucus. Once back onboard the boat, samples were transferred to 15 ml sterile conical tubes and placed in a 4°C cooler. Upon returning to the lab, samples were frozen to -20°C until analysis.

To capture the surrounding seawater physical and chemical environment an Exo2 multiparameter sonde (YSI, Yellow Springs, OH, USA) was used with probes for temperature, salinity, dissolved oxygen, pH, and turbidity (Table S1). On the day before sampling (February 10, 2020), each sonde probe was calibrated following the manufacturer’s protocols, and was not calibrated between sampling dates.

DNA extraction, PCR, and sequencing

DNA extraction, PCR, and sequencing preparation protocols were specifically designed for the Illumina iSeq 100 System (Illumina Inc., San Diego, CA, USA), a portable, high-quality sequencing technology. In an approximately 1 cu. ft. size, the Illumina iSeq 100 System produces 4 million paired-end sequence reads of high quality (<1% error rate), making it an attractive technology to adapt for field-based microbiome studies. Sequencing runs took approximately 17 hours, and generated about 2 GB of data, which could easily be offloaded and processed on a standard laptop without the need for Wi-Fi. In the USVI, we brought the iSeq 100 System to a home rental, which was transformed into a remote laboratory where we successfully conducted all DNA extractions, PCR and subsequent sequencing.

DNA extraction on both seawater and tissue slurry samples and associated extraction controls was conducted using the DNeasy PowerBiofilm Kit (Qiagen, Germantown, MD, USA). Modifications at the beginning of the extraction were applied based on the sample type. For seawater filters, the 0.22 mm filter was placed directly into the bead tube, and then manufacturer instructions were followed. Three seawater microbiome DNA extraction controls were included (named D4-D6). For those, an unused 0.22 mm filter was used in the DNA extraction method. For tissue slurry samples and syringe method control samples, the slurries were thawed at room temperature, then immediately transferred to 4°C prior to extraction. Slurries were vortexed for 10 seconds, then 1.8 ml of each slurry was transferred to an emptied bead tube. Samples were vortexed at 12,045 rcf (maximum rcf available on centrifuge) for 10 min to concentrate tissue, mucus, and the associated microorganisms at the bottom of the tube, and supernatant was removed. For samples that were very clear (very little tissue collected via syringe), a second aliquot of 1.8 ml of tissue slurry was centrifuged on top of the existing pellet in order to capture more tissue- and mucus-associated microorganisms. Beads were then returned to the tube now containing the pellet and the extraction continued following manufacturer protocols. Three DNA extraction controls were included for the slurry extractions (named D1-D3), and included DNA extractions proceeding using an empty bead tube.

PCR was used to amplify the V4 region of the small sub-unit ribosomal RNA (SSU rRNA) gene of bacteria and archaea. This amplification occurred in two stages. In stage one, 2 ml of tissue slurry template DNA was added to a 50 ml PCR reaction. One ml template and a 25 ml total reaction volume was used for seawater samples. For negative PCR controls, 1 or 2 ml of sterile PCR-grade water was used. One Human Microbiome Project mock community, Genomic DNA from Microbial Mock Community B (even, low concentration), v5.1L, for 16S rRNA Gene Sequencing, HM-782D was included as a sequencing control, and 1 μl was used in a PCR reaction with seawater samples. Fifty ml PCR reactions contained 0.5 ml polymerase (GoTaq, Promega, Madison, WI, USA), 1 ml each of 10 mM forward and reverse primers, 1 ml of 10 mM dNTPs (Promega), 5 ml MgCl2 (GoTaq), 10 ml 5X colorless buffer (GoTaq), and 29.5 ml UV-sterilized and PCR-grade water. PCR reactions with 25 ml proceeded with the exact proportions of reagents as 50 ml reactions. Earth microbiome project primers, 515F and 806R, targeted bacteria and archaea and were used with Illumina-specific adapters (Apprill et al., 2015; Parada et al., 2016). The PCR reaction continued using two small, portable PCR thermocyclers: mini8 (miniPCR, Cambridge, MA, USA), which contained 8 wells and connected to a laptop for programming and initiation of the run, and the BentoLab (Bento Bioworks Ltd, London, UK), which contained 32 wells and was programmable as a unit. The usage of both machines was ideal because the targeted number of samples per iSeq run was 40. PCR in preparation for iSeq sequencing occurred in two separate stages. Stage one PCR proceeded as follows: 2 min at 95°C, then 35 cycles (coral slurry) or 28 cycles (seawater) of 20 sec at 95°C, 20 sec at 55°C, and 5 min at 72°C, followed by 10 min at 72°C and a final hold at 12°C. The final hold at 12°C was used due to the limitations of the BentoLab thermocycler, and samples were generally removed within an hour of the completed PCR program. Thirty ml of the resulting slurry PCR products were mixed with 6 ml 5X loading dye (Bioline, London, UK), then visualized on a 1.5% agarose gel stained with SYBR Safe DNA gel stain (Invitrogen, Thermo Fisher Scientific, Waltham, MA, USA). Bands from tissue samples that were approximately 350 bp long were excised for purification using the MinElute Gel Extraction Kit (Qiagen) following manufacturer protocols. For seawater PCR products, 5 ml of product was visualized on a 1% agarose gel to verify successful amplification, and the rest of the PCR product was purified with the MinElute PCR Purification Kit (Qiagen).

The stage two PCR procedure attached unique index primers to each sample using the Nextera XT v2 set A kit (Illumina). Purified DNA (5 ml) from stage one PCR products was added to a 50 ml reaction with the following: 5 ml Nextera index primer 1, 5 ml Nextera index primer 2, 5 ml MgCl2 (GoTaq), 10 ml 5X colorless buffer (GoTaq), 0.5 ml Taq polymerase (GoTaq), 1 ml of 10 mM dNTPs (Promega), and 18.5 ml UV-sterilized and PCR-grade water. The PCR stripettes were placed in the BentoLab or mini8 thermocyclers for the following program: 3 min at 95°C, 8 cycles of 30 sec at 95°C, 30 sec at 55°C, 30 sec at 72°C, followed by 5 min at 72°C and a hold at 12°C. A subset of PCR products was visualized on a 1% agarose gel stained with SYBR Safe DNA gel stain (Invitrogen) using 5 ml product with 1 ml 5X loading dye (BioLine) to verify bands of approximately 450 bp, indicating successful attachment of sample-specific indexes. The PCR products were purified with the MinElute PCR purification kit (Qiagen) following manufacturer protocols. Purified products were quantified using the Qubit 2.0 fluorometer dsDNA high sensitivity (HS) assay (Invitrogen) following manufacturer protocols to obtain stock concentrations in ng/ml. Concentrations were then converted to nM assuming an average amplicon length of 450 bp and average nucleotide mass of 660 g/mol. Samples were diluted to 5 nM and pooled. Pooled samples were quantified via Qubit HS assay as before, and diluted to 1 nM, quantified again, and diluted to a loading concentration of 90 pM. A 10% PhiX spike-in was added to the pooled 90 pM library and 20 ml of the resulting library was run on the iSeq 100 System using paired-end 150 bp sequencing with adapter removal. Samples were sequenced over three sequencing runs.

Data analysis

Forward reads were exclusively used for the downstream processing and data analysis due to minimal overlap between forward and reverse reads. 515F and 806R primers were removed from all sequence reads, and they were filtered for quality and chimeras using the DADA2 pipeline and amplicon sequence variants (ASVs) were generated for each sample (Callahan et al., 2016). The following parameters were used in the filterAndTrim function: trimLeft = 19, truncLen = 145, maxN = 0, maxEE = 1, rm.phix = TRUE, compress = TRUE, multithread = TRUE. This resulted in the production of 17,190 ASVs of all the same length (126 bp) across all samples. ASVs that classified to mitochondria, chloroplast, eukaryote, or an unknown Kingdom were removed from the analysis. This removed many spurious ASVs, resulting in 7,366 remaining ASVs. We further filtered our dataset based on controls, then conducted a sequence-count-based filtering method. The R package decontam (v. 1.6.0) was used to filter out DNA extraction contaminants in all seawater and tissue samples by using a combined frequency and prevalence method (Davis et al., 2018). The method identified 26 ASV contaminants, of which only 11 contained enriched frequency in DNA extraction controls, so those 11 ASVs were removed (Appendix 1). Following this, tissue/mucus slurry samples were subset and ASVs associated with the bulk (off-coral) seawater samples were removed using the prevalence method in decontam (v. 1.6.0) by comparing all slurry samples to the nine seawater-only syringe samples that were collected as methodological controls because the syringe method by nature collects a significant portion of seawater. The contamination method identified 184 ASVs most prevalent in the seawater controls (typically oligotrophic bacteria such as SAR11, Prochlorococcus, OM60 clade, Synechococcus, “Candidatus Actinomarina”, AEGEAN-169 clade, etc.), which were removed from the slurry sample ASV table. These ASVs were generally found at low relative abundance in tissue samples and was at most 0.0074 (Appendix 2). Tissue and near-coral collected seawater samples were re-merged into one large dataset for further filtering. Next, only ASVs with a count greater than 0.5 when averaged across all samples were kept for further analysis to remove sparse ASVs (present at a count of 0 in the majority of samples). This left 2,010 ASVs, which were used for all downstream analyses.

Count data were transformed to relative abundance and coral tissue microbial communities were visualized using stacked bar charts. Data were then further log transformed following addition of a pseudo count of one in preparation for beta diversity analyses. Bray-Curtis dissimilarity between samples was calculated using the R package vegan and the resulting dissimilarities were presented in a Principal Coordinates Analysis (PCoA) (Oksanen et al., 2019). PERMANOVA (Permutational Analysis of Variance) tests using 999 permutations were conducted between healthy and diseased corals to test the hypothesis that coral microbiomes are significantly different between healthy and SCTLD-afflicted tissues using the adonis function in the vegan R package (Oksanen et al., 2019). We also tested the hypotheses that species, reef location, and health state nested within species would significantly structure the microbial community. We tested the same hypotheses on the near-coral seawater directly overlying the coral colony. Dispersion of beta diversity within coral tissue samples was calculated by measuring the distance to centroid within the PCoA as grouped by health state (HH & HD compared to DD) by implementing the betadisper function in R package, vegan (Oksanen et al., 2019). Significant difference in dispersion by health state was measured with an independent Mann-Whitney U test. Variability of beta diversity was additionally measured by extracting the Bray-Curtis dissimilarity values calculated within a tissue condition (diseased or healthy).

To detect ASVs enriched in diseased coral tissue compared to healthy tissue, the R package, corncob, was employed, which modeled the relative abundance of each ASVs and tested for differential abundance between healthy and diseased coral tissue (Martin et al., 2019). Following analysis of significantly differentially abundant ASVs in the coral tissue, we hypothesized that disease-associated ASVs would be recoverable in the near-coral seawater and graphed relative abundances of each disease-associated ASV in the near-coral seawater. We additionally employed corncob to test each identified disease-associated ASV to see if it was detectable at significantly higher abundances in seawater over diseased tissue compared to healthy tissue or apparently healthy colonies. Furthermore, we compared the ASV sequences of disease-associated ASVs to existing literature on SCTLD to understand if identical taxa were associated with SCTLD in other studies.

Sequences of putative pathogens were identified to the species level, when possible, as part of the DADA2 pipeline. To obtain better genus and species-level identification of putative pathogen ASVs, as well as to relate these ASVs to other studies of coral disease-associated bacteria, we constructed phylogenetic trees for disease-associated ASVs classifying to Vibrio, Arcobacter, Rhizobiaceae, and Rhodobacteraceae. Vibrio and Arcobacter were targeted due to their increased representation in SCTLD-associated ASVs in this study as well as their previous association with SCTLD (Meyer et al., 2019) and coral disease in general (Ben-Haim et al., 2003; Ushijima et al., 2012). Rhizobiaceae and Rhodobacteraceae were targeted for phylogenetic tree analysis given their previous association with SCTLD (Rosales et al., 2020). Phylogenetic trees of coral-associated Vibrio and Rhodobacteraceae bacteria that were previously constructed from the Coral Microbiome Database (Huggett and Apprill, 2019) were used as reference trees for the insertion of SCTLD-associated ASVs that classified as Vibrio or Rhodobacteraceae. Insertion of short SCTLD-associated sequence reads was achieved using the ‘quick add marked’ tool in ARB (version 6.0.6.rev15220). Trees produced from ARB were exported using xFig. Phylogenetic trees for Arcobacter and Rhizobiaceae were constructed de novo using tools from the CIPRES Science Gateway (Miller et al., 2010). For each tree, long-read (~1,200 bp) 16S rRNA gene sequences from closely-related type strains, strains isolated from the marine environment, or clone sequences from corals were identified via BLAST searches and compiled into a .fasta file and used for a sequence alignment in MAFFT (v7.402). This sequence alignment was then used to generate a reference tree using RAxML-HPC (v.8) with the following commands to produce a bootstrapped maximum-likelihood best tree: raxmlHPC-HYBRID -T 4 -f a -N autoMRE -n [output_name] -s [input_alignment] -m GTRGAMMA -p 12345 -x 12345. Next, SCTLD-associated short sequence reads were compiled into a .fasta file and added to the long-read sequence alignment in MAFFT (v7.402) using the “--addfragments” parameter. The sequence alignment with both short and long reads and the reference tree were then used as inputs for the Evolutionary Placement Algorithm, implemented in RAxML (Berger et al., 2011). RAxML was called as: raxmlHPC-PTHREADS -T 12 -f v -n [output_name] -s [long_and_short_read_alignment] -m GTRGAMMA -p 12345 -t [reference_tree]. The output tree including short read sequences (RAxML_labelledTree.[output_name]) was visualized and saved using the interactive tree of life (iTOL v5.6.3) (Letunic and Bork, 2016).