Coral microbiome diversity reflects mass coral bleaching susceptibility during the 2016 El Niño heat wave

Abstract

Repeat marine heat wave‐induced mass coral bleaching has decimated reefs in Seychelles for 35 years, but how coral‐associated microbial diversity (microalgal endosymbionts of the family Symbiodiniaceae and bacterial communities) potentially underpins broad‐scale bleaching dynamics remains unknown. We assessed microbiome composition during the 2016 heat wave peak at two contrasting reef sites (clear vs. turbid) in Seychelles, for key coral species considered bleaching sensitive (Acropora muricata, Acropora gemmifera) or tolerant (Porites lutea, Coelastrea aspera). For all species and sites, we sampled bleached versus unbleached colonies to examine how microbiomes align with heat stress susceptibility. Over 30% of all corals bleached in 2016, half of which were from Acropora sp. and Pocillopora sp. mass bleaching that largely transitioned to mortality by 2017. Symbiodiniaceae ITS2‐sequencing revealed that the two Acropora sp. and P. lutea generally associated with C3z/C3 and C15 types, respectively, whereas C. aspera exhibited a plastic association with multiple D types and two C3z types. 16S rRNA gene sequencing revealed that bacterial communities were coral host‐specific, largely through differences in the most abundant families, Hahellaceae (comprising Endozoicomonas), Rhodospirillaceae, and Rhodobacteraceae. Both Acropora sp. exhibited lower bacterial diversity, species richness, and community evenness compared to more bleaching‐resistant P. lutea and C. aspera. Different bleaching susceptibility among coral species was thus consistent with distinct microbiome community profiles. These profiles were conserved across bleached and unbleached colonies of all coral species. As this pattern could also reflect a parallel response of the microbiome to environmental changes, the detailed functional associations will need to be determined in future studies. Further understanding such microbiome‐environmental interactions is likely critical to target more effective management within oceanically isolated reefs of Seychelles.

INTRODUCTION

Coral reef ecosystems are exceptionally vulnerable to anthropogenic disturbance and have been decimated by climate change‐driven marine heat waves during 2015–2017, with >30% of all corals lost at many locations worldwide through bleaching (Hughes et al., ). Western Indian Ocean (WIO) reefs were particularly affected by Sea Surface Temperature (SST) anomalies throughout 2016 broadly exceeding 10–15 maximum degree heating weeks (DHW), driving severe bleaching and mortality throughout this region (Hughes et al., ). Reefs within the WIO have in fact been repeatedly impacted by heat waves throughout the last 20 years (Graham, Jennings, MacNeil, Mouillot, & Wilson, ; Hughes et al., ; McClanahan, Ateweberhan, Darling, Graham, & Muthiga, ), and the combination of smaller scale thermal anomalies with other stressors has increasingly limited long‐term coral recovery (e.g., Graham et al., ; Zinke et al., ).

Seychelles coral reefs were among the most impacted globally during the 1998 mass bleaching, with coral cover reduced by >90% across the inner islands (Graham et al., ; Wilson et al., ). Recovery has been limited by strong phase shifts toward algal dominance in several reefs (Graham et al., ; Wilson et al., ), and recruitment bottlenecks (Chong‐Seng, Graham, & Pratchett, ), such that coral cover and diversity in Seychelles is now lower than for many other regions in the WIO (Harris, Wilson, Graham, & Sheppard, ). Reefs are characterized by patch, granitic, and carbonaceous habitats across coastal fringing and oceanic environments (Graham et al., ; Jennings, Grandcourt, & Polunin, ). However, loss of coral cover and diversity since 1998 has been greatest for the carbonaceous reefs (Graham et al., ; Wilson et al., ), resulting in reduced carbonate budgets, accretion potential, and structural maintenance (Januchowski‐Hartley, Graham, Wilson, Jennings, & Perry, ) that underpin the critical ecosystem service value of Seychelles’ reefs (Clifton et al., ). Existence of functional diversity at macroecological scales of lower latitude reefs in the WIO, including Seychelles, has recently been identified as a key determinant of longer‐term resilience to environmental stress (Zinke et al., ). However, such a role for microecological processes in the WIO remains largely unknown, in particular microbial community composition and functioning that can be critical in determining coral stress resilience in other reef regions worldwide (Putnam, Barott, Ainsworth, & Gates, ; Suggett, Warner, & Leggat, ).

Various components of the coral holobiont, that is, the cnidarian host and its microbial community (the “microbiome”), play key roles in regulating corals’ resistance to heat stress (Morrow, Muller, & Lesser, ; Sweet & Bulling, ). Changes in phylogenetic diversity of the endosymbiont (family Symbiodiniaceae, formerly genus Symbiodinium; LaJeunesse et al., ) are well described to correspond with bleaching susceptibility (Suggett et al., ), whereby corals often associate with known heat‐tolerant Symbiodiniaceae taxa (e.g., Howells et al., ; Hume et al., ; Pettay, Wham, Smith, Iglesias‐Prieto, & LaJeunesse, ; Ziegler, Eguíluz, Duarte, & Voolstra, ). However, the coral holobiont also hosts a highly diverse bacterial community (Bourne & Munn, ; Frias‐Lopez, Zerkle, Bonheyo, & Fouke, ; Rohwer, Breitbart, Jara, Azam, & Knowlton, ; Rohwer, Seguritan, Azam, & Knowlton, ), which in some cases form species‐specific associations with corals (Neave, Rachmawati, et al., ; Rohwer et al., ). Microbial communities may facilitate acclimatization of the coral holobiont to changes in the environment through rapid restructuring of the community (Reshef, Koren, Loya, Zilber‐Rosenberg, & Rosenberg, ; Torda et al., ; Ziegler, Seneca, Yum, Palumbi, & Voolstra, ), and studies indicate that an intact (Krediet, Ritchie, Paul, & Teplitski, ; Roder, Bayer, Aranda, Kruse, & Voolstra, ; Rosenberg, Koren, Reshef, Efrony, & Zilber‐Rosenberg, ) and diverse (Hadaidi et al., ) coral microbiome may be essential to coral immunity and health. Responses of coral‐associated bacterial communities to shifts in coral health (Bourne, Iida, Uthicke, & Smith‐Keune, ; Cárdenas, Rodriguez‐R, Pizarro, Cadavid, & Arévalo‐Ferro, ; Glasl, Herndl, & Frade, ; Guest et al., ; Jones, Berkelmans, Oppen, Mieog, & Sinclair, ; Ziegler, Seneca et al., ) and environmental stressors (e.g., Garren, Son, Tout, Seymour, & Stocker, ; Jessen et al., ; Kwiatkowski, Cox, Halloran, Mumby, & Wiltshire, ; Vega Thurber et al., ; Vega Thurber et al., ) have been extensively explored and reported, with recent evidence suggesting flexibility of these communities may determine holobiont resistance to environmental stress (Pogoreutz et al., ).

Microbiome composition is important in determining coral health over space and time, yet is completely undescribed for the coral communities in Seychelles. Therefore, as part of a long‐term program examining carbonaceous coral communities within Curieuse Marine National Park (Seychelles), we evaluated Symbiodiniaceae and bacterial diversity and community composition among key reef‐building coral species (Acropora gemmifera, Acropora muricata, Coelastrea aspera, and Porites lutea) and for populations within a clear water versus turbid reef environment. Importantly, this combination of species included those previously described as either heat stress sensitive (Acropora spp.) or tolerant (C. aspera, P. lutea) in the Seychelles (Harris et al., ; Smith, Wirshing, Baker, & Birkeland, ) and wider WIO (McClanahan et al., ). Microbial sampling coincided with the 2016 mass bleaching marine heat wave in April 2016 (Figure ) enabling us to uniquely differentiate microbiomes for bleached versus unbleached colonies within these two environments. In doing so, we show for the first time how microbial signatures vary during chronic heat exposure across complex species‐environment interactions in the WIO. Microbiomes were generally conserved among bleached versus unbleached colonies for the two species of Acropora sp. suggesting small‐scale environmental variability (e.g., shading from surrounding substrates, sensu Hoogenboom et al., ) are likely critical in ensuring persistence of heat stress‐sensitive corals in oceanically isolated reef systems such as Seychelles.

