Environmental Variability and Threshold Model’s Predictions for Coral Reefs

Current models of the future of coral reefs rely on threshold (TM) and multivariate environmental variability models (VM) that vary in how they account for spatial and temporal environmental heterogeneity. Here, a VM based on General Additive Model (GAM) methods evaluated the empirical relationships between coral cover (n = 905 sites pooled to 318 reef cells of the Western and Central Indian Ocean Provinces) and 15 potentially influential variables. Six environmental and one fisheries management variables were selected as significant including SST shape distributions, dissolved oxygen, calcite, and fisheries management. Common predictive variables, including cumulative degree-heating weeks (DHW), pH, maximum light, SST bimodality and rate of rise, and two multivariate metrics were either weak or not significant predictors of coral cover. A spatially-resolved 2020 baseline for future predictions of coral cover within 11,678 reef ∼6.25 km2 cells within 13 ecoregions and 4 fisheries management categories using the 7 top VM variables was established for comparing VM and TM coral cover prediction for the year 2050. We compared the two model’s predictions for high and low Relative Concentration Pathway (CMIP5; RCP8.5 and 2.6) scenarios using the four available future-cast SST variables. The excess heat (DHW)-coral mortality relationship of the TM predicted considerably lower coral cover in 2050 than the VM. For example, for the RCP8.5 and RCP2.6 scenarios, the decline in coral for the TM predicted was 81 and 58% compared to a 29 and 20% for the VM among reef cells with >25% coral cover in 2020, if a proposed optimal fisheries management was achieved. Despite differences, coral cover predictions for the VM and TM overlapped in two environmental regions located in the southern equatorial current region of the Indian Ocean. Historical and future patterns of acute and chronic stresses are expected to be more influential than cumulative heat stress in predicting coral cover, which is better accounted for by the VM than the TM.


INTRODUCTION
The ecosystem services provided by coral reefs are under threat from a number of forces of which climate warming, eutrophication, and fishing are among the most wide-spread concerns (Cornwall et al., 2021;Donovan et al., 2021). The two main threatened ecological services are the organic production that produce fisheries and the inorganic production that produces the calcium carbonate reef framework. Net reef growth is particularly critical for providing shallow-water habitat and shoreline protection, which are threatened by the current rates of sea level rise (Perry et al., 2018). Therefore, the cover of calcifying coral and algal functional groups and the factors that influence them are of concern for societies dependent on the above reef services. This is particularly true in the many tropical regions that are dependent on fisheries and coastal tourism for food and economic security (Spalding et al., 2017;Selig et al., 2019).
Current predictions of generic coral health and cover are based on thermal exposure predictions (van Hooidonk et al., 2020). Predictions are frequently based on climate model's predicted temperature deviations from historical warmest seasonal months from satellite temperature histories and rate of rise in seasurface temperatures and the production of future "excess heat" (Cornwall et al., 2021). This accumulated excess heat above +1 • C of warm-season temperatures is the current threshold model (TM) best-estimate of both the historical and future exposure of corals to thermal stress. This excess-heat model can be modified to include and account for other environmental factors. Common modifications include between-year thermal thresholds, coral acclimation and genetic adaptation rates, light attenuation, and future climate model scenarios (Donner et al., 2005;McManus et al., 2020;Gonzalez-Espinosa and Donner, 2021). Most of this class of models can generally be described as a TM because their core predictions rely on a threshold to provoke the coral stress responses that often but not always coincident with +1 • C above summer mean temperatures (Donner, 2011).
Criticisms of threshold models are (1) they poorly account for local variability and environmental contexts (Chollett et al., 2014;McClanahan et al., 2020b), (2) overly reliant in treating exposure as a single latent surface heat stress variable , (3) do not account well for the variable responses in corals to the 1 • C hot season threshold or other heat stress metrics (DeCarlo, 2020;McClanahan et al., 2020a), (4) lack the multivariate aspects of stress that include both thermal and nonthermal stress elements (Maina et al., 2011), and 5) rely solely on coral cover responses to heat stress and treat cover as a single pooled entity without taxa-specific niches, life history, or regional characteristics (Cacciapaglia and van Woesik, 2016;Darling et al., 2019;McClanahan et al., 2020a,b). These critiques often arise from models that can generally be described as environmental variability models (VMs).
Defense of the threshold models are that they were created to illuminate the general problem of global climate change on coral reefs and less to make specific local reef-scale predictions. TM should, therefore, help to promote policies and actions to reduce future impacts expected from different human responses and climate scenarios. Yet, modifications of these models are currently being used to identifying climate sanctuaries and prioritize conservation spending (Beyer et al., 2018). There are, however, many cases where these models fail to predict climatewarming disturbances and changes on a large scale, which can undermine the confidence of threshold metrics and usefulness for influencing management policies (McClanahan et al., 2015b;DeCarlo and Harrison, 2019;Kim et al., 2019;Mollica et al., 2019). Moreover, over-reliance on TMs prompts concerns about policy foci, as per the differential efficacy of managing local fisheries, eutrophication, and the global atmospheric commons (Abelson, 2020). If models poorly account for impacts of local human action, the agency for change is centered solely on the global and not local coral reef commons.
Consequently, comparing and considering a new generation of climate change models is needed to address the policy debate and scales of action. Specifically, building and testing models that account for the multivariate nature of chronic and acute stress and coral sensitivity relationships in specific contexts (Mumby et al., 2011). This requires less reliance on core theoretical concepts and greater empiricism in building and testing models that use the power of modern statistical algorithms. A greater diversity of models needs to be created, examined, and competed against each other in terms of their ability to predict localscale environmental and coral responses. For example, are there missing but potentially important variables that are not being considered by the TMs that might greatly modify the outcomes of climate change? What affects do the five TM simplifications above have on the ability to predict coral cover? In order to begin this process, we explored seldom considered variables that are currently spatially resolved and available to produce a VM option. Our VM examined coral responses in the context of acute and chronic thermal stress and fisheries management in two faunal provinces of the western Indian Ocean (WIO). Thereafter, a common version of the TM that related excess heat and coral mortality was compared with this new VM in terms of their predictions of coral cover in 2050 using Coupled Model Intercomparison Project Phase 5 (CMIP5) global climate model projections.

