Contrasting heat stress response patterns of coral holobionts across the Red Sea suggest distinct mechanisms of thermal tolerance

doi:10.21203/rs.3.rs-117181/v1

Download PDF

Article

Contrasting heat stress response patterns of coral holobionts across the Red Sea suggest distinct mechanisms of thermal tolerance

https://doi.org/10.21203/rs.3.rs-117181/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 03 Aug, 2021

Read the published version in Molecular Ecology →

Version 1

posted

You are reading this latest preprint version

Corals from the northern Red Sea, in particular the Gulf of Aqaba (GoA), have exceptionally high bleaching thresholds approaching >5°C above their maximum monthly mean (MMM) temperatures. These elevated thresholds are thought to be due to historical selection, as corals passed through the warmer Southern Red Sea during re-colonization from the Arabian Sea. To test this hypothesis, we determined thermal tolerance thresholds of GoA versus Central Red Sea (CRS) Stylophora pistillata corals using the Coral Bleaching Automated Stress System (CBASS) to run a series of standardized acute thermal stress assays. Relative thermal thresholds of GoA and CRS corals were indeed similar and exceptionally high (~7°C above MMM). However, absolute thermal thresholds of CRS corals were on average 3°C above those of GoA corals. To explore the mechanistic underpinnings, we determined gene expression response and microbiome dynamics of coral holobiont compartments. Transcriptomic responses differed markedly, with a strong response to the thermal stress in GoA corals versus a remarkably muted response in corals from the CRS. This pattern was recapitulated in the algal symbionts that showed site-specific genetic differentiation. Concomitant to this, a subset of coral and algal genes showed temperature-induced expression in GoA corals, while exhibiting fixed high expression, i.e. front-loading, in CRS corals. Bacterial community composition of GoA corals changed dramatically under heat stress, whereas CRS corals displayed consistent assemblages, indicating distinct microbial response patterns. Our work demonstrates distinct patterns underlying thermal tolerance across spatial scales, even for the same species and ocean basin. We interpret the response of GoA corals as that of a resilient population approaching a tipping point in contrast to a pattern of consistently elevated thermal resistance in CRS corals that cannot further attune. Such response differences suggest distinct thermal tolerance mechanisms that affect the response of coral populations to ocean warming.

Oceanography

General Microbiology

coral bleaching

ocean warming

thermal tolerance

thermal resilience

coral holobiont

coral metaorganism

heat stress

short-term acute heat stress assays

CBASS Coral Bleaching Automated Stress System

Coral reefs are at the brink of ecological collapse due to climate change-driven mass coral bleaching. A concern is that without human-assisted intervention, coral reef ecosystems are globally lost in the coming decades. However, we are still far away from a criterion framework to provide information regarding tolerant coral species, populations, or genotypes for incorporation into efforts that mitigate anthropogenic impact. Here we determined thermal tolerance thresholds of corals across a large geographical range with subsequent elucidation of the underlying molecular signatures. Our analysis demonstrates that molecular response patterns can be indicative of differences in thermal tolerance phenotypes and might prove useful in prospecting the fate of coral populations at scale and their prioritization for conservation and restoration efforts.

Coral reefs have undergone drastic declines in recent years due to increasing anthropogenic pressure. In particular, coral bleaching, i.e. the loss of Symbiodiniaceae, due to ocean warming is now among the main drivers of reef degradation (Hughes et al. 2018). Without their algal symbionts, corals lose their primary source of nutrition, often leading to widespread coral mortality. Hence, a better understanding of coral bleaching and the mechanistic underpinnings that influence susceptibility or resilience to thermal stress are of major importance.

The common notion is that corals can bleach and suffer mortality at just 1–2 °C above their maximum monthly mean (MMM) temperatures (Jokiel and Coles 1977, 1990; Glynn and D’croz 1990), but we actually know little about empirical temperature thresholds, as most of our knowledge comes from predictive modeling, observational studies, or long-term thermal exposure experiments. Similarly, susceptibility of corals to bleaching has been attributed to genomic differences of corals and algal symbionts (DeSalvo et al. 2008, 2010; Barshis et al. 2013; Aranda et al. 2016; Bhattacharya et al. 2016; Voolstra et al. 2017), association with different microalgal symbiont species (Hume et al. 2016, 2020; Ziegler et al. 2017a), or changes in the associated bacteria (Ziegler et al. 2017b; Osman et al. 2020).However, it is unknown whether the same mechanisms or factors contribute to thermal tolerance of corals at large, i.e., across species and sites. Corals from the Northern Red Sea, for instance, harbor exceptionally high thermal tolerance limits approaching > 5 °C above their summer maxima, which are lethal for corals elsewhere (Fine et al. 2013; Krueger et al. 2017; Osman et al. 2018). The Red Sea at large is one of the warmest ocean basins where corals thrive despite water surface temperatures commonly exceeding 30 °C in the summer (Roik et al. 2016, 2018; Ziegler et al. 2019b; Voolstra et al. 2020). Nevertheless, we are uncertain about the consistency of this thermal tolerance across the gradient of the Red Sea, and whether the same holobiont assemblages and molecular mechanisms are at play across different spatial scales.

