High coral heat tolerance at local-scale thermal refugia

Marine heatwaves and mass bleaching have devastated coral populations globally, yet bleaching severity often varies among reefs. To what extent a reef’s past exposure to heat stress influences coral bleaching and mortality remains uncertain. Here we identify persistent local-scale hotspots and thermal refugia among the reefs of Palau, Micronesia, based on 36 years of satellite-derived cumulative heat stress (degree heating weeks–DHW, units: ˚C-weeks). One possibility is that hotspots may harbour more heat tolerant corals due to acclimatisation, directional selection, and/or loss of tolerant genotypes. Historic patterns of assemblage-wide mass bleaching and marine heatwaves align with this hypothesis, with DHW-bleaching responses of hotspots occurring at 1.7˚C-weeks greater heat stress than thermal refugia. This trend was consistent yet weaker for Acropora and corymbose Acropora , with severe bleaching risk reduced by 4–10% at hotspots. However, we find a contrasting pattern for Acropora digitifera exposed to a simulated marine heatwave. Fragments of 174 colonies were collected from replicate hotspot and thermal refugium outer reefs with comparable wave exposure and depth. Higher heat tolerance at thermal refugia (+0.7˚C-weeks) and a correlation with tissue biomass suggests that factors other than DHW may overwhelm any spatially varying effects of past DHW exposure. Further, we found considerable A . digitifera heat tolerance variability across sites; compared to the least-tolerant 10% of colonies, the most-tolerant 10% could withstand additional heat stresses of 5.2 and 4.1˚C-weeks for thermal refugia and hotspots, respectively. Our study demonstrates that hotspot reefs do not necessarily harbour more heat tolerant corals than nearby thermal refugia, and that mass bleaching patterns do not necessarily predict species responses. This nuance has important implications for designing climate-smart initiatives; for instance, in the search for heat tolerant corals, our results suggest that investing effort into identifying the most tolerant colonies within individual reefs may be warranted.


Introduction
Climatic stressors are causing profound impacts on organisms and ecological systems globally [1,2].This is especially true for sessile organisms which cannot move to escape their environment.To persist under climate change, many species will need to adapt, and this will require sufficient standing genetic variability of climate resilience traits within populations [3].Over large geographic scales (> several degrees of latitude, and between populations), gradients in climatic exposure are known to shape organism tolerance to climatic stressors [4][5][6].Higher exposure to stress has been shown to drive higher tolerance through local adaptation (e.g., drought tolerance in European trees [5], thermal tolerance in corals [4], and upper/lower thermal limits for both terrestrial and marine ectotherms [6]).However, this pattern has also been shown to break down over the smaller spatial scales at which organisms and populations operate (< several degrees of latitude, or within populations) due to the influence of different abiotic selective forces (e.g., upper thermal limits and cold tolerance not explained by mean temperature [7,8]) and abiotic factors (e.g., genetic correlations between different thermal tolerance traits [9]).Furthermore, there is an emerging appreciation that microgeographic adaptation-adaptive divergence that occurs within a single population even amid high gene flowcan play an important role in ecological and evolutionary dynamics in nature [10].To effectively manage populations under climate change, there is a crucial need for more information collected at local spatial scales (i.e., within populations, which could vary in terms of absolute distance depending on species range and dispersal capacity [10]).Despite this, the extent to which organism tolerance is influenced by local scale gradients in climatic exposure remains poorly understood.
Reef-building corals epitomise this problem because of their acute sensitivity to marine heatwave stress, which disrupts their symbiosis with phototrophic microalgae and can lead to mass bleaching and mortality [11][12][13].The ability of corals to tolerate heat stress is influenced by numerous other thermal drivers (e.g., temperature variability [14], heating rate [15], and peak temperature [16]), non-thermal environmental variables (e.g., light stress [17], turbidity as a mediator of irradiance exposure [18], water flow [19], and water quality [20]), and biological factors (e.g., energy reserves [21] and Symbiodiniaceae taxa [22]).However, it is the accumulation of heat stress over the warm season-a variable that is mapped at daily resolution in near real-time by satellites [23]-that is the ultimate trigger for mass coral bleaching.Therefore, leveraging spatial variations in heat stress exposure to identify reefs with high heat tolerant corals would be a useful tool for coral reef management.
Across seascapes, some coral reefs are routinely exposed to higher-than-average levels of heat stress while other reefs are consistently exposed to less-intense levels of heat stress, herein referred to as hotspots and thermal refugia, respectively.These distinct thermal stress regimes have been identified across broad geographic regions (>1000 km, e.g., the Great Barrier Reef [24], Red Sea [25], and the Caribbean [26]).As for many other marine and terrestrial species, the thermal limit of corals largely tracks broad latitudinal climatological gradients [4,12,27].When exposed to marine heatwaves, corals at hotspots may be expected to fare better than those at thermal refugia, through means of acclimatisation [28,29], adaptation [30,31], and/or loss of heat-sensitive coral genotypes [32].However, the opposite trend has been noted for the Red Sea, where corals in the cooler north have a higher heat tolerance than those in the south, possibly because they are theorised to have experienced a historic selective filter (upon recolonising the Red Sea after the last ice age), and now have a thermal threshold above the ambient water temperature [25].From a perspective of managing species populations, it is desirable if heterogeneity in recent (i.e., in the satellite record) exposure of reefs to marine heatwave stress leads to predictable differences in coral population responses across seascapes, as this could faac387cb999778055cc to be completed).All datasets analysed are publicly available as of the date of publication.Any additional information required to reanalyse the data reported in this paper is available from the lead contact upon request.allow managers to prioritise population interventions (e.g., protection, outplanting, translocations) in locations where they are likely to be most effective.This principle also holds true from the wider coral reef management perspective, where target coral assemblages comprise numerous species with different bleaching susceptibilities [33], and the ability to predict locations of tolerant assemblages could support effective management decision-making (e.g., which reefs to protect).However, there is growing recognition that the heat stress responses of individual conspecific corals are highly variable even within single populations [34], which could influence the feasibility of spatial seascape management plans if, for instance, speciesspecific heat tolerance variability were to be greater within than between reefs.Because coral reef management typically takes place over local within-population scales (i.e., 10s to 100s km for corals given their dispersal capacity [35]), there is a growing need to study smaller environmental gradients.Specifically in this paper we test whether the large-scale among-population expectation-that hotspots should harbour higher levels of heat tolerance than thermal refugiaholds true over smaller, more manageable distances within coral populations, or whether this effect is overwhelmed by drivers other than past heat stress exposure.
Here, we combine geographic-scale concepts of local adaptation of different populations, with a more detailed exploration of within-population variability of coral heat tolerance, given new research on the role of microgeographic adaptation.These smaller within-population local spatial scales are highly important for management actions (e.g., patrolling marine protected areas, or collecting corals for restoration) as logistics become far more challenging when managing reefs over larger distances (i.e., multiple populations).Our study takes a local-scale seascape perspective, drawing on long-term historic spatial datasets available from this region on mass coral bleaching and marine heatwave intensity, and a species-specific experimental manipulation.Because the mass bleaching dataset is given at the assemblage level [36], possible adaptive or acclimatory responses to past heat stress could be masked by numerous environmental variables (e.g., temperature variability, irradiance) and biotic factors (e.g., coral community composition, Symbiodiniaceae community).For this reason, we also performed a reanalysis of Acropora-specific bleaching response for the only marine heatwave from which data were available (2010), and conducted a 5-week simulated marine heatwave experiment under controlled conditions to test the heat tolerance of a common coral species collected from both hotspots and thermal refugia.If we find that the broad geographic-scale patterns (i.e., higher heat tolerance at warmer reefs [4]) are also seen at the within-population scale, then it simplifies the provision of advice to ecosystem managers, that hotspot reefs might harbour corals with predictably higher levels of heat tolerance than thermal refugia.