MATERIALS AND METHODS
This study evaluated the empirical relationships between coral cover and fisheries management and 15 geospatially available and resolved SST and seawater chemistry variables. Model results were then used to predict coral cover within 13 ecoregions of the Western and Central Indian Ocean in 2020 (WIO). The inputs were the significant SST, chemistry variables and fisheries management categories selected by a machine-learning algorithm described below. Thereafter, using the modeled coral cover initial condition in 2020, cover predictions were made for 2050 based on the four future-cast SST variables available for the CMIP5 for low (RCP2.6) and high Relative Concentration Pathway (RCP8.5) scenarios. The two models' predictions in 2050 were then compared to a DHW-coral mortality relationship or threshold model starting with the same 2020 modelled benchmark (Hughes et al., 2018;Cornwall et al., 2021). Finally, we evaluated the potential impacts of a common fisheries sustainability proposals on coral cover by establishing 30% of the cells as fisheries closures and 70% as restricted fishing in each country, which is generally considered a best-practices fisheries policy (O'Leary et al., 2016).

Study Region and Management
The evaluation included the western Indian Ocean coral reefs or western and central Indian Ocean faunal provinces with a spatial extent of 32 • E, 27 • S to 73 • E −7 • N. The United Nations Environmental Program -World Conservation Monitoring Center (UNEP-WCMC) coral reef distribution data were used to create a ∼6.25 km 2 grid of reef cells for the entire region. The map resulted in 11,678 reef cells that were subsequently classified by ecoregion, nation, and the national and local fisheries management status (Supplementary Table 1). The cell grid was appropriate for the local nature of reef fishing, management, and encompasses the home range size of most large coral reef fishes and also close to the extent of small-scale fishing movements and gear used in this region (McClanahan and Abunge, 2018). Fishery management classification were either restricted (limited gear often excluding spear guns, explosives, and small-meshed or drag nets) or unrestricted (all gears including spearguns and small-meshed nets used) based on national and local laws but also modified by observations of compliances with restriction (McClanahan et al., 2015a). Legally established closures were classified as low and high compliance based on published literature and personal and stakeholder observations of compliance. The management categories for each cell was based on a combination of legal records and maps, personal observations, and communications with colleagues. The dominant management coverage was recorded in the few locations where the reef area was within more than one management system.

