High-resolution modeling of thermal thresholds and environmental influences on coral bleaching for local and regional reef management

Coral reefs are one of the world’s most threatened ecosystems, with global and local stressors contributing to their decline. Excessive sea-surface temperatures (SSTs) can cause coral bleaching, resulting in coral death and decreases in coral cover. A SST threshold of 1 °C over the climatological maximum is widely used to predict coral bleaching. In this study, we refined thermal indices predicting coral bleaching at high-spatial resolution (1 km) by statistically optimizing thermal thresholds, as well as considering other environmental influences on bleaching such as ultraviolet (UV) radiation, water turbidity, and cooling effects. We used a coral bleaching dataset derived from the web-based monitoring system Sango Map Project, at scales appropriate for the local and regional conservation of Japanese coral reefs. We recorded coral bleaching events in the years 2004–2016 in Japan. We revealed the influence of multiple factors on the ability to predict coral bleaching, including selection of thermal indices, statistical optimization of thermal thresholds, quantification of multiple environmental influences, and use of multiple modeling methods (generalized linear models and random forests). After optimization, differences in predictive ability among thermal indices were negligible. Thermal index, UV radiation, water turbidity, and cooling effects were important predictors of the occurrence of coral bleaching. Predictions based on the best model revealed that coral reefs in Japan have experienced recent and widespread bleaching. A practical method to reduce bleaching frequency by screening UV radiation was also demonstrated in this paper.