Methods
Here, we evaluated local-scale thermal stress regimes among the reefs of a small remote archipelago using a 36-year satellite record of sea surface temperatures (SST).The Republic of Palau serves as an ideal case study as the spatial scale of its reef system is an order of magnitude smaller (<150 km long, within one coral population-local) than reef systems previously linked to gradients of thermal exposure and coral heat tolerance (1000s km [24][25][26], between multiple coral populations-not local).To test whether persistent hotspots or thermal refugia are present in Palau, spatial patterns in thermal history and marine heatwaves were assessed based on accumulated warm season heat stress.Firstly, the spatial variability of bleaching risk across the entire coral assemblage was tested using a Bayesian statistical modelling approach linking heat stress data with historic bleaching survey observations from across Palau.Then we restricted our spatial comparison to a single species to avoid any possible confounding effects of species compositional changes on coral responses to heat stress.We achieved this by conducting a long-term (sensu [37]) 5-week marine heatwave experiment with assays of bleaching and mortality every 1-3 days.

Historic heat stress data
Heat stress on Palauan coral reefs was calculated from CoralTemp version 3.1, a daily 0.05˚x 0.05˚latitude-longitude (~5 km) resolution satellite-based Sea Surface Temperature (SST) dataset (1985 to 2020) available from the National Oceanic and Atmospheric Administration's Coral Reef Watch (NOAA CRW) [38].Although coral bleaching mortality can be influenced by light intensity, cold spells, nutrient enrichment, and sedimentation, among other factors, chronic levels of accumulated heat stress is the primary driver of mass coral bleaching and mortality [11].NOAA CRW uses the CoralTemp satellite dataset to measure accumulated heat stress in terms of Degree Heating Weeks (DHW, units: ˚C-weeks), a metric combining both the intensity and duration of marine heatwaves.Coral reefs exposed to higher levels of DHW have a greater risk of mass bleaching and mortality [39].As such, DHW is used by NOAA CRW to provide real-time bleaching risk forecasts, whereby DHW of 4-8˚C-weeks corresponds to the expectation of significant bleaching, and DHW > 8˚C-weeks corresponds to the expectation of significant bleaching and mortality.Here we use the same standard NOAA DHW dataset as presented in Lachs et al. (2023a).Briefly, we follow the NOAA CRW methodology [23], DHW on a given day (i) is computed as the sum of the last 12 weeks (84 days) of daily temperature anomalies (HotSpots, a standard parameter used by NOAA, not the thermal regime) relative to a standard coral stress threshold baseline (MMM, maximum of monthly means) which is held constant through time, and only HotSpots > 1˚C are accumulated and then divided by 7 to make the standard DHW a weekly (rather than daily) metric.

Identifying persistent hotspots and thermal refugia
Whether or not hotspots and thermal refugia were present in Palau was determined using annual maps (n = 36, 1985-2020) of maximum DHW (the maximum DHW in the year reflects peak bleaching heat stress, herein the maximum annual DHW is referred to as DHW for simplicity), combined with a thermal consistency biplot analysis (Fig 1 ), which computes the typical heat stress received by a given reef and compares this to the consistency of the reef's heat stress performance relative to other reefs.For each year independently, maps of DHW (S1A Fig) were transformed to anomalies (based on the Palau mean DHW for that year) and percentile ranks (S1B Fig) to capture relative spatial differences in heat stress performance for all grid cells that intersect coral reefs in Palau [40].Then for each reef, we calculated the mean DHW anomaly through time to describe the reef's heat stress exposure relative to other reefs, and the standard deviation of DHW percentile ranks through time to measure the inter-annual consistency with which the reef performs relative to other reefs (different than sub-daily thermal variability).As some reefs are always more or less heat stressed in comparison to others, the mostand least-heat stressed reefs in Palau could be determined using the mean DHW anomaly based on the upper and lower terciles.However, to be considered hotspots or thermal refugia, reefs must have such thermal conditions consistently through time.Therefore, the reefs were also categorised based on DHW percentile SD (standard deviation), with the most-consistent thermal regimes occurring in the lower tercile (DHW percentile SD cutoff of 23.4), and the  We also report on the median daily SST of Palauan reefs and the variability of daily SSTs computed as the 99% range (0.5 th percentile SST to 99.5 th percentile SST).

Assemblage-wide coral bleaching data
The historic coral bleaching survey data used in this study is the same subset of the publicly available global dataset [36,42] presented in Lachs et al. (2023a) (1998-2017, n = 239).Notably, bleaching observations (primarily reported as percentage of corals bleached) are given as severity scores, ranging from 0 to 3, namely: no bleaching (0%), mild bleaching (1-10%), moderate bleaching (11-50%), and severe bleaching (>50%).We filtered out any survey records from the southern Palauan islands which are over 500 km from the main island of Palau (n = 2, Helen Reef and Tobi Island in 1998) and any records with suspected incorrect coordinates, where the location did not intersect with reef area (n = 3 in 1998).

Spatial variability in assemblage-wide coral bleaching
To quantify the spatial variability of DHW-bleaching responses across the coral assemblage (i.e., all hard coral taxa) in Palau, we combined CoralTemp with historic coral bleaching survey observations.The effect of DHW on bleaching was fitted using a spatial beta GLM via Integrated Nested Laplace Approximation (INLA) in R-INLA.While this model was also used as part of a separate project focussing on temporal variability [43], the present study provides extended analysis of the spatial component of the model and additional sensitivity analyses.
The Bayesian statistical approach we apply provides an intuitive spatially explicit estimation of model uncertainty [44,45].Essentially, this is a map of spatially correlated uncertainty that shows where bleaching is under-or over-predicted for a given heat stress dosage (DHW).To fit the beta distribution, bleaching severity scores were converted to proportions (divide by 3) and transformed to remove extremes of 0 and 1, which in practice is ðy � ðn À 1Þ þ 0:5Þ=n, where y is the bleaching severity score scaled from 0 to 1, and n is the sample size [45,46].This approach preserves differences among bleaching categories that would be lost using other transformations (e.g., binary transformation for a binomial GLM [47]).The full model specifications are given in [43], but briefly, spatial variation in the uncertainty of bleaching-DHW responses was estimated across a high-resolution Delaunay triangulation mesh of the study area (5, the occurrence or not of severe bleaching-bleaching in >50% of colonies) as a function of heat stress (DHW i ) whilst accounting for additional underlying spatial correlation among bleaching observations (random effect: u i ) and spatially uncorrelated error (ε i ) where β 0 is the intercept, β 1 is the slope for DHW, u i represents the smoothed spatial random effect.Notably, the probability of severe bleaching pertains to the occurrence, or not, of bleaching in > 50% of colonies [36,42] and so cannot resolve spatial patterns in the magnitude of severe bleaching (i.e., where 50-100% of colonies are bleached).Whether DHW-based predictions of bleaching severity (u values) were significantly different between hotspots and thermal refugia was tested using a Wilcoxon sum rank test on the subset of nodes that intersected with the area of these thermal regimes identified through the satellite data analysis described above.
We also ran five sensitivity analyses to address assumptions of our modelling approach.4) Transformation to fit beta distribution: we reran the iterative test above using gaussian R-INLA models with raw severity scores (0,1,2,3) as the response variable (S5 Fig) .(5) Equidistance between ordinal bleaching severity score: we tested for spatial effects in DHW-bleaching based on raw ordinal data using a cumulative link model (bleaching_severity ~DHW + Latitude), where latitude is modelled as a proxy of gradient between hotspots in the southwest and thermal refugia in the north (S6 Fig) .Together these analyses provide support the findings of the model in our main result.