Here we sought to investigate the consistency of thermal tolerance limits across spatial scales using the Coral Bleaching Automated Stress System (CBASS) (Voolstra et al. 2020) to run a series of identical short-term thermal stress assays on colonies of the coral Stylophora pistillata from three warmer Central Red Sea (CRS) sites (0.3–10 km overwater distance) and a cooler Northern Red Sea site in the Gulf of Aqaba (GoA) (~ 1,000 km overwater distance). Subsequent molecular examination allowed for a direct comparison of the response patterns underlying the determined temperature thresholds and revealed differences in the gene expression response and holobiont configurations that suggest distinct thermal tolerance mechanisms are at play across spatial scales. Our work demonstrates the utility of running standardized short-term heat stress assays for determination of thermal thresholds with subsequent molecular interrogation, which provide a powerful framework to identify differences in thermal tolerance response patterns and associated holobiont signatures (e.g., alleles, genes, microbes).

Study sites and sample collection

In August 2018, fragments of Stylophora pistillata coral colonies were collected from one site in the northern Red Sea in the Gulf of Aqaba (GoA), the Interuniversity Institute for  Marine Science (IUI) coral nursery (ICN) (N 29.50000° E 34.93333°), and three sites in the central Red Sea, Al Fahal reef (AF) (N 22.29634° E38.95914°), the ocean-facing exposed site of Tahala reef (ExT) (N 22.26189°, E 39.04878°), and the land-facing protected site of Tahala reef (PrT) (N 22.26302° E 39.05165°). At each of the four sites, 7 coral colonies were sampled, with 7 ramets collected per colony, using SCUBA at depths of 2-8 m. Colonies were sampled at least 5 m apart from each other to minimize the potential of sampling clonal genotypes. Sampled specimens were stored in Ziploc plastic bags upon underwater collection and transported to shore in a cooler filled with seawater where they were subjected to short-term thermal stress assays (see below).

Short-term thermal stress assays

Ramets of 7 colonies of Stylophora pistillata from each of the four sites were subjected to short-term thermal stress assays. For the GoA site, experiments were run using a set of small containers with manual heating and chilling controls. For the central Red Sea sites, we employed the Coral Bleaching Automated Stress System (CBASS) (Voolstra et al. 2020). For each site, coral ramets from all 7 colonies were distributed across four tanks, so that ramets from each colony (i.e., genotype) were exposed to each one of four distinct temperature treatment conditions (baseline = 30°C, medium = 33°C, high = 36°C, extreme = 39°C), totaling 112 fragments (7 colonies x 4 temperatures x 4 sites) (Table S1). The four temperature treatment conditions were as follows: Start of the experiment at noon (12.00h). The baseline (control) tank was maintained at 30°C for the entire duration of the experiment. The three other tanks were heated to 33°C, 36°C, and 39°C over a period of three hours. The respective temperatures were held for three hours, and then decreased to 30°C over the course of one hour (19.00h). All corals were kept at 30°C overnight until sampling the following morning, completing the 18-hour short-term thermal stress assay (Voolstra et al. 2020). All tanks were continuously supplied with Red Sea seawater and light intensity was set to photosynthetic active radiation of 600 µmol photons m⁻²s⁻¹ to match in situ light fields (Grottoli et al. 2020).

Photosynthetic efficiency and visual coral bleaching assessment

