Local dynamics of a white syndrome outbreak and changes in the microbial community associated with colonies of the scleractinian brain coral Pseudodiploria strigosa

WS outbreak follow-up

Pseudodiploria strigosa is an abundant reef building coral species at our study site, located near Puerto Morelos in the NE Mexican Caribbean (20°53′04″N, 86°50′50″W). The study site is characterized by a dense coral assemblage with more than 14 scleractinian species spatially distributed so that few of them are in direct contact with each other. We selected P. strigosa colonies given their characteristic high abundance and also because it seemed to be one of the most severely affected species (Alvarez-Filip et al., 2019). Ninety-six colonies of P. strigosa were monitored during a WS outbreak in this reef site. Each colony condition was followed through direct observation and photographic records in surveys carried out at non regular intervals (dates of surveys are given in Fig. S1), until most of the colonies showed signs of the WS disease or died (306 days in total, initiating on August 22nd, 2018). Initial colony’s photographs included a 5 cm scale to allow for colony size estimations as colony projected surface area, using the free program ImageJ. The progression of lesions was assessed by measuring tissue loss on consecutive images.

Initial WS prevalence (number of diseased colonies/total number of colonies × 100), incidence rates (IR), and survival probability were estimated for P. strigosa at the study site. The incidence rate describes how quickly an event (developing signs or dying) occurs in the sampled population at risk; the index calculation accepts colonies leaving or entering the study by incorporating time directly into the denominator of the equation. In our study, some colonies died before the study ended, others were lost for a given survey, and some others entered the study later, therefore each colony has its own tracking time. The average IR for the whole study time for both sets of colonies at risk (becoming diseased or dying if already diseased) were calculated as the total number of new cases divided by the total colony-time observed (Olsen et al., 2010). A more detailed description of the WS dynamics was accomplished by estimating the IR at each survey period for both sets at risk; these incident rates for each survey were calculated as the number of new cases per time interval between surveys, divided by the total colony-days observed during the interval. The number of colony-days were adjusted for mid interval times whenever a colony developed signs or died, since it is not known when the event occurred exactly within each interval (actuarial method; Benichou & Palta, 2005) and thus, colony time contribution is half of the corresponding interval. The data from the first survey were not considered in the analyses, unless the colony condition was apparently healthy.

To complete the description of the WS dynamics at the study site, survival probability of diseased colonies since the onset of disease (emergence of WS signs) was estimated with a non-parametric maximum likelihood estimator for interval censored data (R survival package; Therneau, 2015) and expressed as Kaplan–Meier survival curves. Comparisons of survival probability as a function of colony size were assessed by means of log-rank non-parametric tests for interval data (R interval package) that assumes disease and mortality onset times as the start and endpoint of the interval for each colony at risk of dying (Fay & Shaw, 2010), and the step curve reflects the fact that it is not known when along the interval the death event occurred. We also described gross lesions for the surveyed colonies through direct observations and photographic detailed analysis following the protocol by Work & Aeby (2006).

Coral tissue and mucus sampling

Coral samples were collected from 12 colonies of P. strigosa between September and October 2018, under permit number PPF/DGOPA-033/19 from CONAPESCA. The sampled colonies were scattered throughout the same coral reef site and their condition is shown in Table S1. We compromised for a low number of colonies to be sampled, to minimize the impact of sampling in the face of the unknown outbreak extent. Six out of the 12 colonies had no visible signs of white syndrome (named healthy), and six colonies showed visible signs of the white syndrome (named unhealthy). Two core samples (2 cm diam. to a depth of 2–3 cm into the skeleton of colonies) were collected from unhealthy colonies, one close to the edge of the lesion (UCL), and another sample 15 cm away from the lesion over healthy-looking tissue (UAL) (Fig. S2), to make a total of 18 samples (six healthy, 12 unhealthy). Also, one liter of seawater next to one healthy and one unhealthy colony was collected as background reference. After collection, each sample core was individually stored in a hermetic container and transferred to the lab in a cooler. At the lab, mucus (ca. 150 µl per core) from each sample was collected with a sterile micropipette after an air exposure (20 min) in a sterilized chamber. We processed 10 µL of mucus for DNA extraction. The rest of the mucus was either used fresh for antibacterial assays or refrigerated (7 °C) until immunological assays (2–3 days). After mucus collection, the cores were individually stored with 10% buffered formalin in filtered seawater and mailed to the Histology Lab at UABC for processing.

