High frequency temperature variability reduces the risk of coral bleaching

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.

C oral reef ecosystems provide subsistence nutrition, coastal protection, and revenue from tourism to hundreds of millions of people globally 1,2 , and are valued at trillions of dollars annually 3 . Especially during recent years, coral reefs are increasingly threatened by accelerated rises in ocean temperatures owing to global warming [4][5][6] . Elevated seawater temperatures are the primary cause of mass coral bleaching 5,7 , or the loss of pigmentation due to the collapse of the symbiotic relationship between the coral host and its endodermal dinoflagellate algae (zooxanthellae) 7,8 . Bleached corals are susceptible to disease 9 and reduced carbonate accretion 9,10 , and prolonged bleaching will lead to mortality 5,11,12 .
Thermal stress on corals and regional bleaching events are most often predicted by the magnitude and duration of remotely sensed sea surface temperatures (SSTs) above a fixed, locally defined average summer threshold temperature 5,8,13 . A conventionally used metric for quantifying these temperature anomalies is provided by the National Oceanic and Atmospheric Administration's (NOAA) Coral Reef Watch program, which has reported cumulative thermal stress on reefs twice a week since 1997 14 . Furthermore, bleaching predictions from remotely sensed temperatures can be improved by including SST-based calculations of interannual temperature variability 15,16 and coral sensitivity to thermal stress exposure 17 . However, the relatively coarse spatiotemporal resolution of the remotely sensed data prevents ensuing thermal stress quantifications from identifying the often observed significant spatial heterogeneity in bleaching that occurs within reef regions and individual reefs [18][19][20] . The response of reefs to temperature at these smaller spatial scales is complex and putatively depends on a combination of organism-level and reefscale factors such as coral life-history strategies and stressor cotolerances 21 , the history and duration of thermal stress exposure 22,23 , the rate of change in seawater temperature 24,25 , flow conditions 26 , heterotrophic feeding 6 , turbidity 27 , and the intensity and history of exposure to solar radiation 28,29 . In turn, many of these environmental conditions are mediated by reef-scale factors such as waves 30 , winds, tides 31 , and daily heating and cooling 32 .
Site-specific studies suggest that historical temperature variability within diurnal time scales affects corals' physiological tolerance 19,26,33,34 and performance 35 under thermal stress. For example, it has been theorized that corals located in areas characterized by large temperature fluctuations, such as reef flats or shallow lagoons, may be better acclimatized or adapted to thermal stress, and therefore more resistant to anomalous temperatures and bleaching, than corals in areas where temperatures are more stable, such as on reef crests or reef slopes [36][37][38] . Other studies have suggested that water temperatures in the weeks or months leading up to peak temperatures are critical in determining the coral physiological response. A recent analysis of experimentally heated corals from the Great Barrier Reef showed that bleaching and cell death responses were indeed lower when the thermal exposure included a moderate pre-stress followed by a short recovery period (i.e., a "protective temperature trajectory") 39 . Depending on intrinsic properties of coral physiology 40 such as energy reserves and algal phenotypic plasticity 41 , pre-peak temperatures may either protect against or exacerbate bleaching at peak temperatures 41 . Taken together, a growing body of evidence thus suggests that historical temperature variability, and particularly, "high-frequency" temperature variability, which we define as occurring within diurnal or shorter periods, may play an important role in determining corals' physiological responses to thermal stress and thereby reef-scale vulnerability to bleaching. In turn, a better understanding of reef-scale bleaching risk factors could help coastal management efforts to identify natural refugia and may be important for the recovery of coral communities following a bleaching event 42 .
Here, using a global suite of in situ data, we compare and assess the ability of 20 commonly used environmental variables and 7 remotely sensed variables to explain observed bleaching prevalence, testing the hypothesis that including high-frequency temperature variability as one of these model variables will lead to more accurate predictions. Analyzed data include records of in situ temperature time series at 118 reef locations from five major reef regions with sampling intervals of ≤3 h and continuous measurements of ≥1 year, as well as precise information on habitats and depths (Supplementary Data 1), along with 81 spatially and temporally coincident, quantitative coral bleaching observations (Supplementary Data 2). Each of the 81 bleaching observations was matched to its own spatiotemporally coincident temperature time series data, such that 46 of the 118 temperature time series were used in the subsequent bleaching analysis. Bleaching observations, which are most often reported as the average percent of colony or transect area bleached, were standardized to ordinal-valued "bleaching prevalence scores" (1: ≤10%; 2: 10−25%; 3: 25−50%; 4: >50% of reef area bleached), representing mild to pervasive bleaching, respectively (Methods). The influence of different factors on bleaching prevalence scores are evaluated by selecting covariates from a pool of 20 explanatory variables (depth, latitude, and 18 thermal metrics) grouped into 8 categories of metrics often used to predict bleaching (Table 1). In addition to these in situ variables, we also include 7 analogous and conventional remotely sensed SST thermal stress metrics (Table 1). After standardizing all covariates and fitting them to ordinal-valued bleaching prevalence scores using ordinal logistic regression (OLR) models (Methods), we conclude that high-frequency temperature variability, specifically the average daily temperature range (DTR) of the 30 days preceding a bleaching observation, is the most influential covariate in predicting the bleaching response, and serves to attenuate the prevalence of bleaching.