Coral Cover Response Variable
Live coral cover was the biological response used in this model and was retrieved from a large database of coral cover collected in 905 sites from 2005 to 2020 compiled by three experienced observers (T.R. McClanahan, N.A. Muthiga, and N.A.J Graham; Supplementary Table 2) and often collected across regional thermal disturbances, such as 2010(McClanahan et al., 2007, 2020a. Coral cover data from multiple survey sites within the grid cells were averaged. Where multiple temporal records existed for a site, only the three most recent records were retained, generating a final cell-based pooled dataset of n = 318 percell means.

Sea-Surface Temperature Variables
We accessed daily sea-surface temperature (SST) data for the 1985-2020 study period from the U.S. National Oceanic and Atmospheric Administration's Coraltemp v3.1 SST website 1 . These SST data were available at 5-km resolution. From these data, a number of SST metrics were calculated for the period prior to the empirical coral cover sampling and used to develop the empirical coral cover-environmental metrics model. The metrics calculated included the mean, standard deviation, cumulative degree heating weeks (cumDHW), skewness, kurtosis, and a bimodality coefficient found to be strongly related to bleaching .
These variables characterize the temporal structure of temperatures experienced by corals and can act as proxies for 1 https://coralreefwatch.noaa.gov/product/5km/index_5km_sst.php either chronic or acute thermal stress. For example, kurtosis is a measure of the thickness of the tails of the SST distribution with negative values indicating increasing thickness. Thus, it represents a measure of the frequency of exposure to SSTs away from the mean and negative values should be associated with some degree of acclimation of corals to unusual temperatures (McClanahan et al., 2020a). Skewness is a measure of the frequency and direction of rare SST events, with high negative and positive values indicating more extreme rare cold and warm SST values. Thus, increased skewness implies more acute stress with more positive values indicating warm and negative values cold extremes. The bimodality coefficient measures the degree of bimodality in SST, and ranges from 0-1, with a value equal to or greater than 0.55 indicating a bimodal or multimodal distribution. Increases in bimodality may represent corals experiencing a stressful 'jump' to extreme temperatures, rather than a more continuous acclimation to thermal variability .

Cumulative Degree Heating Weeks
Degree heating weeks (DHW) were calculated as the sum of positive thermal anomalies (hotspots) greater than 1 • C above the historically warmest month (maximum monthly mean) over a 12-week period (Liu et al., 2014)   , which is the accumulation of daily hotspots > 1 0 C over a 12-week period.
Equation 3 | Calculation of cumulative degree heating weeks, which is an aggregated sum of maximum yearly (t) DHW at a reef location from 1985 to 2020.
To understand the geographical variation in the thermal regime across the WIO, we performed a principal components analysis (PCA) and cluster analysis based on thermal stress metrics. We used the 1985-2020 NOAA CRW data to calculate 7 thermal stress metrics for each 6.25 km 2 reef cell, including the mean, standard deviation, cumulative degree heating weeks (cumDHW), skewness, kurtosis, bimodality coefficient, and the SST rate of rise. The SST rate of rise is the linear regression slope of annual mean SST from 1985 -2020. We conducted a PCA on the 7 thermal stress metrics for 11,678 reef cells using the Factominer package in R (Lê et al., 2008). Subsequently, we conducted a cluster analysis using the clustering large applications algorithm (CLARA) with Euclidean dissimilarities and 50 samples to group coral reef locations based on thermal stress metrics (Kaufman and Rousseeuw, 1990).

Additional Environmental Variables
To broaden the scope of the potential impacts of the environment on coral cover beyond SST, we extracted a number of abiotic variables previously found to influence the distribution and abundance of hard scleractinian corals (Vercammen et al., 2019). Variables extracted from Bio-ORACLE 2 (Tyberghein et al., 2012;Assis et al., 2018) included surface photosynthetically active radiation (PAR), diffuse water-penetration attenuation coefficient at 490nm, calcite concentration, dissolved oxygen, and pH (Supplementary Table 2). Bio-ORACLE datasets were obtained from satellite and in situ observation collected over various time periods from 1898 to 2009. We extracted values for each cell using the cell's centroid and the data available prior to the field sampling of coral cover. Calcite was estimated from backscattering light reflectance bands near 443 and 555 nm supported by field measurements and dissolved oxygen from compiled field measurements and interpolation methods.
An additional source of data included values from a multivariate coral climate exposure stress index (Maina et al., 2011). The methods provided both a Climate Stress Model (CSM) including radiation, sediment and also a second option or Global Stress model (GSM) that included reinforcing and reducing variables, such as eutrophication. The CSM and GSM were compared against alternative models that included the above thermal history variables but did not improve the deviance explained and were therefore excluded from the final model.
Variable selection procedures were applied to the above variables that included evaluating the correlations among variables (Supplementary Figure 1). Highly correlated variables (Pearson r > 0.70) were extracted for further investigation and excluded from the final model based on previous evaluations of potential causation. For example, the variables of SST-SD, bimodality, and kurtosis were highly correlated and, based on previous variable-selected evaluations and relationships with coral community metrics, kurtosis was selected (Ateweberhan et al., 2018;Zinke et al., 2018;McClanahan et al., 2020a,b). Consequently, we considered fifteen environmental and one fisheries management variable but based on the above selection and preliminary models using AIC criteria, a final model of six environmental and one fisheries management variables is presented.

Modelling Coral Cover Benchmarks
We used a generalized additive model (GAM) as a specific case of a VM, to model the in-situ percent coral cover against the above predictors using a beta distribution and logit link. A beta distribution is appropriate for variables bound between 0-1, thus we divided our percent coral cover by 100 prior to analysis. GAMs are useful for modelling non-linear responses between responses and predictor variables. Therefore, we included smoothing terms to our predictor variables to capture the expected non-linear relationships but limited the knots for the smoothed predictors to five to prevent overfitting. Country was included as a random effect capture to address potential bias in survey effort.
Selection of predictor variables was performed to choose the most parsimonious models using the Akaike Information Criteria (AIC), where models with <2AIC were further compared on basis of the variation explained by the addition of variables. We randomly selected ∼70% of the cell-pooled data (n = 223) to develop the model and reserved the remaining ∼30% (n = 95) to test the accuracy of the model predictions. This procedure was repeated 10 times, where we calculated the average Pearson correlation coefficient, r squared and the mean squared error for all 10 runs (Supplementary Table 3). The model was then used to extrapolate hard coral cover predictions in 2020 to the entire planning grid and to 2050 as described below.
Future sea surface temperature data were obtained from the CMIP5 (Taylor et al., 2012) and included simulations of relative concentration pathways (RCP) consistent with the BAU or high (RCP 8.5) and reduced carbon emission (RCP 2.6) scenarios. We extracted daily SST for the 2006 -2050 period from two global climate models available from NOAA Geophysical Fluid Dynamics Laboratory-CM3 and MIROC5. We averaged the daily SST observations for each cell location from the two CMIP5 models for each RCP scenario to reduce individual model bias. The average daily SST were used to calculate four thermal metrics (mean, skewness, kurtosis, and cumDHW) for the 2020-2050 period for low and high Relative Concentration Pathway (RCP) scenarios. We were unable to forecast the future conditions of dissolved oxygen and calcite because they are not predicted or considered in current RCP scenarios. In present-day models, cumulative DHW was the weakest among the four variables and had a p value of 0.08, or short of formal significance (Figure 3). It was, however, retained because it was used in most previous future-cast evaluations, historically found to be statistically significant, and to have similar responses pattern in some other large-sample studies McClanahan et al., 2020a). Maps of the distribution of key variables as well as a comparison of the satellite data with selected in situ temperature data are presented (Supplementary Figures 2, 3).