(a) Coral cover for East Bay (black symbols) and Praslin (light‐gray symbols), (b) sea surface temperatures (SST) between 1982 and 2017 for the Seychelles (dark blue symbols) showing the annual maximum degree heating weeks (DHW; red triangles). Dotted line indicates the average sea surface temperature recorded for Seychelles. Averages ±SE shown for coral cover (n = 3 for East Bay and 6 for Praslin). (c) Photographs showing coral cover and condition between 2015 and 2017 for Acropora sp.

MATERIALS AND METHODS

Site description and benthic sampling

Coral communities have been examined from two fringing carbonaceous reef sites within Curieuse Marine National Park (CMNP) since 2009 (Supporting Information Figure S1): A relatively clear water reef along the southern coast of Curieuse Island (East Bay; 4°16′55″S 55°44′32″E) and a relatively turbid reef on the northern coast of Praslin Island (Praslin; 4°18′35″S 55°43′28″E). Both sites are subject to a semidiurnal tidal cycle, with mean ± SE tidal range of 1.3 ± 0.2 m and assessed annually (2009–2017) via 2‐week sampling campaigns to coincide with the end of the northwest monsoon wind season when SST is warmest (April–May). In 2016, this campaign coincided with the peak of the prolonged global heatwave, where cumulative heat stress for Seychelles exceeded and 10 DHW (^oC‐weeks; https://coralreefwatch.noaa.gov/satellite/index.php; Mahe). Reefs within this region were on “Bleaching Watch” from November 2015 and “Alert Level 2” by April 2016.

Coral cover and species identification were determined for all years, and sampling for Symbiodiniaceae and bacterial community diversity was conducted in 2016. Coral cover data are presented as means ± standard error (SE).

Abiotic measurements

To assess the historical thermal stress specifically around CMNP, remotely sensed Sea Surface Temperatures (SST) and Degree Heating Weeks (DHWs) were extracted for Curieuse Island (lat −4.297764, lon 55.727171) for the last 35 years (1982–2017) from two high‐resolution sources: (a) weekly CoRTAD data (Pathfinder 5.2, AVHAR, 4 km resolution) (https://data.nodc.noaa.gov/cgi-bin/iso?xml:id=gov.noaa.nodc:0126774) for the period between 1982 and 2012 (Casey et al., ); and (b) daily Coral Reef Watch (CRW; 5 km resolution) data for the period from 2013 to 2017 (Liu et al., ). Data were extracted from single grid (pixel) avoiding land interference using ArcGIS 10.3.1 software. In addition, temperature was measured at each location over the two‐week sampling period using calibrated HOBO® pendant loggers set at 20‐s measuring intervals, to verify local conditions relative to the global SST retrieval, as well as contrast SSTs on the reef flat (1–2 m) versus upper slope (5–8 m). Consistent with the remotely sensed SSTs (Figure b), maximum daily SSTs throughout our 2016 sampling period exceeded 32ºC, ranging 32.3–32.9ºC and 31.8–32.3ºC at 1–2 m and 5–8 m, respectively. Light was measured at two depths (6 and 1 m) using intercalibrated (no breaking) HOBO® pendant light loggers in 2014–2015. Data were converted from Lux to µmol photons m⁻² s⁻¹ (as per Long, Rheuban, Berg, & Zieman, ) to return light attenuation coefficients (K_D (PAR), m⁻¹) for each site. As expected, K_D (PAR) was highly variable at both sites as a result of tidal resuspension of sediment, but the hourly integrated range for clear water East Bay (0.14–0.21 m⁻¹; n = 83) generally was half that for turbid Praslin (0.23–0.38 m⁻¹; n = 59). Previous work at these sites has demonstrated that pH, A_T, salinity, PCO₂, and nitrate concentrations are similar between locations for this time of year, ca. 8.12 (total scale), 2.360 μmol kg/SW, 35.5 ppm, 323 μatm, and 1.1 μM, respectively (for details, see Camp et al., ).

Characterizing benthic habitat

Standard continuous line intercept transects were conducted using SCUBA diving and were used to quantify coral community structure at these sites. Data were recorded by high‐definition video (Sony HDR‐SR1E, 2009; Canon G10, 2011–2013; GoPro Hero series from 2014–2017) and later analyzed to determine habitat type as percentage live coral cover to species level. Both sites were characterized by a shallow (1–3 m) reef flat that transitioned to a slope of semicontinuous consolidated carbonaceous reef to depths of 12–13 m (Praslin) or 10–14 m (East Bay). Importantly, these reef flats were very different habitat types at the two sites, specifically patch reef and sand at Praslin versus rubble at East Bay; thus, we restricted our analysis to the upper reef slope to standardize habitat type, performed as 3 × 30 m (2009, 2012, 2014–2017) continuous line intercept transects (5–8 m depth) at each site per year. All coral data from 2016 were further categorized to account for bleaching intensity. Our initial categorization adopted that of McClanahan et al. (), which assigns a score based on % live coral surface area fully bleached. However, since all observations during the mass bleaching were 80%–100% of the live coral surface area fully bleached, we eventually defined categories as unbleached, bleached (pale and with host tissue remaining), and recently dead (host tissue lost and/or first signs of algal overgrowth).

Microbiome sampling and DNA extraction

All colonies were sampled 6–7 days after strong bleaching was first observed. Triplicate fragments were taken from independent colonies for each of the four key coral species found throughout the upper reef slope at both turbid (Praslin) and clear water (East Bay) sites: A. gemmifera, A. muricata, P. lutea, and C. aspera (three fragments × four species × two sites × two coral condition = 48 fragments total). Replicate colonies of the same species were generally sampled >20 m apart, although unbleached colonies were often within close proximity (5–10 m) of neighboring bleached colonies. A small portion of each fragment was removed for Symbiodiniaceae cell density assessment using a haemocytometer and normalized to surface area (as per Camp et al., ). The remaining fragment was immediately preserved in RNAlater for bacterial 16S and Symbiodiniaceae ITS2 rRNA genotyping and stored at −80°C until processing.

For DNA extraction, coral fragments were thawed slowly on ice and removed from the RNAlater solution using sterile forceps and kimwipes to remove excess RNAlater (Tout et al., ; Vega Thurber et al., ). Coral fragments were transferred into sterile ziplock bags and mucus and tissues subsequently air‐blasted using airflow from a sterile pipette tip (1,000 µl filter barrier tips; Neptune, USA) into 5 ml PBS‐EDTA. Care was taken to exclude any skeletal fragments (to reduce contamination by the endolithic microbiome). DNA was extracted using the Qiagen DNeasy Plant Mini Kit (Qiagen, Hilden, Germany) according to manufacturer's instructions. To disrupt Symbiodiniaceae cells, 200 µl tissue slurry in 200 µl AP‐1, 4 µl RNase A stock solution (100 g/ml), and 200 µl 0.5 mm sterile glass beads (BioSpec, Bartlesville, OK, USA) were bead‐beaten at 30 Hz for 90 s with a Tissue Lyser II (Qiagen). An additional step using phenol:chloroform:isoamyl alcohol (25:24:1) prebuffered to pH 8 was used to purify the DNA. Mock samples (i.e., empty) were used for DNA extraction to account for kit contaminants (Pogoreutz et al., ; Salter et al. ). Extracted DNA was quantified and quality checked using a NanoDrop 2000C spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). PCR amplifications (including mock [i.e., no template] PCRs) to account for contaminants (Pogoreutz et al., ; Salter et al. ) were performed in triplicate reactions with Qiagen Multiplex PCR Kit (Qiagen; see below).

PCR amplification and sequencing

For Symbiodiniaceae typing, amplification of the ribosomal Internal Transcribed Spacer 2 (ITS2) region, a multicopy genetic marker commonly used to assess Symbiodiniaceae diversity was used (Arif et al., ; Smith, Ketchum, et al., ; Smith, Vaughan, Ketchum, McParland, & Burt, ). For amplicon‐specific PCRs, we used the primers ITSintfor2 5'‐TCG TCG GCA GCG TCA GAT GTG TAT AAG AGA CAG GAA TTG CAG AAC TCC GTG‐3' and ITS2‐reverse 5'‐GTC TCG TGG GCT CGG AGA TGT GTA TAA GAG ACA GGG GAT CCA TAT GCT TAA GTT CAG CGG GT‐3' with concatenated Illumina sequencing adapters (underlined). Initial PCR amplification was achieved with the following thermal cycles: 94°C for 15 min, then 35 cycles of 94°C for 30 s, 51°C for 30 s, 72°C for 30 s, followed by one cycle of 72°C for 10 min and 4°C hold (Coleman, Suarez, & Goff, ; LaJeunesse, ).

To amplify the bacterial 16S rRNA gene, we used the primers 16SMiSeqF‐Andersson 5‐TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGAGGATTAGATACCCTGGTA‐3′ and 16SMiSeqR‐Andersson 5′‐GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGCRRCACGAGCTGACGAC‐3´ (Illumina sequencing adapters underlined) that target the variable regions five and six of the 16S gene (Andersson et al., ). These primers have previously been reported to amplify well with coral DNA (Bayer et al., ; Röthig, Ochsenkühn, Roik, Merwe, & Voolstra, ). The thermal conditions for 16S amplicon PCRs were as follows: 95°C for 15 min, followed by 27 cycles of 95°C for 40 s, 55°C for 40 s, 72°C for 40 s, and a final extension cycle of 72°C at 10 min (Ziegler et al., ).

For individual PCR reactions, DNA was aliquoted to 12–50 ng/L, with 10 µl Qiagen Mix, 0.5 µl of each 10 M primer mix, 1 µl of DNA template, and RNAse‐free water to adjust the reaction volume to 20 µl. 10 µl of each PCR product was run on an 1% agarose gel to visualize successful amplification.

Sample triplicates were subsequently pooled and then purified using the Agencourt AMPure XP magnetic bead system (Beckman Coulter, Brea, CA, USA). Purified PCR products were subjected to an indexing PCR (eight cycles) to add Nextera XT indexing and sequencing adapters (Illumina) according to the manufacturer's protocol. Indexed amplicons were again purified, quantified on the QuBit (Quant‐IT dsDNA Broad Range Assay Kit; Invitrogen, Carlsbad, CA, USA) and pooled in equimolar ratios on the BioAnalyzer (Agilent Technologies, Santa Clara, CA, USA). The final pooled library was purified on a 2% agarose gel to remove excess primer dimer. The library was sequenced at 8 p.m. with 10% phiX on the Illumina Miseq, 2× 300 bp end version 3 chemistry according to the manufacturer's specifications at the Bioscience Core Lab (KAUST, Saudi Arabia).

Microbial analysis—Symbiodiniaceae

The SymPortal analytical framework (symportal.org, github.com/SymPortal) was used to predict putative Symbiodiniaceae taxa. Briefly, the multicopy nature of the rRNA gene means that every Symbiodiniaceae genome contains hundreds to thousands of copies of it. Each of these gene copies is able to accrue mutations somewhat independently. As such, considerable intragenomic sequence diversity that may be leveraged for purposes of taxonomic delineation is found within every Symbiodiniaceae cell (Hume, D'Angelo, Burt, & Wiedenmann, ). SymPortal aims to make use of this diversity, using next‐generation amplicon sequencing data, to resolve between genetically differentiated taxa. SymPortal works by identifying specific sets of defining intragenomic ITS2 sequence variants (DIVs) that are used to define the taxonomic unit of SymPortal, the ITS2 type profile, indicative of genetically differentiated Symbiodiniaceae taxa. Demultiplexed and paired forward and reverse fastq.gz files outputted from the Illumina sequencing were submitted directly to SymPortal. Sequence quality control was conducted as part of the SymPortal pipeline using Mothur 1.39.5 (Schloss et al., ), the BLAST + suite of executables (Camacho et al., ), and minimum entropy decomposition (MED; Eren et al., ).

Over the years, development and variation in the range of techniques employed to genetically differentiate within the Symbiodiniaceae have led to a range of different terms being used to describe the genotypic units of resolution. For example, “type,” “ITS2 type,” “ITS2 profile,” “ITS2 fingerprint,” “clade,” “subclade,” “subtype,” and most recently “ITS2 type profile” are among the most common. Some of these are used interchangeably in one setting, while representing different entities in others. To clarify and limit ambiguity, for the purposes of this study, we will restrict our use to “type” and “ITS2 type profile.” A type refers to Symbiodiniacea taxa that have a specific sequence or set of sequences as their most abundant sequence. An ITS2 type profile is a set of sequences that are used to define either a putative or defined taxa. For example, Durusdinium trenchii is a D1 type and has an ITS2 type profile of D1–D4.

Microbial analysis—bacteria

Sequencing reads for 16S bacterial analysis were processed as outlined in Kahlke (). In summary, sequences were joined using FLASH (Magoč & Salzberg, ) and subsequently trimmed using MOTHUR (Parameters: minlength = 294, maxlength = 301, maxhomop = 6, maxambig = 0). The resulting fragments were clustered into operational taxonomic units (OTUs) and identity threshold of 97% and chimeric sequences were identified using VSEARCH (Rognes, Flouri, Nichols, Quince, & Mahé, ) and the SILVA v128 database. Unwanted sequences related to Archaea, mitochondria, chloroplasts, and kit contaminants (confirmed with sequencing of mock samples and PCRs; OTUs 2, 4, and 7; Brevibacter casei, Shigella flexneri, and Brachybacterium sp., respectively) were removed. To assign taxonomy, QIIME (Caporaso et al., ) was used with the BLAST algorithm against the SILVA v128 database at a 97% similarity cutoff. Data were rarefied to 1,000 sequences per sample and the subsequent.biom file was used for all downstream applications. Alpha diversity indices (Chao1 index, Simpson's evenness, Shannon diversity index, and Phylogenetic diversity) were computed using QIIME.

Data analysis

To initially examine trends of total live coral cover, we pooled data according to species of Acropora and Pocillopora versus “other species” since, as in previous bleaching episodes in Seychelles (e.g., see Wilson et al., ), these two genera are typically most susceptible to bleaching. Unpaired t‐tests were used to analyze the percent of coral cover data, using Welch's correction for variable standard deviations (GraphPad Prism v.6).

A univariate general linear model (GLM) for bacterial alpha diversity indices (including number of reads, OTUs per sample, phylogenetic diversity, Chao1, Simpson's evenness, and Shannon diversity index) was used to compare interactive effects of species, site, and coral condition and a one‐way ANOVA used to detect significant differences within species using IBM SPSS Statistics (v.21; IBM Corporation, New York, NY, USA). Coral‐associated bacterial assemblages were tested for differences between species (A. gemmifera, A. muricata, C. aspera, and P. lutea), sites (East Bay and Praslin), and coral condition (bleached and unbleached) using permutation multivariate analysis of variance (PERMANOVA) at family level. Here, all fixed factors (sites and coral condition) were nested according to hierarchy, and 9,999 permutations of residuals were conducted based on Bray–Curtis distances between fourth‐root transformed samples using the PRIMER‐E software with the PERMANOVA+ add‐on package v1.0.6 (Anderson, Gorley, & Clarke, ; Clarke & Gorley, ). Analysis of Similarity (ANOSIM) was conducted for species, site, and coral condition, and a significant difference was only detected for species. Similarity percentage analysis (SIMPER; using species as the factor) revealed the main contributing bacterial families responsible for differences using PRIMER v6.1.16 (Clarke & Gorley, ). Analyses were conducted at family level (i.e., at a taxonomic level where the majority of OTUs are assigned) to observe differences that were not be resolved at the strain level. Beta diversity differences for bacterial community composition were visualized in a principal coordinate analysis (PCoA) based on a Bray‐Curtis dissimilarity matrix and Pearson's correlation. In addition, analysis of bacterial community data excluding the most abundant bacterial families (Hahellaceae, Rhodospirillaceae, and Rhodobacteraceae) were run using PERMANOVA to determine whether any community composition changes were masked by the dominant taxa (see Supporting Information).

Symbiodiniaceae diversity analysis was conducted using the ITS2 type profile data output by SymPortal. We tested for significant individual and interactive effects of species, site, and coral condition using permutational multivariate analysis of variance (PERMANOVA using the PRIMER‐E software with the PERMANOVA + add‐on package v1.0.6 (Anderson et al., ; Clarke & Gorley, ). ANOSIM was performed for species, site, and coral condition and a significant difference was only detected between species. SIMPER revealed the main contributing Symbiodiniaceae ITS2 type profile. Beta diversity was visualized in a principal coordinate analysis (PCoA) based on a Bray‐Curtis dissimilarity matrix and Pearson's correlation.

RESULTS

Impact of the 2016 heat wave event on live coral cover

Prior to 2016, Seychelles reefs experienced few years with prominent annual maximum degree heating weeks DHWs (>4°C‐weeks), notably, 1998, 2003, 2010 (see Figure b). In 2016, DHWs were the highest recorded for the last 35 years, where the annual maximum DHW in CMNP reached 8.2°C‐weeks (Figure b). By 2017, much of the live coral cover that had bleached in 2016 had transitioned to dead eroding coral covered by algal mats (Figure c).

Total live coral cover has been steadily increasing within CMNP at both sites between 2009 and 2016 (from 35% ± 5% to 49% ± 9% and from 30% ± 2% to 42% ± 7%, in East Bay and Praslin, respectively; Figure a). However, following the 2016 heat wave, total live coral cover at East Bay significantly declined from 49% ± 9% to 22% ± 2% (t = 2.91, df = 4, p < 0.05) and at Praslin from 42% ± 7% to 13% ± 3% (t = 3.61, df = 10, p < 0.005) between 2016 and 2017 (Figure a and a). Approximately one‐third of the coral cover identified at each site was bleached (39% ± 7% at East Bay and 32% ± 6% at Praslin, Figure a) at the peak of the heat wave during our sampling campaign in April 2016, resulting in a loss in coral cover of 27% ± 7% (East Bay) and 29% ± 4% (Praslin; Figure b).

(a) Total coral cover for 2016 (dark blue bars) and 2017 (light blue bars; left y‐axis) and total percent of bleached coral in 2016 (red symbols; right y‐axis) and (b) total percent of coral lost between 2016 and 2017 (yellow bars) for East Bay (EB) and Praslin (P), (c) total coral cover for Acropora and Pocillopora and Other groups in 2016 (dark blue bars) and 2017 (light blue bars; left y‐axis) and (d) total percent of bleached coral in 2016 classified into Acropora and Pocillopora (dark green bars) and Other (light green bars; right y‐axis) groups for East Bay (EB) and Praslin (P). Averages ± SE are shown (n = 3 for East Bay and n = 6 for Praslin). Symbiodiniaceae cell density for unbleached (dark gray bars) and bleached (light gray bars) colonies of (e) Acropora gemmifera, (f) Acropora muricata, (g) Coelastrea aspera, and (h) Porites lutea at East Bay and Praslin. Asterisks (*) indicate significant differences between health, arrows (^) indicate a significant difference between site for the bleached colonies of A. muricata where p < 0.05. Data were log10 transformed for A. muricata and P. lutea and square‐root transformed for C. aspera for homogeneity of variance. Averages ± SE shown (n = 3, n = 2 for bleached A. muricata at Praslin)

Of the 49% total live cover at East Bay before the heat wave in 2016, almost half (21% ± 4%) was Acropora sp. and Pocillopora sp. with the remaining 28% ± 7% as other (Figure c; particularly Porites sp. but also Montipora sp. and Favia sp.). The large loss of coral cover by 2017 was attributed largely to mortality of Acropora sp. and Pocillopora sp., where live cover had declined to only 1% ± 1%, whereas all other taxa had declined to only 21% ± 3% (Figure c), despite no difference in the percent of bleached corals in 2016 for both sites (Figure d). Similar trends were observed at the more turbid Praslin site, where total live cover for Acropora sp. and Pocillopora sp. versus all other taxa was 16% ± 5% versus 26% ± 3% in 2016, declining to 0.9% ± 0.5% versus 12% ± 3% in 2017 (Figure c). As such, loss of taxa other than species of Acropora and Pocillopora was also important at this turbid site (and notably Montipora sp., Lobophyllia sp.).

Symbiodiniaceae communities

Symbiodiniaceae cell densities were variable across samples but overall lower for all bleached colonies (2 × 10⁵−1 × 10⁴ cm⁻²) compared to unbleached (8 × 10⁷−5 × 10⁵ cm⁻²) at both locations (see Figure e–h). Thus, substantial numbers of Symbiodiniaceae were still present in highly visibly bleached corals at the time of sampling in April 2016.

Symbiodiniaceae genera detected in 48 ITS2 gene libraries in all samples from the East Bay and Praslin sites were Cladocopium (formerly clade C; LaJeunesse et al., ) (13 different ITS2 type profiles, indicative of genetically differentiated Cladocopium spp. taxa) and Durusdinium (formerly clade D; 12 ITS2 type profiles, Figure a–d, Supporting Information Table S6). Among all hosts, the greatest Symbiodiniaceae diversity was observed for C. aspera (16 ITS2 type profiles identified) and the least for A. gemmifera (three identified). PERMANOVA revealed no significant interaction between species, site, and coral condition for all samples (pseudo F_3,47 = 1.46, p = 0.164), although a significant difference between species was detected (pseudo F_3,47 = 6.33, p < 0.001; Supporting Information Table S1a). Further analysis confirmed that no significant difference in Symbiodiniaceae ITS2 type profiles for site (ANOSIM; Global R = 0.01, p = 0.261) or coral condition (ANOSIM; Global R = 0.041, p = 0.073), and only significant differences between coral species detected (ANOSIM; Global R = 0.25, p < 0.0001; Supporting Information Table S1b), corresponding to distinct coral taxa specific responses: (a) A. gemmifera, generally showed little variance across sites or between coral condition (contribution from ITS2 type profiles of 51.96% from C3z‐C3‐C3.10‐C3bq and of 47.22% from C3z/C3‐C3.10‐C3an; Figure a, Supporting Information Table S2); (b) A. muricata, had different communities between unbleached sites (C3z/C3‐C3.10‐C3an for East Bay vs. D1‐D2.2‐D1m‐D4‐D2c dominant for Praslin) that then converged to a single “bleached” community of C3z/C3‐C3.10‐C3an and C3z‐C3‐C3.10‐C3bq for both sites (Figure b); (c) C. aspera, which had very different communities across site and bleached versus unbleached colonies, with the highest contribution from D1‐D4‐D4c‐D2‐D4f (28.62%) and 28.24% contribution from C3z/C3‐C3.10‐C3an (Figure c, Supporting Information Table S2) and finally (d) P. lutea, a single community of C15‐C15ad (91.61%) and C15 (5.75%) for both sites unbleached, but additions of D4 and D9 in the bleached samples (Figure d, Supporting Information Table S2). The defining intragenomic variant (DIV) counts for Symbiodiniaceae profiles show the breakdown of raw data in high resolution to predict taxa at each site and condition (Supporting Information Figure S2e–h). Principle Coordinate Analysis (PCoA) shows clear separation for the C3z/C3‐C3.10‐C3an, C3z‐C3‐C3.10‐C3bq and C15‐C15ad ITS2 type profiles along the axes where PCO1 explained 32.5% of the variation while Pco₂ explained 25.1% (Figure e).

Average relative abundance (%) showing the Symbiodiniaceae ITS2 type profile for bleached and unbleached corals at East Bay and Praslin for (a) Acropora gemmifera, (b) Acropora muricata, (c) Coelastrea aspera and (d) Porites lutea. Colors represent different Symbiodiniaceae ITS2 type profiles. (e) Principle Coordinate Analysis (PCoA) for dominant Symbiodiniaceae ITS2 type profiles found in four coral species at East Bay and Praslin in Seychelles. Data was fourth‐root transformed, and a Bray‐Curtis similarity matrix was used with a correlation of 0.2. Ellipses denote similarity clusters of 20% (green dashed line). Percentages on axes indicate variation explained by the two coordinates

Bacterial community structure

Overall, the data set comprised 43 16S rRNA gene libraries (three replicates × four coral species × two sites × two coral conditions) totaling 1,412,100 sequences with a mean length of 294 bp. After quality filtering and exclusion of chimeras, 1,065,414 sequences were annotated to bacteria. Clustering of these sequences at the 97% similarity level resulted in 2,362 OTUs (Supporting Information Table S7), presented as a taxonomy stacked column plot to the phylogenetic level of family (Supporting Information Figure S3). Significant interactions in bacterial diversity were found between species, site, and colony condition for the number of OTUs per sample, phylogenetic diversity and Chao1 (Univariate GLM; see Table ). There was a significant difference between coral species for the number of OTUs (One‐way ANOVA; F_3,42 = 21.91, p < 0.001), phylogenetic diversity (F_3,42 = 22.16, p < 0.001), Chao1 (F_3,42 = 23.76, p < 0.001), Simpson's diversity (F_3,42 = 6.70, p < 0.001), and Shannon's diversity (F_3,42 = 14.74, p < 0.001). Also, we found a significant difference in Chao1 between sites (F_1,42 = 22.26, p < 0.0001), but no differences detected for site or coral condition for the remaining parameters. Highest number of OTUs per sample was observed for P. lutea across all sites (average 208.46 ± 5.8), but with the highest number recorded for the bleached C. aspera samples from Praslin (360 ± 27) and lowest for the unbleached A. gemmifera samples at East Bay (55 ± 23; Table ). Total species richness was generally highest for C. aspera with the maximum recorded for unbleached samples at East Bay (Chao1 = 662.96 ± 6.05) indicating a more diverse and heterogeneous bacterial community, and lowest for the bleached A. muricata samples from East Bay (60.93 ± 31.90; Table ). Species evenness (Simpson's diversity) was highest for P. lutea ranging between 0.81 ± 0.16 and 0.97 ± 0.01, and lowest for A. muricata (ranging 0.41 ± 0.13 to 0.60 ± 0.09; Table ). Bacterial diversity was highest for P. lutea, where Shannon's diversity values ranged from 4.58 ± 1.64 to 6.36 ± 0.21, and lowest for A. muricata (Shannon's diversity between 1.35 ± 0.39 and 2.71 ± 0.55; Table ).

Summary statistics for alpha diversity indices of 16S rRNA gene amplicon sequences for microbial communities associated with Acropora gemmifera, Acropora muricata, Coelastrea aspera and Porites lutea at East Bay and Praslin for each category of coral condition (i.e., bleached and unbleached)

Site	Condition	n	# reads	OTU/Sample	PD_whole_tree	Chao1	Simpsons	Shannons
A. gemmifera
East Bay	Bleached	3	16,970.67 ± 848.99	67.33 ± 18.90	4.59 ± 0.99	130.76 ± 59.35	0.64 ± 0.14	2.73 ± 0.67
Unbleached	3	14,076.33 ± 3,515.56	55.67 ± 23.80	4.16 ± 1.46	85.70 ± 22.37	0.65 ± 0.16	2.75 ± 1.00
Praslin	Bleached	3	19,063.67 ± 4,324.15	83.33 ± 14.97	5.48 ± 0.65	154.16 ± 28.44	0.78 ± 0.09	3.57 ± 0.66
Unbleached	3	14,725.33 ± 2044.84	105.33 ± 18.76	6.67 ± 0.95	168.70 ± 34.66	0.82 ± 0.04	3.86 ± 0.49
A. muricata
East Bay	Bleached	3	22,968.33 ± 4,380.83	30.00 ± 12.50	2.46 ± 0.84	60.93 ± 31.90	0.41 ± 0.13	1.35 ± 0.39
Unbleached	3	15,612.67 ± 2,796.26	71.00 ± 23.71	5.03 ± 1.40	109.05 ± 26.55	0.60 ± 0.06	2.71 ± 0.55
Praslin	Bleached	3	16,103.67 ± 3,189.49	40.00 ± 8.88	3.08 ± 0.57	75.38 ± 16.22	0.60 ± 0.09	2.16 ± 0.43
Unbleached	3	10,968.00 ± 1954.59	41.00 ± 13.01	3.06 ± 0.89	75.17 ± 33.10	0.64 ± 0.12	2.33 ± 0.63
C. aspera
East Bay	Bleached	2	13,306.50 ± 5,128.50	184.50 ± 106.50	9.62 ± 3.95	306.88 ± 162.88	0.70 ± 0.28	4.61 ± 2.61
Unbleached	2	12,907.00 ± 6,179.00	360.50 ± 27.50	15.80 ± 1.33	662.96 ± 6.05	0.99 ± 0.01	7.71 ± 0.24
Praslin	Bleached	3	16,738.33 ± 696.74	253.33 ± 23.56	12.29 ± 0.99	510.21 ± 40.09	0.93 ± 0.03	6.01 ± 0.35
Unbleached	1	13,853.00 (n/a)	99.00 (n/a)	7.37 (n/a)	251.50 (n/a)	0.75 (n/a)	3.44 (n/a)
P. lutea
East Bay	Bleached	3	24,632.67 ± 113.81	183.00 ± 49.21	10.27 ± 2.26	317.75 ± 90.47	0.94 ± 0.02	5.70 ± 0.54
Unbleached	3	19,577.00 ± 6,361.46	277.67 ± 41.59	13.58 ± 1.65	470.00 ± 60.58	0.92 ± 0.05	6.24 ± 0.95
Praslin	Bleached	3	32,955.67 ± 5,782.41	141.67 ± 64.01	8.19 ± 2.96	269.59 ± 105.85	0.81 ± 0.16	4.58 ± 1.64
Unbleached	2	25,426.50 ± 8,188.50	231.50 ± 65.50	12.58 ± 2.51	438.13 ± 145.99	0.97 ± 0.005	6.36 ± 0.21
F(df) statistic	F3,42 = 0.58	F3,42 = 3.76	F3,42 = 2.92	F3,42 = 4.10	F3,42 = 1.41	F3,42 = 1.87
p value	0.634	0.020*	0.049*	0.014*	0.259	0.152

Note:

Statistical analysis (Univariate general linear model [GLM] for interactions between species ×site × coral condition are shown and significance is shown with an asterisk (*). Average ± SE are shown. n/a = no SE calculated for C. aspera unbleached sample from Praslin (n = 1).

Although there was no significant interaction for bacterial community composition between species, site, and coral condition, a significant difference was detected between species (PERMANOVA; pseudo F_3,47 = 4.84, p < 0.0001, Supporting Information Table S3a). This was further confirmed with ANOSIM (Global R = 0.376, p < 0.0001; Supporting Information Table S3b) with significant differences between all species except C. aspera and P. lutea. No differences (Global R = −0.012, p = 0.568) were detected between bleached and unbleached coral samples, and once the most abundant (top three ranked by relative abundance) families were excluded from analysis (Hahellaceae, Rhodospirillaceae, and Rhodobacteraceae), an interaction between species (PERMANOVA; pseudo F_3,47 = 5.20, p < 0.001; Supporting Information Table S4a) and species × site (pseudo F_3,47 = 1.63, p < 0.020; Supporting Information Table S4a) became apparent. The largest difference detected was between A. muricata and P. lutea (ANOSIM; R statistic = 0.675, p = 0.0002; Supporting Information Table S4b) while the most similar species were C. aspera and A. muricata (ANOSIM; R statistic = 0.07, p = 0.102; Supporting Information Table S4b) when the most abundant families were excluded from analysis. This was not detected in the full dataset, highlighting the importance of additional analyses excluding the most abundant taxa to detect differences in the less abundant ones that can be masked by those in higher abundance.

Similarity percentage analysis was used to identify the main contributing bacterial families within each of the coral species. Hahellaceae were the top contributing bacterial family in all coral species, contributing 31.25% in A. gemmifera, 44.87% in A. muricata, 15.78% in C. aspera, and 10.45% in P. lutea (Figure a–d, Supporting Information Table S5). Alteromonadaceae were the second highest contributing bacterial family for both Acropora sp. (4.81% A. gemmifera and 5.05% in A. muricata) while Rhodospirillaceae and Rhodobacteraceae comprised the top 3 contributing taxa for C. aspera and P. lutea (Supporting Information Table S5). PCoA of the dominant bacterial communities showed a distinct community composition for A. gemmifera and A. muricata, clustering separately from the majority of P. lutea and C. aspera (Figure e) along the primary axis (PCO1), explaining 38.8% of the total variation.

Average relative abundance (%) of bacterial community composition for (a) Acropora gemmifera, (b) Acropora muricata, (c) Coelastrea aspera and (d) Porites lutea classified as bleached or unbleached from East Bay and Praslin in Seychelles as a taxonomy stacked column plot to family level. Remaining taxa are grouped as “other”. Values displayed are mean relative abundances (n = 2–3). Each color represents one of the 10 most abundant bacterial families. (e) Principle Coordinate Analysis (PCoA) for the dominant bacterial taxa found in four coral species at East Bay and Praslin in Seychelles. Data was fourth‐root transformed, and a Bray‐Curtis similarity matrix was used with a correlation of 0.2. Ellipses denote similarity clusters of 20% (green dashed line). Percentages on axes indicate variation explained by the two coordinates

DISCUSSION

Microbiomes play a key role in contributing to coral fitness over space and time (Putnam et al., ; Suggett et al., ) and are known to exhibit broad changes across reefs persisting under different environmental conditions (Roder et al., ) and when subjected to atypical stress (Grottoli et al., ; Röthig et al., ; Ziegler, Seneca et al., ). Here, we provide the first characterization of the microbial community composition (i.e., Symbiodiniaceae and bacteria) for key reef‐building coral taxa of Seychelles across two different environments. In addition, we characterize the microbial communities associated with states of coral health (i.e., bleached and unbleached), collected during the 2016 marine heat wave that induced mass coral bleaching and mortality. In examining coral species that have previously been shown to be broadly heat stress sensitive (A. muricata, A. gemmifera) and tolerant (P. lutea, C. aspera) in Seychelles (Harris et al., ) and WIO (McClanahan et al., ), we have shown that bleaching susceptibility among coral species is indeed broadly consistent with differences in microbiomes. However, conserved microbiome signatures observed for bleached and unbleached colonies of all coral species suggest complex regulation of bleaching severity by the coral holobiont and genotype.

Decline in coral cover during 2016 mass bleaching

Heat stress is recognized as the most common cause of coral bleaching, and high record temperatures between 2015 and 2017 triggered the third global mass bleaching event, the most damaging to date (Hughes et al., ). At CMNP, about 40% of all corals bleached in April 2016 (Figure a), in particular species of Acropora and Pocillopora that resulted in large declines in coral cover recorded in April 2017. Specifically, decline of Acropora sp. and Pocillopora sp. from ca. 15%–20% (2016) to 1%–2% of total benthic cover was recorded (2017; Figure c). Such dramatic loss of these same taxa was similarly observed during previous recent marine heat wave events in Seychelles (Graham et al., ; Wilson et al., ) and other sites in the WIO (Baker, McClanahan, Starger, & Boonstra, ; McClanahan et al., ), reflecting the typically stress‐sensitive “boom and bust” nature of commonly fast‐growing branching taxa (Darling, Alvarez‐Filip, Oliver, McClanahan, & Côté, ; Zinke et al., ). The significant decline in abundance of Acropora sp. and Pocillopora sp. in 2016 suggests that these populations in Seychelles remain inherently susceptible to heat stress, and currently show no evidence for acclimatization to recent repeat heat wave events in this region unlike other WIO reef locations (Kenya; McClanahan, ).

Species‐specific Symbiodiniaceae composition

As expected, the major symbiont communities for corals sampled in CMNP broadly reflected those sampled from other Indo‐Pacific regions. For example, ITS2 type C3 (and additional types characterized by sequences from the C3 radiation; Thornhill, Lewis, Wham, & LaJeunesse, ) for Acropora sp. in Eastern Africa (Chauka, ), the Chagos Archipelago (Yang et al., ), Red Sea (Ziegler, Eguíluz et al., ) and Persian‐Arabian Gulf (Hume et al., ; Smith, Ketchum et al., ; Smith, Vaughan et al., ). The C3 sequence represents one of several major radiations within the genus Cladocopium (Thornhill et al., ) and our C3z type profile among Acropora sp. sampled here is 2 bp different from the C3 basal sequence. As such, it is currently unclear how genetically and/or phenotypically comparable our C3 type(s) are to these reports from elsewhere in the WIO previously (but see Ziegler, Eguíluz et al., for the Red Sea). Heat stress tolerance is clearly highly variable among the C3 radiation (Hume et al., ). However, the mass bleaching response for Acropora sp. within the CMNP commonly hosting C3z‐C3 would suggest that this Cladocopium type is inherently heat stress sensitive, and a highly conserved Acropora sp. host‐symbiont association. While we did not assess symbiont types associated with Pocillopora sp. that also experienced mass bleaching in CMNP, previous observations from the Chagos Archipelago (Yang et al., ) and Tanzania (Chauka, ) have identified almost synonymous associations with the C1 group (Pocillopora damicornis, Pocillopora verrucosa, Pocillopora eydouxi), another major radiation within Cladocopium (Thornhill et al., ).

Alternate host‐symbiont associations were observed among the coral taxa that exhibited comparatively little mass bleaching within CMNP, P. lutea and C. aspera. As with elsewhere in Eastern Africa (Chauka, ), the Red Sea and Persian‐Arabian Gulf (Smith, Ketchum et al., ; Smith, Vaughan et al., ; Ziegler, Eguíluz et al., ) and Pacific (e.g., LaJeunesse et al., ; LaJeunesse et al., ), P. lutea in the CMNP almost exclusively associated with ITS2 type C15. Perhaps most intriguingly was the high diversity of Symbiodiniaceae in C. aspera (formerly known as Goniastrea aspera—see Huang et al., ) comprising predominantly D1–D4 (–D6) ITS2 type profiles, in comparison to previously being associated with C3 types in the Indian Ocean (LaJeunesse et al., ), C1 types in Western Australia (Silverstein, Correa, LaJeunesse, & Baker, ), and D1a types in Thailand (Brown, Dunne, Edwards, Sweet, & Phongsuwan, ). Thermally tolerant D. trenchii (ITS2 type profile D1–4 is routinely observed in the WIO with a wide host range and where more than one Symbiodiniaceae type is detected within a coral species (LaJeunesse et al., ; see alsoSmith, Ketchum et al., ; Smith, Vaughan et al., ).

Bleaching observed for the heat stress‐sensitive taxa, A. muricata and A. gemmifera, in both reef environments was reflected by a loss of symbiont density (Figure e–h) and shifts in dominant C3 radiation types (Figure a,b; Supporting Information Figure S2), whereas the few unbleached colonies sampled at Praslin contained D1 types (A. muricata) or C3z/C3 types (A. gemmifera). Similarly, bleached colonies of the heat stress‐tolerant taxa were characterized by a loss of Symbiodiniaceae cells (but maintaining type C15 dominance, P. lutea) or alternate symbiont types (changing from D1 to predominantly C3 types, C. aspera). Together, these highlight important localized regulatory processes of bleaching susceptibility within these reefs. First, given the apparent stress‐sensitive nature of our Acropora sp. C3z/C3 associations, lack of bleaching for A. gemmifera at Praslin may indicate small‐scale environmental variability, as similarly observed for in‐shore Acropora sp. colonies on the Great Barrier Reef also in 2016 (Hoogenboom et al., ). Turbidity has been suggested to be critical in providing refuge against bleaching intensity (LaJeunesse et al., ; Oliver & Palumbi, ; Ulstrup & Van Oppen ). However, the mass bleaching observed at both our relatively turbid and clear sites (of similar extent, Figure ) would thus suggest small‐scale environmental variability afforded through complex reef habitats (e.g., shading from overhangs, Cacciapaglia & Woesik, ) are more important in providing refuge from heat stress. Bleaching of P. lutea and loss of C15 type cells at both sites similarly suggests localized small‐scale amplification of heat stress, for example, by high light (Hoogenboom et al., ). Second, alternate symbiont types of unbleached typically heat‐sensitive A. muricata (D1 type, Praslin) and bleached typically heat‐tolerant C. aspera (D1 type) reflects some capacity for Symbiodiniaceae community re‐organization needed for thermal acclimatization. D. trenchii (ITS2 type profile D1–D4) is a known coral host‐generalist (LaJeunesse et al., ) and heat stress‐tolerant taxon (LaJeunesse, Smith, Finney, & Oxenford, ). While it is unclear why unbleached A. muricata only associated with D1 types at Praslin, enhanced prevalence of Durusdinium spp. is consistent with increased turbidity (LaJeunesse et al., ; Oliver & Palumbi, ) that is typical of Praslin.

Stability in bacterial communities under thermal stress

Differences in bleaching susceptibility among taxa that are not easily resolvable through Symbiodiniaceae identity or host phenotype alone may result from differences in the associated bacterial communities (Morrow et al., ; Sweet & Bulling, ). Given the putative critical contribution of bacteria to coral holobiont functioning, changes in the identity and abundance of associated bacteria may allow for rapid acclimatization or adaptation to environmental change (Bourne et al., ; Cárdenas et al., ; Roder et al., ; Roder, Arif, Daniels, Weil, & Voolstra, ; Vega Thurber et al., ). Community shifts (such as a decrease in Endozoicomonas sp.) can occur in visibly healthy corals in degraded ecosystems (Ziegler et al., ) and changes to the bacterial community prior to visual signs bleaching (e.g., increase in Vibrio related sequences) have been shown to occur (Bourne et al., ). In contrast, we did not observe such obvious bacterial community shifts between bleached and unbleached conspecifics for heat stress susceptible or resistant coral species in Seychelles. This outcome is surprising since bacterial communities are expected to exhibit more rapid responses to stressors than Symbiodinium due to faster metabolism and generation times (Pogoreutz et al., ). As such, the observed bleaching susceptibility in our study may be more likely driven by differences in Symbiodiniaceae (via stress susceptibility) rather than by the hosts’ bacterial assemblages. That said, we acknowledge that our opportunistic sampling during the 2016 mass bleaching event meant that we were not able to compare samples before and after the beaching event. While more resolute temporal sampling throughout bleaching is needed to resolve for potentially different dynamics (and the role) of Symbiodiniaceae versus bacterial communities during heat stress, the direct comparison of bleached and unbleached colonies for all coral taxa during the heat stress allowed us to uniquely assess a putative contribution of bacterial taxa to bleaching susceptibility.

We identified three dominant bacterial families across all factors (species, site, and coral condition), suggesting a putatively important potential role in coral health and function due to their widespread prevalence and high abundance. Bacterial community stability was largely attributed to the abundance of Hahellaceae (comprising Endozoicomonas—ranging between 10% and 44% for all coral species; Supporting Information Table S5), congruent with recent work (Bayer et al., ; Neave, Apprill, Ferrier‐Pagès, & Voolstra, ; Neave, Rachmawati, et al., ; Pogoreutz et al., ; Pootakham et al., ). While the function of Endozoicomonas has not yet been defined, their genomes are significantly enriched in genes for carbohydrate transport and recycling as well as for protein and amino acid provision (Neave, Michell, Apprill, & Voolstra, ) and phenotypic assays confirm a high metabolic versatility in vitro (Yang et al., ). Despite the high abundance of Endozoicomonas in apparently healthy corals (Apprill, Hughen, & Mincer, ; Roder et al., ), strongly reduced abundances have been reported for stressed, diseased or bleached corals, suggesting they may be an indicator of coral health or habitat suitability (Bourne et al., ; Cárdenas et al., ; Meyer, Paul, & Teplitski, ; Röthig et al., ; Ziegler et al., ). However, similar to previous findings in bleached Red Sea P. verrucosa (Pogoreutz et al., ), we did not observe reduced abundances of Endozoicomonas in bleached colonies, even though colonies were visibly strongly bleached 1 week prior to sampling (see methods).

Alteromonadales and Rhodobacteraceae are commonly identified as members of the coral microbiome including larvae and juvenile early stage colonies (Apprill, Weber, & Santoro, ), some of which are able to degrade and assimilate dimethylsulfoniopropionate (DMSP; Reisch, Moran, & Whitman, ), an abundant carbon source in corals (Raina et al., ). Members of both bacterial families have previously been observed in visibly healthy and stressed corals (Li et al., ; Pantos, Bongaerts, Dennis, Tyson, & Hoegh‐Guldberg, ) and were fairly abundant in some of the investigated coral species in the present study (Alteromonadaceae were 3%–4% relative abundance and Rhodobacteraceae were up to 5% relative abundance; for details refer to Supporting Information Table S5). Given their presence in a range of hosts and environments (including seawater), it suggests this taxon is metabolically flexible and may provide important functions to the coral holobiont (Röthig et al., ).

Vibrionaceae are an opportunistic and potentially pathogenic bacterial family commonly associated with coral disease and have previously been linked to bleaching (Bourne et al., ; Garren et al., ; Tout et al., ). We found Vibrionaceae in all corals, albeit in lower abundance, ranging between a contribution of 1.12% (P. lutea) to 2.61% (A. muricata). The reported association of Vibrionaceae associated with unimpaired, apparently healthy corals has been reported (Bourne & Munn, ). Nevertheless, the consistent association of these bacteria suggests that corals in Seychelles harbor stable, and presumably locally adjusted microbiomes, allowing the corals to cope well in the ambient environment (sensu Hernandez‐Agreda, Leggat, Bongaerts, & Ainsworth, ).

Bacterial diversity aligns with bleaching susceptibility

Healthy corals generally comprise specific, stable, and uneven microbial assemblages indicating host‐selected microbiomes (Bayer et al., ; Bourne et al., ). We observed the highest species richness (Chao1), evenness (Simpson's), and bacterial diversity (Shannon's) in the more heat‐tolerant massive corals (P. lutea, followed by C. aspera) compared with stress‐sensitive branching corals (A. muricata higher than A. gemmifera), consistent with findings from Liang et al. (). Acropora muricata also exhibited a highly uneven bacterial community that was dominated by Endozoicomonas (up to 44% contribution), compared with P. lutea that had the highest evenness and diversity, with only a 10% contribution of Endozoicomonas (Table , Supporting Information Table S5). The higher bacterial diversity found at Praslin (the turbid site) compared to East Bay (clear water site) is also supported by work showing corals on anthropogenically impacted reefs harbor higher bacterial diversity than those further from the disturbance (Morrow, Moss, Chadwick, & Liles, ). Praslin is within a bay bordered by human populations and thus more prone to anthropogenic influence from run‐off and pollution, compared with East Bay (which is located further from human influence off Curieuse Island). Furthermore, bacterial community composition in corals from more impacted shallow sites can display higher intraspecies variability with thermal stress (Littman, Willis, & Bourne, ). Consequently, it is plausible that the changes in bacterial communities we observed between sites may in part reflect differences in overall reef health through proximity to more localized human populations at Praslin compared to East Bay.

Responses of less abundant or rare members of the bacterial microbiome may not become apparent until the dominant bacterial members are excluded from analysis. Indeed, we identified a significant interaction between species and site (Supporting Information Table S4) when the most abundant families were excluded from analysis. Here, the bacterial family contributing the highest percentage for both sites was Enterobacteriaceae, with Vibrionaceae also in the top three. In contrast, East Bay had a 2.94% contribution from Pseudomonadaceae while Praslin had Alteromonadaceae (7.7%) in the top three dominant taxa (once the most abundant families were excluded). This difference was only detected upon exclusion of the most abundant families, and a comparison of bacterial taxa for coral condition revealed no difference in the order of the most abundant taxa (Enterobacteriaceae, Alteromonadaceae and Vibrionaceae for unbleached and bleached). Such site‐specific differences in dominant bacterial taxa have similarly been reported for the heat‐sensitive species (Acropora sp.) exposed to a range of stressors (Apprill et al., ; Littman, Willis, Pfeffer, & Bourne, ; McDevitt‐Irwin, Baum, Garren, & Vega Thurber, ; Pantos et al., ). As such, the low flexibility of bacterial communities associated with heat‐sensitive corals in Seychelles may reduce their ability to rapidly respond to environmental stress (Pogoreutz et al., ), as highlighted by the significant loss in coral cover reported during the 2016–2017 mass bleaching event for Acropora sp. (and Pocillopora sp.; Figure a).

Finally, it is possible that suboptimal partnerships with Symbiodiniaceae algal symbionts could enhance vulnerability to opportunistic bacterial infection. As such, changing Symbiodiniaceae taxa following bleaching events may not only provide immediate benefits to the coral holobiont in terms of thermal tolerance, but could also result in a longer‐term trade‐off with disease resistance (Littman, Bourne, & Willis, ). Recent work has shown coral host intraspecific differences in Symbiodiniaceae composition correlated with disease susceptibility (Rouzé, Lecellier, Saulnier, & Berteaux‐Lecellier, ). Specifically, the predisposition to disease and infection by Vibrio spp. was positively correlated with Symbiodiniaceae genus Symbiodinium (formerly clade A), but negatively correlated with Durusdinium in Acropora cytherea (Rouzé et al., ). In our study, we found Durusdinium types D4 and D9 unique to P. lutea (and bleached samples of C. aspera) and previous work has shown increased occurrence of Durusdinium types is consistent with increased turbidity. Furthermore, Durusdinium‐infected corals have shown shifts in the associated bacterial community under heat stress, while no shifts were reported for Cladocopium within the same coral host (Littman et al., ). It might therefore be informative to examine for bacterial associates that co‐occur with Symbiodiniaceae types in future studies.

In characterizing the microbiome composition of four species of Seychelles’ corals during the most severe mass bleaching event on record, we have shown that susceptibility to stressors is reflected by underlying microbiome community structures. Specifically, bleaching susceptibility among coral taxa corresponds largely to differences in specific host‐Symbiodiniaceae associations, while the bacterial microbiome community remains largely stable. As such, unbleached colonies of bleaching‐susceptible corals likely persist through availability of small‐scale (“micro”) environmental refuges that can dampen the effect of heat stress, (e.g., shading from surrounding substrate). While bacterial communities were highly similar between bleached and unbleached corals, differences were observed between species, adding to previous evidence for species‐specificity. Thus, microbiome profiling of both Symbiodiniaceae and bacterial communities may provide new capacity to more broadly identify stress susceptible versus tolerant coral populations, which is needed to aid targeted management in reef systems such as Seychelles.

ACKNOWLEDGMENTS

We wish to extend our greatest thanks to the Seychelles National Parks Authority (SNPA) for unparalleled support across all years in generation of the annual ecological assessments, and sample collection in 2016. Also, we are indebted to the Seychelles Ministry for Environment, Energy and Climate for permitting non‐commercial transfer of genetic material of analysis for this publication. We also thank Craig Michell for sequencing library preparation and the KAUST Bioscience Core Lab (BCL) for sequencing. Data collection, analysis and write up was supported by The Earthwatch Institute and the Mitsubishi Corporation (DJSm), an Endeavour Fellowship (EFC), an ARC Discovery Grant (DJSu; DP160100271) as well as a KAUST baseline research funds to CRV.

CONFLICT OF INTEREST

None declared.

AUTHOR CONTRIBUTION

EFC, DJSm, and DJSu designed and conceived the experiment. BCCH, CP, and CRV contributed new reagents or analytical tools. SGG, EFC, DJSm, TK, EOO, GG, BCCH, CP, CRV, and DJSu performed research. All authors analyzed and interpreted data. SGG, DJSu and EFC wrote the first draft of the manuscript. All authors contributed to and approved the manuscript.

DATA ACCESSIBILITY

Raw sequencing data determined in this study are available under NCBI BioProject ID PRJNA477731 (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA477731) for 16S rRNA gene data and under NCBI BioProject ID PRJNA477732 (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA477732) for ITS2 sequence data. Other data are available in the Supporting Information.