Machine learning reveals distinct gene expression signatures across tissue states in stony coral tissue loss disease

Abstract

Stony coral tissue loss disease (SCTLD) has rapidly degraded Caribbean reefs, compounding climate-related stressors and threatening ecosystem stability. Effective intervention requires understanding the mechanisms driving disease progression and resistance. Here, we apply a supervised machine learning approach-support vector machine recursive feature elimination-combined with differential gene expression analysis to describe SCTLD in the reef-building coral Montastraea cavernosa and its dominant algal endosymbiont, Cladocopium goreaui. We analyse three tissue types: apparently healthy tissue on apparently healthy colonies, apparently healthy tissue on SCTLD-affected colonies and lesion tissue on SCTLD-affected colonies. This approach identifies genes with high classification accuracy and reveals processes associated with SCTLD resistance, such as immune regulation and lipid biosynthesis, as well as processes involved in disease progression, such as inflammation, cytoskeletal disruption and symbiosis breakdown. Our findings support evidence that SCTLD induces dysbiosis between the coral host and Symbiodiniaceae and describe the metabolic and immune shifts that occur as the holobiont transitions from healthy to diseased. This supervised machine learning methodology offers a novel approach to accurately assess the health states of endangered coral species, with potential applications in guiding targeted restoration efforts and informing early disease intervention strategies.

Keywords: coral; gene expression; machine learning; stony coral tissue loss disease; symbiosis; transcriptomics.

Figure 1.

Experimental design and analysis workflow. (a,b) Photographs showing the three tissue health states collected: HH, HD and LD. (c) Sampling data. (d) Feature selection pipeline using a combination of SVM-RFE and DE t-statistic with the sigFeature R package [52]. (e) External stratified k-fold cross-validation (CV) results on the top 400 M. cavernosa features from each tissue health state. (f) External stratified k-fold CV results on the top 400 C. goreaui features from each tissue health state (SVM-RFE, support vector machine recursive feature elimination; DE, differential expression). (a,b) Photo credit: Amy Apprill.

Figure 2.

Functional enrichment of upregulated and downregulated genes within the top 500 M. cavernosa features from each tissue health state. The number of genes (features) in each list of upregulated and downregulated genes from each tissue health state is shown beneath each bubble plot. Upregulated enrichments are shown on the top panel, and downregulated enrichments are shown on the bottom panel. The size of bubbles represents the log-transformed false discovery rate (FDR) for the enrichment, corrected for multiple testing using the Benjamini–Hochberg procedure. Colour indicates the enrichment strength, calculated as log₁₀(observed/expected), where ‘observed’ is the number of genes associated with a function and ‘expected’ is the number expected by chance. Darker colours reflect functions that exceed random expectations. Strength values from the downregulated enrichments were multiplied by −1 for visualization purposes.

Figure 3.

Functional enrichment of upregulated and downregulated genes within the top 500 C. goreaui features from each tissue health state. The number of genes (features) in each list of upregulated and downregulated genes from each tissue health state is shown beneath each bubble plot. Upregulated enrichments are shown on the top panel, and downregulated enrichments are shown on the bottom panel. The size of bubbles represents the log-transformed FDR for the enrichment, corrected for multiple testing using the Benjamini–Hochberg procedure. Colour indicates the enrichment strength, calculated as log₁₀(observed/expected), where ‘observed’ is the number of genes associated with a function and ‘expected’ is the number expected by chance. Darker colours reflect functions that exceed random expectations. Strength values from the downregulated enrichments were multiplied by −1 for visualization purposes.

Figure 4.

Top features in HH M. cavernosa. (a) Relative expression heatmap of the top 15 HH features from M. cavernosa. Red boxes signify elevated expression relative to the row mean, and blue boxes signify lowered expression relative to the row mean. The top feature is shown in bold, and genes plotted in (b) end in an asterisk. (b) Boxplots showing the rlog-transformed expression of five selected features from (a), organized by tissue health state. p-values represent two-sample Wilcoxon test results. The colour of boxplots corresponds to tissue health state. Tmem145: TGFβ signalling; PROKR1: G protein-coupled receptor signalling; INHBB: TGFβ signalling; NAIP: apoptosis regulation; Dmbt1: microbial homeostasis. Boxplot elements: centre line, median; box limits, upper and lower quartiles; whiskers, 1.5× interquartile range; points beyond whiskers, outliers.

Figure 5.

Top features in HD M. cavernosa. (a) Relative expression heatmap of the top 15 HD features from M. cavernosa. Red boxes signify elevated expression relative to the row mean, and blue boxes signify lowered expression relative to the row mean. The top feature is shown in bold, and genes plotted in (b) end in an asterisk. (b) Boxplots showing the rlog-transformed expression of five selected features from (a), organized by tissue health state. p-values represent two-sample Wilcoxon test results. The colour of boxplots corresponds to tissue health state. Ucp2: antioxidant; Rab5a: symbiont uptake; cfap298: cilia motility; fam166b: cilia motility; dnaJ: chaperone protein. Boxplot elements: centre line, median; box limits, upper and lower quartiles; whiskers, 1.5× interquartile range; points beyond whiskers, outliers.

Figure 6.

Top features in LD M. cavernosa. (a) Relative expression heatmap of the top 15 LD features from M. cavernosa. Red boxes signify elevated expression relative to the row mean, and blue boxes signify lowered expression relative to the row mean. The top feature is shown in bold, and genes plotted in (b) end in an asterisk. (b) Boxplots showing the rlog-transformed expression of five selected features from (a), organized by tissue health state. p-values represent two-sample Wilcoxon test results. The colour of boxplots corresponds to tissue health state. TPM1: cytoskeleton stabilization; col2a1 [a]: ECM structural component; col1a2 [a]: ECM structural component; GFPL: photoprotection; sodA: antioxidant. Boxplot elements: centre line, median; box limits, upper and lower quartiles; whiskers, 1.5× interquartile range; points beyond whiskers, outliers.

Figure 7.

Selected top features in C. goreaui. Boxplots show the rlog-transformed expression of three selected features from each tissue health state. HH features are shown in the top panel, HD features in the middle panel and LD features in the bottom panel. The top feature from each tissue health state is shown in the first column in bold. P-values represent the two-sample Wilcoxon test results. The colour of boxplots corresponds to tissue health state. Boxplot elements: centre line, median; box limits, upper and lower quartiles; whiskers, 1.5× interquartile range; points beyond whiskers, outliers. Gene functions are listed in electronic supplementary material, table S20.

Figure 8.

Characterization of SCTLD tissue health states in M. cavernosa and C. goreaui. Processes are shown in regular font, and selected features are shown in italic font. Processes are inferred from both functional enrichment of the top 500 features and analysis of selected features (genes listed). Upregulated processes/features are denoted with a red arrow, and downregulated processes/features are denoted with a blue arrow. (FA, fatty acid; ROS, reactive oxygen species).