High frequency temperature variability reduces the risk of coral bleaching

ARTICLE Corrected: Author correction DOI: 10.1038/s41467-018-04074-2 OPEN High frequency temperature variability reduces the risk of coral bleaching 1 2,3 4 5 6 Aryan Safaie , Nyssa J. Silbiger , Timothy R. McClanahan , Geno Pawlak , Daniel J. Barshis , 7 8 9 1 James L. Hench , Justin S. Rogers , Gareth J. Williams & Kristen A. Davis Coral bleaching is the detrimental expulsion of algal symbionts from their cnidarian hosts, and predominantly occurs when corals are exposed to thermal stress. The incidence and severity of bleaching is often spatially heterogeneous within reef-scales (<1 km), and is therefore not predictable using conventional remote sensing products. Here, we systematically assess the relationship between in situ measurements of 20 environmental variables, along with seven remotely sensed SST thermal stress metrics, and 81 observed bleaching events at coral reef locations spanning five major reef regions globally. We find that high-frequency temperature variability (i.e., daily temperature range) was the most influential factor in predicting bleaching prevalence and had a mitigating effect, such that a 1 °C increase in daily temperature range would reduce the odds of more severe bleaching by a factor of 33. Our findings suggest that reefs with greater high-frequency temperature variability may represent particularly important opportunities to conserve coral ecosystems against the major threat posed by warming ocean temperatures. 1 2 Department of Civil and Environmental Engineering, University of California, Irvine, CA 92697, USA. Department of Biology, California State University, Northridge, 18111 Nordhoff Street, Northridge, CA 91330-8303, USA. Department of Ecology and Evolutionary Biology, University of California, Irvine, CA 4 5 92697, USA. Marine Programs, Wildlife Conservation Society, 2300 Southern Boulevard, Bronx, NY 10460, USA. Department of Mechanical and Aerospace Engineering, University of California San Diego, 9500 Gilman Drive, MC0411, La Jolla, CA 92093, USA. Department of Biological Sciences, Old Dominion University, Mills Godwin Building 110, Norfolk, VA 23529, USA. Nicholas School of the Environment, Duke University, 135 Duke Marine Lab Road, Beaufort, NC 28516, USA. Department of Civil and Environmental Engineering, Stanford University, 473 Via Ortega, Y2E2 Rm 126, Stanford, CA 94305, USA. School of Ocean Sciences, Bangor University, Anglesey LL59 5AB, UK. Correspondence and requests for materials should be addressed to A.S. (email: safaiea@uci.edu) or to K.A.D. (email: davis@uci.edu) NATURE COMMUNICATIONS (2018) 9:1671 DOI: 10.1038/s41467-018-04074-2 www.nature.com/naturecommunications 1 | | | 1234567890():,; ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04074-2 oral reef ecosystems provide subsistence nutrition, coastal Here, using a global suite of in situ data, we compare and assess protection, and revenue from tourism to hundreds of the ability of 20 commonly used environmental variables and 7 1,2 Cmillions of people globally , and are valued at trillions of remotely sensed variables to explain observed bleaching pre- dollars annually . Especially during recent years, coral reefs are valence, testing the hypothesis that including high-frequency increasingly threatened by accelerated rises in ocean temperatures temperature variability as one of these model variables will lead to 4–6 owing to global warming . Elevated seawater temperatures are more accurate predictions. Analyzed data include records of 5,7 the primary cause of mass coral bleaching , or the loss of pig- in situ temperature time series at 118 reef locations from five mentation due to the collapse of the symbiotic relationship major reef regions with sampling intervals of ≤3 h and continuous between the coral host and its endodermal dinoflagellate algae measurements of ≥1 year, as well as precise information on 7,8 9 (zooxanthellae) . Bleached corals are susceptible to disease and habitats and depths (Supplementary Data 1), along with 81 spa- 9,10 reduced carbonate accretion , and prolonged bleaching will tially and temporally coincident, quantitative coral bleaching 5,11,12 lead to mortality . observations (Supplementary Data 2). Each of the 81 bleaching Thermal stress on corals and regional bleaching events are observations was matched to its own spatiotemporally coincident most often predicted by the magnitude and duration of remotely temperature time series data, such that 46 of the 118 temperature sensed sea surface temperatures (SSTs) above a fixed, locally time series were used in the subsequent bleaching analysis. 5,8,13 defined average summer threshold temperature . A con- Bleaching observations, which are most often reported as the ventionally used metric for quantifying these temperature average percent of colony or transect area bleached, were stan- anomalies is provided by the National Oceanic and Atmospheric dardized to ordinal-valued “bleaching prevalence scores” (1: Administration’s (NOAA) Coral Reef Watch program, which has ≤10%; 2: 10−25%; 3: 25−50%; 4: >50% of reef area bleached), reported cumulative thermal stress on reefs twice a week since representing mild to pervasive bleaching, respectively (Methods). 1997 . Furthermore, bleaching predictions from remotely sensed The influence of different factors on bleaching prevalence scores temperatures can be improved by including SST-based calcula- are evaluated by selecting covariates from a pool of 20 explana- 15,16 tions of interannual temperature variability and coral sensi- tory variables (depth, latitude, and 18 thermal metrics) grouped tivity to thermal stress exposure . However, the relatively coarse into 8 categories of metrics often used to predict bleaching spatiotemporal resolution of the remotely sensed data prevents (Table 1). In addition to these in situ variables, we also include 7 ensuing thermal stress quantifications from identifying the often analogous and conventional remotely sensed SST thermal stress observed significant spatial heterogeneity in bleaching that occurs metrics (Table 1). After standardizing all covariates and fitting 18–20 within reef regions and individual reefs . The response of them to ordinal-valued bleaching prevalence scores using ordinal reefs to temperature at these smaller spatial scales is complex and logistic regression (OLR) models (Methods), we conclude that putatively depends on a combination of organism-level and reef- high-frequency temperature variability, specifically the average scale factors such as coral life-history strategies and stressor daily temperature range (DTR) of the 30 days preceding a cotolerances , the history and duration of thermal stress expo- bleaching observation, is the most influential covariate in pre- 22,23 24,25 sure , the rate of change in seawater temperature , flow dicting the bleaching response, and serves to attenuate the pre- 26 6 27 conditions , heterotrophic feeding , turbidity , and the intensity valence of bleaching. 28,29 and history of exposure to solar radiation . In turn, many of these environmental conditions are mediated by reef-scale factors Results 30 31 32 such as waves , winds, tides , and daily heating and cooling . Variation among in situ explanatory variables. A principal Site-specific studies suggest that historical temperature varia- components analysis (PCA) displays the projection for each site bility within diurnal time scales affects corals’ physiological tol- onto the 2D plane that accounts for the most variance in the 20 19,26,33,34 35 erance and performance under thermal stress. For in situ explanatory variables (Fig. 1), and the locations of the example, it has been theorized that corals located in areas char- loading vectors reveal how these explanatory variables relate to acterized by large temperature fluctuations, such as reef flats or their respective groupings. The first principal component shallow lagoons, may be better acclimatized or adapted to thermal accounts for 44.2% of the variation in the explanatory variables, stress, and therefore more resistant to anomalous temperatures and is largely driven by high-frequency temperature variability and bleaching, than corals in areas where temperatures are more and cumulative thermal stress. 36–38 stable, such as on reef crests or reef slopes . Other studies have suggested that water temperatures in the weeks or months leading up to peak temperatures are critical in determining the Spatiotemporal dependence of diurnal temperature variability. coral physiological response. A recent analysis of experimentally The thermal metrics computed from temperature time series were heated corals from the Great Barrier Reef showed that bleaching highly variable across sites, but regardless of location and depth, and cell death responses were indeed lower when the thermal all 118 time series show significant temperature variations in the exposure included a moderate pre-stress followed by a short high-frequency band (Supplementary Note 1; Supplementary recovery period (i.e., a “protective temperature trajectory”) . Fig. 1), which we define as 0.727−4 cycles per day (cpd). Power Depending on intrinsic properties of coral physiology such as spectra of temperature variations were calculated for each loca- energy reserves and algal phenotypic plasticity , pre-peak tem- tion, and the ratios of high-frequency band to seasonal band peratures may either protect against or exacerbate bleaching at (0.012−0.143 cpd, or 1/7 to 1/84 days) variance in these spectra peak temperatures . Taken together, a growing body of evidence were used to characterize the relative importance of variance thus suggests that historical temperature variability, and parti- within the high-frequency band. This ratio correlates with the cularly, “high-frequency” temperature variability, which we define inverse of depth (r = 0.381, Student's t-test p < 0.05), indicating as occurring within diurnal or shorter periods, may play an that the relative contribution of high-frequency variability to the important role in determining corals’ physiological responses to variance within a temperature time series is stronger at shallower thermal stress and thereby reef-scale vulnerability to bleaching. In sites (Supplementary Fig. 2a). At back reef, reef flat, and reef slope turn, a better understanding of reef-scale bleaching risk factors habitats, these ratios were on average 1.83, 0.68, and 0.44, could help coastal management efforts to identify natural refugia respectively, while across all locations, this ratio was 1.02 (Sup- and may be important for the recovery of coral communities plementary Fig. 2b). Furthermore, these ratios differed sig- 42 2 following a bleaching event . nificantly among the three habitats (Kruskal-Wallis, χ = 24.66, 2 NATURE COMMUNICATIONS (2018) 9:1671 DOI: 10.1038/s41467-018-04074-2 www.nature.com/naturecommunications | | | NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04074-2 ARTICLE Table 1 List of explanatory variables used in the ordinal logistic regression analysis Category Variable [Units] Identifier Description Ref. 1. Depth Instrument depth [m] depth In situ water depth of instrument 2. Background Latitude [DD] lats Latitude of instrument Conditions Maximum Monthly Mean (MMM) [°C] MMM Maximum of monthly mean climatology from entire time Total series MMM Maximum of monthly mean climatology using data only before and during bleaching event MMM Maximum of monthly mean climatology using 4 km weekly 4km CoRTAD SST data MMM Mean of maximum monthly SST from each year in Max climatological time period 3. Cumulative Degree Heating Weeks (DHW) [°C-weeks] DHW Trapezoidal integration of temperatures in excess of MMM+ Thermal Stress 1 °C during 90 days preceding a bleaching event DHW Trapezoidal integration of temperatures in excess of MMM+ 1 °C during 30 days preceding a bleaching event DHW Degree heating week product from 4 km weekly CoRTAD SST 4km data Cumulative Summer Anomaly (CSA) [°C-days] CSA Trapezoidal integration of temperatures in excess of MMM+ Total 1 °C during all summer periods through entire time series CSA Trapezoidal integration of temperatures in excess of MMM+ Before 1 °C during summer periods before and during a bleaching event CSA Trapezoidal integration of temperatures in excess of MMM+ During 1 °C during summer of bleaching event 4. Acute Presence/absence of acute temperature Acute1 Binary value indicating whether any of the daily mean Thermal Stress anomaly [binary] temperatures within 90 days preceeding a bleaching event exceeded MMM+ 1°C Acute1 Acute1 computed using 4 km weekly CoRTAD SST data 4km Acute2 Binary value indicating whether any of the daily mean temperatures within 90 days preceeding a bleaching event exceeded MMM+ 2°C Acute2 Acute2 computed using 4 km weekly CoRTAD SST data 4km 5. Thermal Type of induced thermal tolerance prior to acute TT 0: No thermal stress (temperatures do not exceed MMM+ 2° Trajectory thermal stress, using twice-weekly averaged C within 90 days prior to survey date) temperatures [ordinal] 1: Protective Trajectory (temperatures exceed MMM, then have a recovery period below MMM for at least 10 days prior to exceeding MMM+ 2 °C) 2: Single Bleaching Trajectory (temperatures exceed both MMM and MMM+ 2 °C without a 10-day recovery period in between) 3: Repetitive Bleaching Trajectory (temperatures exceed MMM+ 2 °C in two peaks separated by 9 days) 6. Heating Rate Rate of spring-summer temperature change [° ROTC Mean rate of temperature change during spring and summer SS C/day] of all years ROTC Mean rate of temperature change during 90 days preceding a 90- bleaching event using CoRTAD SST data 4km ROTC Mean rate of temperature change during spring and summer SS-4 of all years using CoRTAD SST data km 7. High- Daily Temperature Range (DTR) [°C] DTR Mean DTR over entire time series Total Frequency Temperature Variability DTR Mean DTR of all spring and summer periods SS DTR Mean DTR of all fall and winter periods FW DTR Mean DTR over 90 days preceding a bleaching event DTR Mean DTR over 30 days preceding a bleaching event 8. DTR Measure of shape of distribution of all DTR kurtosis Kurtosis of full time series of DTR values Distribution values w/in a time series [−] Shape skewness Skewness of full time series of DTR values Variables are grouped according to eight categories representing different aspects of ecologically relevant environmental and temperature factors. Seasons were defined such that each season spanned three complete months, and austral and boreal summers were December through February and June through August, respectively NATURE COMMUNICATIONS (2018) 9:1671 DOI: 10.1038/s41467-018-04074-2 www.nature.com/naturecommunications 3 | | | ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04074-2 30°N Background Caribbean Tahala 0.08 Red conditions 15°N Sea 0° Western Pacific Indian Heating Ocean 0.06 Nelly Bay 15°S Great rate Barrier Reef 30°S 120°W60°W0° 60°E120°E 0.04 0.02 High-frequency Tahala temperature variability Thermal trajectory Nelly Bay −0.02 Distribution shape of Cumulative daily temperature range thermal stress −0.04 Acute thermal stress Depth −0.02 0 0.02 0.04 0.06 0.08 Component 1 (44.2% of variation) Fig. 1 First two axes of variation of site-specific explanatory variables. Biplot of principal components analysis (PCA) showing the first two components (44.2% and 18.8%, respectively) that explain the majority of the variance in the matrix of 20 in situ explanatory variables (Table 1) used to explain bleaching prevalence. The light gray dots (“scores”) each represent temperature time series associated with a distinct bleaching event at a given reef site. Gray dots that are close to each other have more similar temperature environments than dots further apart. The vectors are colored according to the categories described in Table 1. The time series inspected later in Fig. 4 are also indicated by red squares (Tahala and Nelly Bay shoreward habitats) and blue triangles (Tahala and Nelly Bay seaward habitats). The magenta circles in the inset map indicate the locations of all 118 in situ time series, with their associated reef regions labeled. The map was created using the MATLAB package “M_Map”, created by Rich Pawlowicz under the license Copyright (c) 2014, Chad Greene. All rights reserved df = 117, p< 0.05; Supplementary Fig. 2b). Although the mag- shoreward locations or thermally stable seaward ones) are nitude of diurnal temperature fluctuations varies by location, the reflected in the mean, skewness, and kurtosis of DTR values ubiquity and prominence of temperature variance in this fre- (Fig. 2d). Shallower and more shoreward sites have a peak in their quency band indicated by these average ratios reflects the DTR distributions corresponding to a larger DTR value, and importance of some common physical forcing processes gov- furthermore, their distributions take on more extreme values than erning the flow and heating of reef waters, such as diurnal solar those from sites in deeper and more seaward locations. For 31,32,43,44 heating, tides, winds, and waves . example, at Heron Island in the Great Barrier Reef, the mean Power spectra of six representative time series from different DTR of 4.23 °C on the reef flat was over three times as large as reef regions (Fig. 2a) show a broad range of temperature that of the reef slope (Fig. 2d). The implications of these different variability from annual to hourly periods (see Supplementary thermal microclimates for resistance to thermal stress and Fig. 1 for other spectra). Yearly composites of mean water resilience to bleaching are discussed below. temperature and DTR (Fig. 2b) both show prominent seasonal cycles (Supplementary Note 2): the magnitude of daily tempera- ture fluctuations was seasonally dependent (Kruskal-Wallis, p < The effect of diurnal temperature variability on bleaching. 0.01) for 96% of reefs in our study (113 of 118 time series), with Ordinal logistic regression (“logit”) models were computed for all maximum DTRs occurring most often in spring and summer permutations of selecting at most one variable from each of the months (74% of time series, Supplementary Fig. 3), and eight categories in Table 1 (a total of 10,367 models), with minimum DTRs occurring most often in fall and winter months bleaching prevalence scores as the response variable. Corrected (also 74% of time series, Supplementary Fig. 3). On global scales Akaike’s Information Criterion (AIC ) values were used to rank (~10 km), latitudinal gradients in solar forcing drive variations the logit models, where the model with the lowest AIC value was in seasonal temperature patterns on reefs (Supplementary Fig. 4), ranked the highest (Fig. 3b). The model coefficients indicate the but there is also considerable heterogeneity in thermal environ- association of tested variables with bleaching prevalence score, ments at reef-scales (~10 m) due to variation in depth and such that positive coefficients indicate a “mitigating” effect on 32,42,45 circulation . The differences in thermal environments at bleaching prevalence, and negative coefficients an “exacerbating” reef-scales are often greatest in the high-frequency band (daily effect on bleaching prevalence. and tidal timescales; Fig. 2c). Dramatically different thermal “High-Frequency Temperature Variability” (Table 1) was used environments can be found at locations separated by 10s or 100s to capture temperature variability on diurnal and shorter time of meters on a reef, as illustrated by 7-day temperature time series periods, a metric that is important for characterizing differential 19,32,46 from various locations on the same island, or different habitats reef- and habitat-scale microclimates . In the best model within a given reef (Fig. 2c). For example, during a week in (Fig. 3a), high-frequency temperature variability, specifically the November 2009, two sites in American Samoa that are separated average DTR over the 30 days preceding a bleaching event by <2 km and at similar water depths experienced average DTRs (DTR , Table 1) was the most influential metric for predicting of 1.78 and 0.51 °C (Fig. 2c, sites OF3 and OF5 respectively). bleaching prevalence score, with greater daily temperature Differences in the distributions of DTRs that distinguish variability serving as a mitigating factor (Fig. 3b). Furthermore, microclimates within a reef system (e.g., thermally variable among all models within 2 AIC units of the highest ranked 4 NATURE COMMUNICATIONS (2018) 9:1671 DOI: 10.1038/s41467-018-04074-2 www.nature.com/naturecommunications | | | Component 2 (18.8% of variation) NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04074-2 ARTICLE a b c d 37 200 200 m apart Depths: 34.5 100 0.9 m –3 0.9 m 1.2 years 32 0 19 20 21 22 23 24 25 26 0 2 4 6 8 Aug-2010 30 30 500 604 m apart 29 Depths: 28.5 4.3 m –3 28 5.2 m 1.8 years 27 0 5 6 7 8 9 10 11 12 0 0.5 1 1.5 2 2.5 May-2013 31 30.5 400 1641 m apart 10 30 Depths: 28.75 200 0.8 m 0.85 m –3 7.7 years 27 0 21 22 23 24 25 26 27 024 68 Nov-2009 31.5 400 9 m apart Depths: 30.5 200 2 m 10 m –3 28 7.1 years 29.5 0 123 4 5678 0 0.5 1 1.5 2 Apr-2010 31.5 1000 125 m apart 1 Depths: 30.25 500 2.4 m 5.6 m –3 17.7 years 29 0 45678 910 11 0246 Jan-2006 35 400 1173 m apart Depths: 29.75 200 1 m 7.9 m 16.3 years 24.5 0 18 19 20 21 22 23 24 25 0 246 8 10 1/365 1 2 90 180 270 360 Feb-1998 Daily temperature range (°C) Cycles per day Yearday Fig. 2 Temperature variability of six reef records. a Power spectra of temperature for TA3, P21, OF3, VT1, HW1, and HR1, with asterisks marking significant peaks, b yearly composites of mean daily temperatures and temperature ranges (red and pink shading respectively) for the same six time series in a, c 7-day trends in temperatures at two different habitats on the reef, and d histograms of daily temperature range at the same two habitats on each reef. In each case, reef locations are shown in maps on the left (for site information see Supplementary Data 1), the full duration of temperature records are indicated in a, and the great-circle distances between same-reef sites are indicated in d. The maps were created using the MATLAB package “M_Map”, created by Rich Pawlowicz under the license Copyright (c) 2014, Chad Greene. All rights reserved 45 19,48 model (i.e. ΔAIC = AIC – min(AIC ) ≤ 2, Fig. 3a), high- depths , may develop greater thermal tolerance , while C C C frequency temperature variability was both the greatest mitigating deeper coral habitats, despite their propensity for milder diurnal factor of bleaching prevalence score and the most influential temperature variability (outside of internal wave-influenced 49–51 covariate—more influential than widely used metrics of acute and regions ), may serve as refuge areas resistant to the intrusion of hot water , perhaps facilitating recovery of coral cover cumulative thermal stress by a factor of 2 and 3 times, respectively (Fig. 3c). Using globally averaged values of following bleaching events . explanatory variables, our highest-ranked logit model (Fig. 3b) “Background Conditions”, “Cumulative Thermal Stress”, and implies that, in native units, a 1 °C increase from the mean DTR “Acute Thermal Stress” were the three explanatory variable value would decrease the odds of more prevalent bleaching by a categories largely suspected of exacerbating bleaching. “Back- factor of 33. To standardize this, each unit increase in high- ground Conditions” (Table 1) consisted of the average summer- frequency temperature variability (i.e., DTR ) would reduce the time, or maximum monthly mean (MMM), temperature, but 2.66 odds of more prevalent bleaching by a factor of e = 14.3. computed from our in situ time series data, as opposed to Contrasting this against a unit increase in cumulative thermal conventional remotely sensed SST data . “Background Condi- stress (i.e., DHW ), which would only increase the odds of more tions” also included the latitude of the temperature logger, a prevalent bleaching by a factor of 2.6, highlights the dominant variable that served as a proxy for unresolved oceanographic influence of diurnal temperature variability on reef-scale bleach- factors related to the large-scale processes that influence ing prevalence. climatologies. The “Cumulative Thermal Stress” category “Depth” (Table 1) was taken as the mean depth of the water (Table 1) encompassed various methods for the computation of temperature measurement, in meters below the surface, for each the magnitude and duration of acute in situ thermal stress site, and is also representative of local water column depth as exposure on reefs. Similar to the MMM, cumulative thermal sensors were placed near the bed. Depth was the second-most stress is traditionally derived from remotely sensed SSTs and is effective predictor of bleaching prevalence (Fig. 3c), with deeper among the most common metrics used to predict coral 25,33 reefs less likely to experience pervasive bleaching. However, bleaching . The “Acute Thermal Stress” category (Table 1) “depth” is also a proxy for other characteristics of the reef sites was included as a safeguard to differentiate sites with tempera- such as habitat (e.g., deeper forereefs and lagoons, shallow reef tures that may not have exceeded MMM+ 1 °C (i.e. no thermal flats) and light intensity, which decays exponentially with depth. stress) yet still experienced bleaching. Consistent with the well- Although the logit models preclude significant collinearity of established perspective that anomalously high temperatures are tested variables (Methods), corals at shallow depths may the primary cause of coral bleaching , among our highest ranked experience greater high-frequency temperature variability , models, bleaching was most exacerbated by greater cumulative although accounting for water flow can complicate this and acute thermal stress, and also, to a lesser degree, by increases 26,47 25 interpretation as it pertains to bleaching . High-frequency in MMM temperature and heating rate . “Heating Rate” temperature variability and depth may mitigate bleaching in (Table 1) was the average rate of change in spring to summer complementary ways: habitats with greater high-frequency temperatures, which is believed to have a positive relationship temperature variability, which are likely to be found at shallower with bleaching-induced tissue damage, and this time period has NATURE COMMUNICATIONS (2018) 9:1671 DOI: 10.1038/s41467-018-04074-2 www.nature.com/naturecommunications 5 | | | Power spectral density, (°C /cyclers per day) Daily temperature range (°C) Temperature (°C) Histogram counts 29.2 ± 0.4 46.0 ± 0.4 –3.3 ± 0.1 –4.2 ± 0.2 –5.0 ± 0.4 –13.7 ± 0.8 –20.2 ± 0.2 ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04074-2 acute stress exposure. Although not as influential as the above variables, a no-stress or protective thermal trajectory (i.e., a pre- stress, sub-bleaching warming period, followed by a cooler recovery period) is more likely to result in lower bleaching prevalence than a single bleaching trajectory (temperatures that cross the bleaching threshold without a prior protective event) or a repetitive bleaching trajectory (Fig. 3b, c and Table 1). Finally, the “Shape of DTR Distribution” category (Table 1) was used to 25 capture the skewness and kurtosis of DTR values derived from each time series to represent the symmetry and tail-density of DTR distributions. While these variables were not present in any of the highest ranked models, kurtosis and skewness of n = 1 temperature time series have been associated with site-specific n = 19 n = 10,347 increased thermal tolerance . To summarize the results of our highest-ranked logit model, we ΔAIC = 2 can examine how manipulating each covariate, while holding all others at their mean values, will change the probability of bleaching (Fig. 4). For example, a 0.88 °C decrease in high- Number of model parameters frequency temperature variability (DTR ) from its mean value would increase the probability of Category 4 bleaching from 12% to 75%, for a change of 63% (Fig. 4a), and a depth decrease of 5 m 2.66 ± 0.10 would increase this probability by 41% (Fig. 4b). Similarly, a 0.03 °C/day increase in ROTC from its mean value would SS increase the probability of Category 4 bleaching by 34% (Fig. 4c), 1.85 ± 0.09 and a 1 °C-weeks increase in DHW would increase this probability by 44% (Fig. 4f). To broaden the applicability of our conclusions, we repeated the OLR analysis, with the addition of remotely sensed SST- derived covariates, to determine how the results would differ from the in situ driven models. We obtained weekly 4 km resolution CoRTAD SST data , using the data pixels closest to the coordinates of our in situ loggers, and used this SST dataset to –0.66 ± 0.03 compute covariates within the Background Conditions, Acute –0.74 ± 0.05 –1 –0.90 ± 0.13 and Cumulative Thermal Stress, and Heating Rate categories –0.96 ± 0.13 (Table 1). This resulted in an improved highest ranked model (Methods), that included six covariates, three of which (MMM, DHW, Rate of Temperature Change) were computed using the SST, as opposed to in situ, data (Fig. 5a). However, similar to the highest-ranked model fit to exclusively in situ data, covariates representing High-Frequency Temperature Variability, specifi- cally DTR , and depth were again the dominant drivers of bleaching, and served as mitigating factors (Fig. 5a, b). Similarly, among covariates that exacerbated bleaching, Acute and Cumu- lative Thermal Stress provided the strongest influence (Fig. 5a, b), while Background Conditions (MMM , Table 1) represented a 4km mild exacerbating effect. A notable difference occurring in these new models was the opposite effect Heating Rate had from before; whereas in the exclusively in situ models, Heating Rate exacerbated bleaching, these SST-based models imply stronger heating rates serve to mitigate bleaching. Ultimately, these –10 Fig. 3 In situ explanatory variables of bleaching and their standardized logit –20 coefficients with greatest predictive power. a ΔAIC , computed as AIC – C C min(AIC ), values of all 10,367 runs of an ordinal logistic regression model, –30 where models within ΔAIC ≤ 2 (dashed line and gray shaded region) are statistically indistinguishable, of which there were 20. b The best model been shown to be crucial for determining the fate of corals to (i.e. ΔAIC = 0) included six variables, of which high-frequency summertime bleaching susceptibility . The “Thermal Trajectory” temperature variability was the absolute most influential and also greatest (Table 1) category followed the methodology of a previous study mitigating factor to bleaching prevalence. c Summing across 20 that highlighted the role of protective warm, pre-stress tempera- indistinguishably good models (i.e. within ΔAIC ≤ 2), high-frequency tures as being important for resilience to bleaching from intense temperature variability was consistently most influential. Variable acute stress temperature events . Our results reinforce recent categories are shown in Table 1. Delete-1 jackknife standard error bars are findings that a reef’s thermal trajectory is a significant predictor of shown in (b), while the standard error bars shown in (c) were obtained by bleaching prevalence (Fig. 3c), with thermal tolerance conferred summing in quadrature the individual standard errors from each of the 20 by exposure to a protective, sub-lethal bleaching stress prior to models computed after delete-1 jackknife resampling 6 NATURE COMMUNICATIONS (2018) 9:1671 DOI: 10.1038/s41467-018-04074-2 www.nature.com/naturecommunications | | | Cumulative effect Standardized logit model coefficients ΔAIC High-frequency High-frequency temperature temperature variability variability Depth Depth Heating Rate Heating Rate Background Conditions Acute Thermal Stress Thermal Trajectory Thermal Trajectory Cumulative Thermal Stress Cumulative Thermal Stress Acute Thermal Stress Exacerbates Mitigates Exacerbates Mitigates bleaching bleaching bleaching bleaching j = 3 j = 3 j = 2 j = 2 j = 1 j = 1 NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04074-2 ARTICLE Specific reef cases. Our results reveal the importance of high- a 1 = 0.76 °C frequency temperature variability at locations worldwide, but include reef-scale observations where such variability influences 0.5 More severe bleaching prevalence of corals in different locations of the same bleaching High-frequency 32,46 temp. variability reef during the same bleaching event . Here, we present two such case studies: one from Tahala Reef, a platform reef in the –2.0 –1.5 –1.0 –0.5 0 0.5 1.0 1.5 2.0 central Red Sea, and another from Nelly Bay, a fringing reef in the Δ DTR (°C) Great Barrier Reef in Australia (Fig. 6). These sites were chosen 32,56 due to the availability of additional meteorological data at 1 = 4.36m these reefs. At each location, temperature time series (Fig. 6a, b) for both a seaward and a shoreward location show that, whereas 0.5 low-frequency variations in water temperature are often very Depth similar over reef-scales (Fig. 6c, d), high-frequency variations may be quite distinct (Fig. 6e, f). In these cases, bleaching events were –10 –5 0 5 10 more widespread and severe at the seaward locations where DTRs Δ Depth (m) were smaller (Fig. 6a, b), consistent with our best logit models. 1 = 0.01 °C/day Discussion 0.5 More severe For corals, a shift in thermal tolerance can occur due to adap- bleaching Heating rate tation of the coral animal or algal symbionts through natural 57,58 selection of heat-tolerant lineages , or physiological acclima- –0.06 –0.04 –0.02 0 0.02 0.04 0.06 tion through the expression of heat shock proteins and regulation 23,36 Δ Rate of temperature change, ROTC ( °C/day) of apoptosis (i.e., programmed cell death) . As discussed, SS recent work highlights the importance of short-term temperature history (daily-weekly periods) for coral acclimatization to higher 1 = 0.49 temperatures , such that corals subject to warmer than average 0.5 temperatures prior to thermal stress may exhibit a greater toler- Acute 23 ance to acute temperature stress . In the context of these studies thermal stress and in keeping with other site-specific and experimental stu- 19,26,45,47,59 dies , our results suggest that temperature fluctuations –2.0 –1.5 –1.0 –0.5 0 0.5 1.0 1.5 2.0 on daily or tidal timescales are often sufficient to expose corals to Δ Acute1 , [–] temperatures high enough to encourage greater tolerance (via acclimation or adaptation) to thermal stress, but for time periods 1 = 0.64 19,48 short enough to avoid mortality . Further, our results establish 0.5 that the resistance of corals located in areas of high-frequency Thermal temperature variability to bleaching occurs in reef regions trajectory throughout the world. While we lack sufficient species-level data, we fully acknowledge that intrinsic coral properties , differences –3.0 –2.0 –1.0 0 1.0 2.0 3.0 in reef-scale community compositions , and taxonomic sus- Δ Thermal trajectory, TT [–] ceptibility are likely to influence heterogeneous reef-scale 1 bleaching responses, and an improvement to our model frame- 1 = 0.43 °C-weeks work would include species-level covariates. Our results also demonstrate the potential to both improve 0.5 predictions of bleaching and anticipate heterogeneous patterns in Cumulative bleaching prevalence at reef-scales, considering that beyond mean thermal stress values, accounting for variability in temperature regimes yields –2.5 –2.0 –1.5 –1.0 –0.5 0 0.5 1.0 1.5 2.0 2.5 better predictions of organismal responses to anomalous envir- Δ DHW (°C-weeks) 60 onmental events . Although SSTs from satellite remote sensing are not yet available at the spatiotemporal resolution required to Fig. 4 Influence of each in situ covariate on bleaching. Using the covariates calculate reef-scale high-frequency temperature variability, from the highest-ranked logit model, the probability of observing bleaching observational work on a range of reef structures suggests that it prevalence greater than the jth category is plotted against changes in each may be possible to predict reef-scale thermal environments using covariate from their respective mean values (where 0 corresponds to the relatively simple hydrodynamic models given readily available mean value), while keeping all other covariates at their mean values. bathymetry and basic hydrographic data such as tidal range, wave Bleaching prevalence categories are defined as 1: ≤10%; 2: 10−25%; 3: 25 30,32,43 height, and offshore mean SST . While we did not assess −50%; 4: >50% of reef area bleached. Highest-ranked model covariates other biogeochemical parameters, it is worth noting that the same include: (a) High-frequency temperature variability (DTR ), (b) Depth, physical circulation that drives spatially variable thermal envir- (c) Heating Rate (ROTC ), (d) Acute Thermal Stress (Acute1), (e) SS onments, in locations with an active benthic community, can also Thermal Trajectory (TT), and (f) Cumulative Thermal Stress (DHW ). 61,62 create dynamic oxygen, pH, and nutrient environments . Standard deviations for each covariate within our data set are also indicated Urgent global efforts at reducing anthropogenic greenhouse gas emissions must remain a priority for reef preservation, due to the SST-based OLR models indicate that upon consideration of the acute thermal stress that is now arising from global warming 5,63 consistent importance of DTR to bleaching, a globally available projections on reefs . However, combating the effects of local remotely sensed metric for diurnal temperature variability would stressors on reefs through conservation tools such as marine be valuable for improved bleaching predictions. protected areas is likely to increase the chances of reef persistence NATURE COMMUNICATIONS (2018) 9:1671 DOI: 10.1038/s41467-018-04074-2 www.nature.com/naturecommunications 7 | | | j = 1 j = 1 j = 2 j = 1 j = 1 j = 2 j = 3 j = 3 j = 2 j = 2 j = 3 j = 3 Probability ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04074-2 2.82 ± 0.10 2.77 ± 0.13 1.21 ± 0.05 1 Cumulative Background Acute Thermal Conditions Thermal Stress Stress (MMM ) (Acute1) (DHW ) 4 km 4 km High-frequency Depth Heating temperature (in situ) Rate variability (ROTC ) 90–4 km (DTR ) −1 −0.80 ± 0.05 −1.57 ± 0.05 −2 −1.89 ± 0.07 −3 (In situ): 37.29 (In situ): 35.02 (SST): 0.00 (SST): 0.00 Total: 37.29 ± 3.49 Total: 35.02 ± 4.06 (In situ): 0.00 (SST): 14.45 Total: 14.45 ± 1.68 (In situ): −0.34 (SST): 0.00 Total: −0.34 ± 0.59 −10 (In situ): −1.83 (SST): −8.75 −20 Total: −10.58 ± 2.09 (In situ): −14.61 (SST): −2.10 Total: −16.71 ± 1.84 −30 (In situ): −10.11 (SST): −14.00 Total: −24.11 ± 3.55 −40 Fig. 5 Remotely sensed SST OLR results. a Parameter estimates for standardized model coefficients of the covariates used in the highest-ranked OLR model when weekly 4 km CoRTAD SST-based variables are added to the pool of possible covariates; standard error bars were computed from delete-1 jackknife resampling. b The summation of the standardized covariate coefficients grouped by category from the highest-ranked models when including CoRTAD SST-based covariates; the standard error bars shown are obtained by summing in quadrature the individual standard errors computed after delete-1 jackknife resampling. The remotely sensed SST-based covariate contribution to each Cumulative Effect is colored gray through future warming, as well as to facilitate recovery following additional information for each time series are listed in Supplementary Data 1, 2,52,64,65 which also lists the source for each time series. Water temperature records included a bleaching event . Warming ocean temperatures are in this analysis spanned at least 12 months in duration, with a sampling interval projected to result in annual severe bleaching regimes by the less than or equal to 3 h. In cases where instrument substitution resulted in varying middle of this century, with spatial variability in the onset of these sampling intervals, time series were sub-sampled to the largest of these intervals, or, events on the order of ±10 years . Considering our results in the in rare cases, interpolated to remain below a 3 h sampling interval. The tempera- ture time series data used in this study originate from instruments that were context of this inevitable and persistent acute thermal stress, calibrated using varying methodologies, such as by placing loggers together at one focusing management efforts on the more resistant reef locations location and comparing recorded temperatures with a reference temperature that also experience delayed onsets in annual severe bleaching dataset, ice bath calibration, or multiplying raw field recorded temperatures by would maximize the likelihood that at least some healthy reefs normalized logger calibration coefficients. Our analysis is largely based on relative temperature variations, however, so that absolute temperature accuracy will not will exist in the future. affect the results presented here. In an effort to examine how representative our sample temperature time series data was, we used one-sample t tests to compare the overall means and extremes of our data to a global temperature data set taken Methods from nearly 1000 reef locations . From these tests, we cannot conclude any sig- Data synthesis. Water temperature time series from 118 locations, representing nificant differences (α = 0.05) between our time series data and this larger global five major ocean basins where tropical coral reefs are found (Western Indian data set. All data analyses were done using MATLAB 7.14 (The Mathworks, Natick, Ocean, Pacific Ocean, Caribbean Sea, Great Barrier Reef, and Red Sea), were MA, USA). obtained from existing records of in situ temperature data. Many of these time series were obtained directly from the researchers (Supplementary Data 1), or from publicly available databases including the Australian Institute of Marine Science Spectral analysis. Power spectral density (PSD) estimates were computed for each (AIMS), the National Data Buoy Center (NDBC), and the Florida Institute of temperature time series. First, if necessary, temperature time series were resampled Oceanography (FIO). However, time series within our dataset were also selected to or linearly interpolated to maintain a constant sampling interval, chosen to be 3 h match precise reef locations that had sufficiently documented bleaching events, so as to resolve spectral frequencies of up to 4 cpd. In order to examine tem- while also containing as long and consistently sampled records as possible. Site perature variability for a broad range of frequencies spanning annual to diurnal names and three-letter codes, locations, depths, instrument descriptions, and and shorter periods, a PSD was calculated as follows: time series greater than or 8 NATURE COMMUNICATIONS (2018) 9:1671 DOI: 10.1038/s41467-018-04074-2 www.nature.com/naturecommunications | | | Cumulative effect Standardized model coefficients High-frequency temperature variability Depth Heating Rate DTR Distribution Shape Background Condition Acute Thermal Stress Cumulative Thermal Stress NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04074-2 ARTICLE a b Tahala Nelly Bay % Bleaching prevalence Shoreward, 0.9 m Shoreward, 2.4 m Seaward, 0.9 m Seaward, 5.1 m 11% 50−60% % Mortality 33% 99% c d 1000 10,000 Local MMM 35 0 20 25 30 24 26 28 30 32 34 Biweekly average, °C 30 Local Biweekly average, °C 30 MMM e f 6 6 4 4 0 0 0 2 4 6 0 1 2 3 4 Daily temperature range, °C Daily temperature range, °C 2 2 0 0 −2 −2 Fig. 6 Same-reef case studies. a Temperature time series taken from the wave-exposed (blue) and wave-protected (red) edges of the Tahala reef platform in the Red Sea, which are separated by ~200 m. The percentage of observed mortality, which was associated with the bleaching event in September 46,84 2010 , is indicated for each corresponding platform edge. c 2-week low-pass filtered time series of the raw Tahala temperature data, with the Maximum Monthly Mean (MMM) temperature calculated using the in situ data for each time series. e 33-h high-pass filtered time series of the raw Tahala temperature data, with a histogram of the Daily Temperature Range (DTR) values for each time series. b, d, f Analogous versions of a, c, and e, respectively, but for the Nelly Bay reef flat (red) and reef slope (blue) habitats in the Great Barrier Reef, separated by ~122 m. Bleaching prevalence as proportions of belt transects are also indicated in b. The gray bars highlight the approximate periods of reported bleaching events equal to 10 years in duration were divided into 3-year sections, while all others Bleaching data synthesis. Coral bleaching observations (81 events) that corre- were divided into 4-month sections. Sections were overlapped by 50% and wind- sponded in time and location with temperature data (Supplementary Data 2) owed with a Hamming function before spectra for each section were calculated, were obtained from a variety of sources, but primarily from peer-reviewed pub- which were then ensemble averaged to obtain the PSD estimate. The statistical lications. Bleaching reports were based upon various quantification schemes, significance (α = 0.05) of observed spectral peaks was ascertained by comparison but some common methods included recording bleaching as: (i) a proportion of 67 69 with the upper confidence level of a background red noise fit to the spectrum .We transect area , (ii) severity among different colonies of different coral empirically defined annual, seasonal, and diurnal frequency bands as 0.00185 to species , and (iii) prevalence categories based on aerial (and ground-verified) 0.0111, 0.0119 to 0.143, and 0.727 to 1.333 cycles per day (cpd), respectively. Band surveys . In publications that reported percentages of colony bleaching within variance was computed using trapezoidal integration of PSD values within each different coloration and paling categories, we used the weighted average of the respective frequency range, and ratios of high-frequency to seasonal band variances different percentages within each category, though the majority of bleaching among habitat types (“Back Reef”, “Reef Flat”, “Reef Slope”) were compared using a records naturally translated into our four bleaching prevalence scores. To aggregate Kruskal−Wallis test (Supplementary Fig. 2b). or standardize the bleaching reports for use in this study, we defined a bleaching response variable in terms of percentage of spatial area bleached, and we therefore assigned the following categories as ordinal values of bleaching pre- valence score: 1: “bleaching prevalence” ≤ 10% (n = 48); 2: 10% < “bleaching pre- Spatiotemporal variability in water temperature. To quantify the magnitude of valence” ≤ 25% (n = 5); 3: 25% < “bleaching prevalence” ≤ 50% (n = 6); 4: diurnal patterns of heating and cooling, a DTR was calculated as the difference “bleaching prevalence” > 50% (n = 22). Note that bleaching events of 0 (reports of between maximum and minimum temperatures for each day of each time series. no bleaching observed or negligibly mild paling) are binned into bleaching pre- Temporal variations in DTR values were examined in multiple ways. First, DTRs valence score 1, to provide a conservative grouping for mild bleaching events, were composite-averaged for yeardays 1−366, and based on evident seasonal DTR 71 seasonal patterns of discoloration and variation in zooxanthellae densities , variability from panel b of Fig. 2, a non-parametric Kruskal−Wallis test was used 72 unresolved “background bleaching” levels , and observation errors. As opposed to to assess the seasonal dependence of DTR distributions . Seasons were defined continuous interval variables, ordinal variables represent categorical values that can such that each season spanned 3 complete months, and austral and boreal summers be ranked and have a natural ordering to them, offering the advantage of creating were December through February and June through August, respectively. Tem- bins for ranges of values. perature metrics are summarized in Table 1. Large-scale spatial patterns in water temperature variability relative to latitude were also characterized; annual tem- perature ranges, calculated as the range of monthly mean temperatures, were Explanatory variables. To assess the influence of explanatory environmental compared against latitude for all sites. Variance in high-frequency (33 to 6 h variables, including high-frequency temperature variability, on bleaching response, periods) and seasonal (7 to 84-day periods) spectral bands were computed via we performed a multivariate statistical analysis of the observed bleaching events. As integration of their respective power spectral densities, and the ratios of high- there are multiple aspects of thermal stress and environmental conditions that may frequency to seasonal variance for all time series were compared by the habitat explain the bleaching response, we selected 20 experimental variables for the from which each time series was recorded. Habitats were divided into three groups: ordinal regression (defined in Table 1), organized into eight broad categories: (1) (i) back reefs and back reef lagoons (labeled as “BR”), (ii) reef flats and non-back Depth, (2) Background Oceanographic Conditions, (3) Cumulative Thermal Stress, reef lagoons (“RF”), and (iii) reef crests, reef slopes, forereefs, and anything further (4) Acute Thermal Stress, (5) Thermal Trajectory, (6) Heating Rate, (7) High- offshore (“RS”). Frequency Temperature Variability, and (8) Shape of the Distribution of HF NATURE COMMUNICATIONS (2018) 9:1671 DOI: 10.1038/s41467-018-04074-2 www.nature.com/naturecommunications 9 | | | Feb10 Mar10 Apr10 May10 Jun10 Jul10 Aug10 Sep10 Oct10 Nov10 Dec10 Jan11 Feb11 Mar11 Apr11 May11 Jan96 Jan97 Jan98 Jan99 Jan00 Jan01 Jan02 Jan03 Jan04 Jan05 Temperature, °C Temperature, °C Temperature, °C Count Count Temperature, °C Temperature, °C Temperature, °C Count Count ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04074-2 Here i indexes each of N observations, with bleaching observation y , and the left- Temperature Variability. We lacked sufficient habitat (i.e., reef flat, reef crest, reef i hand side quantity is referred to as the logit of the probability of observing slope, etc.) information at which many loggers were placed, and therefore did not bleaching prevalence score j or lower, for j = 1, 2, or 3 (observations with bleaching include habitat as an explanatory variable. Environmental variables in the prevalence scores of 4 contribute to the regression through calculation of the log- Cumulative Thermal Stress category were calculated from in situ temperature odds). Note that the odds are defined as the ratio of the probability of an event measurements rather than the remotely sensed National Oceanic and Atmospheric occurring to the probability of the event not occurring, which is exactly the ratio Administration’s Coral Reef Watch (NOAA CRW) products because many of our inside the natural logarithm. Each C is an MLE-computed model intercept, and logger locations were outside of the coverage areas of the CRW Virtual Stations, each B is the MLE coefficient corresponding to each standardized independent and furthermore, as the first suite of CRW products was released in 2000, they do variable z , for k = 1, …, p, where p is the number of independent variables used in ik not include temperature data corresponding to the 1998 global bleaching event. a given model. A fundamental component of this model is the assumption of The explanatory variables we use here are commonly appended with a subscript to proportional odds, or parallel regression, which implies B values are independent indicate the period of time, relative to the bleaching observation date, used to of the logit level j. The validity of this parallel regression assumption was calculate the metric. For example, the subscript “SS” denotes the spring to summer ascertained using Brant’s Wald test , as well as a likelihood ratio test (α = 0.05). period, and the subscript “30” denotes the 30 days preceding (and including) a As each time series was the sole source of its explanatory variables, we can bleaching observation. expect many of these variables to be correlated with each other (multicollinear). If left unaccounted for, multicollinearity obscures the interpretation of the Principal components analysis. Principal components analysis was performed to explanatory variables and their coefficients, and may decrease the statistical power examine the spatial structure of the environmental forcings and determine the 78 of the logit analysis . The degree of multicollinearity among the 20 explanatory association between independent variables within each of the eight categories variables used in the logistic regression was assessed by calculating the condition (Fig. 1). The first two PC axes accounted for 44.2 and 18.8% of the variance within 79 indices and variance-decomposition proportions of the matrix of explanatory the matrix of independent variables. The magnitude and orientation of the loading variables. This revealed the following nontrivial multicollinearities: (1) DHW vectors indicate the importance of each parameter in describing the variance of the with CSA ; (2) DTR with DTR ; (3) MMM with MMM ; and (4) During 90 30 Total PCA components. DTR with DTR and DTR (Table 1). Note that within each of these four Total SS FW multicollinear groupings, the multicollinear variables come from the same explanatory variable category (for example, both DHW and CSA are both Computation of thermal trajectory and acute stress variables. An ordinal- 90 During from the Cumulative Thermal Stress category). Therefore, for each logit model, we valued Thermal Trajectory variable was included as an independent variable to selected at most one variable from a category without multicollinearity erroneously assess the degree to which environmental conditions confer thermal tolerance. The influencing our results. Spatial autocorrelation within each covariate was calculation of this was as follows: first the MMM and MMM+ 2 °C (the latter determined by calculating Moran’s I and examining correlograms, from which we quantity is referred to as the “local bleaching threshold”) were computed for a determined significant spatial autocorrelation typically up to distances of 500 km given time series. Then the 33-h low pass filtered time series for the 90 days (Supplementary Fig. 5). To account for this, we also added a random effect to the preceding a bleaching event was inspected to determine the type of Thermal highest-ranked logit models, from which we failed to conclude a marginal model Trajectory. If temperatures exceeded the MMM, then fell below the MMM for a improvement (see below). required 10-day “recovery period” before proceeding to exceed the local bleaching All permutations of all possible explanatory variables were used to compute a threshold, a Protective Trajectory with an ordinal value of 1 was recorded. If total of 10,367 logit models, where all logit models were computed using a temperatures increased from below the MMM to above the local bleaching multinomial logistic regression function (“mnrfit”) in MATLAB. Model threshold, without a 10-day recovery period, a Single Bleaching Trajectory with a comparison was performed using a bias correction for small sample sizes to value of 2 was recorded. If temperatures exceeded the local bleaching threshold at Akaike’s Information Criterion, AIC , and all models within ΔAIC ≤ 2 of the least twice, with a required 9-day recovery period between threshold exceedances, a C C best model (ΔAIC = 0), which have statistically indistinguishable performances , Repetitive Bleaching Trajectory with an ordinal value of 3 was assigned. Finally, if C are presented in Supplementary Fig. 6. McFadden’s pseudo-R was also computed temperatures did not exceed the local bleaching threshold, an ordinal value of 0 for the highest ranked models, and ranged from 0.26 to 0.30, with that of the corresponding to no thermal stress was assigned. Justification of our ordinal-value highest-ranked model equal to 0.30. While the logit model with the lowest AIC ,as scheme comes from an analysis of experimentally heated corals from the Great C well as all models within ΔAIC ≤ 2, provide a general outline for the coefficients of Barrier Reef, whereby in the face of thermal stress, corals with a Protective Tra- C the critical independent variables in explaining bleaching prevalence, these jectory experienced localized cell death of approximately 30%, while that of corals parameter estimates have errors of unknown distribution. Additionally, the under Single and Repetitive Bleaching Trajectories was approximately 60% and possible existence of outliers in the high-frequency temperature variability data 70%, respectively . Furthermore, experimental results showed that corals under a may influence the results of the logit parameter estimates, and therefore, delete-1 Protective Trajectory maintained significantly greater symbiont density than those jackknife resampling was used to compute estimates of bias and standard errors. under Single and Repetitive Bleaching Trajectories, and hence the ordinal scores We found a slight positive bias that did not significantly alter the influence of high- from 0 to 3 aptly account for the monotonic nature of coral tissue detriment frequency temperature variability relative to the other covariates. Estimates of associated with no heat stress, a Protective Trajectory, a Single Bleaching Trajec- standard errors for all logit parameters of all models within ΔAIC ≤ 2 can be seen tory, or a Repetitive Bleaching Trajectory, respectively. The “No Thermal Stress” C in standard error bars (Fig. 3b and Supplementary Fig. 6). Specifically, the error trajectory corresponded to ~0% cell death and greater symbiont density than the bars seen in Fig. 3c are the standard errors from each contribution summed in Protective Trajectory. quadrature. A modified jackknife resampling scheme was also performed, in which The Acute Thermal Stress category was composed of two binary variables to instead of leaving out one site at a time, sites within 10 km of each other were indicate the presence/absence of acute thermal stress. The calculation of this was as grouped together and each of these proximity groups was left out incrementally follows: first the MMM+ 1 °C and MMM+ 2 °C were computed for a given time before fitting OLR models to the remaining data. This spatial resampling analysis series, then the daily mean temperatures within 90 days before a bleaching event (Supplementary Fig. 7) did not result in significantly different parameter estimates were inspected to determine if temperatures exceeded MMM+ 1 °C (in which case than the full model presented in Fig. 3b. Estimates for the intercept terms C , C , Acute1 would equal 1) and MMM+ 2 °C (in which case Acute2 would equal 1). 1 2 and C (±non-resampled standard errors) were found to be 0.72 ± 0.35, 1.41 ± 0.38, and 2.00 ± 0.42, respectively, indicating no significant difference between bleaching Ordinal logistic regression. Using the eight explanatory variable categories prevalence categories within our dataset. Using the statistical computing software described above, and bleaching prevalence scores from 1 to 4 as the response R (https://www.r-project.org), a random effect grouping variable was added to variable, we performed an ordinal logistic regression analysis to determine how the each of the highest-ranked models, which grouped reefs within 5 km of each other. relative log odds of a given bleaching prevalence score depends on the interactions These resulting mixed effects models were compared to their fixed effects among the explanatory variables. Ordinal regression models have been adeptly equivalents to determine model fit improvement, and the inclusion of a random used in ecological studies where data are often present as semi-quantitative effect did not improve the fit of any of the 20 highest-ranked models 73,74 variables in which relative differences between values are of importance . (Supplementary Data 3). Furthermore, to account for possible nonlinear Furthermore, logit functions have been previously implemented to predict the interactions between covariates, such as Depth × DTR or Acute1 × DHW ,we 30 30 72,75 presence/absence of bleaching using gridded remote-sensed data , or to explain also included a nonlinear interaction term to the highest-ranked model and the influence of a range of environmental and coral physiological factors on reef determined model fit improvement. The AIC of these nonlinear models ecosystem response following a disturbance . These logit models are multivariate (Supplementary Table 1) indicated that the addition of a nonlinear interaction term extensions of generalized linear regression models , providing parameter estimates did not significantly improve the fit of the original highest-ranked model displayed via maximum likelihood estimation (MLE) to model the relative log odds of, for in Fig. 3b. Therefore, as our main result, we ultimately report the fixed effects ORL our purposes, observing one bleaching prevalence score or less versus observing the model parameter estimates with no nonlinear or interaction terms (Fig. 3b). remaining greater bleaching prevalence scores: The SST-based OLR analysis summarized in Fig. 5 was performed using the CoRTAD 4 km weekly SST data pixels that were closest to our in situ bleaching PyðÞ
j 15 observations. The quantities MMM, MMM , Acute1, Acute2, DHW, ROTC , Max SS ln ¼ C þ B z þ þ B z : ð1Þ j 1 i1 p ip PyðÞ >j and ROTC (Table 1) were computed from the CoRTAD temperature data. This i 90 resulted in a total of 27 covariates (20 as described above, and 7 new ones computed from the CoRTAD data). We then proceeded to fit OLR models using all 10 NATURE COMMUNICATIONS (2018) 9:1671 DOI: 10.1038/s41467-018-04074-2 www.nature.com/naturecommunications | | | NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04074-2 ARTICLE permutations of covariates, with the constraint that within each model, we only 18. DeCarlo, T. M. et al. Mass coral mortality under local amplification of 2 include ≤1 covariate from each category listed in Table 1. This produced 60,479 degrees C ocean warming. Sci. Rep. 7, 44586 (2017). models, from which we performed model comparison using AIC , resulting in 12 19. Oliver, T. & Palumbi, S. Do fluctuating temperature environments elevate dominant models being identified as statistically indistinguishable. From here, we coral thermal tolerance? Coral Reefs 30, 429–440 (2011). concluded that when 4 km SST metrics are incorporated into our model 20. Guadayol, O., Silbiger, N. J., Donahue, M. J. & Thomas, F. I. Patterns in framework, high-frequency temperature variability and depth remain the most temporal variability of temperature, oxygen and pH along an environmental influential covariates on bleaching prevalence, acting to attenuate bleaching gradient in a coral reef. PLoS ONE 9, e85213 (2014). (Fig. 5). Furthermore, acute and cumulative thermal stress exacerbate bleaching 21. Darling, E. S., McClanahan, T. R. & Côté, I. M. Life histories predict coral prevalence the most, and the SST-derived version of the latter contributes community disassembly under multiple stressors. Glob. Change Biol. 19, substantially to this effect. The results presented in Fig. 5 ultimately mirror those 1930–1940 (2013). from Fig. 3b, c, indicating that DTR is the main driver of bleaching prevalence, 22. Berkelmans, R. Time-integrated thermal bleaching thresholds of reefs and regardless of whether SST or in situ quantities of the other predictors are used. their variation on the Great Barrier Reef. Mar. Ecol. Progress Ser. 229,73–82 Common covariates in these models were MMM, MMM , MMM , DHW , Max 4km 30 (2002). DHW , Acute1, Acute2 , and ROTC , and each of these 12 highest- 4km 4km 90-4 km 23. Middlebrook, R., Hoegh-Guldberg, O. & Leggat, W. The effect of thermal ranked models included at least one SST-based covariate. The AIC values of these history on the susceptibility of reef-building corals to thermal stress. J. Exp. models ranged from 129.19 to 131.12, with McFadden’s pseudo-R values ranging Biol. 211, 1050–1056 (2008). from 0.35 to 0.38, which presents a considerable improvement to the highest- 24. Middlebrook, R., Anthony, K. R., Hoegh-Guldberg, O. & Dove, S. Heating ranked model that was fit to only in situ values, which had an AIC equal to 143.75 rate and symbiont productivity are key factors determining thermal stress in and a McFadden’s pseudo-R of 0.30. To highlight the importance of High- the reef-building coral Acropora formosa. J. Exp. Biol. 213, 1026–1034 Frequency Temperature Variability in model improvement, we also fit these (2010). models to the data and excluded the Daily Temperature Range covariate, which 25. Chollett, I., Enriquez, S. & Mumby, P. J. Redefining thermal regimes to design was either DTR or DTR . This significantly decreased model performance, with 30 90 reserves for coral reefs in the face of climate change. PLoS ONE 9, e110634 AIC values ranging from 153.46 to 169.57, and McFadden’s pseudo-R values (2014). ranging from 0.12 to 0.23, which supports the idea that DTR is the dominant 26. McClanahan, T., Maina, J., Moothien-Pillay, R. & Baker, A. Effects of driving variable for the bleaching response. geography, taxa, water flow, and temperature variation on coral bleaching intensity in Mauritius. Mar. Ecol. Progress Ser. 298, 131–142 (2005). Data availability. All 118 temperature time series used in this study are archived 27. Williams, G. J., Knapp, I. S., Maragos, J. E. & Davy, S. K. Modeling patterns of with the NOAA National Centers for Environmental (NCEI), with accession coral bleaching at a remote Central Pacific atoll. Mar. Pollut. Bull. 60, number 0170826 . 1467–1476 (2010). 28. Lesser, M. P. & Farrell, J. H. Exposure to solar radiation increases damage to Received: 29 June 2017 Accepted: 2 April 2018 both host tissues and algal symbionts of corals during thermal stress. Coral Reefs 23, 367–377 (2004). 29. Brown, B. E., Dunne, R. P., Goodson, M. & Douglas, A. Marine ecology: bleaching patterns in reef corals. Nature 404, 142 (2000). 30. Lentz, S. J. et al. The characteristics and dynamics of wave-driven flow across a References platform coral reef in the Red Sea. J. Geophys. Res.: Oceans 121, 1360–1376 1. West, J. M. & Salm, R. V. Resistance and resilience to coral bleaching: (2016). implications for coral reef conservation and management. Conserv. Biol. 17, 31. Leichter, J. J. & Miller, S. L. Predicting high-frequency upwelling: spatial and 956–967 (2003). temporal patterns of temperature anomalies on a Florida coral reef. Cont. Shelf 2. Hoegh-Guldberg, O. et al. Coral reefs under rapidclimate change and ocean Res. 19, 911–928 (1999). acidification. Science 318, 1737–1742 (2007). 32. Davis, K. A. et al. Observations of the thermal environment on Red Sea 3. Costanza, R. et al. Changes in the global value of ecosystem services. Glob. platform reefs: a heat budget analysis. Coral Reefs 30,25–36 (2011). Environ. Change 26, 152–158 (2014). 33. Carilli, J., Donner, S. D. & Hartmann, A. C. Historical temperature variability 4. van Hooidonk, R. et al. Local-scale projections of coral reef futures and affects coral response to heat stress. PLoS ONE 7, e34418 (2012). implications of the Paris Agreement. Sci. Rep. 6, 39666 (2016). 34. Palumbi, S. R., Barshis, D. J., Traylor-Knowles, N. & Bay, R. A. Mechanisms 5. Hughes, T. P. et al. Global warming and recurrent mass bleaching of corals. of reef coral resistance to future climate change. Science 344, 895–898 Nature 543, 373–377 (2017). (2014). 6. Grottoli, A. G., Rodrigues, L. J. & Palardy, J. E. Heterotrophic plasticity and 35. Castillo, K. D., Ries, J. B., Weiss, J. M. & Lima, F. P. Decline of forereef corals resilience in bleached corals. Nature 440, 1186–1189 (2006). in response to recent warming linked to history of thermal exposure. Nat. 7. Glynn, P. Coral reef bleaching: ecological perspectives. Coral Reefs 12,1–17 Clim. Change 2, 756–760 (2012). (1993). 36. Barshis, D. et al. Genomic basis for coral resilience to climate change. Proc. 8. Goreau, T. J. & Hayes, R. L. Coral bleaching and ocean “hot spots”. Ambio-J. Natl. Acad. Sci. USA 110, 1387–1392 (2013). Human Environ. Res. Manag. 23, 176–180 (1994). 37. Rogers, J. S., Monismith, S. G., Koweek, D. A., Torres, W. I. & Dunbar, R. B. 9. Jokiel, P. L. in Coral Health and Disease (ed. Rosenberg, E.) 401−425 Thermodynamics and hydrodynamics in an atoll reef system and their (Springer, New York, 2004). influence on coral cover. Limnol. Oceanogr. 61, 2191–2206 (2016). 10. Cantin, N. E. & Lough, J. M. Surviving coral bleaching events: porites growth 38. Schoepf, V., Stat, M., Falter, J. L. & McCulloch, M. T. Limits to the thermal anomalies on the Great Barrier Reef. PLoS ONE 9, e88720 (2014). tolerance of corals adapted to a highly fluctuating, naturally extreme 11. Dunn, S. R., Thomason, J. C., Le Tissier, M. D. & Bythell, J. C. Heat stress temperature environment. Sci. Rep. 5, 17639 (2015). induces different forms of cell death in sea anemones and their endosymbiotic 39. Ainsworth, T. D. et al. Climate change disables coral bleaching protection on algae depending on temperature and duration. Cell Death Differ. 11, the Great Barrier Reef. Science 352, 338–342 (2016). 1213–1222 (2004). 40. Swain, T. D. et al. Coral bleaching response index: a new tool to standardize 12. Berkelmans, R., De’ath, G., Kininmonth, S. & Skirving, W. J. A comparison of and compare susceptibility to thermal bleaching. Glob. Change Biol. 22, the 1998 and 2002 coral bleaching events on the Great Barrier Reef: spatial 2475–2488 (2016). correlation, patterns, and predictions. Coral Reefs 23,74–83 (2004). 41. Grottoli, A. G. et al. The cumulative impact of annual coral bleaching can turn 13. Goreau, T., McClanahan, T., Hayes, R. & Strong, A. Conservation of coral some coral species winners into losers. Glob. Change Biol. 20, 3823–3833 reefs after the 1998 global bleaching event. Conserv. Biol. 14,5–15 (2000). (2014). 14. Liu, G. et al. Reef-scale thermal stress monitoring of coral ecosystems: new 5- 42. Riegl, B. & Piller, W. E. Possible refugia for reefs in times of environmental km global products from NOAA Coral Reef Watch. Remote Sens. 6, stress. Int. J. Earth Sci. 92, 520–531 (2003). 11579–11606 (2014). 43. Lowe, R. J., Pivan, X., Falter, J., Symonds, G. & Gruber, R. Rising sea levels will 15. Donner, S. D. An evaluation of the effect of recent temperature variability on reduce extreme temperature variations in tide-dominated reef habitats. Sci. the prediction of coral bleaching events. Ecol. Appl. 21, 1718–1730 (2011). Adv. 2, e1600825 (2016). 16. Logan, C. A., Dunne, J., Eakin, C. & Donner, S. A framework for comparing 44. Herdman, L. M., Hench, J. L. & Monismith, S. G. Heat balances and thermally coral bleaching thresholds. In Proc. of the 12th International Coral Reef driven lagoon ocean exchangeson a tropical coral reef system (Moorea, French Symposium (eds Yellowlees D., & Hughes T. P.) (James Cook University, Polynesia). J. Geophys. Res.: Oceans 120, 1233–1252 (2015). 2012). 45. Craig, P., Birkeland, C. & Belliveau, S. High temperatures tolerated by a 17. Boylan, P. & Kleypas, J. New insights into the exposure and sensitivity of coral diverse assemblage of shallow-water corals in American Samoa. Coral Reefs reefs to ocean warming. In Proc. of the 11th International Coral Reef 20, 185–189 (2001). Symposium (eds Riegl, B. & Dodge, R. E.) 854−858 (CNSO, 2008). NATURE COMMUNICATIONS (2018) 9:1671 DOI: 10.1038/s41467-018-04074-2 www.nature.com/naturecommunications 11 | | | ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04074-2 46. Pineda, J. et al. Two spatial scales in a bleaching event: corals from the mildest 75. Yee, S. H. & Barron, M. G. Predicting coral bleaching in response to and the most extreme thermal environments escape mortality. Limnol. environmental stressors using 8 years of global-scale data. Environ. Monit. Oceanogr. 58, 1531–1545 (2013). Assess. 161, 423–438 (2010). 47. McClanahan, T. R., Ateweberhan, M., Muhando, C. A., Maina, J. & 76. McCullagh, P. Regression models for ordinal data. J. R. Stat. Soc. Ser. B Mohammed, M. S. Effects of climate and seawater temperature variation on Methodol. 42, 109–142 (1980). coral bleaching and mortality. Ecol. Monogr. 77, 503–525 (2007). 77. Brant, R. Assessing proportionality in the proportional odds model for ordinal 48. Castillo, K. D. & Helmuth, B. S. T. Influence of thermal history on the logistic regression. Biometrics 46, 1171–1178 (1990). response of Montastraea annularis to short-term temperature exposure. Mar. 78. Graham, M. H. Confronting multicollinearity in ecological multiple Biol. 148, 261–270 (2005). regression. Ecology 84, 2809–2815 (2003). 49. Davis, K. A., Leichter, J. J., Hench, J. L. & Monismith, S. G. Effects of western 79. Belsley, D. A., Kuh, E. & Welsch, R. E. Regression Diagnostics: Identifying boundary current dynamics on the internal wave field of the Southeast Florida Influential Data and Sources of Collinearity Vol. 571 (John Wiley & Sons, New shelf. J. Geophys. Res. 113, C09010 (2008). York, 2005). 50. Davis, K. A. & Monismith, S. G. The modification of bottom boundary layer 80. Hurvich, C. M. & Tsai, C. -L. Regression and time series model selection in turbulence and mixing by internal waves shoaling on a barrier reef. J. Phys. small samples. Biometrika 76, 297–307 (1989). Oceanogr. 41, 2223–2241 (2011). 81. Burnham, K. P. & Anderson, D. R. Model Selection and Multimodel Inference: 51. Leichter, J. J., Deane, G. & Stokes, M. Spatial and temporal variability of A Practical Information-Theoretic Approach (Springer, New York, 2002). internal wave forcing on a coral reef. J. Phys. Oceanogr. 35, 1945–1962 (2005). 82. R Core Team. R: A Language and Environment for Statistical Computing (R 52. Graham, N. A., Jennings, S., MacNeil, M. A., Mouillot, D. & Wilson, S. K. Foundation for Statistical Computing, Vienna, Austria, 2017). Predicting climate-driven regime shifts versus rebound potential in coral reefs. 83. Barshis, D. J. et al. Water temperature from thermistors deployed at coral reef Nature 518,94–97 (2015). sites in the Indian and Atlantic Oceans, and other locations from 1987-11-02 53. McClanahan, T. R. et al. Predictability of coral bleaching from synoptic to 2016-06-2. NOAA https://data.nodc.noaa.gov/cgi-bin/iso?id=gov.noaa. satellite and in situ temperature observations. Coral Reefs 26, 695–701 (2007). nodc:0170826 (2018). 54. Baker, A. C., McClanahan, T. R., Starger, C. J. & Boonstra, R. K. Long-term 84. Furby, K. A., Bouwmeester, J. & Berumen, M. L. Susceptibility of central Red monitoring of algal symbiont communities in corals reveals stability is taxon Sea corals during a major bleaching event. Coral Reefs 32, 505–513 (2013). dependent and driven by site-specific thermal regime. Mar. Ecol. Progress Ser. 85. Liu, G., Strong, A. E., Skirving, W. & Arzayus, L. F. Overview of NOAA 479,85–97 (2013). coral reef watch program’s near-real time satellite global coral bleaching 55. Casey, K. S. S. et al. The Coral Reef Temperature Anomaly Database monitoring activities. In Proc. of the 10th International Coral Reef Symposium (CoRTAD) Version 5—Global, 4 km sea surface temperature and related (ed. Suzuki, Y.) 1783-1793 (Japanese Coral Reef Society, 2006). thermal stress metrics for 1982–2012. NOAA https://data.nodc.noaa.gov/cgi- bin/iso?id=gov.noaa.nodc:0126774 (2015). Acknowledgements 56. Marshall, P. & Baird, A. Bleaching of corals on the Great Barrier Reef: The authors thank the National Data Buoy Center and the Florida Institute of Ocea- differential susceptibilities among taxa. Coral Reefs 19, 155–163 (2000). nography for the use of data from the Florida Keys, and the Moorea Coral Reef Long- 57. Berkelmans, R. & van Oppen, M. J. The role of zooxanthellae in the thermal Term Ecological Research site for South Pacific data. The authors thank the UC Irvine tolerance of corals: a ‘nugget of hope’ for coral reefs in an era of climate Data Science Initiative, as well as the UCI OCEANS Initiative, for funding to support this change. Proc. Biol. Sci. 273, 2305–2312 (2006). work. The authors also thank Professor Peter J. Edmunds and Professor Steven J. Davis 58. Rowan, R. Coral bleaching: thermal adaptation in reef coral symbionts. Nature for their assistance in writing this manuscript. This is CSUN Marine Biology contribu- 430, 742–742 (2004). tion #267. The authors acknowledge that A. Safaie and K. Davis were supported by 59. McClanahan, T. R., Ateweberhan, M., Omukoto, J. & Pearson, L. Recent National Science Foundation Award No. 1436254 and G. Pawlak was supported by seawater temperature histories, status, and predictions for Madagascar’s coral Award No. 1436522. reefs. Mar. Ecol. Progress Ser. 380, 117–128 (2009). 60. Sheldon, K. S. & Dillon, M. E. Beyond the mean: biological impacts of cryptic temperature change. Integr. Comp. Biol. 56, 110–119 (2016). Author contributions 61. Falter, J. L., Lowe, R. J., Zhang, Z. & McCulloch, M. Physical and biological K.A.D. and A.S. designed the study. T.R.M., G.P., D.J.B., J.L.H., J.S.R., G.J.W., and K.A.D. controls on the carbonate chemistry of coral reef waters: effects of metabolism, collected the data. A.S., K.A.D., and N.J.S. analyzed the data. A.S., K.A.D., N.J.S wrote the wave forcing, sea level, and geomorphology. PLoS ONE 8, e53303 (2013). paper, and all authors contributed significantly to the interpretation and editing of the 62. Gruber, R. K., Lowe, R. J. & Falter, J. L. Metabolism of a tide-dominated reef manuscript. platform subject to extreme diel temperature and oxygen variations. Limnol. Oceanogr. 62, 1701–1717 (2017). 63. Selig, E. R., Casey, K. S. & Bruno, J. F. Temperature-driven coral decline: the Additional information Supplementary Information accompanies this paper at https://doi.org/10.1038/s41467- role of marine protected areas. Glob. Change Biol. 18, 1561–1570 (2012). 64. Baker, A. C., Glynn, P. W. & Riegl, B. Climate change and coral reef bleaching: 018-04074-2. an ecological assessment of long-term impacts, recovery trends and future outlook. Estuar., Coast. Shelf Sci. 80, 435–471 (2008). Competing interests: The authors declare no competing interests. 65. Mellin, C., Aaron MacNeil, M., Cheal, A. J., Emslie, M. J. & Julian Caley, M. Marine protected areas increase resilience among coral reef communities. Reprints and permission information is available online at http://npg.nature.com/ Ecol. Lett. 19, 629–637 (2016). reprintsandpermissions/ 66. Kleypas, J. A., McManus, J. W. & Meñez, L. A. Environmental limits to coral reef development: where do we draw the line? Am. Zool. 39, 146–159 (1999). Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in 67. Warner, R. M. Spectral Analysis of Time-Series Data (Guilford Press, New published maps and institutional affiliations. York, 1998). 68. Green, D. H., Edmunds, P. J., Pochon, X. & Gates, R. D. The effects of substratum type on the growth, mortality, and photophysiology of juvenile corals in St. John, US Virgin Islands. J. Exp. Mar. Biol. Ecol. 384,18–29 (2010). Open Access This article is licensed under a Creative Commons 69. Elvidge, C. D. et al. Satellite observation of Keppel Islands (Great Barrier Reef) Attribution 4.0 International License, which permits use, sharing, 2002 coral bleaching using IKONOS data. Coral Reefs 23, 123–132 (2004). adaptation, distribution and reproduction in any medium or format, as long as you give 70. Baird, A. & Marshall, P. Mortality, growth and reproduction in scleractinian appropriate credit to the original author(s) and the source, provide a link to the Creative corals following bleaching on the Great Barrier Reef. Mar. Ecol. Progress Ser. Commons license, and indicate if changes were made. The images or other third party 237, 133–141 (2002). material in this article are included in the article’s Creative Commons license, unless 71. Fitt, W., McFarland, F., Warner, M. & Chilcoat, G. Seasonal patterns of tissue indicated otherwise in a credit line to the material. If material is not included in the biomass and densities of symbiotic dinoflagellates in reef corals and relation to article’s Creative Commons license and your intended use is not permitted by statutory coral bleaching. Limnol. Oceanogr. 45, 677–685 (2000). regulation or exceeds the permitted use, you will need to obtain permission directly from 72. Yee, S. H., Santavy, D. L. & Barron, M. G. Comparing environmental the copyright holder. To view a copy of this license, visit http://creativecommons.org/ influences on coral bleaching across and within species using clustered licenses/by/4.0/. binomial regression. Ecol. Model. 218, 162–174 (2008). 73. Schabenberger, O. The use of ordinal response methodology in forestry. For. Sci. 41, 321–336 (1995). © The Author(s) 2018 74. Guisan, A. & Harrell, F. E. Ordinal response regression models in ecology. J. Veg. Sci. 11, 617–626 (2000). 12 NATURE COMMUNICATIONS (2018) 9:1671 DOI: 10.1038/s41467-018-04074-2 www.nature.com/naturecommunications | | |