Experimental tanks were covered with a tarp at 19.00h (following the heat stress and subsequent ramping down to the control temperature) to ensure complete darkness. After one hour (at 20.00h), we measured dark-adapted maximum quantum yield (F_v/F_m) of photosystem II of all coral fragments in all temperature treatments using a Pulse Amplitude Modulated (PAM) fluorometer (WALZ, Effeltrich, Germany). Temperature tolerance thresholds were determined for each site as the mean (across all genets) temperature setpoint at which photosynthetic efficiency dropped to 50% of the value at baseline temperatures, defined as the Effective Dose 50 or ED50 following (Evensen et al. 2020) using the DRC package in R (Ritz et al. 2015). Statistical differences among site-specific ED50s were assessed via a one-way ANOVA (aov function) with individual genet ED50s as the response variable and site as the factor, with subsequent Tukey’s post-hoc testing via the TukeyHSD function in R. R script and data are available at GitHub (https://github.com/reefgenomics/CBASS84). For visual bleaching assessment, coral fragments were photographed before being snap-frozen in liquid nitrogen and stored at -80°C until further processing (Data S1). All coral fragments showed increased paling/whiteness at higher heat stress temperatures.

RNA-Seq library generation

For nucleic acid extraction, coral tissue was sprayed off from frozen fragments using airflow from a sterile, 1000 µl pipette tip connected via a rubber hose to a benchtop air pressure valve and 1 ml of ice‐cold 0.22 µm filtered seawater (FSW) for a maximum of 3 min. Following this, an aliquot of 100 µl of the tissue slurry was added to each of 400 µl of Buffer RLT (Qiagen) and 400 µl of Buffer ATL (Qiagen) for subsequent RNA and DNA isolation, respectively. Total RNA was isolated from corals exposed to 30°C, 33°C, and 36°C to make a total of 84 fragments (7 colonies x 3 temperatures x 4 sites). RNA was isolated using the Qiagen RNeasy 96 kit according to the manufacturer’s recommendations. Total RNA was quantified using the Qubit RNA HS Assay Kit on a Qubit 2.0 Fluorometer (Thermo Fisher Scientific, Waltham, US) and quality-checked using the Agilent RNA 6000 Nano kit on the Agilent 2100 bioanalyzer (Agilent Technologies, Palo Alto, US). RNA-Seq libraries were prepared using the TruSeq Stranded Total RNA kit (Illumina, San Diego, US) and sequenced on 8 lanes of an Illumina HiSeq 4000 using 2 × 151bp paired-end reads in groups of 8-11 samples/lane. Samples from the 39°C heat stress were not processed, given that we observed visible tissue loss that potentially compromised sample integrity (Data S1).

Gene expression analysis

Paired‐end Illumina reads were checked for technical artifacts using TrimGalore version 0.4.3 (Krueger 2012) following Illumina default quality filtering steps. Reads were further trimmed for low-quality ends (Phred score < 20) and cleaned up for adapter contamination with TrimGalore. For sequence alignment, reference genomic gene sets for Stylophora pistillata (v1.0) (Voolstra et al. 2017) and Symbiodinium microadriaticum (v1.0) (Aranda et al. 2016) were obtained from spis.reefgenomics.org (n = 25,769 genes) and smic.reefgenomics.org (n = 49,109 genes), respectively (Liew et al. 2016), and gene sets were combined to create a merged reference. Transcript abundance estimation was performed by using kallisto v0.44.0 (Bray et al. 2016). Differential gene expression analysis was performed using DESeq2 package v1.22.2 (Love et al. 2014) in R after importing kallisto transcript abundance estimates with tximport package v1.14.2 (Soneson et al. 2015). To account for large dispersion with low read counts and create more accurate log2 fold change (LFC) estimates, lfcShrink function for shrinking LFC estimates was applied. Transcripts with s-values (Stephens 2017) smaller than 0.005 were defined as significantly differentially expressed (Data S2). In addition, transcript quantifications in transcripts per million (TPM) were obtained using kallisto v0.44.0 (Bray et al. 2016) (Data S2). The analysis workflow, implemented as custom python and R scripts, is available at GitHub (https://github.com/reefgenomics/CBASS84). Sequence data determined in this study are available under NCBI BioProject PRJNA681108 (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA681108). To determine front- and back-loaded gene candidates (Barshis et al. 2013), we filtered for genes that were at least 4-fold differentially expressed at 30°C between ICN and AF corals or their symbionts and that had expression values above an accumulative sum of 30 (TPM). This set of genes was subsequently selected for those that showed significant differential expression between 30°C and 33°C and between 30°C and 36°C in ICN or AF corals and their symbionts, respectively. We also conducted the afore-mentioned analyses between corals from the ExT and PrT sites to test for regional front-/back-loading. Moreover, we identified the most variably expressed genes across all temperatures and reef sites to elucidate the cohort of genes that are most responsive to thermal stress.

Symbiodiniaceae ITS2 and bacterial 16S marker gene library generation

DNA isolation was performed using the Qiagen DNeasy 96 Blood & Tissue kit (Qiagen, Hilden, Germany) following the manufacturer's instructions with minor adjustments. Briefly, coral tissue samples aliquoted for DNA isolation (see above) were thawed and equilibrated to room temperature. The slurry was vortexed and 180 µl of each sample were transferred to a microtube with 20 µl of proteinase K for incubation at 56°C for 1 hour. DNA extractions were then continued according to manufacturer’s instructions. In addition to the coral samples selected for 16S marker gene sequencing, 3 negative controls testing for ‘contaminants’ from (i) the FSW used for spraying off the coral tissue, (ii) the DNA extraction procedure, and (iii) PCR consumables were included to account for putative lab and kit contaminants introduced during sample preparation. DNA concentrations were quantified by Qubit (Qubit dsDNA High Sensitivity Assay Kit, Invitrogen, Carlsbad, USA). ITS2 (Symbiodiniaceae) and 16S (bacteria) amplicon libraries were prepared for sequencing on the Illumina MiSeq platform. To amplify the ITS2 region, the primers SYM_VAR_5.8S2 [5’-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGGAATTGCAGAACTCCGTGAACC-3’] and SYM_VAR_REV [5’-GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGCGGGTTCWCTTGTYTGACTTCATGC-3’] (Hume et al. 2013, 2015, 2018) were used (Illumina adaptor overhangs underlined). To amplify the variable regions 5 and 6 of the 16S rRNA gene, the primers 784F [5′-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGAGGATTAGATACCCTGGTA-3′] and 1061R [5′-GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGCRRCACGAGCTGACGAC-3′] (Andersson et al. 2008; Bayer et al. 2013) were used (Illumina adaptor overhangs underlined). Triplicate PCRs were performed (10-50 ng of DNA from each coral sample per PCR reaction) using the Qiagen Multiplex PCR kit and a final primer concentration of 0.5 μM in a reaction volume of 10 μl. Thermal cycling conditions for ITS2 PCR amplifications were: 95°C for 15 min, 30 cycles of 95°C for 30s, 56°C for 90s, and 72°C for 30s, followed by a final extension step of 72°C for 10min. Thermal cycling conditions for 16S PCR amplifications were: 95°C for 15min, 27 cycles of 95°C for 30s, 55°C for 90s, and 72°C for 30s, followed by a final extension step of 72°C for 10 min. After the PCRs, 5 µl of the PCR products were run on a 1% agarose gel to confirm successful amplification. Triplicate PCRs for each sample were pooled, and samples were cleaned using ExoProStar 1-step (GE Healthcare, Chicago, USA). Samples were indexed using the Nextera XT Index Kit v2. Successful addition of indexes was confirmed by comparing the length of the initial PCR product to the corresponding indexed sample on a 1% agarose gel. Samples were then cleaned and normalized using the SequalPrep Normalization Plate Kit (Invitrogen, Carlsbad, USA). The ITS2 and 16S amplicon libraries were pooled separately (4μl per sample) and concentrated using a CentriVap Benchtop Vacuum Concentrator (Labconco, Kansas City, USA). Quality of the libraries was assessed using the Agilent High Sensitivity DNA Kit on the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, USA), quantification was done using the Qubit dsDNA High Sensitivity Assay Kit (Invitrogen, Carlsbad, USA). The ITS2 and 16S amplicon libraries were sequenced on a single MiSeq run. Libraries were pooled in a ⅓ ITS2 to ⅔ 16S ratio to achieve higher read numbers for 16S amplicons. Sequencing was performed at 6 pM with 20% phiX on the Illumina MiSeq platform at 2 x 301bp paired-end V3 chemistry according to the manufacturer’s specifications. Sequence data determined in this study are available under NCBI BioProject PRJNA681108 (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA681108).

Symbiodiniaceae Community analysis

The SymPortal (symportal.org) analytical framework was used to analyze Symbiodiniaceae ITS2 sequence data (Hume et al. 2019). Briefly, demultiplexed and paired forward and reverse fastq.gz files outputted from the Illumina sequencing were submitted directly to SymPortal. Firstly, sequence quality control was conducted as part of the SymPortal pipeline using mothur 1.43.0 (Schloss et al. 2009), the BLAST+ suite of executables (Camacho et al. 2009), and minimum entropy decomposition (MED) (Eren et al. 2015). Then, ITS2 type profiles (representative of putative Symbiodiniaceae taxa or genotypes) were predicted and characterized by specific sets of defining intragenomic ITS2 sequence variants (DIVs). Finally, the ITS2 sequence and ITS2 type profile abundance count tables (Data S3), as well as the Bray Curtis-based between-sample and between-ITS2 type profile dissimilarities, as output by the analysis, were directly used to plot data. Scripts used for data curation and plotting are available at GitHub (https://github.com/reefgenomics/CBASS84).

Bacterial Community analysis

Sequence data fastq files were analyzed in mothur 1.39.5 (Schloss et al. 2009). Briefly, remaining adaptors and primers were trimmed and read pairs were concatenated. Unique sequences were removed using split.abund, followed by alignment to the 16 rRNA SILVA database, release 132 (Quast et al. 2013). Chimeras were identified using VSEARCH and removed for subsequent analysis. Clustering was done using the OptiClust algorithm and taxonomic annotation was done using Greengenes (release May 2013) and SILVA (release 138) 16S rRNA databases. Sequences annotated as chloroplast, mitochondria, unknown, or archaea were discarded. OTU abundance table, taxonomy, and fasta files (Data S4) generated by the make.shared, classify.otu, and get.oturep commands in mothur respectively, were used as input files for statistical analyses in R. OTUs were considered putative contaminants if their relative abundance in negative controls was higher than 10% in relation to coral samples. Barplots showing the relative abundances of the most abundant 20 families were plotted using ggplot2 v2.3.1.0 (Wickham 2011). Ordination plots and Bray-Curtis dissimilarities were calculated in Phyloseq v1.26.1 (McMurdie and Holmes 2013). The same analyses were carried out in ASV space with resembling results, however, due to backward compatibility to reference OTU sequences from previous work (in particular from the family Endozoicomonadaceae, (Neave et al. 2017)), we opted for OTU-based analysis.

Spearman rank correlation between coral holobiont compartments

To quantify the similarity across holobiont compartments of S. pistillata, we determined Spearman’s rank correlation coefficients across samples and temperatures for a given site (+1 = perfect correlation, 0 = no correlation). For the transcriptomes of coral host and microalgal symbionts, transcripts with low TPM expression and a variance of zero across all samples were removed from the analysis to reduce noise, yielding 19,154 coral genes and 29,412 microalgal symbiont genes. Similarly, to quantify the similarity between bacterial microbiomes, we calculated Spearman’s rank correlation coefficients based on log10(x+1)-converted OTU relative abundances. To reduce noise, all low abundant bacteria (25% percentile) were removed, yielding 1,660 out of 2,264 OTUs considered.

Differential thermal tolerance thresholds between corals from the northern and central Red Sea

To assess the consistency of the exceptional thermal tolerance of corals from the Red Sea and the underlying holobiont patterns, we investigated the response to heat stress in Red Sea Stylophora pistillata corals. Using a set of small aquaria and the Coral Bleaching Automated Stress System (CBASS) (Voolstra et al. 2020), we determined thermal tolerance thresholds by running a series of identical short-term acute heat stress assays (Voolstra et al. 2020). Thermal thresholds were calculated as the temperature point at which a loss of 50% photosynthetic efficiency was inferred, henceforth referred to as the effective dose 50 (ED50) (Evensen et al. 2020) to provide an empirical, comparable, and standardized proxy for coral thermotolerance (similar to the LD50 used in toxicological studies), based on the notion that increased/decreased loss of photosynthetic efficiency with temperatures reflects differences in coral bleaching susceptibility (Voolstra et al. 2020). We determined thermal thresholds at two spatial scales: large-scale (~ 1,000 km overwater distance) between corals from the northern Red Sea at IUI Coral Nursery (ICN) in the Gulf of Aqaba (GoA) and corals from the central Red Sea (CRS) at Al Fahal reef (AF) and the exposed (ExT) and protected (PrT) sites of Tahala reef, and small-scale (0.3–10 km overwater distance) between corals from the AF, ExT, and PrT CRS sites (Fig. 1A). The thermal tolerance threshold of GoA corals was significantly different from all three CRS sites (all P < 0.0001, Fig. 1B). Notably, the ED50 thermal limit for GoA corals, deemed to be among the most temperature resilient corals in the world (Fine et al. 2013; Osman et al. 2018), was on average > 3 °C below that of CRS corals (ED50 GoA = 34.08 °C ± 0.30 °C vs. ED50 CRS = 37.54 °C ± 0.30 °C (mean ± SE)). However, corals from both the GoA and the CRS exhibited a temperature tolerance threshold of about + 7 °C above their regional maximum monthly mean (MMM) temperature (MMM GoA: 27.01 °C; MMM AF: 30.76 °C; MMM ExT: 30.74 °C; MMM PrT 30.75 °C). Among the CRS reef sites, corals from the protected site of Tahala reef exhibited the highest ED50 thermal threshold (ED50 PrT 38.27 °C vs. ED50 ExT 37.56 °C vs. ED50 AF 36.80 °C) and were significantly different from corals at the AF site (P < 0.05, Fig. 1B). Apart from the differences in ED50 thermal thresholds, and conversely, the similarities found with respect to relative thermal limits, the distinct shapes of the temperature response curves of corals from the GoA and CRS are noteworthy (Fig. 1B). While GoA corals exhibit a continuous loss over increasing temperatures (Fig. 1B blue curve), all CRS corals retained their F_v/F_m values up to a critical threshold above which they dropped drastically (Fig. 1 yellow, pink, red curves).

Disparate patterns of gene expression in coral hosts and algal symbionts from the Gulf of Aqaba and central Red Sea under heat stress

To elucidate the molecular response patterns underlying the differences in temperature tolerance thresholds, we analyzed gene expression of the coral host (n = 24,971 genes) and algal symbionts (n = 39,795 genes) at the control/baseline (30 °C) and 33 °C and 36 °C heat stress temperatures (i.e., 7 colonies/site across 4 sites across 3 temperatures = 84 samples) (Data S1, Data S2). Coral samples from the GoA and CRS sites were separated by region. Closer inspection revealed that corals from both regions differed with regard to response to temperature treatments: corals from the GoA displayed discrete expression responses over temperatures, as shown by the distinct colony clustering (Fig. 2A) and the high number of differentially expressed genes (DEGs) between temperature comparisons (GoA site ICN, 30 °C vs 33 °C: 105 DEGs; 30 °C vs 36 °C: 2,766 DEGs; 33 °C vs 36 °C: 1,903 DEGs; Fig. 2C). In contrast, corals from the CRS did not exhibit a clear temperature response, as indicated by the clustering of all CRS corals across temperatures (Fig. 2A) and the almost complete absence of DEGs between temperatures for any of the CRS sites (CRS sites AF/ExT/PrT, 30 °C vs 33 °C: 0/0/0 DEGs; 30 °C vs 36 °C: 14/16/4 DEGs; 33 °C vs 36 °C: 22/12/4 DEGs; Fig. 2C). Similarly, algal symbiont gene expression clustered foremost by region with some differentiation across temperatures for algal symbionts from GoA corals and an absence thereof in CRS corals (Fig. 2B). As observed for the coral host gene expression, we found a higher number of DEGs across temperatures for GoA algal symbionts versus a largely mute response in CRS Symbiodiniaceae (Fig. 2D).

To further explore the pattern of gene expression responsiveness in GoA corals vs. the transcriptional stasis of CRS corals, we looked for genes that exhibited front-/back-loading, suggested as a mechanism that conveys increased thermal tolerance (Barshis et al. 2013). Characterization of DEGs that had at least 4-fold higher expression at the control/baseline temperature (30 °C) in corals from the AF site in comparison to the ICN site and that also exhibited differential expression between 30 °C and 33 °C and between 30 °C and 36 °C within the ICN site were considered front-/back-loaded in CRS corals and inducible in GoA corals. We determined 26 genes exhibiting front-loading and 1 gene showing back-loading in AF corals, as well as a further 33 and 88 front- and back-loaded genes respectively in AF coral-associated Symbiodiniaceae (Fig. 2E, Table S2). Conversely, and importantly, we could not detect any gene in GoA corals or their algal symbionts that exhibited front-loading (and only 1 gene exhibiting back-loading in AF corals); further, no front-/back-loaded genes were identified between PrT and ExT sites for corals and their symbionts. This supports the notion that only CRS corals exhibit fixed high/low expression of genes. By contrast, these genes show consecutive induction/suppression across temperatures in GoA corals. This, in turn, suggests distinct mechanisms of thermal tolerance in S. pistillata across Red Sea regions, i.e., even for the same ocean basin and species. Conversely, mechanisms appear similar across sites, i.e. within a region.

Closer examination of the front-loaded genes (Table S2) revealed three S. pistillata homologs of matrix metalloproteinases (Spis5851, Spis5854, Spis5860), previously identified among the set of front-loaded genes between temperature-tolerant and -sensitive Acropora hyacinthus in American Samoa (Barshis et al. 2013). Additionally, two Pax3 transcription factors and two homologs of TNF receptor-associated factor 3 putatively involved in apoptosis and activation of the immune response were identified (Margue et al. 2000; Häcker et al. 2011). Notably, we could not detect any heat shock proteins (hsps) among the front-loaded genes, converse to what was previously found (Barshis et al. 2013). Rather, heat shock proteins (Spis4493, Spis4494, Spis4499, Spis19922, Spis19923, Spis19925) were among a list of genes that were most temperature-responsive across sites in line with their ubiquitous induction upon exposure to stressful conditions found across animal and plant species (Whitley et al. 1999) (Fig. S1, Table S3). Among the front-loaded genes in algal symbionts, we found three genes with homology to a high-affinity nitrate transporter (Smic4635, Smic4659, Smic19959) and two oxidoreductases (Smic8517, Smic40364), implicated to play a role in thermal stress susceptibility (Baumgarten et al. 2013), in addition to a suite of metabolic genes (Table S3). Further, a heat shock protein homolog exhibited back-loading in the algal symbiont, i.e., lower expression at higher temperatures, for reasons that are currently unknown.

Thermally stable site-specific algal symbiont communities

We elucidated the fine-scale structure of the algal symbiont community of coral hosts from the GoA and CRS sites using the analytical framework SymPortal (Hume et al. 2019). Our analysis revealed unanimous association with microalgae of the genus Symbiodinium (LaJeunesse et al. 2018). Despite this broad-scale similarity, we found fine-scale genetic differentiation of associated Symbiodinium symbionts between sites that was stable across temperatures, as evidenced by the ITS2 type profiles (Fig. 3A). Clustering by sample similarity further revealed that algal symbionts from all CRS sites were close to each other (with some outliers) and distinct from GoA algal symbionts (Fig. 3B). The observed pattern of symbiont association indicates high host fidelity of algal symbionts across sites, as reported previously (Terraneo et al. 2019; Howells et al. 2020; Hume et al. 2020). Further to that, our results support fine-scale differentiation even between geographically close-by sites, as the central Red Sea sites (AF, ExT, ExP) all differed with regard to symbiont identities (Fig. 3). Notably, sample dissimilarity increased with geographic distance, although the distinct clustering of some microalgal symbionts from ExT corals may suggest that divergence aligns with prevailing environmental differences, rather than geographical distance (Fig. 3B).

Bacterial community dynamics resemble gene expression patterns of corals from the northern and central Red Sea

Bacterial community composition (n = 2,264 OTUs, Data S4) profoundly changed in GoA corals in contrast to stable bacterial community association of colonies across temperatures for CRS sites (all P_adj > 0.05, PERMANOVA, Fig. 4AC). This resembles the coral host gene expression pattern, in particular with regard to a ramping response, as evidenced by the increasing number of significant differentially abundant bacterial taxa between temperature comparisons in GoA corals and absence thereof in CRS corals (Fig. 4B). In line with this, all CRS sites were significantly different from GoA corals at all treatment temperatures (all P_adj < 0.05, PERMANOVA). Among CRS sites, bacterial community composition of AF corals was significantly different from ExT and PrT sites across all temperatures (P_adj < 0.05, PERMANOVA), whereas the latter two resembled each other (P_adj > 0.05, PERMANOVA).

Changes in the bacterial community structure in GoA corals followed a distinct pattern over the temperature increase: only 4 OTUs (Fig. 4B) were differentially abundant between 30 °C and 33 °C in GoA corals, whereas 80 and 42 OTUs were differentially abundant in the 30 °C vs. 36 °C and 33 °C vs. 36 °C comparisons, respectively (Fig. 4B). Despite the large number of differentially abundant bacterial taxa, the overall pattern was driven by a relative loss of Endozoicimonaceae (2 OTUs, mean log2 FC = -7.31) and a relative increase of Vibrionaceae (15 OTUs, mean log2 = + 9.77) and Rhodobacteraceae (21 OTUs, mean log2 = + 6.46), suggesting dysbiosis and loss of microbiome structure (Fig. 4C) (Voolstra and Ziegler 2020). S. pistillata colonies from CRS sites by comparison largely retained association with Endozoicomonas. In this context, bacteria of the genus Endozoicomonas are typically ascribed to be abundantly associated with S. pistillata in the Red Sea and globally (Bayer et al. 2013; Neave et al. 2017). However, Endoizoicomonas were only dominantly associated with corals from the ICN and the PrT sites, but comprised a rather small relative portion of the bacterial community in corals from the AF and ExT sites (Data S4). This suggests fine-scale structuring of bacterial community composition across spatially adjacent reefs (Roder et al. 2015; Ziegler et al. 2017b, 2019a). In light of the known thermotolerance differences between the ExT and PrT reef sites (Pineda et al. 2013; Voolstra et al. 2020), it is tempting to speculate that abundance of Endozoicomonas may contribute to the observed differences in bleaching susceptibility, although prior studies found that their relative abundance was not correlated with bleaching (Pogoreutz et al. 2018; Shiu et al. 2020). Corals from the PrT site were further characterized by a relative prevalence of Simkaniaceaea, which comprise obligate intracellular bacteria of eukaryotes (Everett 2014) and were recently shown as a hallmark of corals from pristine environments (Ziegler et al. 2019a), though their function is at present unknown. Notably, the nature of the data does not allow us to draw conclusions on bacterial biomass, i.e. increase/decrease of bacterial cell densities. It is desirable to understand bacterial biomass dynamics in future experiments to better assert whether bacteria are lost from the holobiont with increasing temperature stress (similar to algal symbiont loss in bleaching) or whether we indeed witness bacterial community restructuring.

The need for a standardized experimental framework to further our understanding of the factors underlying thermal tolerance

Despite global coral reef loss driven by ocean warming that triggers coral bleaching (Hughes et al. 2018), coral reef regions (Guest et al. 2012; Fine et al. 2013; Osman et al. 2018), populations (Barshis et al. 2013; Palumbi et al. 2014; Voolstra et al. 2020), and individual coral genotypes (Palumbi et al. 2014; Dixon et al. 2015) with enhanced bleaching resilience exist. However, a detailed understanding of the underlying factors and molecular underpinnings that drive increased thermal tolerance remains elusive. Here, we conducted identical acute thermal stress assays across four reef sites (totaling ~ 1,000 km overwater distance) to assess thermal limits of corals from the northernmost Red Sea in the Gulf of Aqaba and the central Red Sea, employing a standardized experimental framework and subsequent molecular characterization of the coral holobiont heat stress response. Counter to previous work suggesting that the extraordinary thermal tolerance from corals of the Gulf of Aqaba might equal that of southern/central Red Sea populations (Fine et al. 2013; Osman et al. 2018), we could demonstrate that while corals throughout the Red Sea harbor consistent thermal tolerance thresholds exceeding + 7 °C above their regional maximum monthly mean (MMM) temperature, absolute thermal limits differed according to the prevailing MMM temperature. Although it is unclear at present how ED50-derived thermal thresholds relate to natural bleaching thresholds, our determined thermal limits correspond very well with the exceptional thermal tolerance of S. pistillata from the Gulf of Aqaba that was found to sustain temperatures of + 5 °C above MMM for a period of 2 weeks, equivalent to 9.4 DHW (Bellworthy and Fine 2017; Krueger et al. 2017). Irrespective of how exactly ED50s relate to DHWs (or other means of determining bleaching thresholds), our work and previous work supports the notion that higher thermal tolerance in short-term heat stress assays are representative of increased thermal stress resilience in natural settings with ED50s serving a standardized proxy for bleaching susceptibility (Ziegler et al. 2017b; Morikawa and Palumbi 2019; Evensen et al. 2020; Voolstra et al. 2020). Importantly, we found different thermal thresholds across different sites at both large and small scales, which were reflected by the holobiont molecular response differences. This highlights the utility of employing mobile short-term heat stress assays for standardized phenotyping to detect thermal tolerance differences with subsequent investigation of the molecular response to obtain a better understanding of the factors shaping and the mechanisms underlying thermal tolerance.

Signatures of thermal tolerance across regions and coral holobiont compartments

The different thermal thresholds seen in the F_v/F_m data were reflected in the gene expression responses of corals and their symbionts as well as in the composition of the algal symbiont and bacterial communities. Transcriptional response patterns differed substantially between GoA and CRS corals with a strong response to the thermal stress in GoA corals vs. a muted/static response in CRS corals (up to 2,483 vs. 20 differentially expressed genes between temperature treatments, respectively). Notably, transcriptional differences in GoA corals were minor between the 30 °C control and the 33 °C heat stress temperature (95 DEGs), corroborating the notion that 33 °C is below the critical thermal threshold of GoA corals (Bellworthy and Fine 2017; Krueger et al. 2017; Evensen et al. 2020). Conversely, differences were most pronounced between 30 °C and 36 °C, but we also found substantial gene expression differences between 33 °C and 36 °C. Overall, our data indicate the ramping up of a gene expression response with increasing heat stress in GoA corals in stark contrast to the largely muted response in CRS corals. The ramping up vs. mute response of GoA vs. CRS corals is further highlighted by the identification of front- and back-loaded genes in CRS corals and the complete lack thereof in GoA corals. Notably, a subset of the front-loaded genes we identified in CRS corals were previously identified in Acropora hyacinthus corals from American Samoa (Barshis et al. 2013). On the one hand, this suggests the conservation of thermal tolerance genes across large geographic distances and coral species, and consequently, universal molecular signatures of thermal tolerance exist. On the other hand, corals can show strikingly disparate heat stress response patterns, even for coral colonies/populations of the same species of the same ocean basin, as found here; the latter suggests distinct mechanisms of thermal tolerance. To further understand the dynamics of such response patterns and their utility as indicators of differential thermal tolerance, determination of expression patterns over the course of the heat stress (e.g., before, during, after heat stress, and after the recovery period) may be critical. Similarly, reciprocal transplant or common garden follow-up experiments can further evidence whether the differences observed here are genetically fixed or a consequence of acclimation (Palumbi et al. 2014; Thomas et al. 2018).

It is also striking that the coral host gene expression patterns were mirrored in the algal symbiont. At large, the level of responsiveness in algal symbionts from GoA corals was unexpected, given that the common notion is that of a paucity of regulation at the transcriptional level for Symbiodiniaceae (Bayer et al. 2012; Baumgarten et al. 2013), signified by fixed expression differences between symbiont genera and species (Barshis et al. 2014; Parkinson et al. 2016). Notably, corals from all sites associated with distinct and largely consistent microalgal symbiont genotypes, suggesting host fidelity and local adaptation (Howells et al. 2020; Hume et al. 2020). Despite the differences across sites, the consistent gene expression stasis of algal symbionts across CRS corals and the disparate responsiveness in Symbiodiniaceae from GoA corals potentially suggests that geographically more distant symbionts are also more differentiated.

Bacterial community composition of GoA corals changed dramatically under heat stress, in contrast to CRS corals that displayed rather consistent assemblages, corroborating previous work that thermally more susceptible corals have less stable microbiomes (Ziegler et al. 2017b). At present, however, any presumptive functional consequences of these differences for coral holobiont physiology are unknown, despite it becoming broadly accepted that bacteria contribute to metaorganism biology and fitness (McFall-Ngai et al. 2013; Bang et al. 2018; Rosado et al. 2019; Pogoreutz et al. 2020; Voolstra and Ziegler 2020). In cases of rapid environmental change, the role of the microbiome in supporting adaptation of the metaorganism may become even more pronounced, given that microbial-mediated change can be brought about through association with different microbes or exchange/incorporation of new genetic material (e.g., via horizontal gene transfer) and does not rely on mutation/recombination (Cárdenas et al. 2020; Voolstra and Ziegler 2020).

Contrasting heat stress response patterns of coral holobionts suggest distinct mechanisms of thermal tolerance

Our analyses indicate distinct patterns in the response to heat stress of corals from geographically distant sites in the Red Sea: corals from the CRS sites show coral-algal transcriptional and bacterial community stasis vs. the responsiveness across holobiont compartments in corals from the GoA site. To differentiate patterns of resilience and resistance we assessed the change in correlation coefficients among samples for the three holobiont compartments (i.e., host gene expression, algal gene expression, bacterial community assemblage).

All CRS sites (AF, ExT, PrT) showed a high and stable correlation across temperatures for the coral host compartment (Fig. 5). By contrast, we found a clear loss of correlation in corals from the GoA at 36 °C suggesting (systemic) collapse, possibly due to crossing the temperature tipping point of the system (as corroborated by the predicted ED50 of ~ 34 °C). Similarly, algal symbionts from CRS sites showed consistently high correlation, while they exhibited a low correlation in GoA corals, irrespective of the temperature. Thus, besides the difference in symbiont genotypes associated per se, GoA coral algal symbionts may respond to temperature stress ‘early’ (30 °C exceeds the MMM of the GoA), and consequently, exhibit variability across treatment temperatures (i.e., at 30 °C, 33 °C, 36 °C). The bacterial assembly also exhibited a highly consistent, albeit lower, correlation across corals from all CRS sites, reflecting their stable community composition across temperatures. The overall lower correlation may be attributable to site-specific bacterial community composition differences. Corals from the GoA, however, exhibited a low correlation between replicate colonies at 30° and 33 °C that increased considerably at 36 °C, due to the substantial relative increase of Vibrionaceae and Rhodobacteraceae, driving the similarity between samples and indicating dysbiosis/loss of microbiome structure, as alluded earlier (Sweet and Bulling 2017; Boilard et al. 2020; Lima et al. 2020). Thus, the decreasing correlation coefficient(s) with increasing temperature observed in GoA corals is that of a population approaching a tipping point with subsequent systemic collapse: GoA corals ramp up their response to the increasing stressor until they reach a tipping point upon which the system collapses. As such, they presumably have a sensing feedback loop that adjusts the response based on the environmental input. Conversely, corals from the CRS exhibit a fixed response, irrespective of the strength of stress: CRS corals are not able to adjust their response, but rather employ a fixed response, in line with the front-loading of genes, that cannot be modulated (any further).

We posit that such differences may be reflective of the nature of the underlying differences in thermal tolerance: we interpret the GoA coral response as that of a resilient population responding to the strength of a stressor until reaching a tipping point upon which the holobiont collapses (loss of correlation). In contrast, the static pattern of consistently elevated but fixed expression states is suggestive of thermal resistance in CRS corals. Clearly both mechanisms evolved in their respective environments, and it is not clear whether they present different solutions to the same problem, or one may prove more successful than the other in countering ocean warming. Intuitively, a response-resilience mechanism seems more adaptable to changing environments than a static-resistance mechanism that makes coral holobionts operate at maximum capacity irrespective of the external input, as any additional stress cannot lead to a change in response. Further insight may come from the photosynthetic efficiency response curves, which reflect the molecular patterns observed: CRS corals exhibit stable F_v/F_m values over the experimental temperatures studied (i.e., 30 °C, 33 °C, 36 °C). In contrast, GoA corals exhibit a (slow and) increasing F_v/F_m loss with temperature. Pending whether such differences in mechanisms are predictive of susceptibility to future ocean warming, they may help in targeting coral colonies, species, and reefs for conservation and restoration efforts. Importantly, our work supports the existence of common gene ‘signatures’ (e.g., high expression of metalloproteinases under heat stress) irrespective of the underlying thermal tolerance mechanism. Such putatively universal signatures may become important in the screening for and engineering of corals with superior thermal resilience/resistance.

Acknowledgements

Research reported in this publication was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) project number 433042944 to CRV. We further acknowledge funding by the Paul G. Allen Family Foundation (PGAFF, Vulcan) and a BiNational Science Foundation grant (#2016403) to DJB and MF.

Author contributions

Conceived and designed the experiments: CRV, DJB

Collected/processed samples: CRV, GP, AC, CB-L, GB-P, JB, MF, DJB

Generated data: CRV, JJV, ST, AC, BCCH, GP, CB-L, KR, GB-P, JB, MF, DJB

Generated figures: CRV, JJV, AC, BCCH, DJB

Analyzed data: CRV, JJV, ST, AC, BCCH, DJB

Interpreted data: CRV, JJV, ST, DJB

Provided tools/reagents/methods: CRV, JJV, ST, AC, BCCH, MVO, NSB, SP, MF,SF-T, DJB

Wrote the manuscript: CRV, JJV, DJB

All authors approved the final manuscript.

Competing Interests statement

The authors declare no conflict of interest exists.