Acropora mass bleaching and simulated marine heatwave experiment
Mass bleaching data captures trends at the coral assemblage level, missing genus-and specieslevel heat tolerance responses, especially for Acropora, a thermally sensitive yet keystone taxon of reef-building corals [48].To address this, we performed an Acropora-specific spatial analysis of bleaching for the 2010 marine heatwave, and a simulated marine heatwave tank experiment on a common corymbose Acropora species, A. digitifera.
Widespread data on Acropora mass bleaching (including corymbose Acropora) was only available from 2010 [49], and thus spatial trends do not reflect the residual long-term patterns discernible from the 20-year assemblage-level bleaching data used above.Acropora bleaching data were reported as the percentage of colonies bleached.Therefore, in order to keep the comparison with assemblage-level data robust, we converted percentage bleaching into bleaching severity scores following Virgen-Urcelay and Donner (2023) (see above).We then tested for spatial patterns in heat stress responses using the same Bayesian spatial beta GLM approach as mentioned for the assemblage-level data.To test the influence of transforming raw percentages to severity scores, we also ran a second version of the model based on the raw percentage of colonies bleached.As such, this second model predicts the probability of 100% coral bleaching, while the former predicts the probability of severe bleaching (50% of colonies bleached).Finally, we also ran the model on a subset of the data for corymbose Acropora, specifically.
For the aquarium-based controlled manipulation, we quantified how thermal regimes influence heat tolerance by conducting a long-term 5-week simulated marine heatwave experiment, inducing bleaching and mortality in 1020 replicate fragments from 174 colonies of Acropora digitifera, a widely distributed reef-building coral in Palauan reefs (see S1 Table for sample size by site).Coral colonies were sampled haphazardly at 1-5 m depth but had to be large enough to provide at least 10 replicated nubbins (approximately 5 cm in length)-six for heat tolerance, three for tissue biomass (see below), and one for DNA extraction (not presented here).Fragments from 90 colonies were collected from six replicate northwest-facing reef crests located in hotspot grid cells (3 sites) and refugia grid cells (3 sites, with one chosen to be slightly more protected, but still with a northwest-facing aspect for logistical reasons), with at least 4 km between sites.Sites showed a very similar consistency in their thermal regimes with hotspots showing marginally higher DHW percentile SDs (22.8 across-site average) than thermal refugia (22.4 across-site average) (site-specific values shown in S1 Table ).In comparison to short heat-shock experiments that typically last 1-2 days, the experimental temperature profile used here (S1 Table ) was designed to match the intensity and duration of a natural marine heatwave more closely [34,37], with the assumption that the phenotypic bleaching and mortality responses would be ecologically relevant to the heat stresses that wild coral populations face during bleaching events.All collections were conducted with relevant National Marine Research Permits (Ministry of Natural Resources, Environment, and Tourism: RE-22-11) and individual State Permits (Ngarchelong State, Kayangel State, and Koror State).All work was in accordance with the ethical standards of and given ethical approval by Newcastle University.
After collection between 3 rd and 5 th April 2022, fragments were acclimatised for 7-10 days in tanks maintained at ambient sea temperature (29.79 ± 1.18 SD ˚C).Relative to the site-specific adjusted (see below) climatological stress accumulation thresholds (MMM adj + 1˚C, among-site mean of 30.1˚C), the temperature in stress tanks was increased on 13 th April gradually over two weeks (+ 0.4˚C on day 1 and + 0.5˚C on days 3, 5, 7, and 13) to reach a final bleaching-level temperature of approximately + 2.5˚C.This temperature was maintained for the remainder of the experiment (S1 Fig, S1 Table).In stress tanks, night-time cooling was not allowed as this would be below the set temperature, however, the use of flow-through tank systems allowed an element of natural diel temperature variability due to short midday thermal peaks [50] (S7 Fig) .Tanks were illuminated using aquarium lights at a constant light intensity matching the midday values from a Palauan outer reef (400 μmol photons m −2 s −1 , S1 Table) using a 12 h:12 h diurnal cycle.Following our experimental protocol to determine A. digitifera heat tolerance in Palau [34], no ramping of light intensity to mimic sunrise and sunset was used, as this could confound comparison to that work.Calibrated HOBO Pendant loggers were placed in each tank and recorded temperatures at 10-minute intervals.To relate coral bleaching and mortality responses to accumulated heat stress, not instantaneous temperature, we calculated heat stress in the experiment using the DHW metric.To compare our experimental DHWs for each site with the NOAA coral bleaching forecasts (satellite-derived DHW), we used local in situ temperature data (night-time averages) to adjust the satellite-based Coral-Temp baselines (MMM-maximum of monthly means climatology-from the 0.05˚grid cell encompassing a collection site), following Humanes et al. (2022).This involved adjusting the MMM based on the linear regression between daily 0.05˚sea surface temperature (CoralTemp v3.1) and daily averaged in situ temperatures (recorded from additional U22 HOBO Pro V2 loggers deployed at each collection site).This produced six adjusted local climatological baselines (MMM adj −adjusted MMM), one for each collection site.

Heat tolerance
Since corals can recover from bleaching, a pragmatic definition of heat tolerance is a coral's ability to survive levels of heat stress that would be sufficient to cause bleaching mortality in nature [51][52][53].Notably, we have previously shown that the visual bleaching mortality responses of A. digitifera in the same experimental system used in this study correlate strongly to image-based measurements of bleaching (based on the values of red, green and blue channels), the density of zooxanthellae cells, and pigment concentrations (supplementary materials of Humanes et al. ( 2022)).We did not measure the photochemical efficiency of symbiotic microalgae (Fv/Fm), because this metric has uncertain relationships to coral mortality (the primary variable of interest in our study) and cannot be used as an indicator to differentiate bleached from recently dead corals.The six replicate fragments of each colony to be assayed for heat tolerance were dispersed in random locations among the ten heat stress tanks (which received 5 of the fragments per colony) and one of the two procedural control tanks (which received the remaining colony fragment and was used to control for possible colony handling effects during collection, see below).With more stress tanks (n = 10) than replicate stressed nubbins (n = 5), the fragments were dispersed in such a way that each tank had an equal representation of the six different sites, and a maximum of one fragment from any given colony.If a fragment died in a procedural control tank, held at non-stressful ambient temperature conditions, it was an indication of handling effects for that colony, so all remaining fragments from the colony were removed from the experiment and the colony was not assigned a heat tolerance score (only one colony out of 174 tested had a procedural control fragment that died).The health status of each fragment was scored visually at intervals of between 1 and 3 days.We used the bleaching and mortality index (BMI) [34] adapted from McClanahan et al. ( 2004) [54] to categorise coral bleaching and mortality responses, such that c 1 to c 5 are the proportion of replicate fragments (per colony) recorded as healthy (c 1 ), half bleached (c 2 ), bleached (c 3 ), partial mortality (c 4 ), or dead (c 5 ), and N is the total number of categories (here N = 5).For example, a colony whose replicate fragments are either all healthy or all dead at a certain time point, will have a BMI value of 0 or 1, respectively.
The BMI of a particular colony is representative of only a single timepoint and will change throughout the heat stress exposure.Therefore, we calculated DHW 50 , a colony-specific metric of heat tolerance analogous to the effective dosage (EC 50 ) metric used in toxicology [55] to determine the dosage at which 50% of the fragments for each colony experiences an effect.DHW 50 for each colony was calculated by fitting a logistic dose response curve for all colonies, allowing a different intercept and slope per colony, and then computing the DHW value at which the fitted curve reaches 50% of the full bleaching and mortality response (i.e., a BMI value of 0.5) using the 'qlogis' function in R.This builds on coral heat stress experiments using the inverse average BMI through time as a measure of heat tolerance (1 -mean BMI).To allow for comparisons with previous work [34,56], we also computed the regression between the compliment of mean BMI (a unitless measure of heat tolerance) and DHW 50 Intermediate BMI values (i.e., not zero or one) can reflect different combinations of health status across a colony's replicate fragments.For instance, a colony having a BMI of 0.5 could be achieved if all replicate fragments are bleached, if half are healthy and half are dead, or other combinations of health status scores.Given this limitation of BMI metric, we also tested for differences in the effect of DHW (continuous fixed effect) on heat tolerance (health status scores, ordinal response) between hotspots and thermal refugia (2-level categorical fixed effect) using a logistic mixed effect cumulative link model, accounting for site and colony ID with random intercepts, in the form: Status ~DHW + Region + (1|Colony_ID) + (1|Site).The logit cumulative link model (clmm) was set with the 'flexible' threshold parametrisation which does not impose any a priori restrictions on the 'distance' between ordinal scores (S9 Fig).