Results
Variation among in situ explanatory variables. A principal components analysis (PCA) displays the projection for each site onto the 2D plane that accounts for the most variance in the 20 in situ explanatory variables (Fig. 1), and the locations of the loading vectors reveal how these explanatory variables relate to their respective groupings. The first principal component accounts for 44.2% of the variation in the explanatory variables, and is largely driven by high-frequency temperature variability and cumulative thermal stress.
Spatiotemporal dependence of diurnal temperature variability. The thermal metrics computed from temperature time series were highly variable across sites, but regardless of location and depth, all 118 time series show significant temperature variations in the high-frequency band (Supplementary Note 1; Supplementary  Fig. 1), which we define as 0.727−4 cycles per day (cpd). Power spectra of temperature variations were calculated for each location, and the ratios of high-frequency band to seasonal band (0.012−0.143 cpd, or 1/7 to 1/84 days) variance in these spectra were used to characterize the relative importance of variance within the high-frequency band. This ratio correlates with the inverse of depth (r = 0.381, Student's t-test p < 0.05), indicating that the relative contribution of high-frequency variability to the variance within a temperature time series is stronger at shallower sites ( Supplementary Fig. 2a). At back reef, reef flat, and reef slope habitats, these ratios were on average 1.83, 0.68, and 0.44, respectively, while across all locations, this ratio was 1.02 (Supplementary Fig. 2b). Furthermore, these ratios differed significantly among the three habitats (Kruskal-Wallis, χ 2 = 24.66,  Supplementary Fig. 2b). Although the magnitude of diurnal temperature fluctuations varies by location, the ubiquity and prominence of temperature variance in this frequency band indicated by these average ratios reflects the importance of some common physical forcing processes governing the flow and heating of reef waters, such as diurnal solar heating, tides, winds, and waves 31,32,43,44 . Power spectra of six representative time series from different reef regions (Fig. 2a) show a broad range of temperature variability from annual to hourly periods (see Supplementary  Fig. 1 for other spectra). Yearly composites of mean water temperature and DTR ( Fig. 2b) both show prominent seasonal cycles (Supplementary Note 2): the magnitude of daily temperature fluctuations was seasonally dependent (Kruskal-Wallis, p < 0.01) for 96% of reefs in our study (113 of 118 time series), with maximum DTRs occurring most often in spring and summer months (74% of time series, Supplementary Fig. 3), and minimum DTRs occurring most often in fall and winter months (also 74% of time series, Supplementary Fig. 3). On global scales (~10 3 km), latitudinal gradients in solar forcing drive variations in seasonal temperature patterns on reefs ( Supplementary Fig. 4), but there is also considerable heterogeneity in thermal environments at reef-scales (~10 2 m) due to variation in depth and circulation 32,42,45 . The differences in thermal environments at reef-scales are often greatest in the high-frequency band (daily and tidal timescales; Fig. 2c). Dramatically different thermal environments can be found at locations separated by 10s or 100s of meters on a reef, as illustrated by 7-day temperature time series from various locations on the same island, or different habitats within a given reef (Fig. 2c). For example, during a week in November 2009, two sites in American Samoa that are separated by <2 km and at similar water depths experienced average DTRs of 1.78 and 0.51°C (Fig. 2c, sites OF3 and OF5 respectively).
Differences in the distributions of DTRs that distinguish microclimates within a reef system (e.g., thermally variable shoreward locations or thermally stable seaward ones) are reflected in the mean, skewness, and kurtosis of DTR values (Fig. 2d). Shallower and more shoreward sites have a peak in their DTR distributions corresponding to a larger DTR value, and furthermore, their distributions take on more extreme values than those from sites in deeper and more seaward locations. For example, at Heron Island in the Great Barrier Reef, the mean DTR of 4.23°C on the reef flat was over three times as large as that of the reef slope (Fig. 2d). The implications of these different thermal microclimates for resistance to thermal stress and resilience to bleaching are discussed below.
The effect of diurnal temperature variability on bleaching.
Ordinal logistic regression ("logit") models were computed for all permutations of selecting at most one variable from each of the eight categories in Table 1 (a total of 10,367 models), with bleaching prevalence scores as the response variable. Corrected Akaike's Information Criterion (AIC C ) values were used to rank the logit models, where the model with the lowest AIC C value was ranked the highest (Fig. 3b). The model coefficients indicate the association of tested variables with bleaching prevalence score, such that positive coefficients indicate a "mitigating" effect on bleaching prevalence, and negative coefficients an "exacerbating" effect on bleaching prevalence.
"High-Frequency Temperature Variability" (Table 1) was used to capture temperature variability on diurnal and shorter time periods, a metric that is important for characterizing differential reef-and habitat-scale microclimates 19,32,46 . In the best model ( Fig. 3a), high-frequency temperature variability, specifically the average DTR over the 30 days preceding a bleaching event (DTR 30 , Table 1) was the most influential metric for predicting bleaching prevalence score, with greater daily temperature variability serving as a mitigating factor (Fig. 3b) Fig. 3a), highfrequency temperature variability was both the greatest mitigating factor of bleaching prevalence score and the most influential covariate-more influential than widely used metrics of acute and cumulative thermal stress by a factor of 2 and 3 times, respectively (Fig. 3c). Using globally averaged values of explanatory variables, our highest-ranked logit model ( Fig. 3b) implies that, in native units, a 1°C increase from the mean DTR 30 value would decrease the odds of more prevalent bleaching by a factor of 33. To standardize this, each unit increase in highfrequency temperature variability (i.e., DTR 30 ) would reduce the odds of more prevalent bleaching by a factor of e 2.66 = 14.3. Contrasting this against a unit increase in cumulative thermal stress (i.e., DHW 30 ), which would only increase the odds of more prevalent bleaching by a factor of 2.6, highlights the dominant influence of diurnal temperature variability on reef-scale bleaching prevalence. "Depth" (Table 1) was taken as the mean depth of the water temperature measurement, in meters below the surface, for each site, and is also representative of local water column depth as sensors were placed near the bed. Depth was the second-most effective predictor of bleaching prevalence (Fig. 3c), with deeper reefs less likely to experience pervasive bleaching. However, "depth" is also a proxy for other characteristics of the reef sites such as habitat (e.g., deeper forereefs and lagoons, shallow reef flats) and light intensity, which decays exponentially with depth. Although the logit models preclude significant collinearity of tested variables (Methods), corals at shallow depths may experience greater high-frequency temperature variability 45 , although accounting for water flow can complicate this interpretation as it pertains to bleaching 26,47 . High-frequency temperature variability and depth may mitigate bleaching in complementary ways: habitats with greater high-frequency temperature variability, which are likely to be found at shallower depths 45 , may develop greater thermal tolerance 19,48 , while deeper coral habitats, despite their propensity for milder diurnal temperature variability (outside of internal wave-influenced regions [49][50][51] ), may serve as refuge areas resistant to the intrusion of hot water 25 , perhaps facilitating recovery of coral cover following bleaching events 52 .
"Background Conditions", "Cumulative Thermal Stress", and "Acute Thermal Stress" were the three explanatory variable categories largely suspected of exacerbating bleaching. "Background Conditions" (Table 1) consisted of the average summertime, or maximum monthly mean (MMM), temperature, but computed from our in situ time series data, as opposed to conventional remotely sensed SST data 14 . "Background Conditions" also included the latitude of the temperature logger, a variable that served as a proxy for unresolved oceanographic factors 53 related to the large-scale processes that influence climatologies. The "Cumulative Thermal Stress" category (Table 1) encompassed various methods for the computation of the magnitude and duration of acute in situ thermal stress exposure on reefs. Similar to the MMM, cumulative thermal stress is traditionally derived from remotely sensed SSTs and is among the most common metrics used to predict coral bleaching 25,33 . The "Acute Thermal Stress" category (Table 1) was included as a safeguard to differentiate sites with temperatures that may not have exceeded MMM+ 1°C (i.e. no thermal stress) yet still experienced bleaching. Consistent with the wellestablished perspective that anomalously high temperatures are the primary cause of coral bleaching 7 , among our highest ranked models, bleaching was most exacerbated by greater cumulative and acute thermal stress, and also, to a lesser degree, by increases in MMM temperature and heating rate 25 . "Heating Rate" (Table 1)  been shown to be crucial for determining the fate of corals to summertime bleaching susceptibility 25 . The "Thermal Trajectory" (Table 1) category followed the methodology of a previous study that highlighted the role of protective warm, pre-stress temperatures as being important for resilience to bleaching from intense acute stress temperature events 39 . Our results reinforce recent findings that a reef's thermal trajectory is a significant predictor of bleaching prevalence 39 (Fig. 3c), with thermal tolerance conferred by exposure to a protective, sub-lethal bleaching stress prior to acute stress exposure. Although not as influential as the above variables, a no-stress or protective thermal trajectory (i.e., a prestress, 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 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 temperature time series have been associated with site-specific increased thermal tolerance 54 .
To summarize the results of our highest-ranked logit model, we 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 highfrequency temperature variability (DTR 30 ) 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 would increase this probability by 41% (Fig. 4b). Similarly, a 0.03°C/day increase in ROTC SS from its mean value would increase the probability of Category 4 bleaching by 34% (Fig. 4c), and a 1°C-weeks increase in DHW 30 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 SSTderived covariates, to determine how the results would differ from the in situ driven models. We obtained weekly 4 km resolution CoRTAD SST data 55 , using the data pixels closest to the coordinates of our in situ loggers, and used this SST dataset to compute covariates within the Background Conditions, Acute and Cumulative Thermal Stress, and Heating Rate categories (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, specifically DTR 30 , 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 Cumulative Thermal Stress provided the strongest influence (Fig. 5a, b), while Background Conditions (MMM 4 km , Table 1) represented a 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  The best model (i.e. ΔAIC C = 0) included six variables, of which high-frequency temperature variability was the absolute most influential and also greatest mitigating factor to bleaching prevalence. c Summing across 20 indistinguishably good models (i.e. within ΔAIC C ≤ 2), high-frequency temperature variability was consistently most influential. Variable categories are shown in Table 1. Delete-1 jackknife standard error bars are shown in (b), while the standard error bars shown in (c) were obtained by summing in quadrature the individual standard errors from each of the 20 models computed after delete-1 jackknife resampling SST-based OLR models indicate that upon consideration of the consistent importance of DTR to bleaching, a globally available remotely sensed metric for diurnal temperature variability would be valuable for improved bleaching predictions.
Specific reef cases. Our results reveal the importance of highfrequency temperature variability at locations worldwide, but include reef-scale observations where such variability influences bleaching prevalence of corals in different locations of the same reef during the same bleaching event 32,46 . Here, we present two such case studies: one from Tahala Reef, a platform reef in the central Red Sea, and another from Nelly Bay, a fringing reef in the Great Barrier Reef in Australia (Fig. 6). These sites were chosen due to the availability of additional meteorological data 32,56 at these reefs. At each location, temperature time series (Fig. 6a, b) for both a seaward and a shoreward location show that, whereas low-frequency variations in water temperature are often very similar over reef-scales (Fig. 6c, d), high-frequency variations may be quite distinct (Fig. 6e, f). In these cases, bleaching events were more widespread and severe at the seaward locations where DTRs were smaller (Fig. 6a, b), consistent with our best logit models.

Discussion
For corals, a shift in thermal tolerance can occur due to adaptation of the coral animal or algal symbionts through natural selection of heat-tolerant lineages 57,58 , or physiological acclimation through the expression of heat shock proteins and regulation of apoptosis (i.e., programmed cell death) 23,36 . As discussed, recent work highlights the importance of short-term temperature history (daily-weekly periods) for coral acclimatization to higher temperatures 19 , such that corals subject to warmer than average temperatures prior to thermal stress may exhibit a greater tolerance to acute temperature stress 23 . In the context of these studies and in keeping with other site-specific and experimental studies 19,26,45,47,59 , our results suggest that temperature fluctuations on daily or tidal timescales are often sufficient to expose corals to temperatures high enough to encourage greater tolerance (via acclimation or adaptation) to thermal stress, but for time periods short enough to avoid mortality 19,48 . Further, our results establish that the resistance of corals located in areas of high-frequency temperature variability to bleaching occurs in reef regions throughout the world. While we lack sufficient species-level data, we fully acknowledge that intrinsic coral properties 40 , differences in reef-scale community compositions 56 , and taxonomic susceptibility 26 are likely to influence heterogeneous reef-scale bleaching responses, and an improvement to our model framework would include species-level covariates. Our results also demonstrate the potential to both improve predictions of bleaching and anticipate heterogeneous patterns in bleaching prevalence at reef-scales, considering that beyond mean values, accounting for variability in temperature regimes yields better predictions of organismal responses to anomalous environmental events 60 . Although SSTs from satellite remote sensing are not yet available at the spatiotemporal resolution required to calculate reef-scale high-frequency temperature variability, observational work on a range of reef structures suggests that it may be possible to predict reef-scale thermal environments using relatively simple hydrodynamic models given readily available bathymetry and basic hydrographic data such as tidal range, wave height, and offshore mean SST 30,32,43 . While we did not assess other biogeochemical parameters, it is worth noting that the same physical circulation that drives spatially variable thermal environments, in locations with an active benthic community, can also create dynamic oxygen, pH, and nutrient environments 61,62 .
Urgent global efforts at reducing anthropogenic greenhouse gas emissions must remain a priority for reef preservation, due to the acute thermal stress that is now arising from global warming projections on reefs 5,63 . However, combating the effects of local stressors on reefs through conservation tools such as marine protected areas is likely to increase the chances of reef persistence through future warming, as well as to facilitate recovery following a bleaching event 2,52,64,65 . Warming ocean temperatures are projected to result in annual severe bleaching regimes by the middle of this century, with spatial variability in the onset of these events on the order of ±10 years 4 . Considering our results in the context of this inevitable and persistent acute thermal stress, focusing management efforts on the more resistant reef locations that also experience delayed onsets in annual severe bleaching would maximize the likelihood that at least some healthy reefs will exist in the future.

Methods
Data synthesis. Water temperature time series from 118 locations, representing five major ocean basins where tropical coral reefs are found (Western Indian Ocean, Pacific Ocean, Caribbean Sea, Great Barrier Reef, and Red Sea), were 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 (AIMS), the National Data Buoy Center (NDBC), and the Florida Institute of Oceanography (FIO). However, time series within our dataset were also selected to match precise reef locations that had sufficiently documented bleaching events, while also containing as long and consistently sampled records as possible. Site names and three-letter codes, locations, depths, instrument descriptions, and additional information for each time series are listed in Supplementary Data 1, which also lists the source for each time series. Water temperature records included in this analysis spanned at least 12 months in duration, with a sampling interval less than or equal to 3 h. In cases where instrument substitution resulted in varying sampling intervals, time series were sub-sampled to the largest of these intervals, or, in rare cases, interpolated to remain below a 3 h sampling interval. The temperature time series data used in this study originate from instruments that were calibrated using varying methodologies, such as by placing loggers together at one location and comparing recorded temperatures with a reference temperature dataset, ice bath calibration, or multiplying raw field recorded temperatures by normalized logger calibration coefficients. Our analysis is largely based on relative temperature variations, however, so that absolute temperature accuracy will not 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 from nearly 1000 reef locations 66 . From these tests, we cannot conclude any significant differences (α = 0.05) between our time series data and this larger global data set. All data analyses were done using MATLAB 7.14 (The Mathworks, Natick, MA, USA).
Spectral analysis. Power spectral density (PSD) estimates were computed for each temperature time series. First, if necessary, temperature time series were resampled or linearly interpolated to maintain a constant sampling interval, chosen to be 3 h so as to resolve spectral frequencies of up to 4 cpd. In order to examine temperature variability for a broad range of frequencies spanning annual to diurnal and shorter periods, a PSD was calculated as follows: time series greater than or 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 equal to 10 years in duration were divided into 3-year sections, while all others were divided into 4-month sections. Sections were overlapped by 50% and windowed with a Hamming function before spectra for each section were calculated, which were then ensemble averaged to obtain the PSD estimate. The statistical significance (α = 0.05) of observed spectral peaks was ascertained by comparison with the upper confidence level of a background red noise fit to the spectrum 67 . We empirically defined annual, seasonal, and diurnal frequency bands as 0.00185 to 0.0111, 0.0119 to 0.143, and 0.727 to 1.333 cycles per day (cpd), respectively. Band variance was computed using trapezoidal integration of PSD values within each respective frequency range, and ratios of high-frequency to seasonal band variances among habitat types ("Back Reef", "Reef Flat", "Reef Slope") were compared using a Kruskal−Wallis test (Supplementary Fig. 2b).
Spatiotemporal variability in water temperature. To quantify the magnitude of diurnal patterns of heating and cooling, a DTR was calculated as the difference between maximum and minimum temperatures for each day of each time series. Temporal variations in DTR values were examined in multiple ways. First, DTRs were composite-averaged for yeardays 1−366, and based on evident seasonal DTR variability from panel b of Fig. 2, a non-parametric Kruskal−Wallis test was used to assess the seasonal dependence of DTR distributions 68 . Seasons were defined such that each season spanned 3 complete months, and austral and boreal summers were December through February and June through August, respectively. Temperature metrics are summarized in Table 1. Large-scale spatial patterns in water temperature variability relative to latitude were also characterized; annual temperature ranges, calculated as the range of monthly mean temperatures, were compared against latitude for all sites. Variance in high-frequency (33 to 6 h periods) and seasonal (7 to 84-day periods) spectral bands were computed via integration of their respective power spectral densities, and the ratios of highfrequency to seasonal variance for all time series were compared by the habitat from which each time series was recorded. Habitats were divided into three groups: (i) back reefs and back reef lagoons (labeled as "BR"), (ii) reef flats and non-back reef lagoons ("RF"), and (iii) reef crests, reef slopes, forereefs, and anything further offshore ("RS").
Bleaching data synthesis. Coral bleaching observations (81 events) that corresponded in time and location with temperature data (Supplementary Data 2) were obtained from a variety of sources, but primarily from peer-reviewed publications. Bleaching reports were based upon various quantification schemes, but some common methods included recording bleaching as: (i) a proportion of transect area 69 , (ii) severity among different colonies of different coral species 70 , and (iii) prevalence categories based on aerial (and ground-verified) surveys 22 . In publications that reported percentages of colony bleaching within different coloration and paling categories, we used the weighted average of the different percentages within each category, though the majority of bleaching records naturally translated into our four bleaching prevalence scores. To aggregate 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 prevalence score: 1: "bleaching prevalence" ≤ 10% (n = 48); 2: 10% < "bleaching prevalence" ≤ 25% (n = 5); 3: 25% < "bleaching prevalence" ≤ 50% (n = 6); 4: "bleaching prevalence" > 50% (n = 22). Note that bleaching events of 0 (reports of no bleaching observed or negligibly mild paling) are binned into bleaching prevalence score 1, to provide a conservative grouping for mild bleaching events, seasonal patterns of discoloration and variation in zooxanthellae densities 71 , unresolved "background bleaching" levels 72 , and observation errors. As opposed to continuous interval variables, ordinal variables represent categorical values that can be ranked and have a natural ordering to them, offering the advantage of creating bins for ranges of values.
Explanatory variables. To assess the influence of explanatory environmental variables, including high-frequency temperature variability, on bleaching response, we performed a multivariate statistical analysis of the observed bleaching events. As there are multiple aspects of thermal stress and environmental conditions that may explain the bleaching response, we selected 20 experimental variables for the ordinal regression (defined in Table 1  °C J a n 9 6 J a n 9 7 J a n 9 8 J a n 9 9 J a n 0 0 J a n 0 1 J a n 0 2 J a n 0 3 J a n 0 4 J a n 0 5 Temperature Variability. We lacked sufficient habitat (i.e., reef flat, reef crest, reef slope, etc.) information at which many loggers were placed, and therefore did not include habitat as an explanatory variable. Environmental variables in the Cumulative Thermal Stress category were calculated from in situ temperature measurements rather than the remotely sensed National Oceanic and Atmospheric Administration's Coral Reef Watch (NOAA CRW) products because many of our logger locations were outside of the coverage areas of the CRW Virtual Stations, and furthermore, as the first suite of CRW products was released in 2000, they do not include temperature data corresponding to the 1998 global bleaching event.
The explanatory variables we use here are commonly appended with a subscript to indicate the period of time, relative to the bleaching observation date, used to calculate the metric. For example, the subscript "SS" denotes the spring to summer period, and the subscript "30" denotes the 30 days preceding (and including) a bleaching observation.
Principal components analysis. Principal components analysis was performed to examine the spatial structure of the environmental forcings and determine the association between independent variables within each of the eight categories (Fig. 1). The first two PC axes accounted for 44.2 and 18.8% of the variance within the matrix of independent variables. The magnitude and orientation of the loading vectors indicate the importance of each parameter in describing the variance of the PCA components.
Computation of thermal trajectory and acute stress variables. An ordinalvalued Thermal Trajectory 39 variable was included as an independent variable to assess the degree to which environmental conditions confer thermal tolerance. The calculation of this was as follows: first the MMM and MMM+ 2°C (the latter quantity is referred to as the "local bleaching threshold") were computed for a given time series. Then the 33-h low pass filtered time series for the 90 days preceding a bleaching event was inspected to determine the type of Thermal Trajectory. If temperatures exceeded the MMM, then fell below the MMM for a required 10-day "recovery period" before proceeding to exceed the local bleaching threshold, a Protective Trajectory with an ordinal value of 1 was recorded. If temperatures increased from below the MMM to above the local bleaching threshold, without a 10-day recovery period, a Single Bleaching Trajectory with a value of 2 was recorded. If temperatures exceeded the local bleaching threshold at least twice, with a required 9-day recovery period between threshold exceedances, a Repetitive Bleaching Trajectory with an ordinal value of 3 was assigned. Finally, if temperatures did not exceed the local bleaching threshold, an ordinal value of 0 corresponding to no thermal stress was assigned. Justification of our ordinal-value scheme comes from an analysis of experimentally heated corals from the Great Barrier Reef, whereby in the face of thermal stress, corals with a Protective Trajectory experienced localized cell death of approximately 30%, while that of corals under Single and Repetitive Bleaching Trajectories was approximately 60% and 70%, respectively 39 . Furthermore, experimental results showed that corals under a Protective Trajectory maintained significantly greater symbiont density than those under Single and Repetitive Bleaching Trajectories, and hence the ordinal scores from 0 to 3 aptly account for the monotonic nature of coral tissue detriment associated with no heat stress, a Protective Trajectory, a Single Bleaching Trajectory, or a Repetitive Bleaching Trajectory, respectively. The "No Thermal Stress" trajectory corresponded to~0% cell death and greater symbiont density than the Protective Trajectory. The Acute Thermal Stress category was composed of two binary variables to indicate the presence/absence of acute thermal stress. The calculation of this was as follows: first the MMM+ 1°C and MMM+ 2°C were computed for a given time series, then the daily mean temperatures within 90 days before a bleaching event were inspected to determine if temperatures exceeded MMM+ 1°C (in which case Acute1 would equal 1) and MMM+ 2°C (in which case Acute2 would equal 1).
Ordinal logistic regression. Using the eight explanatory variable categories described above, and bleaching prevalence scores from 1 to 4 as the response variable, we performed an ordinal logistic regression analysis to determine how the relative log odds of a given bleaching prevalence score depends on the interactions among the explanatory variables. Ordinal regression models have been adeptly used in ecological studies where data are often present as semi-quantitative variables in which relative differences between values are of importance 73,74 . Furthermore, logit functions have been previously implemented to predict the presence/absence of bleaching using gridded remote-sensed data 72,75 , or to explain the influence of a range of environmental and coral physiological factors on reef ecosystem response following a disturbance 52 . These logit models are multivariate extensions of generalized linear regression models 76 , providing parameter estimates via maximum likelihood estimation (MLE) to model the relative log odds of, for our purposes, observing one bleaching prevalence score or less versus observing the remaining greater bleaching prevalence scores: Here i indexes each of N observations, with bleaching observation y i , and the lefthand side quantity is referred to as the logit of the probability of observing bleaching prevalence score j or lower, for j = 1, 2, or 3 (observations with bleaching prevalence scores of 4 contribute to the regression through calculation of the logodds). Note that the odds are defined as the ratio of the probability of an event occurring to the probability of the event not occurring, which is exactly the ratio inside the natural logarithm. Each C j is an MLE-computed model intercept, and each B k is the MLE coefficient corresponding to each standardized independent variable z ik , for k = 1, …, p, where p is the number of independent variables used in a given model. A fundamental component of this model is the assumption of proportional odds, or parallel regression, which implies B k values are independent of the logit level j. The validity of this parallel regression assumption was ascertained using Brant's Wald test 77 , as well as a likelihood ratio test (α = 0.05).
As each time series was the sole source of its explanatory variables, we can expect many of these variables to be correlated with each other (multicollinear). If left unaccounted for, multicollinearity obscures the interpretation of the explanatory variables and their coefficients, and may decrease the statistical power of the logit analysis 78 . The degree of multicollinearity among the 20 explanatory variables used in the logistic regression was assessed by calculating the condition indices and variance-decomposition proportions 79 Table 1). Note that within each of these four multicollinear groupings, the multicollinear variables come from the same explanatory variable category (for example, both DHW 90 and CSA During are both from the Cumulative Thermal Stress category). Therefore, for each logit model, we selected at most one variable from a category without multicollinearity erroneously influencing our results. Spatial autocorrelation within each covariate was determined by calculating Moran's I and examining correlograms, from which we determined significant spatial autocorrelation typically up to distances of 500 km ( Supplementary Fig. 5). To account for this, we also added a random effect to the highest-ranked logit models, from which we failed to conclude a marginal model improvement (see below).
All permutations of all possible explanatory variables were used to compute a total of 10,367 logit models, where all logit models were computed using a multinomial logistic regression function ("mnrfit") in MATLAB. Model comparison was performed using a bias correction for small sample sizes to Akaike's Information Criterion, AIC C 80 , and all models within ΔAIC C ≤ 2 of the best model (ΔAIC C = 0), which have statistically indistinguishable performances 81 , are presented in Supplementary Fig. 6. McFadden's pseudo-R 2 was also computed for the highest ranked models, and ranged from 0.26 to 0.30, with that of the highest-ranked model equal to 0.30. While the logit model with the lowest AIC C , as well as all models within ΔAIC C ≤ 2, provide a general outline for the coefficients of the critical independent variables in explaining bleaching prevalence, these parameter estimates have errors of unknown distribution. Additionally, the possible existence of outliers in the high-frequency temperature variability data may influence the results of the logit parameter estimates, and therefore, delete-1 jackknife resampling was used to compute estimates of bias and standard errors. We found a slight positive bias that did not significantly alter the influence of highfrequency temperature variability relative to the other covariates. Estimates of standard errors for all logit parameters of all models within ΔAIC C ≤ 2 can be seen in standard error bars ( Fig. 3b and Supplementary Fig. 6). Specifically, the error bars seen in Fig. 3c are the standard errors from each contribution summed in quadrature. A modified jackknife resampling scheme was also performed, in which instead of leaving out one site at a time, sites within 10 km of each other were grouped together and each of these proximity groups was left out incrementally before fitting OLR models to the remaining data. This spatial resampling analysis ( Supplementary Fig. 7) did not result in significantly different parameter estimates than the full model presented in Fig. 3b. Estimates for the intercept terms C 1 , C 2 , and C 3 (±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 prevalence categories within our dataset. Using the statistical computing software R 82 (https://www.r-project.org), a random effect grouping variable was added to each of the highest-ranked models, which grouped reefs within 5 km of each other. These resulting mixed effects models were compared to their fixed effects equivalents to determine model fit improvement, and the inclusion of a random effect did not improve the fit of any of the 20 highest-ranked models (Supplementary Data 3). Furthermore, to account for possible nonlinear interactions between covariates, such as Depth × DTR 30 or Acute1 × DHW 30 , we also included a nonlinear interaction term to the highest-ranked model and determined model fit improvement. The AIC C of these nonlinear models (Supplementary Table 1) indicated that the addition of a nonlinear interaction term did not significantly improve the fit of the original highest-ranked model displayed in Fig. 3b. Therefore, as our main result, we ultimately report the fixed effects ORL model parameter estimates with no nonlinear or interaction terms (Fig. 3b).
The SST-based OLR analysis summarized in Fig. 5 was performed using the CoRTAD 4 km weekly SST data 55 pixels that were closest to our in situ bleaching observations. The quantities MMM, MMM Max 15 , Acute1, Acute2, DHW, ROTC SS , and ROTC 90 (Table 1) were computed from the CoRTAD temperature data. This 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 permutations of covariates, with the constraint that within each model, we only include ≤1 covariate from each category listed in Table 1. This produced 60,479 models, from which we performed model comparison using AIC C , resulting in 12 dominant models being identified as statistically indistinguishable. From here, we concluded that when 4 km SST metrics are incorporated into our model framework, high-frequency temperature variability and depth remain the most influential covariates on bleaching prevalence, acting to attenuate bleaching (Fig. 5). Furthermore, acute and cumulative thermal stress exacerbate bleaching prevalence the most, and the SST-derived version of the latter contributes substantially to this effect. The results presented in Fig. 5 ultimately mirror those from Fig. 3b, c, indicating that DTR is the main driver of bleaching prevalence, regardless of whether SST or in situ quantities of the other predictors are used. Common covariates in these models were MMM, MMM Max , MMM 4 km , DHW 30 , DHW 4 km , Acute1, Acute2 4 km , and ROTC 90-4 km , and each of these 12 highestranked models included at least one SST-based covariate. The AIC C values of these models ranged from 129.19 to 131.12, with McFadden's pseudo-R 2 values ranging from 0.35 to 0.38, which presents a considerable improvement to the highestranked model that was fit to only in situ values, which had an AIC C equal to 143.75 and a McFadden's pseudo-R 2 of 0.30. To highlight the importance of High-Frequency Temperature Variability in model improvement, we also fit these models to the data and excluded the Daily Temperature Range covariate, which was either DTR 30 or DTR 90 . This significantly decreased model performance, with AIC C values ranging from 153.46 to 169.57, and McFadden's pseudo-R 2 values ranging from 0.12 to 0.23, which supports the idea that DTR 30 is the dominant driving variable for the bleaching response.