A meta-analysis of the stony coral tissue loss disease microbiome finds key bacteria in unaffected and lesion tissue in diseased colonies

Summary of SCTLD microbiome studies

Initially, datasets were acquired from 17 SCTLD studies, but one study [24] did not pass quality filtering and was removed from the analysis, resulting in 16 SCTLD studies used in this meta-analysis. In addition, one Acropora spp. rapid tissue loss (RTL) disease study was included for comparison of bacteria which may be associated more generally with coral tissue loss diseases (Supplementary Table 1). The combined dataset included 2425 samples, representing various coral species and environments described below. A total of 63 miscellaneous samples such as lab controls were included in this total (Supplementary Table 1). Samples from the studies were sequenced using five primer pairs: CS1-515F/CS2-806R [31] with additional 5’ linker sequences [32] (n = 79), 515FY [33]/806RB [34] (n = 1219), S-D-Bact-0341-b-S-17/S-D-Bact-0785-a-A-21 [35] (n = 31), 515F/806R [31] (n = 49), and 515F [31]/Arch806R [36] (n = 984; Fig. 1A). Although five primer pairs were used across studies, only the forward reads were evaluated in this analysis (see “Methods”). A description of the differences between 515F primers can be found in detail [34].

Samples were collected throughout Florida and the U.S. Virgin Islands (USVI). Field samples totaled 1274, representing 40 sites, and a further 1088 samples were from aquaria (i.e., laboratory-based experiments; Fig. 1). Thirteen SCTLD-susceptible coral species were included, with Montastraea cavernosa (MCAV; n = 543) and Orbicella faveolata (OFAV; n = 357) most represented and Pseudodiploria clivosa (PCLI; n = 6) and Orbicella franksi (OFRA; n = 7) least represented (Fig. 1). Coral samples (n = 2031) were from three compartments: mucus only (n = 393), mucus and surface tissue (tissue slurry; n = 1585), and skeleton samples with embedded coral tissue (tissue slurry skeleton; n = 53). Seawater (n = 198) and sediment (n = 133) samples from both the field and aquaria experiments also were included to evaluate potential sources of transmission of disease-associated bacteria (Fig. 1B). For seawater from aquaria experiments, 18 L samples were collected [27], while in the field between 60 mL and 1 L samples were collected [11, 25]. In sediment aquaria experiments, 2 mL samples were collected [12], and in the field, approximately 5 mL samples were collected (of the 5 mL, DNA was extracted from 0.25 g sediment [11]). Coral samples represented three SCTLD health states: apparently healthy colonies (AH), which was the most represented (n = 1021), followed by lesions on diseased colonies (DL; n = 661), and unaffected areas on diseased colonies (DU; n = 349; Fig. 1C). AH represents grossly normal tissue, DU grossly normal tissue on diseased colonies, and DL grossly abnormal tissue.

Differences in the microbial composition were found in AH corals among zones (vulnerable, endemic, and epidemic)

Differences in alpha-diversity were tested among three SCTLD zones: vulnerable (i.e., locations where the disease had not been observed/reported), endemic (i.e., locations where a disease outbreak had moved through the reef and no or few colonies had active lesions), and epidemic (i.e., locations where the outbreak was active and prevalent). For alpha-diversity, for AH field-sourced samples, after filtering, 41,504 amplicon sequence variants (ASVs) remained, which were reduced to 15,021 following rarefaction. Among the filtered AH samples, Shannon (alpha) diversity from the vulnerable zone was slightly higher (estimated marginal means (emmean) = 3.95) compared to the epidemic zone (emmean = 3.70), but this was not significant (Supplementary Fig. 1). For beta-diversity, both within and between-group differences were tested using a filtered counts table. Within-group beta-diversity (variation in microbial composition or dispersion) was not different between zones, but was significant for all comparisons between zones (PERMANOVA, P-adjusted (Padj) <0.03; Fig. 2A). Differential abundance analysis found 61 ASVs enriched between vulnerable and endemic sites (Fig. 2B, Supplementary Fig. 2, and Supplementary Table 2). In the endemic zone, the orders Synechococcales (Cyanobium PCC-6307; log-fold = 12.67) and an uncultured Flavobacteriales (log-fold = 9.96) contained ASVs with the highest log-fold change, but the order Flavobacteriales was the group of bacteria with the most enriched ASVs (n = 13), followed by SAR11 clade (n = 4) and Rhodobacterales (n = 3). Fewer ASVs were enriched between the vulnerable and epidemic zones (n = 31; Fig. 2C and Supplementary Table 3), with the highest log-fold ASV changes found in the orders Burkholderiales (Delftia; log-fold = 5.84) and Peptostreptococcales–Tissierellales (Fusibacter; log-fold = 5.65). Like in endemic sites, Flavobacteriales was the group with the most enriched ASVs in the epidemic zone (n = 5) and were detected in the three disease states (AH, DU, and DL; Supplementary Fig. 2).