Colony size and tissue biomass
We also assessed regional differences in coral size and tissue biomass, as well as their influence on heat tolerance, from the same set of aforementioned colonies.Notably, we measured tissue biomass as a measure of energy reserves rather than measuring proteins and lipids directly due to logistical constraints, and the fact that a previous study of A. digitifera in Palau (Platz et al. in preparation) showed a strong correlation between tissue biomass and heat tolerance, supporting previous studies [21,57].For each sampled A. digitifera colony, a downward-facing image was taken prior to sampling with a ruler for scale.These images were then analysed using the SizeExtractR workflow [58] to determine the geometric mean diameter of each colony.This is achieved by manually outlining the irregular-shaped coral colony and allowing the SizeExtractR workflow to extract the maximum and minimum perpendicular diameters from this 2D region of interest and compute the geometric mean diameter to account for variation in shape across corals ( p D max ×D min ).Three additional fragments were sampled for each colony to determine tissue biomass, as this is known to have an influence on energy reserves and potentially the ability to survive extended marine heatwave stress [21].Total tissue biomass for each fragment was determined as ash-free dry weight, or the difference between dry weight and ash weight.To achieve this, fragments were dried in an oven at the Palau International Coral Reef Center at 60˚C until constant weight (approximately 3 days).Branch tips were removed so that fragments were approximately cylindrical and surface area was estimated from total length and midpoint width recorded using vernier callipers (mean of 1.55 cm ± 0.43 SD and 0.88 cm ± 0.18 SD, respectively).Each branch was crushed using a clean hammer-style pollen press and the powdered sample was added to a pre-weighed vial.In Newcastle University, these samples were placed in clean porcelain crucibles, reweighed, and burned in a muffle furnace for 2 hours at 450˚C to remove organic matter.After being transferred to a desiccator to cool down to room temperature, the samples were reweighed to determine the weight of non-organic matter.Ash free dry weight was standardised to tissue surface area of each branch based on the length and midpoint width assuming the area of a cylinder, and colony tissue biomass values were assigned as averages across the three branches per colony.
Differences in colony size and tissue biomass (response variables) between thermal refugia and hotspots (2-level categorical fixed effect) were tested using linear mixed effects models under log-link, accounting for site-level differences by fitting a random intercept by site.The effect on heat tolerance (response variable, DHW 50 ) of the interaction between either (a) colony size (log-transformed) and thermal regime or (b) tissue biomass (log-transformed) and thermal regime was also tested using linear mixed effects models, accounting for site-level differences using random intercepts.For these interaction models, backward stepwise model selection was used to remove non-important predictors with P values > 0.1, leaving the final models as (a) DHW 50 ~log_tissue_biomass × region + (1|site), and (b) DHW 50 ~log_colo-ny_size + (1|site).

Persistent thermal refugia and hotspots
At a first assessment, Palau's coral reefs appear to have broadly similar historic satellite-based thermal environments, with warm season SST climatologies, median daily SST, and daily SST variability (99% range from 0.5 th to 99.The three strongest marine heatwaves on record in Palau occurred in 1998, 2010, and 2017.During these events, thermal refugia were consistently exposed to less-intense heat stress than hotspots (2˚C-weeks less on average), despite absolute DHW values still reaching 4-6˚Cweeks, or NOAA Bleaching Alert 1 [23].While local-scale thermal refugia in Palau are not immune to marine heatwaves, they experience less intense heat stress when such events do occur.Cool years without any accumulated heat stress (i.e., DHW of zero) have been far more frequent at thermal refugia (8-12 years since 1980) than at hotspots (only 3-4 years since 1980).Hotspots have been subjected more frequently to years with low-level heat stress.

Thermal refugia show lower and more variable resistance to mass bleaching
Given historical patterns of mass coral bleaching and marine heatwave stress, we quantified spatial trends in the effect of contrasting thermal regimes on assemblage-wide bleaching severity.As expected [12,47], DHW was an important predictor of coral bleaching in Palau (Bayesian DHW slope: 0.362 [0.287-0.439],mean and 95% credible intervals).For a given heat stress dosage, the probability of severe assemblage-wide bleaching was highest in the north of Palau.This is consistent with observations of more severe bleaching in the north in 2010 [49], where thermal refugia are located.This result aligns closely to the delineation of thermal refugia and hotspots, with significantly less-severe bleaching predictions at hotspots (Fig 2D) , this heat stress exposure would predict only mild bleaching at this thermal refugium site (1-10% bleaching, severity code 1).However, of these four survey records, three reported severe bleaching (11-50% bleaching, severity code 3), and one reported moderate bleaching (11-50% bleaching, severity code 2), much worse than would the DHW prediction.Moreover, the variability of severe bleaching predictions for a given DHW were higher at thermal refugia than hotspots, with 2.7-fold greater standard deviation, 2.3-fold greater range, and 3.1-fold greater interquartile range of spatial random field (SRF) values (Fig 2D).
For the Acropora component of the assemblage during the 2010 marine heatwave, the mass bleaching response followed the same trend as the broader assemblage, with slightly higher resistance to a given DHW at hotspots compared to thermal refugia (S13 Fig).However, this effect was small with on average a 4% lower chance of severe bleaching at hotspots than thermal refugia.Notably, this spatial pattern of Acropora mass bleaching was robust to the transformation of raw bleaching percentages (percentage of colonies bleached) into bleaching severity scores (S14 Fig) , and for the subset of corymbose Acropora (S15 Fig).