Fisheries Management
The impacts of fisheries management were explored by changing the management systems of each cell and estimating cover values under this new fisheries management in 2050. We chose one common optimal fisheries management system that has been proposed and tested (O'Leary et al., 2016), which is to have 30% of reef area in high compliance closures and 70% having restricted fishing. Thus, we changed the management of the cells in 2050 under both the BAU RCP8.5 and the RCP2.6 carbon emissions reduction scenarios to estimate the coral cover under this optimal fishing scenario. The two scenarios of distributing the 30% reserve and 70% restricted fishing selection at either the regional and national level were explored. The differences were minor and therefore we present the more likely by-country scenario. We present differences in the cumulative frequencies of coral cover for the four scenarios (RCP8.5 and RCP2.6 current management and conservation management) against the 2020 benchmark for each SST cluster.

Predictions and Comparisons of Future Coral Cover Between Models
The predictions and differences in prediction of the GAM's multivariate model and the DHW-mortality model were compared. The DHW-mortality decision is a common decision of TMs used to make future predictions for RCP models (Cornwall et al., 2021). We examined the RCP8.5 and 2.6 DHW prediction for the 2020 to 2050 period to identify the decade that had two or more thermal events greater than 4, 6, 8, and 10 DHW. The coral cover was then reduced using the largest DHW by the proportion established by the empirical DHWmorality relationship reported for the Great Barrier Reef in 2016 (Hughes et al., 2018). Therefore, coral cover was reduced by 39, 60, 67, and 90% for any repeated 4, 6, 8, and 10 DHWs increments. The generality of this response beyond the GBR and over time is unknown but these or similar post-bleaching mortality rates have been used to predict future states of reefs on a large scale (Sheppard, 2003;Cornwall et al., 2021) To simplify the presentation of results, we often compare changes of coral cover above and below 10 and 25% to describe some model outcomes, as the 10% cover has been used as a thresholds for net coral reef growth (Perry et al., 2018) and 25% to achieve maximum fish diversity (Wilson et al., 2009).

Sea-Surface Temperature Distribution Patterns
Sea-surface temperature distribution patterns in the WIO were highly variable, differentiated by most SST metrics, and indicated that sites clustered into six SST groups (Table 1 and Figure 1A). The strongest differentiation (60%) was attributable to the measures of mean, SD, kurtosis, and bimodality that may be considered chronic stress metrics that differ across the cooler and less stable southern to warmer and more stable northern reefs. Cumulative DHW was a weaker differentiating variable and stronger among reef cells with cooler water, high SD, bimodality, and thick-tailed distributions (negative kurtosis) in the central to southern part of this region dominated by the westward-flowing equatorial current ( Figure 1B). The weaker variation explained by the second axis (19%) was most influenced by SST rate of rise and skewness, or a measure of the frequency of acute stress.
Clusters 1, 2, and 4 located on the left side of the PCA diagram were distinguished by being located in the central to southern WIO and having low water temperatures and high standard deviations, negative kurtosis, and bimodality (Figure 1). Cluster 1 was the largest cluster with 4794 reef cells and associated with the dominant northern flow of the Southern Equatorial Current (SEC) system. Cluster 2 was the smallest cluster with 621 cells and associated with southwest Madagascar, which was a location with more stable chronic temperatures but higher right skewness than clusters 1 and 4 cells. Cluster 4 was also a moderately large cluster with 1,463 cells that was broadly distributed from the Mauritius to the African coastline, including the southern Seychelles. This cluster included more reefs on the windward than leeward and southern sides of Madagascar than cluster 1.
Clusters 3, 5, and 6 were located on the right sides of the PCA diagram and distinguished by warm temperature with low metrics of variation or SD, bimodality, DHW and positive kurtosis. Cluster 3 was moderate size cluster with 1,548 cells but largely located along a belt from the Chagos Islands to northern Kenya that included reefs in the northern Seychelles. Cluster 3 was a more moderate set of measurements and a transition between the southern and northern WIO. It was distinguished from other southern WIO reefs by having some negative skewness whereas positive skewness distinguished the northern clusters 5 and 6. Cluster 3 also had high cumDHW similar to clusters 2 and 3. Cluster 5 had a moderate number of 1,004 cells with restricted distributions to the northern Maldives and distinguished by warm water, the strongest right skewness, and low rate of SST rise. Cluster 6 was the second most frequent cluster with 2,248 cells but had restricted distributions to the southern Maldives and fewer cells (70) in the northern Chagos. Clusters 5 and 6 both had warm water and low chronic variation in SST but cluster 5 had a thinner-tailed distribution and more episodic warm water than cluster 6.

Model Results
The variable selection procedure for the VM chose 6 of the 14 environmental variables and fisheries restrictions. Variables not selected included the rate of rise, SD, and bimodality of SST, global and climate stress models, maximum PAR, light attenuation, and pH. The full VM (GAM) with the six selected environmental variables found weak and not statistically significant relationships for CumDHW and country. Therefore, the final predictive model included in order of importance: management, dissolved oxygen, calcite, and mean, skewness, and kurtosis of SST distributions (Figure 2 and Supplementary  Table 3). This model explained 40% of the variance, had an adjusted R 2 of 0.35, and was generally supported at this level of prediction by repeated cross validation procedures. All forms of restrictive management had a modest influence of increasing coral cover by 9.3% (SD = 0.6). The observed and predicted patterns were generally consistent across the SST clusters with most of the deviation from expected coral cover among cells occurring at high coral cover levels (i.e., > 60%) among SST clusters 1 and 4 of the southern Indian Ocean (Supplementary Figure 4).
Coral cover responses predicted by the VM suggest that cover declined linearly with an increase in SST skewness and kurtosis and also with increasing dissolved oxygen to between 4.4 and 4.6 ml/l, where it becomes less predictive until 4.9 ml/l. The relationship with calcite was hump-shaped with maximum coral coverage at −7 log mol/m 3 . The best-fit model for CumDHW and coral cover was also a hump-shaped pattern but with high variance above 30 CumDHW that was likely responsible for the lack of significance and the unpredictable nature of this variable. Mean temperature was also predictive but concave shaped with Results based on the multivariate GAM or environmental variability model ( Figure  2) and climate variable predictions in 2020 and 2050. Right side of table presents summaries of SST variables mean (+SD) based on satellite data from 1985 to 2018 from NOAA Coraltemp v3.1 SST product. Organized from left to right based on the strength in ordinating the variables among the 6 SST clusters on the first PCA axis.
the lowest coverage at ∼27 o C and increasing at lower and higher temperatures. Fisheries management was a strong effect but this was largely due to increased coral cover with any form of restriction, whether it be restricted gear use or high or low compliance closures.

Future Coral Cover Predictions
The VM was used to predict coral cover in all cells in 2020 and in 2050 for the BAU (RCP8.5), reduced emissions (RCP2.6), and optimal fisheries management scenarios (Table 1 and Figures 3, 4). Results indicate geographically patchy cover distributions that generally follow the SST cluster groups. In 2020, the predicted mean cover for the whole region was 35.2% (+9.2) in 2020 ( Figure 3A). The lowest coral cover predictions (<25%) were for northern Kenya and Maldives and south-west Madagascar or more generally in clusters 2 in the south and 5 in the north. The highest cover predictions were for clusters 1 and 6 with intermediate values in clusters 3 and 4 with the exception of most of northern Kenya. At the regional level in 2050, coral cover by the VM was predicted to decline to 22.4% (+12.3) and become more spatially variable under the BAU scenario ( Figure 3B). If the low carbon emission plan of RCP2.6 was successfully implemented, an overall higher coral cover values of 27.3% (+13.7) was predicted and some of the high cover reefs were predicted to increase further (Figures 3C, 4). In the BAU scenario, the largest declines by 2050 were predicted for cluster 3 or the Seychelles, Chagos, and northern Kenya but also the northern part of cluster 5 or the Maldives. There was also a high cover loss predicted for most of the southern reefs of the Delagoa and Bight of Sofala/Swamp coast ecoregions of Mozambique from 39.2% (+8.4) to 7.2 (+8.0) and 30.5% (+4.8) to 18.2% (+7.9) cover, respectively. Southeast Madagascar reefs were also predicted to decline considerably to ∼20% for both scenarios. In contrast, the southern region's Mascarene Islands of Mauritius and Reunion, the scenario outcomes were more variable. Reunion Island's coral cover was predicted to increase under both scenarios. Coverage in Mauritius could either decline or increase dependent on the scenario, increasing in the BAU and declining slightly in the reduced carbon emission scenario. Although both climate change scenarios reduced overall coral cover relative to 2020 regional estimates, the reduced emission scenario was predicted to increase coral cover in SST clusters 1, 2, and 3 locations but not in clusters 4, 5, and 6.

Optimal Fisheries Management Scenario
Fishing restrictions consistently increased coral cover predictions relative to the business-as-usual (BAU; RCP8.5) and reduced emissions scenarios (Figures 3, 4). Management induced increases in coral cover were most notable between the 10 and 40% cover predicted in 2020 with maximum increases of around 15% relative to the BAU and 10% by the carbon emissions scenarios. Because of this non-linear initial-condition effect, including fisheries management, had a minor effect on the predicted number of cells having coral cover above 10% relative to the BAU (∼1%) and carbon emission reduction (∼2%) scenarios (Supplementary Table 4). The effect was, however, larger for reef cells with > 25% cover, increasing by 18% for the BAU and 12% for the reduced emission scenario.

Comparing Models
Comparing the TM with the VM indicates considerable variability in regional-level coral cover predictions (Table 1 and Figures 3-5). Initial coral cover in both models used the same VM-predicted 2020 baseline but cover states in 2050 differed considerably dependent on the two climate change scenarios. Overall, the GAM model predicted higher coral cover in the region by both scenarios. For example, for the RCP8.5 BAU scenarios, the number of cells with coral cover >10% in 2050 was 83% by the GAM and 23% by the DHW-mortality model (Supplementary Table 4). Thus, the number of cells with cover <10% increased by 74% for the VM compared to a 14% for the TM.
Reducing carbon emissions in the RCP2.6 scenario increased cover of reefs and cells with >10% cover. The increase was minor (11%) by the TM but more substantial (38.2%) for the VM predictions. If optimal fisheries management were established, the increase in number of reefs predicted to have coral cover <10% was reduced to 8.7% of all WIO reefs. Furthermore, the reefs with initial 2020 coral cover >40% were predicted to have more coral cover in 2050 if both restrictive management and emission reductions were implemented.
Prediction were less positive and more divergent between models for corals if the >25% threshold for biodiversity conservation was used. In particular, there was more divergence between scenarios when initial coral cover ranged between 10 and 30%. The 2020 baseline scenario predicted that 87% of reefs had coral cover >25%. The best-cast reduced emissions and increased management of the VM 2050 scenario predicted only 66% and the worst-case TM predicted only 6% of reefs had >25% cover. In a BAU scenario, 60% of the reefs were predicted to fall below 25% or the level where biodiversity was lost. Consequently, while many reefs may maintain net growth if the VM predictions were upheld, there could still be considerable losses in biodiversity without active local and global management. Where there was overlap between the TM and VM model, it was largely in clusters 4 and 5 in the reefs located from northeastern Madagascar to the Tanzanian-Mozambique border region ( Figure 5A).
The model and 2050 scenario predictions varied considerably within the SST clusters (Table 1; Figure 5 and Supplementary Table 4). The TM predictions were that the percentage of reefs with >10% cover was close to zero in clusters 2, 3, and 5 and only clusters 1 and 4 predicted to have ∼42% of the reefs above this net reef growth threshold. Cover was predicted to decline further by the TM when evaluating reefs by the >25% biodiversity threshold with the persistent clusters 1 and 4 having ∼11% of their reefs above 25% and none in all other clusters.
The 2050 predictions for VM RCP8.5 were considerably more optimistic for corals than TM with >85% of the reefs in clusters 1, 2, 4, and 6 having cover >10%. Clusters 3 and 5 were predicted to be the two most affected clusters with only 18 and 65% of reef cells above 10% cover, respectively. The frequency of reefs having coral cover >25% dropped considerably by the VM for this biodiversity threshold. Clusters 1, 2, 4, and 6 were predicted to have between 35 and 79% while predictions for clusters 3 and 5 were 12% and 25% above 25% cover. Fisheries management increased the frequencies above these thresholds, but mostly for the >25% cover threshold. For example, management in the largest clusters 1, 4, and 6 increased the frequency of reefs above this threshold by 25, 22, and 14%, respectively, compared to without management restrictions ( Figure 5A). Reducing carbon emissions in the RCP2.6 scenario further increased the predicted coral cover, such that only the most disturbed clusters 3 and 5 had significant numbers of reefs with <10% coral cover in 2050.

DISCUSSION
In the WIO region, the empirically-based VM produced a number of outcomes that differed markedly from TM predictions. First, SST heterogeneity in the region was strong and exhibited variable chronic and acute forms of stress that could often be poorly captured solely by a TM mean summer temperature threshold benchmarks. Mean optimal growth and average variances of SSTs, largely the core metrics of the TM, were not significant predictors of coral cover in the WIO, whereas skewness and kurtosis were selected and influential in the final VM. Most of these variables differed among the six large SST clusters dependent on a variety of oceanographic and geographic patterns of current and thermal exposure. Furthermore, there were a number of infrequently considered variables that were shown to be influential. Most notably, the top variables of dissolved oxygen and calcite concentrations, which should be included in future CMIP evaluates to improve future predictions for coral reefs. The importance of these variables was not peculiar to the WIO. For example, a large modelling study of coral cover in 24 ecoregions of the Coral Triangle found dissolved oxygen to be the single strongest environmental predictor of cover and calcite to be 6 th among the 17 variables evaluated (Vercammen et al., 2019). Thus, these variables exhibit similar influences in both regions and are likely to be globally influential rather spatially-restricted drivers of coral cover.
Based on the TM concept and specific application here, coral cover would be expected to decline with cumulative DHW above proposed mortality thresholds. However, large empirical cover compilations have repeatedly found weak, time-since disturbance patterns, or hump-shaped relationships (McClanahan et al., 2015b(McClanahan et al., , 2020aDarling et al., 2019). The weak hump-shaped relationship found here and in previous studies suggests opposing forces that constrain coral cover within an optimal location, which here and elsewhere was reported at around 30-35 DHW. Consequently, better understanding of thermal stress and its local impacts should illuminate why DHW, as currently used, is infrequently a strong long-term predictor of present and future coral cover (McClanahan et al., 2019, 2020a,b).
The proposed +1 • C coral bleaching threshold can appear constant at large scales (Liu et al., 2014). However, this pattern frequently emerges because the +1 • C convention was equivalent to 1.73-2.94 SDs for two-thirds of the world's coral reefs (Donner, 2011). Therefore, DHW values are frequently calculated as temperatures ∼ +2SD SST above the mean temperatures. This deviation from the mean should often but not consistently detect bleaching (van Hooidonk and Huber, 2011). However, using DHW stress thresholds in locations outside of where +1 • C is approximately equivalent to +2SD-SST results in weaker predictions of bleaching . More comprehensive evaluations of bleaching have found that removing the +1 • C threshold and adding additional variables increased the power to predict bleaching Gonzalez-Espinosa and Donner, 2021;Lachs et al., 2021). Moreover, the current DHW metric was expected to become less predictive for bleaching and other coral metrics as corals acclimate, reorganize, and adapt to changing SST conditions (McClanahan, 2017;Hughes et al., 2019). Bleaching represents a stress response to temperatures outside the normal range but is a less accurate measure of mortality (McClanahan, 2004;Welle et al., 2017) -particularly when viewed over time and the subsequent responses in coral mortality and recovery (McClanahan et al., 2015b;Darling et al., 2019). Coral sensitivity to stresses is not constant and appears to be driven by historical differences in SST variability at various temporal scales (Safaie et al., 2018;McClanahan et al., 2020a). Thus, we suggest that the local acute and chronic stress patterns modify the excess thermal stress to produce the hump-shaped DHW-coral cover relationships.
Rare and extreme excess thermal stress, or DHW, will kill corals more than moderate deviations from mean conditions. However, initial warm anomalies will also have a greater impact than subsequent disturbances, as observed in East Africa and the Great Barrier Reef in , 2017(McClanahan, 2017Hughes et al., 2019). Thus, high coral cover may occur in sites where there is a balance between the forces of acute and chronic stress and acclimation, community change, and genetic adaptation (McClanahan et al., 2020a). The lowest coral cover should occur at the two extremes or outside the optimal location between acute and chronic stress. Therefore, the predicted linear decline in coral cover by the TM was not supported here but rather suggests a weak hump-shaped DHW-coral cover relationship that we believe is more influenced by patterns of chronic and acute stress. The distribution and peak of this hump should vary with historical backgrounds and coral life histories and community structure. The considerable spatial variability in SST cluster or ecoregions suggests that coral responses to climate should be better predicted by empirical models specific to each region. An important caveat is that both models used here did not account for potential future adaptation (acclimation and community and genetic change) or differential recovery rates of corals. As empirical data accumulates specific to each region or SST cluster, the accuracy of predictions should increase beyond those provided by the current largerscale options. For example, ecoregion has been shown to be strong predictor of coral sensitivity to exposure, even stronger than dissolved oxygen when examined for the Coral Triangle (Vercammen et al., 2019;McClanahan et al., 2020b). Our study was hampered by incomplete and non-random sampling biases at the ecoregion level, which prevented us from using ecoregion as a predictive variable.
Two other key variables of TMs are the mean summer temperature and rate of SST rise, with pH also being a core part of future predictions (Camp et al., 2018). SST rise and pH did not predict the state of coral cover in our empirical study but pH was a significant variable in the Coral Triangle model (Vercammen et al., 2019). Additionally, while mean SST was the 3rd strongest predictor, the relationship found here was u-shaped rather than a decline in cover as mean temperatures increased. That is, the lowest cover was found near moderate tropical SSTs (i.e., 25-27 o C) associated with the maximum coral growth (Lough and Cantin, 2014). The implication being that in regions with rapidly changing SST and environmental conditions, locations with optimal mean temperatures failed to maintain coral cover FIGURE 5 | (A) Map plot of three future scenarios and management priority options, which include reef cells selected for (1) coral cover >25% in 2050 for reefs where there is<10% differences between the Threshold (TM) and Variability models (VM), (2) where the business-as-usual (BAU) RCP8.5 predicts >25% coral cover by the VM, and (3)>25% cover with reductions in carbon emissions and optimal fisheries management for the VM. Scatterplot predictions for coral cover comparing the VM and TM predictions in all 11,678 cells in the WIO in 2050 for (B) RCP8.5 BAU model in 2050 and (C) RCP2.6 carbon reduction scenario. The specific TM presented here is from the DHW-mortality relationships established from observations of DHW and mortality in the GBR during 2016 (Hughes et al., 2018). Cornwall et al. (2021) used this algorithm to estimate coral cover and net reef calcification in 2050 (see methods).
during this period episodic thermal stress These faunal provinces have been undergoing rapid change in recent years, which may be detrimental to corals with narrow optimal growth conditions (McClanahan et al., 2014;Abram et al., 2020).
Less commonly examined metrics not included in TMs were found to be moderate predictors. For example, SST skewness and kurtosis, which were used here as proxies for chronic and acute stress were the 3rd and 4th strongest variables. When evaluated elsewhere, these variables were often among the strongest predictors of a number of coral metrics (Ateweberhan et al., 2018;Zinke et al., 2018;McClanahan et al., 2020a). Cover declined as SST distribution tails became thinner and right skewed, as would be predicted where there is poor acclimation to moderate chronic stress. In some cases, other related proxies, such as daily, seasonal, and annual temperature ranges were good predictors of bleaching and coral cover (Safaie et al., 2018;Vercammen et al., 2019). The maximum and range of SSTs are among the Bio-Oracle SST data options and found to be the 2nd and 3rd strongest environmental variables in predicting coral cover in the Coral Triangle (Vercammen et al., 2019). Consequently, these alternative proxies support the general importance of acute and chronic stress on very broad scales. Including kurtosis and skewness or other appropriate proxies of chronic and acute stress SST metrics, including dissolved oxygen and calcite, should therefore improve predictions Data distribution variables were shown here to be important but also sensitive to sampling frequency. Therefore, satellite and in situ measurements were likely to differ dependent on the frequency of sampling (McClanahan, 2020; Supplementary  Figure 3). Less frequently sampled and accurate satellite data may not be able to accurately detect chronic and acute stress (McClanahan and Muthiga, 2021). Satellite data has the distinct advantages of wide coverage and usage for RCP-based predictions, which suggests a need to further consider the satellite proxies that are best at measuring chronic and acute stress. Satellite coverage and CMIP usage allowed us to evaluate a reduced-variable model to make predictions aligned with the metrics currently available. This procedure was useful but less convincing for producing the most accurate predictions expected with the future states of all the variables selected by the GAM. Enough is known about some of these variables that including them in future CMIP models and scenarios is possible and their inclusion should be pursued to improve predictions.
Contextual variables of dissolved oxygen and calcite need to be included in future-cast models to improve predictions. Low dissolved oxygen has been considered a potential stress for corals and other reef invertebrates (Hughes et al., 2020). However, the increasing coral cover with declining or fluctuating oxygen found here suggests a more complicated relationship. The high tolerance to low oxygen among some coral species may explain this pattern. Nevertheless, much remains to be known about how corals respond to the direct and indirect effects of oxygen and how they are associated with other oceanographic and physio-chemical drivers of coral acclimation (Camp et al., 2016). The findings presented here and elsewhere should generate hypotheses for future research (Vercammen et al., 2019). An example of potentially complex relationship is the finding that some coral genes provide co-tolerance to both low oxygen and heat (Alderdice et al., 2021). These investigator's experiments suggested that the upregulation of a hypoxia-inducing factor influenced a coral's tolerance to thermal stress. Clearly, addressing the lack of knowledge, covariance, and the possibility of co-tolerance among stressors should improve understanding and predictions. An unexpected future scenario possibility is that the projected declines in dissolved oxygen with ocean warming could increase tolerance of corals to thermal stress and produce greater resistance than expected to climate change.
The role of fisheries management on coral cover is complex when viewed on a global scale Selig et al., 2012;Bates et al., 2019). Nevertheless, previous large compilations of coral cover in the WIO region have shown a weak positive effect of fisheries management on corals but influenced by a number of environmental, ecological, and human use contexts McClanahan et al., 2014;Graham et al., 2020). For example, low compliance reserves and restricted fishing had similar levels of coral cover in the VM. Both fisheries categories have been shown to have similar fish biomass (McClanahan et al., 2015a). Therefore, restrictions on the use of the most destructive gear, such as drag-nets, smallmeshed nets, explosives, and spearguns may be important for maintaining coral cover. Cover is likely to be more resistant to modest human fishing disturbances and therefore not requiring the total elimination of local human impacts. If these restriction patterns remain constant or increase further, there could be potential to maintain or increase future coral, particularly in reef already having moderate to high cover.
While our model had modest predictive power, there is still considerable unexplained variation similar to the Coral Triangle model (Vercammen et al., 2019). Unexplained variation based on recent field data is also expected to increase in the future as corals adapt and other unexamined variables influence outcomes of thermal disturbances (McManus et al., 2021). Missing variables is a concern as shown here but also future changes in the disturbances themselves. Changes include a large variety of possible factors that arise for a number of reasons. Not the least is that coral cover is a composite metric that does not include all of the coral organisms and their various life histories, inter-species relationships, and variable responses to stresses (Thompson et al., 2015;Darling et al., 2019;McClanahan et al., 2020a). While cover is the core metric for evaluating reef status, it lacks a considerable amount of life history and stress-response information. Second, is that the grid scale lacks the ability to predict more localized habitat variability. Moreover, the current environmental data are compiled at a coarser scale of ∼9.2-km and therefore limit the potential to evaluate finer-scale resolutions (Assis et al., 2018), Heterogeneity in smaller-scale stress in reef environments, as well as local organismic sensitivity to stress, is evident. Therefore, any grossscale predictions will fail to capture more resolved variables including habitats, ocean exposure, and depth, among others. Finally, the VM used here reminds us that environmental science can often focus on variables that are assumed to be drivers while failing to examine less considered but potentially important variables. This is evident in the historical focus on water quality parameters that have produced complex, contradictory, and compensatory relationships on coral stress and recovery patterns (MacNeil et al., 2019;Lesser, 2021). Here, eutrophication variables were not examined because they are not widely available on the scale of our study. Nevertheless, future work should consider eutrophication or other proxies to compare their influence relative to calcite and dissolved oxygen. Directed by an incomplete theory, lack of observations, and unobtainable metrics, there is the high possibility for blind spots and unreliable predictions. We see here that the addition of a small number of available variables, such as calcite and dissolved oxygen, had large effects on predictions.
Spatial heterogeneity as uncovered in the SST clustering approach and the consequent future predictions indicate the complexity of potential responses. One notable case of this is the large difference in predictions for clusters 5 and 6 along the Maldives-Chagos islands. Despite their shared oceanic and northern Indian Ocean locations, future predictions were quite sensitive to the patterns of SST skewness along this island chain. Will this observed difference drive changes or will excluded variables of dissolved oxygen, calcite, and others modify outcomes? A large-scale study in 2011 did find lower coral cover in the north compared to the south (Tkachenko, 2015), which is one of the GAM predictions for clusters 5 and 6. These questions are currently unanswerable but could be addressed by further research and monitoring.
This study compared two gross classes of models, the variability and threshold models. It showed how the inclusion of novel environmental variability and empirical relationships can change predictions. The VM predicts less and more spatially variable thermal sensitivity than the TM, which contains fewer influential variables and therefore larger-scale homogeneity in responses to thermal exposure. Predictions overly-reliant on a TM variations and low exposure metrics may have limited spatial applications including identifying the location of sanctuaries. Finding sanctuaries will require expanding and diversifying the suite of ecosystem variables examined for their role in avoidance, resistance, and recovery from climate stresses. The inclusion of five poorly examined variables in the VM produced a more optimistic view of the future than the TM. The heterogeneity of the VM provokes the question of what might be the effect on future predictions for coral reefs by expanding the number of variables and increasing the spatial resolution of models? Moreover, the efficacy of increasing the quantity, quality, and distribution of field data used to calibrate models. Even largescale changes, such as global emissions that are expected improve conditions for corals in most regions, were shown to have quite variable consequences for corals among the various SST clusters and management systems. The current failure to understand and evaluate the many sources of variability suggest a need to adopt risk-spreading policies and management decisions.

DATA AVAILABILITY STATEMENT
Environmental data is publicly available and can be found at: https://www.bio.oracle.org/. Coral cover data are available on a formal request to the authors.

AUTHOR CONTRIBUTIONS
TM: conceived the study, collected the field data, supervised the analysis, interpreted the data, and wrote the manuscript. MA: organized the environmental data, undertook the statistical analyses, wrote the methods section, and edited the manuscript. Both authors contributed to the article and approved the submitted version.