Biome had the highest correlation to bacterial beta-diversity

Microbial dispersion at the ASV level was found to be different across primers, study, biome, year, all coral species, and sample type (Permutest: P < 0.01; Fig. 3). A PERMANOVA test for differences between microbial composition at the ASV level was also significant across all factors, with coral species having the highest correlation (R² = 0.21; Fig. 3E) and disease state showing the lowest correlation (R² = 0.04). Biome (i.e., aquaria and field) had the largest correlation to principal component 1 (PC1, R2 = 0.73; Supplementary Fig. 3) compared to other tested metadata factors, and showed a distinct separation when visualized (Fig. 3C). This was also evident even in sediment and seawater samples that were collected in aquaria studies, which clustered with coral samples from aquaria studies and not with field sediment and seawater samples. Given this pattern, SCTLD-affected corals (with the removal of Acropora spp.) were first combined (i.e., both aquaria and field) and analyzed. In subsequent analyses, the SCTLD-affected corals were divided by biome to identify potential differences between the two.

Bacterial communities differ across disease states, but this may depend upon the biome

When both biomes were combined (Fig. 4A), DL microbial communities were the most highly dispersed compared to both AH and DU (Padj < 0.01 each), but AH and DU were not different. Pairwise PERMANOVA was significant for all comparisons (Padj < 0.001 each; Fig. 4A). Among aquaria samples (Fig. 4B), the dispersion was lower in DU vs both DL (Padj < 0.01) and AH (Padj < 0.005), and was also dissimilar in AH vs DL (Padj = 0.0015). Like the combined samples, all aquaria samples were different in the pairwise PERMANOVA (Padj < 0.001 each; Fig. 4B). In field samples (Fig. 4C), the dispersion was only different between DL and AH. All pairwise PERMANOVA comparisons were significant in the field samples: AH vs DU (Padj < 0.02), AH vs DL (Padj < 0.01), and DU vs DL (Padj < 0.03; Fig. 4C).

Samples were also evaluated for alpha-diversity by disease state in each biome. After quality filtering and rarefaction across disease states, 39,513 ASVs remained. For aquaria and field samples combined, pairwise comparisons showed differences in Shannon diversity for AH vs DU and DL vs DU (Padj < 0.0001 each) but not AH vs DL, with mean alpha-diversity lowest in DL (emmean = 3.42) and highest in DU (emmean = 3.85; Supplementary Fig. 4A). In aquaria samples only, there were no differences in Shannon diversity by disease state, likely due to the low sample size of DU (n = 27, Supplementary Fig. 4B). In field samples, only DU vs DL was different (Padj < 0.01) with DU also showing the highest mean (emmean = 3.90) and DL the lowest (emmean = 3.63; Supplementary Fig. 4C) alpha-diversity.

When comparing differences in mean relative microbial abundances within disease states across biomes, AH samples differed between aquaria and field (Supplementary Fig. 4D): the orders Rhodobacterales (14.20 ± 5.2%) and Cytophagales (9.02 ± 12.32%) were dominant in aquaria samples, but in field samples, the dominant orders were Flavobacteriales (5.75 ± 2.15%) and Synechococcales (3.77 ± 5.88%). Like AH aquaria samples, DU aquaria samples had the highest mean relative abundances in Rhodobacterales, but at a much lower percentage (1.06 ± 3.81%). The DU field samples were also similar to their AH counterparts, showing the highest relative abundances in Flavobacteriales (6.43 ± 1.89%) and Synechococcales (4.45 ± 6.26%). In the DL samples, both aquaria and field samples were dominated by Rhodobacterales, but the aquaria samples had a higher relative abundance of Rhodobacterales (15.34 ± 6.84%) than samples from the field (6.61 ± 4.12%). As with aquaria AH samples, Cytophagales (3.28 ± 11.22%) were also the second most relatively abundant order in DL aquaria samples but were not dominant in field DL samples. Peptostreptococcales–Tissierellales was a dominant DL member at similar mean relative abundances in both aquaria (3.21 ± 6.40%) and field samples (3.79 ± 9.06%; Supplementary Fig. 4D).

Indicator taxa were detected across sample types and zones