Thermal refugia show higher and more diverse species-specific heat tolerance
The long-term 5-week marine heatwave experiment conducted on fragments of Acropora digitifera colonies reached a final DHW of 16.2 ± 0.2˚C-weeks (mean ± SD across tanks).We found marked differences in Acropora digitifera heat tolerance between thermal refugia and hotspots.Although, the progression of bleaching and mortality responses among colonies was highly variable, heat tolerance, measured as DHW 50 (the DHW dosage causing a 50% bleaching and mortality response or BMI = 0.5), was significantly higher at thermal refugia (Fig 3B).
Compared to colonies from hotspots, those from thermal refugia could withstand an additional 0.7˚C-weeks (0.2-1.3 C-weeks 95% confidence interval) of heat stress before the onset of bleaching and mortality (Fig 3A).This colony-level finding was corroborated at the level of individual replicate fragments with an ordinal analysis on raw colony fragment health status categories (S9 Fig) , showing a that the DHW-mediated decline of fragment health status was delayed by 0.8˚C-weeks at thermal refugia compared to hotspots.This contrasts with our finding of lower assemblage-wide heat stress tolerance at thermal refugia.Heat tolerance variability was greater within than between thermal regimes.For hotspots, there was a ΔDHW of 4.1˚C-weeks (2.5-5.8˚C-weeks95% CI) between the bleaching and mortality responses of the least-and most-tolerant 10% of colonies (Fig 4A).This translates to a categorical shift in the NOAA Bleaching Alert system (e.g., Alert Level 1 to 2).However, heat tolerance variability was 25% greater at thermal refugia (Fig 4B ), with a ΔDHW of 5.2˚Cweeks (3.5-7.2˚C-weeks95% CI) between the upper and lower deciles of the population.There was an indication that the greater difference in heat tolerance among colonies from thermal refugia was due to the presence of both the most-sensitive and most-tolerant colonies, yet we detected no significant differences in coral heat tolerance between thermal regimes either for the least-tolerant (GLMM, P = 0.79) or most-tolerant decile subsets (GLMM, P = 0.19).

Phenotypic predictors of species-specific heat tolerance
As an indication of potential phenotypic drivers of heat tolerance, we also assessed colony size and tissue biomass.The A. digitifera colonies sampled for this study-representative of the largest colonies available to sample at these sites-had significantly larger sizes at thermal refugia than at hotspots (Fig 5A).We found no significant relationship between colony size and heat tolerance (Fig 5C ), although the estimated slope value was positive (P = 0.096).For tissue