INTRODUCTION
Biological communities can shift toward alternative stable states in response to changing climate (Parmesan & Yohe, 2003). Coral reefs are one of the most susceptible ecosystems to global warming and local environmental stressors (Hoegh-Guldberg, 1999;. Rising sea-surface temperatures (SST) can cause bleaching in reef-building corals, especially during summer (Hoegh-Guldberg, 1999;Brown et al., 2002;.
Excessive thermal stress leads to expulsion, digestion, or reduced pigmentation of symbiotic dinoflagellate algae in coral cells, resulting in the exposure of white coral skeletons (i.e., bleaching;Hoegh-Guldberg, 1999;Brown et al., 2002). Prolonged warming trends in sea temperature have been predicted to increase the frequency and severity of bleaching in the future, leading to mass mortality of corals (Hoegh-Guldberg, 1999;Donner et al., 2005;Donner, 2009;. Reef management relies on not only global measures to reduce climate warming but also local measures to control environmental influences on coral resilience . Spatial and temporal predictions of coral bleaching under varying environmental conditions could therefore provide valuable information to support local management of coral reefs. The degree heating week (DHW) index of cumulative thermal stress, developed by the National Oceanic and Atmospheric Administration Coral Reef Watch (NOAA CRW), has been widely used to predict coral bleaching. DHW is based on SST derived from satellite images, and is computed as the sum over a period of 12 weeks of thermal stress exceeding 1 C above historical summer monthly SST (Liu, Strong & Skirving, 2003). DHW over 4 C weeks indicate severe coral bleaching and constitute a bleaching alert threshold (Liu, Strong & Skirving, 2003).
Despite its increasing use globally, the predictive performance of DHW may not be sufficient for local reef management, as DHW on average detects only 40% of global coral bleaching events (Donner, 2011). This low predictive performance may be due to the use of a fixed thermal threshold of 1 C above baseline SST. Previous studies have suggested that thermal stress of 1 C or below can induce coral bleaching (Brown, 1997;McWilliams et al., 2005;Kleypas, Danabasoglu & Lough, 2008). In addition, historical temperature variability can affect bleaching and coral resilience . As a consequence, some studies have used modified indices, such as the sum of thermal stress over 0 C above baseline SST Kayanne, 2017). Donner (2011) proposed two-modified DHW indices: an index using historical SST variability as the bleaching alert threshold, and an index using the mean of the warmest monthly SST of each year as the baseline SST.
To evaluate the effects of global and local stressors on corals, a high-performance predictive model operating at high-spatial resolution is required. Global stressors such as thermal stress can vary at a local scale (Strong et al., 2002;Liu et al., 2014). Furthermore, there are potentially interacting environmental stressors such as ultraviolet (UV) radiation (Hoegh-Guldberg, 1999; and variables such as water turbidity Oxenford & Vallés, 2016), topography of the sea floor Oliver, Berkelmans & Eakin, 2009), and exposure to winds  and currents (Nakamura & van Woesik, 2001;) that can affect coral bleaching. For example, increasing the speed of surface currents and winds can reduce bleaching risk by increasing mixing in surface seawater (Nakamura & van Woesik, 2001;. Modeling coral bleaching at a local scale also requires high-resolution observational records, as omission of bleaching events can lead to poor predictive power in models (Oliver, Berkelmans & Eakin, 2009). ReefBase (Tupper et al., 2011) and the Bleaching Database V1.0 (Donner, Rickbeil & Heron, 2017) provide high record coverage in some areas, including in the Great Barrier Reef and the Caribbean (van Hooidonk & Huber, 2009). However, records are still limited for other areas, such as the Pacific islands (Donner, Rickbeil & Heron, 2017). One possible reason for this data gap is language barrier. A considerable amount of data in ReefBase (Tupper et al., 2011) have been provided by nonprofessional (citizen) specialists who are not native English speakers. Few Japanese records (N 64) are found in the global databases, despite the large amounts of research conducted on coral reefs in Japan. To collect and collate observational records of corals throughout Japan, diverse Japanese stakeholders, including professional scientists, government officials, and citizens, constructed a web-based monitoring system for Japanese coral reefs in 2008, the Sango Map Project (Namizaki et al., 2013). Collecting observational records in a web-based database proved to be effective in Japan, as internet service is available to the vast majority of the population. In addition, the use of Japanese language allowed a larger number of stakeholders to contribute, including stakeholders from populated islands where diving services are available. This project contributed key data to the International Year of the Reef Year in Review report (Staub & Chhay, 2009).
In this study, we aimed to improve predictive power in models of coral bleaching at high-spatial resolution, in order to inform local and regional reef management. We used observational records of coral bleaching derived from the Sango Map Project, and we compared the predictive performance of multiple thermal indices and their modifications in models with multiple explanatory variables. We developed a novel derivation of DHW (hereafter "filtering threshold") to compute thermal stress below 1 C in excess of the baseline SST, using historical SST variability as a threshold. We optimized the filtering threshold by statistical estimation of each type of DHW and degree heating month (DHM) index. To maximize predictive performance, we then optimized the combination of multiple explanatory variables while optimizing the filtering threshold. Based on the model with maximum predictive performance, we produced spatial predictions of coral bleaching in the study area, as well as predictions under reduced local environmental stresses. Our results provide a reference for local reef management in Japan, although our methods could be applied for local reef management in other areas.

Observational records of coral bleaching
We used observations from the Japanese coasts submitted to the Sango Map Project (http://www.sangomap.jp/) up to March 2017. Observations were composed of the following information: (1) presence or absence of corals; (2) longitude and latitude of the location, searchable through the Google Maps API (https://developers.google.com/maps/); (3) name of the location; (4) date, month, and year of the observation; (5) method of survey (scuba diving, snorkeling, glass boat, walking, or other); (6) water depth in meters; (7) observer's professional background (professional scientist, nonprofit or nongovernmental organization, tourism, or other); (8) level of severity of coral bleaching (high, medium, low, nonbleaching, or not available) derived from the bleaching dataset in ReefBase. We confirmed or rejected questionable records, such as observations made on land or in the open ocean and observations of doubtful coral species.
After quality control and exclusion of records lacking information on bleaching, we obtained 668 independent records between July 2004 and October 2016. Of these observations, 52 were submitted by professional scientists, 152 by nonprofit or nongovernmental organizations, and 134 by tourists. Fifty-nine observations were conducted as part of CoralWatch (http://www.coralwatch.org/) and 63 as part of ReefCheck Japan (http://www.reefcheck.jp/). The records provided good spatial coverage of coral reefs in Japan (Fig. 1). Most of the records were obtained in the first three years following the launch of the Sango Map Project, including 449 records from 2008 to 2010 alone. In addition, 82 and 111 records were reported in 2013 and 2016, respectively, when mass bleaching events were observed throughout Japan (Kayanne, 2017;Kayanne, Suzuki & Liu, 2017  Records of bleaching not induced by thermal stress were regarded as nonbleaching observations for the purpose of this study. We therefore reclassified 107 bleaching observations as nonbleaching observations (Step 1 in Table 1). Following screening, the prevalence of records was more biased than prior to screening, with 228 bleaching and 440 nonbleaching observations collated (Fig. 1). However, the risk of biased predictions was still deemed low (Step 2 in Table 1). Annual and spatial patterns of bleaching occurrences were consistent with those reported previously in Japan (Kayanne, 2017;Kayanne, Suzuki & Liu, 2017). We assessed spatial autocorrelation in the residuals of the prediction model of the NOAA CRW DHW (Step 2 in Table 1), using the spatial autocorrelation coefficient (Moran's I). We confirmed that there was no significant autocorrelation in the residuals, indicating no significant spatial bias in the data (Dormann et al., 2007).

Thermal indices
To calculate thermal indices, we used daily data at a spatial resolution of 1 km (0.01 ) from the Multi-scale Ultra-high Resolution Sea Surface Temperature (MUR SST) Analysis version 4.1 (JPL MUR MEaSUREs Project, 2015) (http://dx.doi.org/10.5067/GHGMR-4FJ04). The MUR SST product is a blend of SST from six satellites and thus provides higher accuracy than single-satellite products. MUR SST data are only available from 2002 to May 2017, an insufficient period of time to calculate maximum monthly mean (MMM) climatology (Table 2). In addition to data from the MUR SST, we also used data for the years 1985-2002 from the Optimum Interpolation SST (OI SST) version 2 (Reynolds et al., 2007) (http://www.esrl.noaa.gov/psd/data/gridded/data.noaa.oisst.v2.highres.html). To correct the SST bias between MUR SST and OI SST, we added the bias in the 2002-2017 monthly climatology to the OI SST data, after down-scaling to 0.01 using inverse distance weighting interpolation (Tabor & Williams, 2010;Yara et al., 2011).
Using the monthly mean SST from 1985 to 2015, we obtained two types of MMM climatologies. The first MMM climatology follows the NOAA CRW version 3 protocol . The temporal midpoint was recentered to that of the heritage 50 km MMM (1985)(1986)(1987)(1988)(1989)(1990)1993) using the approach of Heron et al. (2014) as follows: where SST i is the SST climatology as obtained above and SST recenteredi is the recentered SST climatology at cell i. The linear trend of monthly mean SST between the center times of the two-time durations (T 1985-2015 , T 1985-1993 ) is represented by slope i at cell i. The down-scaled and recentered MMM correlated significantly with the CRW MMM version 3 (Figs. S1 and S2).
The second climatology, known as MMM max climatology (Donner, 2009(Donner, , 2011, uses the mean of the warmest month of each year instead of the mean of the warmest month in the climatological years in the MMM climatology ( Table 2). The warmest month is not always the same among years and, therefore, MMM max is larger than MMM (Table 2; Figs. S3A-S3C) and represents the seasonal peak in SST more accurately than MMM Table 1 Flowchart summarizing the three steps in our analysis.

Procedure
Approach Reference Step 1

Control of observation errors
Excluding bleaching events not induced by thermal stress Observation records of small bleaching events (e.g., those within microatolls, or caused by disease or predation) and observations made after the small bleaching event were regarded as nonbleaching if the 1 km resolution DHW value at observation site did not exceed zero This study (2018) Step 2 Assumptions for observed data Checking equality in observations of occurrence and absence of bleaching, where higher prevalence (usually biased to occurrences) results in larger predicted probabilities (i.e., biased predictions) Using an evaluation index that is less dependent on prevalence (TSS). The evaluation threshold was also optimized (see Step 4) Allouche, Tsoar & Kadmon (2006) and Liu et al. (2005) Avoiding spatial autocorrelation in the data, which can increase false-positive predictions Evaluating the spatial autocorrelation coefficient (Moran's I) of residuals from a prediction model. If residuals are spatially biased, spatially clustered data should be filtered Dormann et al. (2007) and Boria et al. (2014) Step 3 Assumptions for environmental variables Screening correlated environmental variables If correlations between variables are high (jrj > 0.7), correlated variables should be excluded to reduce multicollinearity, which can affect both GLM and RF Dormann et al. (2013) Step 4 Evaluation and model assessment Multiple performance metrics were used to avoid Type I and Type II errors. Models using standard and optimized thresholds were assessed. A statistical model (GLM) and a machine learning model (RF) were used Optimizing the evaluation threshold Optimizing the threshold to discriminate occurrence and absence from the predicted probability of bleaching. Although statistical models predicting occurrence or absence typically output results as probabilities, using a 0.5 (i.e., midpoint) threshold can yield biased results under unequal class prevalence. To avoid this problem, TPR-TNR sum maximization was used to optimize the threshold (Table 2) Manel, Williams & Ormerod (2001) and Liu et al. (2005) Optimizing the filtering threshold To optimize DHW and DHM, the filtering threshold was adjusted by 0.01 C of precision to maximize predictive power (i.e., TSS) for each combination of explanatory variables This study (2018) Evaluation using 10-fold cross-validation A randomly selected 30% subset of the data were used as testing data, and the remaining data were used as training data. Prediction models were built with the training data and evaluated against the testing data. The test was repeated 10 times for each filtering threshold and combination of explanatory variables Hijmans et al. (2017) climatology. This method is particularly effective in tropical zones with reduced seasonality (Donner, 2011).
We calculated eight types of thermal and cooling indices for each grid cell and observation day, including mean weekly and monthly SST, DHM, DHW (MMM + a C), DHW (MMM + a C) using SST variation as the bleaching alert threshold, DHW (MMM max + a C), DHW (MMM + bs m C), and degree cooling weeks (DCW) (see Table 2 for a detailed derivation of the indices). Historical SST variability (s m ) ( Table 2) was calculated with the monthly mean SST from 1985 to 2015, and ranged from 0.36 to 0.71 with a median of 0.57 (Fig. S2D). Although DCW is calculated with a similar algorithm to that used for DHW, DCW was not significantly correlated with DHW. We therefore included DCW as a covariate in our models. The filtering thresholds (a and b) were fixed to 1 in the standard thermal indices and optimized in our indices.

Additional environmental variables
Monthly UV-B and PAR data were obtained from the Japan Aerospace eXploration Agency Satellite Monitoring for Environmental Studies (JASMES; http://kuroshio.eorc. jaxa.jp/JASMES/; accessed 25 June 2017) and derived from the average of data extracted from the Aqua and Terra sensors of moderate resolution imaging spectroradiometer (MODIS; http://modis.gsfc.nasa.gov/data/dataprod/). Although both UV-B and PAR may affect coral bleaching (Hoegh-Guldberg, 1999), the variables were significantly correlated (r = 0.79) . We excluded PAR from our analysis (Step 3 in Table 1) as parameters may be misestimated in statistical modelings and machine learnings under multicollinearity (Dormann et al., 2013).
To quantify the speed of surface currents, we extracted data from the HYCOM+ NCODA Global 1/12 Analysis GLBu0.08 from 1997 to 2017 (https://hycom.org/ dataserver/gofs-3pt0/analysis/; accessed 22 June 2017). We obtained climatological median from July to September, the months during which most of the recorded bleaching events occurred. To quantify wind speed, typhoon tracking data were obtained from the Regional Specialized Meteorological Center Tokyo (http://www.jma.go.jp/jma/jma-eng/ jma-center/rsmc-hp-pub-eg/trackarchives.html; accessed 22 June 2017). We calculated the wind speed index for each grid cell as the length of time without typhoons, with

Procedure
Approach Reference Step 5 Coral bleaching prediction Prediction under observed environmental conditions Using the best performing model built in each crossvalidation, the probability of coral bleaching was predicted for the study area Hijmans et al. (2017) Prediction under reduced UV radiation due to screening effect Coral bleaching frequency may be reduced by a 40% reduction in UV radiation and a 40% increase in water turbidity due to screening
Step 4: evaluation of predictive models.
Step 5: predictions of coral bleaching. DCW, degree cooling week; DHM, degree heating month; DHW, degree heating week; RF, random forest; TNR, true negative rate; TPR, true positive rate; TSS, true skill statistic; UV, ultraviolet.  typhoons defined as wind speeds over 15 ms -1 . However, the wind speed index was strongly correlated with DHW (r = 0.86) and therefore was excluded from the analysis. We used the diffuse attenuation coefficient (K 490 ) as an index of water turbidity, which can reduce light radiation stress involved in bleaching (Table 2). A monthly composite of K 490 (4 km, Level-3 binned MODIS AQUA products) was obtained from the NOAA OceanColor database (https://oceancolor.gsfc.nasa.gov; accessed 8 September 2017) for the months July-September.
Data on current speed and diffuse attenuation were down-scaled to 1 km using bilinear interpolation. When environmental variables were not available for coastal cells, we used inverse distance weighting interpolation to estimate coastal values.

Model evaluation and optimization
We evaluated coral bleaching models based on the accuracy of both positive (bleaching) and negative (nonbleaching) predictions. Most studies have evaluated models of coral bleaching based only on overall accuracy, such as the proportion of correct predictions and AIC (Maina et al., , 2011Kayanne, 2017;Welle et al., 2017), whereas a few studies have differentiated the accuracy of positive and negative predictions van Hooidonk & Huber, 2009;Donner, 2011). When the number of bleaching and nonbleaching observations is unequal (i.e., under class imbalance), model predictions can be biased. We therefore used four evaluation metrics: overall accuracy, true positive rate (TPR), true negative rate (TNR), and true skill statistic (TSS) (defined in Table 2). TSS quantifies prediction skill as the index weighs positive and negative predictions equally (Allouche, Tsoar & Kadmon, 2006).
To assess the combined effects of thermal stress and multiple environmental influences on coral bleaching, we constructed prediction models of bleaching with two approaches: generalized linear model (GLM) with a binomial error distribution and a logit link function, and random forests (RFs; Breiman, 2001). Although both models compute predictions in the form of probabilities, the models are based on different algorithms. GLM is an extension of regression models, whereas RF is a machine learning method that uses randomly repeating classifications to capture complex interactions among explanatory variables (Table 2). Therefore, GLM has an advantage of which fitted model can be written as a formula that is easy to be used subsequently.
We confirmed that the data met the assumption of binomial GLM (i.e., the residual deviance per degree of freedom was less than 1.5; Zuur et al., 2009). GLM was applied with the "glm" function in base R (R Core Team, 2017), and RF was applied with the "randomForest" function of the randomForest R package. RF was used under standard settings to avoid overfitting the training data. However, we followed the recommendation of Hijmans et al. (2017), by specifying the model as a "regression model," even though the response variable was categorical. The relative importance of explanatory variables was calculated with the "importance" function of the MuMIn package (Barto n, 2015) for GLM and of the randomForest package for RF (Liaw & Wiener, 2002).
Predicted probabilities were transformed into bleaching and nonbleaching categories with the threshold that maximized the sum of TPR and TNR (Liu et al., 2005).
The midpoint (i.e., 0.5) is often used as a threshold (Fig. S1C), although it is sensitive to class imbalance in the training data and may, therefore, lead to inaccurate predictions (Liu et al., 2005). This issue is addressed in studies of species distribution modeling, although it remains poorly addressed in studies of mass bleaching (van Hooidonk & Huber, 2009). We used the "evaluation" function of the dismo R package (Hijmans et al., 2017) for evaluation and optimization of the threshold.
Models were evaluated by 10-fold cross-validations using TSS as the evaluation index. In each repeat, we separated 30% of the data as testing data and used the remaining 70% for constructing GLM and RF (Step 4 in Table 1). We optimized the filtering thresholds (Fig. S1A) for DHM and DHWs by cross-validation, while the filtering thresholds were fixed at 1.0 C in the standard indices (Step 4 in Table 1). We selected the optimum filtering threshold between 0 and 1.5 C for indices using a constant threshold (a), whereas we examined the coefficient of s m (b) between 0.1 and 2.5 for indices based on historical variability (Table 1). We conducted optimizations with 0.01 precision for both types of indices, i.e., with 151 and 241 submodels, respectively.
For models with multiple explanatory variables, we considered DCW, historical SST variability, UV-B, water turbidity, water depth, and current speed, in addition to the thermal index. The two most influential variables (DCW and UV-B) were always included in models with multiple explanatory variables (Step 4 in Table 1). The optimum set of explanatory variables was specified through cross-validation. The set of variables that best explain variation in the testing data was selected among all 15 possible combinations. In total, we evaluated 22,650 and 36,150 models (10 cross-validations Â 15 variable combinations Â 151 or 241 submodels) for each GLM and RF model, respectively.
Finally, we predicted coral bleaching in the warmest month of the main coral habitat in the study area using the best predictive model (Step 5 in Table 1). We also assessed reduction in UV-B as a possible adaptive measure by reducing UV-B radiation by 40% and increasing water turbidity by 40%, thereby simulating the effects of screening with fishnets (Step 5 in Table 1). These levels of changes were consistent with in situ examination in Onna Village in the Ryukyu Islands (Okinawa Prefecture, 2017). Predictions were obtained from each model built in the 10 cross-validations, and subsequently averaged among the models. Spatial data were obtained from the Global Map Japan version 2.1 Vector data, provided by the Geospatial Information Authority of Japan (2015).
All analytical codes (available in Supplemental Information 1) were written in R version 3.4.1 (R Core Team, 2017).

Effects of environmental variables
Predicted probability of bleaching increased with increasing values of thermal indices, including SST, DHM, DHW, and UV-B (Fig. 2). Predicted probability of bleaching decreased with DCW, water turbidity, and water depth (Fig. 2). Relationships between bleaching and historical SST variability and current speed were not significant, with 95% confidence intervals (CIs) ranging from negative to positive. Relationships between bleaching and monthly and weekly SST were positive, although the widths of the 95% CIs suggest these variables are not reliable indices of coral bleaching. Alert thresholds for predicted bleaching were found to be lower than standard thresholds, except for DHM (Table 3). Values of 1 and 0 represent bleaching and nonbleaching, respectively. Solid lines and gray areas indicate mean model fit and 95% confidence intervals, respectively. Dotted lines represent thresholds discriminating bleaching and nonbleaching, which were optimized by true positive rate and true negative rate (TPR-TNR) sum maximization (Table 2). See Table 2  Ranking of important variables was similar between GLM and RF (Fig. 3): the best explanatory variable was DHW (100% in both of GLM and RF), followed by DCW. UV-B, water turbidity, and historical SST variability also explained substantial variation in coral bleaching. The explanatory powers of historical SST variability and current speed were high, despite inconsistent relationships with coral bleaching (Fig. 2). Absolute variable importance differed between GLM and RF. In GLM, most variables explained more than

A B
Relative variable importance (%)

Optimization and assessment of filtering thresholds
Optimization of filtering thresholds improved the predictive performance of DHM and DHWs, although the improvement was small for GLM (Fig. 4). Improvement by the optimization was around 0.01 in TSS in GLM, while the improvement of DHW using the historical SST variability (s m ) as the filtering threshold was 0.02-0.03 in TSS in RF. We compared the predictive performance of all models including standard and optimized thermal indices, and the optimized set of explanatory variables (Fig. 5).  Table 2 for terminology and Tables 3-6  In models based only on a thermal index with the standard threshold, TNR was larger than TPR, indicating that high overall accuracy can result from effective identification of nonbleaching occurrences, despite ineffective identification of bleaching occurrences. TSS was indicative of both TPR and TNR. Among thermal indices with the standard threshold, DHW using historical SST variability (s m ) as the bleaching alert threshold and DHW using s m as the filtering threshold scored the best performances in TSS (0.60 and 0.65, respectively) (Fig. 5). Weekly SST had no prediction skill (TSS = 0) and the skill of monthly SST was low (0.20). TPR of DHW using historical SST variability as the bleaching alert threshold was the highest (0.85) among models, although the false-positive rate (1-TNR = 0.20) was also the highest. Using historical SST variability as the filtering threshold for DHW decreased the false-positive rate (0.09), but also decreased TPR (0.69). DHW using MMM max as the baseline climatology showed the lowest predictive performance (TSS = 0.37) among all DHW indices, with the lowest TPR (0.38). Predictive performance was improved by optimizing the evaluation threshold ( Fig. 5; Table 4). Under optimized evaluation thresholds, DHW with historical SST variability as the filtering threshold performed best (TSS = 0.72 for GLM) among the three types of DHW, although differences were small. Predictive performance of DHM was the lowest (TSS = 0.40). When both evaluation and filtering thresholds were optimized (Table 5), improvements in evaluation thresholds were only small. DHW with historical SST variability as the filtering threshold remained the highest performing index under optimization (TSS = 0.73 in both GLM and RF). Optimized filtering thresholds were less than 1 C in all the models using one thermal index, using the optimized DHM or DHW other than DHW using historical SST variability (s m ) as the filtering threshold (see also Fig. S3D for the distribution of s m ).
The predictive skill of models with multiple explanatory variables was negligibly higher than that of models including only one thermal index. TSS increased by 0.1 at most in all Table 4 Univariate prediction models of coral bleaching using thermal indices with optimized evaluation and filtering thresholds.

Model
Evaluation threshold (Bleaching alert threshold C) In models with multiple explanatory variables, differences in predictive performance between models with optimized evaluation thresholds and models with optimized evaluation and filtering thresholds were smaller than differences between GLM and RF models for most thermal indices. RF always performed better than GLM, with differences in TSS of 0.04 to 0.05. Although the TPR of GLM exceeded that of RF in most cases, the TNR of GLM was lower than those of RF, i.e., the risk of false positives was higher in GLM. The RF model based on DHW with MMM + 0.97 C filtering threshold, DCW, UV-B, water turbidity, historical SST variability, and current speed (Table 6) showed the best predictive performance (TSS = 0.79; TPR = 0.90; TNR = 0.89). Among the GLM, the model consisting of DHW with MMM + 1.83·s m C filtering threshold, DCW, UV-B, and turbidity showed the best predictive skill.  Table 2). The optimized evaluation thresholds (mean ± SE) of the predicted probability of coral bleaching are shown with corresponding bleaching alert thresholds of thermal indices. The optimized formula for predicted probability of bleaching is shown for GLM. logistic(x) = 1/(1 + exp(-x)). SST, sea-surface temperature; DHM, degree heating month; DHW, degree heating week; MMM, maximum of the monthly mean SST climatology; MMM max , mean of the warmest monthly mean SST of each year; GLM, generalized linear model; RF, random forest.  Table 2). The optimized evaluation thresholds (mean ± SE) of the predicted probability of coral bleaching are shown with corresponding bleaching alert thresholds of thermal indices. The optimized formula for predicted probability of bleaching is shown for GLM. logistic(x) = 1/(1 + exp (-x)

Predictions of coral bleaching
We predicted probabilities of coral bleaching in the main coral-habitable areas of Japan with the optimized best multivariate model of RF using DHW with MMM + 0.97 C filtering threshold (Fig. 6). The mean predicted probability of bleaching ranged from 0.46 to 0.74 among areas. Spatial variation in the probability of bleaching was found in both the eastern (Fig. 6A) and the western (Fig. 6C) Ryukyu Islands. Hotspots with higher bleaching probabilities were found in the southeastern part of Okinawa Island, the eastern part of the Kerama Islands (Fig. 6A), and the northern part of Ishigaki Island (Fig. 6C). This resulted in bleaching in 2.5-5 (44-100%) of five years (2008-2010Figs. 6B and 6D  bleaching of up to 0.24 were observed, resulting in a significant decrease in bleaching frequency of up to 56% (Figs. 7B and 7D). Bleaching in fewer than three out of five years occurred in most areas.

DISCUSSION
Optimizing coral bleaching models The ability to predict coral bleaching was improved by optimizing thermal indices, particularly SST and DHWs. As in a previous study (Donner, 2011), we generally found lower TPR than TNR, but both TPR and TNR were improved by optimization. The sole contribution of optimizing the filtering threshold was small. However, optimizing the filtering threshold and the evaluation threshold while combining multiple environmental variables achieved large improvements in TPR and TNR (reaching ∼0.9). We also found that cooling (DCW), UV-B, and screening (water turbidity) were important predictors of bleaching, particularly in RF models.
Our results are mostly consistent with those of Donner (2011), who showed that DHW using historical SST variability as bleaching alert threshold had a higher TPR but a higher false-positive rate than NOAA CRW DHW, and that DHW using MMM max did not predict bleaching accurately but was suitable in equatorial zones. Because our study was conducted at a higher latitude, using MMM max resulted in lowest performance among the DHW indices, as expected. We also used historical SST variability as a filtering threshold, and this method showed the highest prediction skill among the models with only one thermal index. Historical SST variation (s m ) may be a particularly effective predictor of bleaching in Japan, as variation was larger in our study area (ca. 0.56) than in the study area (ca. 0.25) of Donner (2011). Indeed, southern islands in Japan are encircled by the strong Kuroshio boundary current flowing poleward, and tropical waters brought by the current can cause faster warming than the global average (Wu et al., 2012), leading to high levels of historical SST variation in the area.
Optimized filtering thresholds were smaller than 1 C in GLM and RF models using only one thermal index. Previous studies have suggested that thermal stress not exceeding 1 C can induce coral bleaching (Brown, 1997;McWilliams et al., 2005;Kleypas, Danabasoglu & Lough, 2008). Nevertheless, bleaching thresholds had not been statistically optimized before our study.
Random forest was an excellent method for predicting coral bleaching and could be used more widely in studies of coral ecology. The use of RF has much increased in ecological studies in the 10 years since the introduction of Cutler et al. (2007). RF shows considerable potential for ecological analyses including classification, regression, and survival, due to its high predictive accuracy and its ability to model complex interactions among explanatory variables (Cutler et al., 2007).
Models with multiple environmental variables are becoming more popular, and show high explanatory power when modeling bleaching Welle et al., 2017). We found UV radiation to be an important explanatory factor for coral bleaching, consistent with previous studies (Hoegh-Guldberg, 1999;Maina et al., , 2011 The high performance of our bleaching model at 1 km resolution has practical implications for the local and regional management of coral reefs. Our predictions revealed high frequencies of coral bleaching in many parts of the Ryukyu Islands. However, predictions of bleaching frequency were based on the lowest levels of bleaching severity; hence, further analyses may be required to establish full distributions of bleaching frequencies according to levels of severity. Practical management to reduce the risk of coral bleaching should include control of coastal water turbidity (Fabricius, 2005). Increases in water turbidity by terrestrial runoff may decrease the resistance (Wooldridge & Done, 2009) or resilience (Hongo & Yamano, 2013) of corals to bleaching. Turbid coastal regions may provide refuges from climate warming due to limited increases in temperature and solar radiation (Cacciapaglia & van Woesik, 2016). However, coastal turbidity may increase the incidence of coral diseases and promote the growth of competing algae (Fabricius, 2005). Consequently, coastal turbidity should be carefully managed.
Reducing UV radiation may reduce bleaching risk and may constitute a powerful adaptive measure against climate warming. In the Onna Village of the Ryukyu Islands, in situ reduction of UV radiation with no increase in water turbidity has already been tested (Okinawa Prefecture, 2017). Reduction of UV radiation by 30-44% with large fishery nets resulted in a survival rate of 80% in cultured coral colonies in the summer of 2016, when the most severe thermal stress was recorded in the 2004-2016 study period (Kayanne, Suzuki & Liu, 2017). Reduction in UV radiation was similar to that used in our study (40%), so our predictions could provide a quantitative basis for future reef management in this area.

CONCLUSION
Predictive performance of coral bleaching models can be improved by the use of optimized thresholds, multiple environmental influences, and multiple modeling methods. Both high-resolution modeling and observational records (i.e., the Sango Map Project) enabled high performance of bleaching predictions (Oliver, Berkelmans & Eakin, 2009). We provide a template for selecting appropriate indices to predict bleaching, and our research methods could be applied to coral-habitable areas globally. Our highresolution predictions also provide a quantitative basis for the local and regional management of coral reefs . Although corals are suffering from high risks of bleaching globally, our study suggests that reducing UV radiation may be a key tool to improve coral resilience in the coming decades. Holistic bleaching models operating at finer spatial resolutions and incorporating variations in intrinsic thermal tolerance, historical effects of previous thermal impacts, and local environmental conditions should be the focus of future research. Such models will become indispensable as the effects of local and global stressors on corals continue to increase.