The combined three coral compartments (mucus, tissue slurry, and tissue slurry skeleton), from both field and aquaria, yielded a total of 109 differentially abundant ASVs between AH vs DU (Fig. 5A, Supplementary Fig. 5A, and Supplementary Table 4). DU mucus samples showed the highest log-fold change compared to AH in the orders Flavobacteriales (NS5 marine group; log-fold = 6.33) and Synechococcales (Cyanobium PCC-6307; log-fold = 6.19), with Flavobacteriales having the most enriched ASVs (n = 9). Similarly, DU tissue slurry samples were most enriched in Synechococcales (Synechococcus CC9902; log-fold = 20.04) and Flavobacteriales (NS5 marine group; log-fold = 12.71), with Flavobacteriales having the most enriched ASVs (n = 9). Tissue slurry skeleton sample comparisons of AH vs DU identified no ASVs enriched in DU. In addition to coral compartment samples, ASVs enriched in AH and DU samples were also present within sediment and seawater samples. The enriched taxa were also detected across the three zones with Flavobacteriales and Synechococcales found at higher relative abundances in sediment and seawater of endemic and epidemic zones compared to the vulnerable zone (Supplementary Fig. 5). However, some taxa such as Burkholderiales and Staphylococcales were also present at high relative abundances in lab control samples compared to other taxa and thus could be artifacts of contamination (Fig. 5B) [37].

The three combined coral compartments yielded fewer differentially abundant ASVs in AH vs DL (n = 79; Fig. 6A, Supplementary Fig. 6, and Supplementary Table 5) compared to AH vs DU (Fig. 5A). In DL mucus samples, ASVs from the orders Desulfovibrionales (Halodesulfovibrio; log-fold = 13.96) and Rhodobacterales (Shimia; log-fold = 13.18) were the most enriched, and Rhodobacterales had the most enriched ASVs overall (n = 8). In DL tissue slurries, the ASVs with the highest enrichment were two Rhodobacterales from an uncharacterized genus (log-fold = 15.77) and one from the genus Tropicibacter (log-fold = 13.46). Rhodobacterales were also the order with the most enriched ASVs in DL compared with AH tissue slurries (n = 14), followed by Peptostreptococcales–Tissierellales (n = 6). Among tissue slurry skeleton samples, only one ASV was enriched in DL (Burkholderiales, Achromobacter; log-fold = 1.49), but just like in DU samples, Burkholderiales was found at high relative abundances in lab controls and therefore could be a laboratory artifact (Fig. 6B). ASVs enriched in DL were also found in sediment and seawater (Fig. 6B and Supplementary Fig. 6); however, Rhodobacterales was commonly and abundantly found in the vulnerable zone, while Peptostreptococcales–Tissierellales was absent or found at low relative abundances (Supplementary Fig. 6).

When OFAV and MCAV (the most-sampled coral species) were removed from the analysis, similar patterns were still identified in beta-diversity (Supplementary Fig. 7A) and differential abundance when compared to the analysis of all SCTLD-susceptible species. For DU, 18 (35.3%) ASVs were shared between the two analyses (i.e., with vs without OFAV and MCAV), but more unique ASVs were found enriched in the analysis without OFAV and MCAV compared to the analysis that included all SCTLD-susceptible corals (Supplementary Fig. 7B). Still, the two analyses shared more enriched bacterial families/orders compared to the number that was enriched only within each individual analysis. In DL, the differential abundance analysis without OFAV and MCAV compared to that with all SCTLD-susceptible coral species showed that the majority of enriched ASVs were shared (n = 25; 39.1%) between the two analyses (Supplementary Fig. 7C).

Indicator taxa presence varied across coral species and studies

Six coral species were represented by a high number of samples (n > 76 samples each), and all ASVs only enriched in DU were found within all of those species. The seven coral species with lower sampling frequencies (n < 76 each) varied in the numbers of DU-enriched ASVs present (Supplementary Fig. 8A). For example, Dichocoenia stokesii (DSTO) contained all DU-enriched taxa, and Stephanocoenia intersepta (SINT) had all genera present but one, which belonged to Flavobacteriales. In comparison, Pseudodiploria clivosa (PCLI) had the fewest DU-enriched taxa (n = 3) among the coral species. Four orders were not present in Acropora spp. samples and included: Blastocatellales, Pirellulales, Sphingobacteriales, and Peptostreptococcales–Tissierellales. Across studies, the order Sphingobacteriales was not found in any aquaria study but was found in 50% of field studies (Supplementary Fig. 8B). In addition, no aquaria study had representatives from all DU-enriched taxa, likely because of low DU samples in aquaria, but four field studies had all taxa. The two studies with the fewest representatives were studies that used V3–V4 primers (Supplementary Table 1).