Discussion
Large geographic-scale gradients in thermal exposure play a major role in determining the heat tolerance of marine and terrestrial organisms across latitudinal gradients [4][5][6].However, at the scale of individual populations, patterns of organism tolerance do not necessarily follow large-scale thermal trends due to other environmental drivers [7,8].Here, we identified persistent local-scale thermal refugia and hotspot reefs within 60 km of one another in Palau based on their past exposure to marine heatwave stress.Given the higher exposure, hotspots might be expected to be associated with higher coral heat tolerance based on natural selection and acclimatisation alone, unless other biotic or abiotic factors affecting coral heat tolerance [14,15,17,20] overwhelm the effect of past DHW exposure.At the level of the coral assemblage, patterns in historic mass bleaching severity supported these expectations, as we found that refugium coral assemblages are relatively naive to marine heatwave stress ( refugia and hotspots, we found the opposite trend, with higher heat tolerance at thermal refugia.Ultimately, our study demonstrates that, compared to thermal refugia, hotspot reefs do not necessarily harbour corals with higher levels of heat tolerance, although further research such as common garden experiments would be needed to confirm whether trait differences between regions have a genetic basis.We also show that spatial patterns in the assemblagelevel probability of severe mass bleaching does not necessarily predict species responses, which is not wholly surprising given the complexity of biological systems.Hypotheses on coral heat tolerance and local adaptation have been developed and tested across large spatial thermal gradients [30], between reef habitats with extremely different temperature regimes (e.g., outer reefs versus lagoons [59][60][61]), and in natural field experiments such as neighbouring backreef tidal pools (<500 m apart) [62].These studies show that higher temperatures and sub-daily thermal variability can lead to acclimatisation, selective pressure, and ultimately higher heat tolerance [59,60].However, habitat-specific variability in heat tolerance (e.g., at outer reefs only, as in our study) remains poorly understood across spatiotemporal scales that can be mapped by satellite-borne thermal sensors (>1-10 km, daily).In our study, hotspots were associated with higher assemblage-wide tolerance, higher Acropora (and  Each test was performed using separate linear mixed effect models (LMMs), accounting for site-level differences by fitting a random intercept for each site.LMMs with hypothesised interactions were reduced using stepwise backward selection.(a) Coral colonies were significantly larger at thermal refugia (P < 0.001).(c) There was a positive association between colony size and heat tolerance associated with a moderate level of uncertainty (P = 0.096), and with no interaction between thermal regimes (P > 0.05).(b) No differences in mean tissue biomass were detected between thermal regimes (P = 0.84).(d) Tissue biomass was positively associated with colony heat tolerance at thermal refugia (refugia slope, P = 0.015), but this effect was not present at hotspots (interaction, P = 0.037) where the trend did not differ significantly from the null expectation of a slope equal to zero (P > 0.05).P value significance is given by black symbols for P < 0.1 & > 0.05 (.), P < 0.05 (*), P < 0.01 (**), P < 0.001 (***). https://doi.org/10.1371/journal.pclm.0000453.g005 corymbose Acropora) tolerance to the 2010 marine heatwave, but lower heat tolerance for the species investigated.The assemblage-level trend could be due to local differences in coral species composition (i.e., reefs with more abundant heat-sensitive taxa like Acropora are expected to suffer more severe bleaching impacts [63]).However, Palauan reefs have had widespread temporal stability of Acropora populations [64], and Acropora showed overall slightly higher DHW tolerance at hotspots during the 2010 marine heatwave.Therefore, these results suggest that spatially varying acclimatisation or adaptation may be occurring for some members of the assemblage, although we cannot discount the potential role of spatial differences in overall genus composition, species composition within the Acropora assemblage, or even presence of cryptic Acropora species between thermal regimes [65].However, for the study species A. digitifera the pattern between thermal refugia and hotspots was reversed, indicating that other drivers of heat tolerance likely overwhelm any effect of past DHW exposure, such as differences in Symbiodiniaceae community composition (see more discussion below), or non-thermal environmental variables that affect heat tolerance directly (e.g., via phenotypic plasticity) or indirectly (e.g., via genetic trade-offs with other adaptive traits such as skeletal density or reproduction).Variability in A. digitifera heat tolerance was 25% higher at thermal refugia than hotspots (ΔDHW of 5.2 and 4.1˚C-weeks, respectively).Previous work on A. digitifera from a thermally moderate eastern outer reef in Palau detected a ΔDHW of 4.8˚C-weeks [34] between upper and lower deciles of heat tolerance.Combining this with the results presented here, heat tolerance variability may be negatively correlated to historic heat stress exposure in our study system.The most-sensitive A. digitifera colonies were found at thermal refugia but were absent at hotspots, suggesting that thermal selective pressures may be more prevalent at hotspot reefs.Yet the most-tolerant A. digitifera colonies were also absent at hotspots, suggesting a potential role of non-thermal selective pressures influencing heat tolerance (e.g., lower water quality [20,66]).Notably, we also found that heat tolerance variability was widespread.Provided coral heat tolerance has a genetic basis and is heritable [67][68][69][70][71], this indicates a high potential for coral adaptation to climate change across both thermal regimes.
Coral heat stress responses are complex and influenced by biological and environmental factors [12,22,72].Our results show that colonies with thicker tissue, a proxy of energy reserves, could withstand an additional 1˚C-week of heat stress before bleaching and dying, supporting previous findings [21].The larger colony sizes achieved by A. digitifera corals at thermal refugia compared to hotspots suggests that there could be differences in the population age structure, growth rates (in terms of linear extension), and/or calcification rates.Notably, the removal of colony fragments was easier at thermal refugia, indicating lower skeletal density which has been suggested as a potential trade-off with A. digitifera heat tolerance in Palau [56].This lower skeletal density at thermal refugia than at hotspots could be explained by faster linear extension and/or lower calcification rates.All organisms must balance energy allocation between resource-intensive traits, and the energetic cost of physiological stress from being in high energy environments (e.g., tissue repair [73]) or high investment in reproduction [74] could compromise the tolerance of corals to marine heatwaves.While these may be factors in the observed trends in our study, further experimental evidence would be required.
Certain symbiont taxa (e.g., Durusdinium spp.) can provide corals with higher heat tolerance [75,76].However, previous studies of Palau's outer reefs have found A. digitifera is ubiquitously dominated by Cladocopium C40 symbionts with no link to heat tolerance [34,71], and four other coral genera also host predominantly Cladocopium spp.symbionts [59].These studies suggest that Durusdinium spp. is rare on Palauan outer reefs despite prevalence in lagoons [77], reflecting findings from other Pacific islands [78] and suggesting that major differences in Symbiodiniaceae communities are unlikely to be a major driver of our observed heat tolerance trends between hotspots and thermal refugia (although this cannot be ruled out).Thus, other mechanisms may provide more parsimonious explanations of heat tolerance, including host genetics or phenotypic plasticity driven by environmental parameters other than past marine heatwave stress.
Despite the sampled areas being in the same geomorphological zones of the reef with the same depth range, numerous environmental parameters are known to influence Palau's coral reefs, including wave exposure [79] and water quality [80].Such factors can influence coral bleaching [20] and can be correlated to thermal environments.The high A. digitifera heat tolerance at thermal refugia (given that acclimatisation and natural selection should lead to higher tolerance at hotspots in absence of pressures other than DHW), together with other phenotypic and morphological differences between thermal stress regimes, suggest that other environmental factors have overwhelmed the influence of past heat stress exposure on A. digitifera heat tolerance in our study.Indeed, the altered colony morphology observed here between thermal regimes suggests that differences in currents or nutritional opportunities exist, which could explain the patterns in coral heat tolerance we detected.It could prove insightful to categorise physical reef environments-as for the Caribbean [81]-and then revisit the role of thermal environments on coral heat tolerance.We may yet discover that the pattern of higher tolerance existing in more heat stress exposed regions holds true in environments that are truly physically comparable.
Thermal refugia have previously been identified across global [82] and regional [24] scales (>100-1000 km).Yet here we find persistent thermal refugia within 60 km of hotspots.Notably, the methods of identifying thermal regimes are not always consistent, with studies basing their definitions on DHW thresholds [24], the probability of DHW events [82], DHW percentile thresholds [83], or a combination of DHW exposure and consistency of DHW percentiles as in our study.Nonetheless, local variability in thermal regimes at the same spatial scales at which coral populations operate (e.g., self-seeding [84] and dispersal within ~100 km [85]) will likely be important for coral reef persistence and the feasibility of novel reef management approaches.For instance, climate-smart marine protected area networks that aim to boost adaptation rates [30,86,87] could benefit from local-scale variability in thermal regimes detectable from satellites, by protecting key dispersal links between reefs.Similar benefits could apply to selective breeding and assisted gene flow interventions that aim to boost natural rates of climate adaptation by propagating more tolerant individuals.For instance, translocating coral colonies from reefs with different characteristics may be more feasible if donor reefs are closer to recipient reefs (e.g., lower costs, less risks of contamination, fewer permitting issues).However, the science underpinning these novel management approaches is still in its infancy, and they have yet to be implemented in practice.
A dilemma for designing assisted gene flow or selective breeding interventions emerges in how to source heat tolerant genotypes, which could be found either at specific reefs with desirable environmental characteristics (as in the case of assisted gene flow [88]) or by performing assays on numerous colonies at any given reef (as in the case of choosing parental colonies for selective breeding [89]).For Palau, the fact that Acropora heat tolerance was spatially widespread (this study and [90]), and that the variability of A. digitifera heat tolerance was greater within than between sites, suggests that the search for heat tolerant genotypes within the archipelago may benefit from assaying high numbers of colonies within available sites, and not relying solely on the specific hotspots or thermal refugia as sources of colonies.However, this advice may only be applicable in the context of remote coral reef systems with limited habitatspecific environmental gradients, and populations with high genetic diversity (e.g., many remote Pacific reefs).Approaches for finding heat tolerant genotypes may differ for larger reef systems that span considerable thermal gradients (see [24][25][26]), or for populations that have already suffered substantial declines with low remaining genetic diversity (e.g., Acropora palmata in the Caribbean [88]).

Conclusion
Locations with consistent exposure to specific climatic stressors are often thought to harbour individual organisms with higher tolerance to such stressors through acclimatisation and selection.While this pattern holds true across broad climatic gradients [5,25], we find that it can break down at local within-population scales.While marine heatwave stress is the trigger for mass coral bleaching, and past exposure seems to predict local-scale assemblage-level patterns in the occurrence of severe bleaching in Palau, we demonstrate that this pattern does not hold for the coral species studied.This species-level result reflects findings from other ecological systems [8], where the latitudinal-scale thermal tolerance theories were shown to give wrong predictions when applied at local scales.Together, our study highlights the importance of scale when understanding heterogeneity in biological responses to climate change.It follows that a better understanding of local-scale ecological phenomena is needed to aid in the design of climate-smart initiatives capable of addressing global-scale adaptive management problems, with critical developments already underway (e.g., [81]).bleaching (e-f, 11-50%), severe bleaching (g-h, >50%).Because it is not possible to explicitly account for spatial correlation with an ordinal analysis, the results shown are for a cumulative link model in the form: bleaching_severity ~DHW + Latitude.Here we use latitude as a proxy of the change in thermal regimes, showing hotspots (red) and thermal refugia (blue) with representative latitudes of 7 and 8˚N, respectively.This result reflects that in the main study that mass bleaching is predicted to be lesser for a given DHW in hotspots rather than thermal refugia.showing probability of the five different health status categories: healthy (a), partial bleaching (b), bleached (c), partial mortality (d), dead (e).These predictions are from the cumulative link model in form: health_status ~DHW + Region, with a random intercept for colony ID (to account for multiple fragments per colony), and site (multiple sites per region).There is a significant effect of DHW (z = 97.01,P < 0.001) on the declining health status of coral fragments, and significantly more DHW resistance at thermal refugia than hotspots (z = 2.54, P < 0.05), which equates to an additional 0.8˚C-weeks for a probability of 0.5 (dashed horizontal line).highlighting the values of the six sampling sites which vary between 28.98˚C and 29.03˚C but with slightly higher temperatures at hotspots (red points, H in legend) than thermal refugia (blue points, R in legend).(c) The variation in daily SST data for all Palau coral reef pixels shown as the difference between the 0.5 th and 99.5 th percentile SSTs, highlighting the values of the six sampling sites which vary between 3.40˚C and 3.49˚C but with slightly higher variability at thermal refugia (blue points, R in legend) than hotspots (red points, H in legend).Compared to what would have been predicted based on DHW alone (i.e., without a spatial random effect), areas of higher bleaching susceptibility are shown in red (underpredictions), and areas of higher bleaching resistance are shown in blue (overpredictions).(d) Thermal refugia have marginally higher probability of Acropora mass bleaching for a given heat stress dosage (u values) than hotspots (Wilcoxon sum rank test).(e) The difference in Acropora bleaching susceptibility to the 2010 marine heatwave at thermal refugia (blue) versus hotspots (red) translates to only a small yet significant 4% increase in the probability of total Acropora mass bleaching (P(TB)).

Fig 1 .
Fig 1. Characterisation of thermal regimes across Palau showing conceptual diagrams (a-b) and empirical results (c-e) based on the typical heat stress a reef receives and its consistency in performing this way relative to other reefs.(a)The typical heat stress that reefs receive relative to one another, or degree heating week (DHW) anomalies, were calculated from annual maps of peak DHW heat stress.The consistency of a reefs performance relative to other reefs was calculated by first converting annual DHW maps to percentiles, and then computing the standard deviation (SD) of DHW percentiles through time for each reef.(b) A thermal consistency biplot is used to compare these two metrics: typical heat stress (x-axis) and inconsistency (y-axis).Equal-sized tercile subsets of each variable demarcate distinct thermal regimes,

( 1 )
Potential outliers: we reran the beta R-INLA model but without the four most northern survey records (all from 1998) which showed moderate (n = 1) or severe bleaching (n = 3) and could have biased our results (S3 Fig).(2-3) Balanced design and mesh resolution: as the three main survey years had unbalanced sampling effort (1998: n = 29, 2010: n = 80, and 2016: n = 63), we randomly subsampled the more-surveyed years (2010 and 2016) 100 times, resulting in 100 balanced datasets, and reran a beta R-INLA model for each with a lower mesh resolution where maximum triangle edge length was set to 8 km instead of 2 km (S4 Fig).( 5 th percentiles) all differing by < 0.1˚C (S10 Fig).However, distinct spatial trends emerge in the historic exposure of reefs to DHW and their performance relative to one another (i.e., 36 maps of spatial DHW anomalies and percentile rankings, Fig 1A), shown on the thermal consistency biplot (Fig 1B).The average DHW anomaly of a reef through time (Fig 1B x-axis) captures the typical heat stress it receives relative to other reefs, varying by almost 1˚C-week across reefs (Fig 1C).The temporal inconsistency of a reef's relative thermal environment (Fig 1B y-axis) is captured by the standard deviation of spatial DHW percentiles through time, and shows a hump-shaped pattern (Fig 1D).There is a clear gradient of thermal regimes characterised by thermal refugia, hotspots, and thermally moderate reefs (Fig 1D).Of all the reef grid cells in Palau, 13% emerge as persistent thermal refugia distributed sparsely in the north, an area more exposed to the open ocean, while 18% emerge as persistent hotspots located in a southwestern area associated with high water retention of the shallow inner lagoon (Fig 1E).Notably, a small proportion of reefs located at the southern tip of Palau showed the highest interannual inconsistency of DHW percentiles (Fig 1D, yellow reef cells) but with a moderate mean DHW anomaly through time (Fig 1E, grey points and red points at the peak of the hump-shaped pattern).However, the highest and lowest DHW percentiles (leading to high SD) occurred in low heat stress years when DHW was very stable across Palau, demonstrative of moderate thermal regimes.
, characterised by lower model residuals (i.e., vertical distance between points and dose-response curve in Fig 2A) (S11 Fig), and a higher dose response of approximately 1.7˚C-weeks (S12 Fig).This result was robust under sensitivity analyses of outliers (S3 Fig), interannual differences in survey effort (S4 and S5 Figs), mesh resolution (S4 and S5 Figs), and data transformation (gaussian models in S5 Fig; ordinal models in S6 Fig).To illustrate the higher residual bleaching susceptibility at thermal refugia (Fig 2D), consider the thermal refugium reef where surveys were conducted in 1998 (i.e., the four most northerly bleaching survey records).The DHW exposure was 4.4˚C-weeks (S1a Fig).Based on the overall dose-response curve (Fig 2A)

Fig 2 .
Fig 2. Assemblage-wide bleaching risk across Palau based on historic bleaching survey data and satellite-derived heat stress.(a) The predicted effect of DHW on severe coral bleaching is shown as the mean and 95% credible intervals with notable overprediction and underprediction for some bleaching survey records (points shown with vertical jitter).(b) Overlay of historic bleaching survey records on top of thermal regime classifications.Spatial correlated uncertainty in the bleaching predictions shown in (a) is calculated across a high-resolution Delaunay triangulation mesh of the study area comprising 5,710 nodes (c).(d) The Gaussian Markov spatial random field reports the spatial correlated uncertainty (residual bleaching susceptibility) as u values (see linear combination equation).Compared to what would have been predicted based on DHW alone (i.e., without a spatial random effect), areas of higher bleaching susceptibility are shown in red (underpredictions) and areas of higher bleaching resistance are shown in blue (overpredictions).Notably, the blue-white-red colour scheme here (d) is distinct from the hotspots/refugia colour scheme in Fig 1E.(e) Thermal refugia have higher mass bleaching susceptibility for a given heat stress dosage (u values) than hotspots (Wilcoxon sum rank test) and a higher diversity of responses.The Palau shoreline shapefile was from the NOAA National Centre for Coastal Ocean Science [41] (https://products.coastalscience.noaa.gov/collections/benthic/e102palau/).https://doi.org/10.1371/journal.pclm.0000453.g002

Fig 1 )Fig 4 .
Fig 4.Within-population variability in heat tolerance for Acropora digitifera corals at hotspots (a) and those at thermal refugia (b), expressed as differences between the least-tolerant decile (paler shading) and most-tolerant decile (darker shading) of the population (n = 9 colonies for each subsample in each thermal regime) in terms of DHW tolerance (ΔDHW, mean and upper and lower 95% confidence limits) evaluated at BMI = 0.5 (grey dashed line).https://doi.org/10.1371/journal.pclm.0000453.g004

Fig 5 .
Fig 5. Phenotypic drivers of heat tolerance between thermal regimes.(a,b)Differences in trait values between thermal refugia (blue) and hotspots (red), and influence of traits on heat tolerance (c,d).Each test was performed using separate linear mixed effect models (LMMs), accounting for site-level differences by fitting a random intercept for each site.LMMs with hypothesised interactions were reduced using stepwise backward selection.(a) Coral colonies were significantly larger at thermal refugia (P < 0.001).(c) There was a positive association between colony size and heat tolerance associated with a moderate level of uncertainty (P = 0.096), and with no interaction between thermal regimes (P > 0.05).(b) No differences in mean tissue biomass were detected between thermal regimes (P = 0.84).(d) Tissue biomass was positively associated with colony heat tolerance at thermal refugia (refugia slope, P = 0.015), but this effect was not present at hotspots (interaction, P = 0.037) where the trend did not differ significantly from the null expectation of a slope equal to zero (P > 0.05).P value significance is given by black symbols for P < 0.1 & > 0.05 (.), P < 0.05 (*), P < 0.01 (**), P < 0.001 (***).
represent locations where observed bleaching severity was worse than would have been expected based on DHW alone, and negative SRF values show locations where bleaching was less severe than DHW-based predictions.Results are based on 100 independent runs of INLA models set up as follows: (i) Imbalances in sampling effort across bleaching surveys were accounted for by removing years with less than 30 records.This left 1998, 2010, and 2017 with 29, 80 and 63 records, respectively.Then for each replicate INLA model, 29 random records for 2010 and 29 random records for 2017 were chosen to keep and the rest were discarded to achieve a fully balanced design of N = 29 per year over 3 years.(ii) Each replicate INLA model was given a unique low-resolution Delaunay triangulation mesh with a maximum triangle edge length of 8 km achieving 465-502 nodes across replicate meshes, compared with the high-resolution mesh used in the main manuscript with max edge length of 8km achieving 5,710 nodes overall.Results presented here show: (a) the mean spatial random field across all replicate INLA runs (gridded at 5km CoralTemp resolution before averaging to account for mismatch in meshes between replicate runs); (b) the difference in spatial random field values between hotspot and thermal refugia reef cells in terms of the median (bold horizontal line), interquartile range (box) and minimum-to-maximum (vertical feather) computed first for each INLA model and then averaged across all replicate runs for the figure; and (c) Histogram of statistically significant differences in SRF values between hotspots and refugia (as in b) in terms of a P value for each replicate run based on Wilcoxon sum rank tests with an alpha level of 0.05 with Bonferroni correction (red dashed line) and log 10 -scale x-axis.(EPS) S5 Fig. Spatial trends in the DHW-bleaching relationship, and the sensitivity of these results to sampling effort and mesh resolution, as in S3 Fig. S3 Fig figure description also applies to this figure, except bleaching severity scores were used as a response variable (0, 1, 2, 3) with a Gaussian error distribution rather than the beta distribution used for the main manuscript and S3 Fig. (EPS) S6 Fig. Ordinal analysis of mass bleaching data, showing probability of the four different mass bleaching severity categories: no bleaching (a-b, 0%), mild bleaching (c-d, 1-10%), moderate marine heatwave experiment.(A) Temperature profiles for the ten replicate heat stress tanks (Tanks 1-8, 11, and 12) and procedural control tanks (Tanks 9 and 10), sharing the colour legend with panel B. The local climatological baseline (MMM adj ) and the corresponding stress accumulation threshold (MMM adj + 1˚C) are shown for each of the 6 coral collection sites.Specific values for each site are shown in S1 Table.(B) The average accumulated heat stress profile (DHW-degree heating weeks) is shown as the average across all sites in each experimental tank (colour legend).The shaded region shows the total range of absolute DHW across all sites and tanks.(EPS) S8 Fig. Relationship between the compliment of mean BMI (i.e., a measure of heat tolerance) and DHW 50 , the heat stress dosage at which a colonies response is predicted to pass a BMI of 0.5.(EPS) S9 Fig. Ordinal analysis of experimental bleaching and mortality data for Acropora digitifera, of spatial differences in daily sea surface temperature (SST) data from 1985 to 2020.(a) The overall distribution of SSTs for all six sampling sites ordered from bottom to top by increasing latitude.(b) Median temperature for all Palau coral reef pixels, Model residuals from the DHW dose response shown in Fig 2A, for all bleaching survey records in the three most-extensive bleaching survey years conducted in (a) 1998, (b) 2010, and (c) 2016.Residuals are calculated as the vertical distance between the observed bleaching severity (0-1, transformed from severity scores) and the predicted DHW dose response curve in Fig 2A.Colours for each year are shown with a blue-white-red scale, centred on the mean residual for that year.Satellite 5km grid cells identified as hotspots and thermal refugia are shown in dark grey and black, respectively.Comparison of mean residuals of hotspots versus thermal refugia for each year are given in white text on each map.(EPS) S12 Fig.The heat stress (degree heating week DHW) dose response curve for severe mass bleaching for all of Palau (black), only hotspots (red) and only thermal refugia (blue).Predictions are based on the logistic model presented in Fig 2, with an intercept value of -2.08, a slope value of 0.36, and an intercept offset for hotspots and thermal refugia of -0.24 and 0.36, respectively.The DHW required to elicit a 50% chance of severe mass bleaching at hotspots (DHW 50 = 6.4˚C-weeks) is 1.7 DHWs greater than that of thermal refugia (DHW 50 = 4.7˚C-weeks).(EPS) S13 Fig. Risk of severe Acropora bleaching (i.e., >50% of corals bleached) during the 2010 marine heatwave to be compared to the assemblage-level result in the main manuscript.(a) Bleaching survey records from 2010 are overlayed on thermal regime classifications, with (b) the same points also shown on the high-resolution Delaunay triangulation mesh (comprising 4,325 nodes), used to calculate spatial correlated uncertainty in predictions of severe Acropora bleaching.(c) The Gaussian Markov spatial random field reports the spatial correlated uncertainty (residual bleaching susceptibility) as u values (see linear combination equation).Compared to what would have been predicted based on DHW alone (i.e., without a spatial random effect), areas of higher bleaching susceptibility are shown in red (underpredictions), and areas of higher bleaching resistance are shown in blue (overpredictions).(d) Thermal refugia have marginally higher Acropora mass bleaching susceptibility for a given heat stress dosage (u values) than hotspots (Wilcoxon sum rank test).(e) The difference in Acropora bleaching susceptibility to the 2010 marine heatwave at thermal refugia (blue) versus hotspots (red) translates to only a small yet significant 4% increase in the chances of severe Acropora mass bleaching (P(SB)).(EPS) S14 Fig. Probability of total Acropora bleaching (i.e., percentage of corals bleached, reaching up to 100% of colonies) during the 2010 marine heatwave to be compared to S13 Fig.(a) Bleaching survey records from 2010 are overlayed on thermal regime classifications, with (b) the same points also shown on the high-resolution Delaunay triangulation mesh (comprising 4,325 nodes), used to calculate spatial correlated uncertainty in predictions of severe Acropora bleaching.(c) The Gaussian Markov spatial random field reports the spatial correlated uncertainty (residual bleaching susceptibility) as u values (see linear combination equation).
of total corymbose Acropora bleaching (i.e., percentage of colonies bleached, reaching up to 100% of colonies) during the 2010 marine heatwave.(a) Bleaching survey records from 2010 are overlayed on thermal regime classifications, with (b) the same points also shown on the high-resolution Delaunay triangulation mesh (comprising 2,595 nodes), used to calculate spatial correlated uncertainty in predictions of severe corymbose Acropora bleaching.(c) The Gaussian Markov spatial random field reports the spatial correlated uncertainty (residual bleaching susceptibility) as u values (see linear combination equation).Compared to what would have been predicted based on DHW alone (i.e., without a spatial random effect), areas of higher bleaching susceptibility are shown in red (underpredictions), and areas of higher bleaching resistance are shown in blue (overpredictions).(d) Thermal refugia have higher probability of corymbose Acropora mass bleaching for a given heat stress dosage (u values) than hotspots (Wilcoxon sum rank test).(e) The difference in corymbose Acropora bleaching susceptibility to the 2010 marine heatwave at thermal refugia (blue) versus hotspots (red) translates to a 10% increase in the probability of total corymbose Acropora mass bleaching (P(TB)).(EPS) S1 Table.Simulated marine heatwave experiment.Associated metadata following [91].(XLSX) ). Equal-sized tercile subsets of each variable demarcate distinct thermal regimes, [41]nodes).The model estimates the probability of severe coral bleaching (P(SB i ), i.e., where persistent thermal refugia and persistent hotspots are in the lower left and right corners, respectively.(c)Maps of these two metrics for all Palauan reefs.(d)Eachreefs position in the thermal consistency biplot, with an overall hump-shaped pattern across all reefs.(e)Persistentthermal refugia are located among the northern reefs, while persistent hotspots dominate in the southwest.The Palau shoreline shapefile was from the NOAA National Centre for Coastal Ocean Science[41](https://products.coastalscience.noaa.gov/collections/benthic/e102palau/).