Histological observations

Fixed samples were rinsed in freshwater and preserved in 70% ethanol until decalcification with HCL 10% buffered solution (0.7 g EDTA, 0.14 g sodium tartrate, 0.008 g potassium sodium tetrahydrate). Tissues were then embedded in paraffin and sectioned (5 mm thick). Samples were stained with Harris’s hematoxylin and eosin routine procedure (Humason, 1967) for histopathological analysis. Microscopic morphological alterations were assessed in comparison to normal tissues from healthy and unhealthy colonies. Observations were done at different magnifications (x4, x10, x40) using a Leica DM500 light microscope (Leica, Germany) equipped with a Leica ICC50 HD camera. Images were processed using both LAS imaging software (Leica, Germany) and Adobe R Photoshop CC 2015.5®. Tissues were divided for histological interpretation into surface body wall (SBW) comprising the coenenchyme and polyp, and the basal body wall (BBW) including mesenterial filaments. We recorded the observed changes per colony as atrophy, hypertrophy, sloughing of epithelia from mesoglea, zooxanthellae depletion, or malformation. The presence of ciliates, bacterial aggregates, and suspect endolithic algae or fungi were also recorded on the datasheets to relate gross and microscopic pathology findings. Results are shown as the percentage of occurrences, considering six colonies as the total sample size for each condition.

DNA extraction for bacterial analysis

Mucus samples (10 µL) were centrifuged and then transferred to lysis tubes for genomic DNA extraction, adding extraction buffer. Seawater samples were filtered through 2.5 µm filters (Whatman™) followed by the collection of microbes by filtering with Durapore® membranes (0.45 µm and 0.22 µm). DNA from the centrifuged mucus and the microbes collected on the 0.45 and 0.22 µm filters, was extracted with DNeasy® PowerBiofilm Kit. DNA integrity was evaluated by electrophoresis in 1% agarose gels. DNA concentration and quality were assessed and only samples with adequate parameters (DNA concentration ≥ 100 ng, A_260∕280 of 1.6 to 2.0) were processed further.

16S rRNA Illumina MiSeq libraries preparation

PCR amplification of the hypervariable V3 and V4 regions of the 16S rRNA gene was performed using the primers and the conditions suggested by Klindworth et al. (2013). We used primers 341F 5′-CCTACGGGNGGCWGCAG-3′ and 805R 5′-GACTACHVGGGTATCTAATCC-3′ with Illumina overhanging adapters attached. Indexed PCR products were purified and quantified with a Qubit® 3.0 Fluorometer (Life Technologies, USA). Amplicon size was verified by capillary electrophoresis and the sequencing was carried out in CINVESTAV-Mérida using an Illumina-MiSeq platform (Illumina, USA), with the MiSeq reagent kit V3 (2 × 300), following the manufacturer’s recommendations.

Analysis of amplicon libraries