The ASVs enriched only in DL were also present in all high-frequency coral species, while none of the low-frequency coral species had all of the DL-enriched taxa (Supplementary Fig. 9A). PCLI possessed the fewest DL-enriched genera (n = 9) followed by Orbicella franksi (OFRA; n = 15). More DL-enriched orders (n = 11) were absent from Acropora spp. corals than DU-enriched orders (n = 4); the DL orders not present in Acropora were: Bacteroidales, Beggiatoales, Burkholderiales, Cellvibrionales, Clostridiales, Desulfovibrionales, Oligoflexales, Peptostreptococcales–Tissierellales, Rhizobiales, Thiotrichales, and Verrucomicrobiales. Across studies, three had all the DL-enriched orders (all aquaria studies), and the fewest orders were present in those which used V3–V4 primers (Supplementary Fig. 9B), as with the DU-enriched orders.

Alphaproteobacteria and Clostridia were found associated with SCTLD bacterial community interactions

In a network analysis of co-associated ASVs, a total of nine modules were identified, with two that were significantly and positively correlated to AH (R² = 0.1 and 0.26), three to DU (R² = 0.12, 0.31, and 0.46), and four to DL (R² = 0.17, 0.22, 0.46, and 0.47; Supplementary Fig. 10). The modules with the highest positive correlation to each disease state had 134 (AH; blue), 158 (DU; green), and 146 (DL; pink) co-abundant ASVs (Supplementary Fig. 10) and were used for undirected network analysis (Fig. 7). Although AH had the second largest module, the network was smaller than both DU and DL, with only 56 ASV nodes and 59 edges (connections between nodes). DU had the largest network, with 138 nodes and 293 edges, followed by DL, with 123 nodes and 204 edges. In AH, the node with the most neighbors (n = 7) was from the class Polyangia (order Polyangiales), which was also considered a key player (i.e., provides cohesiveness, connectedness, and is embedded in a network [38]; Fig. 7) in the AH network. The two nodes with the highest correlation to the blue weighted correlation network analysis (WGCNA) module (Supplementary Fig. 10) were from the class Bacteroidia (Chitinophagales; R² = 0.88 and 0.87).

In the DU network, highly connected nodes included three orders from the class Alphaproteobacteria (SAR11 clade (n = 16), Rhodobacterales (n = 11), and Rhodospirillales (n = 11); Fig. 7). Alphaproteobacteria were among the classes assigned as key players, but additional key players included: Cyanobacteria, Bacteroidia, and Polyangia. The nodes most highly correlated to their respective WGCNA modules were SAR86 clade (R² = 0.88) and Rhodospirillales (R² = 0.88).

The DL network had nodes with the most neighbors compared to AH and DU and was driven by Alphaproteobacteria (two Rhodobacterales nodes (n = 22 and n = 16), and Rhizobiales (n = 9)), and Bacteroidia (Flavobacteriales (n = 12); Fig. 7). While Alphaproteobacteria (Rhodobacterales and Rhizobiales) were found as key players in DL, Flavobacteriales were not. Additional key players in DL included Clostridia, Chlamydiae, and Campylobacteria. The class Clostridia had the highest correlations to the DL pink WGCNA module (Peptostreptococcales–Tissierellales; R² = 0.77 and Lachnospirales R² = 0.76; Supplementary Fig. 10). The most prevalent classes in DL networks were Alphaproteobacteria (n = 39; mainly Rhodobacterales, n = 29) and Clostridia (n = 23; mainly Peptostreptococcales–Tissierellales, n = 13).

The top microbial functional pathways were more enriched in DL compared to AH and DU

To identify differences in the potential microbial function between disease states, we used the SSU 16S rRNA gene for functional predictions. There was a total of 6307 differently abundant (Padj < 0.05) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways identified across AH (n = 2482), DU (n = 1403), and DL (n = 2422). Of the top ten KEGG pathways, three were enriched in DU and six in DL (Supplementary Fig. 11A). The most enriched pathway in DU was 4-hydroxybutyrate dehydrogenase (effect size = 0.25), and in DL was phospholipase C/alpha-toxin (effect size = 0.97). A total of 392 differentially abundant MetaCyc pathways were found across AH (n = 148), DU (n = 104), and DL (n = 139). Out of the top ten pathways, nine were enriched in DL and one in DU (Supplementary Fig. 11B). Biotin biosynthesis II was the most enriched pathway in DL (effect size = 0.80), and ADP-L-glycero-β-D-manno-heptose biosynthesis in DU (effect size = 0.05).

Source: Ecology - nature.com