Paired end 2 × 250 reads were processed with the QIIME 2 pipeline (https://docs.qiime2.org/, see Supplementary Reference). After manual inspection, both forward and reverse reads were trimmed in position 40 from the 5′ end and truncated in position 250 from the 3′ end, to filter out low quality positions (bases with Q-score <20). DADA2 plugin (Callahan et al., 2016) was used for denoising, error correction, and removing of chimeras, to finally resolve ASVs (Amplicon Sequences Variants, unique sequences indicative of distinct taxa). The taxonomic assignment of the representative sequences of the ASVs was done with the classify-consensus-vsearch plugin (Rognes et al., 2016), using SILVA 132 database as reference. Representative sequences were aligned and masked with MAFFT (Katoh & Standley, 2013), a phylogenetic tree was built with FasterTree 2 (Price, Paramvir & Arkin, 2010). The feature table was filtered for chloroplast and mitochondrial ASVs and imported into the R environment. The feature table was rarefied to 16,600 counts per sample, and the statistical analysis performed with phyloseq (McMurdie & Holmes, 2013), Vegan (Oksanen et al., 2009), and ggplot2 v3.1.0 (Wickham, 2016) libraries. Alpha diversity indexes (Shannon, Simpson, Chao1, and observed ASVs) were calculated. An analysis to identify the ASVs with differential abundance among conditions was performed using the DESeq2 analysis with internal normalization (Love, Huber & Anders, 2014). A phylogenetic analysis was applied to the most abundant ASVs of unhealthy samples. The ASVs were searched for their closest homologs with the blastn program in the RefSeq NCBI database. The ASVs sequences and their best hits were aligned with MUSCLE (Edgar, 2004). The alignment was trimmed, and a phylogenetic tree generated with PhyML with the GTR+I substitution model, supporting branches with the aLRT method (Guindon et al., 2010).

Mucus treatment and immunological assessment

The mucus was processed in two different treatments, the Surface Mucus Layer (SML or mucus without bacteria) and the Mucus Complex (MC or the mucus with its microbiome, without treatment). For the SML treatment, the mucus was sterilized for 20 min under UV light (254 nm) and its effect corroborated by incubating an irradiated sample on marine agar (Sobel Marine Agar, Difco) at 27 °C overnight and checking for the absence of bacterial growth. For antibacterial activity assessments, we employed a swarming assay combined with a double agar overlay assay, as previously described (Rivera-Ortega & Thomé, 2018), challenging the SML and MC samples against two potential coral pathogens, Aurantimonas sp. (strain BAA-667, ATCC) and Serratia marcescens (strain BAA-632, ATCC) (Patterson et al., 2001; Denner et al., 2003). The presence of an inhibition zone was identified under a stereoscopic microscope.

We looked for the presence of melanin in tissues, evidenced as golden deposits in histological slides (Palmer, Mydlarz & Willis, 2008), as well as by the activity of phenoloxidase (PO) in mucus samples. For the quantification of the specific activity of PO in mucus, we followed a protocol described in Mydlarz & Palmer (2011), using 10 µL of mucus from each sample, with three technical replicates. Results are presented as a change in absorbance A₄₉₀ mg protein⁻¹ min⁻¹. Total protein in mucus samples was measured by the Bradford method against a standard curve prepared with BSA.

Statistical analyses

Survival probabilities of diseased colonies at risk of dying were estimated by means of a non-parametric maximum likelihood estimator (NPML) for interval censored data, to take into account the unequal interval monitoring (Therneau, 2015). All statistical analyses were performed in R v.3.6.1. For the antibacterial assays, we applied the exact Fisher test to compare differences among treatments for each tested bacterium. This test is adequate for contingency table analyses with categorical data and small sample sizes. To analyze the statistical differences of the PO activity in the mucus among treatments, we applied a Kruskal–Wallis test since the data was not normally distributed. In all tests, statistical significance was considered at p < 0.05.

White syndrome (WS) outbreak at the study site

At the beginning of the study (August 2018), WS prevalence in P. strigosa colonies was relatively high (5.2%; CI [2.2%–112.9%]). At the end of the study (306 days of observation), out of 96 tracked colonies, eight remained apparently healthy and seven were diseased. The WS outbreak on the monitored colonies showed a progressive pattern of cumulative appearance of signs and colony mortality (Fig. S1). Average incidence rates (IR) for the whole study of healthy colonies at risk of developing WS signs was 0.007/colony-days (for 100 healthy colonies at risk, 0.7 colonies will develop WS signs every day). Incidence risk per survey (IRs) showed two stages: (1) a sudden increase and maintenance of a relatively high IRs in the summer-fall of 2018, when apparently the most susceptible colonies at risk developed WS signs, with an average of 0.012/colonies per day; and (2) a relatively low IRs phase starting in winter (0.005/colonies per day), gradually affecting some of the remaining susceptible colonies of the population at risk to develop WS signs (Fig. 1). For diseased colonies at risk of dying, the average IR was 0.009/colony-days. The isolated peak early in the study (Fig. 1, IR for diseased colonies at risk of dying) is associated with a relatively large number of dead colonies (n = 6) for the observation period (298 colony-days), doubling the average IRs which, except for the spike, correspond with the average IR value. The spike does not seem associated with specific environmental conditions, but with the occurrence of multiple WS lesions in three of the six dead colonies, which rapidly succumbed to the disease. There was a significant correlation (Spearman rho = 0.89, p = 0.003) between the number of new dead colonies and time at risk in days per survey, supporting the constancy of the risk of dying once a colony developed WS signs.

Survival probabilities of diseased colonies at risk of dying were consistent with IR of colonies at risk of dying per survey estimates, where the probability of colony survival drops more or less constantly from the first survey to 0.09 (95% CI [.04–.21]) by the end of the study, as shown in the overall Kaplan–Meier survival curves (Fig. 2). To assess if colony size was related to colony survival probability, we arbitrarily divided the colonies at risk of dying into groups by quantiles (0.2, 0.4, 0.6, 0.8), as to categorize size values into small, medium, large, and very large colony sizes, as the log rank test for interval surveys does not perform well with continuous data. Contrary to expectations, the test found no significant differences (Chi Square = 3.73, p = 0.29) among the survival curves of the different groups, which in general showed the same strong decreasing trend indicated in Fig. 2. Tissue loss speed due to WS varied among colonies, occurring very fast in two cases (41 and 72 cm² per day), but much slower in the remaining colonies with an average of 10.8 ± 7.8 cm² per day.

Macroscopic signs presented high morphological variations (Table 1). According to Work & Aeby (2006) categorization, colonies showed focal, multifocal, and coalescing lesions that were subacute to acute, usually starting as a small spot anywhere on a colony (Table 1). Some colonies developed several lesions, but in all cases, the exposed intact bare skeleton was rapidly covered by turf, with remnant tissue normally pigmented (Fig. 3).

Histological observations

The surface body wall (SBW) in healthy colonies consisted of an epithelium of pseudostratified columnar cells, clearly separated from the gastrodermis by mesoglea (Fig. 4). The gastrodermis comprised a simple layer of cuboidal cells hosting symbionts, with distinct nuclei and pyrenoids. Tissues from UAL and UCL samples, manifested alterations associated with tissue damage such as hypertrophy (66%, n = 4 in UCL, 17%, n = 1 in UAL) and liquefactive necrosis (33%, n = 2 in UCL, 66%, n = 4 in UAL). These changes were evident within polyps from the basal body wall (BBW, Figs. 4B and 4D) (calicodermis and gastrodermis) to the SBW (Fig. 4A), that in general, presented minor signals of damage. In some cases, in samples from healthy (33%, n = 2) and UAL colonies (50%, n = 3), we observed tissue fragmentation from the SBW of the polyps, but with the epidermis of the SBW in good condition (Figs. 4C and 4F). Ablation of tissue from mesoglea and hypertrophy of the epithelia was frequently observed in samples from affected (UCL, 85%, n = 5; UAL, 50%, n = 3) and unaffected colonies (H, 33%, n = 2) (Figs. 4D, 4E and 4F). We did not find evidence of bacteria or crystalline inclusion bodies near the affected tissue; however, in the surrounding areas of the deteriorated calicodermis, we usually found microalgae inclusions (UCL, 50%, n = 3; UAL, 17%, n = 1), cyanobacterial mats (17%, n = 1), and ciliates (UCL, 17%, n = 1; UAL, 17%, n = 1) in association with necrosis of the tissues (Fig. S3).

Bacterial community analysis

We successfully sequenced three healthy, three UCL, and two UAL samples from the original 18 samples collected (6 of each). We obtained a total of 261,176 clean reads which were normalized to 16,600 per sample (see species richness curves, Fig. S4) and finally grouped in 1,167 ASVs (Table S2). Proteobacteria was the dominant Phylum in all samples, except in the mucus layer of UCL samples, where Bacteroidetes and Firmicutes exhibited high relative abundances (Fig. S5). At lower taxonomic level, Gamma- and Alpha-proteobacteria were dominant classes in healthy and UAL samples (Fig. S5) while Bacteroidia, Gamma-proteobacteria, Clostridia, and Bacilli were well represented for the UCL samples. The number of shared ASVs between healthy and UCL samples, and those between UAL and UCL samples, was the same (Fig. S6). Alpha diversity varied in all sampled colonies without a clear pattern (Table S3).

Even though there were a low number of samples, it was possible to apply a differential abundance analysis (DESeq2) to ASVs from healthy and UCL colonies (Fig. 5). This analysis showed the enrichment of three genera in healthy colonies, but 17 genera in unhealthy colonies (Fig. 5A). In healthy colonies, Chlorobia, Iainarchaeia, and Woesearchaeia were significantly different in comparison with the unhealthy colonies near the lesions, where six genera of the Bacteroidia class were found with specific changes (Fig. 5A). Moreover, other abundant genera such as Vibrio, Enterococcus, and Arcobacter were observed in unhealthy colonies (Fig. 5B). The classes Bacteroidia (seven genera), Gamma-proteobacteria (five genera), and Clostridia (three genera) had the most genera that changed in unhealthy colonies. The differential abundance among bacteria from the three colony conditions consisted of 53 ASVs (Fig. 5B). Within the class Bacteroidia, the genus Roseimarinus was particularly enriched. This genus had a relative abundance of 39.27 ± 23.37% on average for UCL samples. Seven different ASVs in the genus Roseimarinus were detected in unhealthy colonies (Table S2). To gain insight into the phylogenetic relationship of Roseimarinus ASVs we did a phylogenetic analysis, including also other abundant ASVs identified in unhealthy colonies (Fig. S7). The most abundant Roseimarinus ASV (tag = ebecc93ba) was clustered with Marinilabilia nitratireducens at 89.3% similarity, the closest relative known; also, for the ASV assigned with JTB215 Clostridia Firmicutes (tag = 5bb6936d9b4) the closest relative known is Peptoclostridium litorale with 89.5% similarity; while all the other ASVs were clustered with the same genus as assigned by the SILVA database. Finally, we show the relative abundance of the orders Rhodobacterales and Rhizobiales in sampled colonies, given their probable importance for white syndromes (Table S4), although our results indicate a reduction in the relative abundance of members in the class Alpha-proteobacteria (Fig. S5).

Immunological responses

The SML of healthy and unhealthy P. strigosa samples inhibited the growth of both tester strains (Table 2). Similarly, the MC of healthy P. strigosa samples inhibited the growth of S. marcescens; however, the growth inhibition of Aurantimonas sp. for the MC treatment was lower. The growth inhibition for S. marcescens by MC was also lower in the UCL and UAL samples compared with that caused by SML (Table 2). None of the treatments had significant differences among the condition of the colonies (p = 1 for Aurantimonas sp.; p = 1 for SML and p = 0.471 for MC, against S. marcescens). We did not find melanin in the tissues. We evidenced the lowest PO activity in mucus collected from healthy compared to unhealthy colonies. Conversely, the mucus of UCL samples showed the highest PO activity, while UAL samples showed intermediate levels of PO activity (Fig. 6). No significant differences were found among the three samples (H = 2, df = 2, p = 0.3679), probably due to the wide dispersion of the data for UCL samples and limited sample size. However, the data in Fig. 6 suggest that most of the UCL samples had 2 to 3 times higher PO activity than healthy and UAL samples.