The dataset utilized in this study was part of the long-term, daytime monitoring efforts of coral reefs in West Hawai‘i conducted by the Hawai‘i Division of Aquatic Resources (DAR) at 23–25 permanent sites from 1999 to 2019 spanning over 140 km of the West Hawai‘i coastline (Supplementary Figure 1). Surveys were repeated 3–6 times annually, where at each site, four replicate transects were surveyed and averaged for the final data point (individuals m^–2) with annual, site-level averages used for the analysis. Each transect was 25 m long by 4 m wide, with the first two transects placed in opposing directions from a central pin, and the second two transects placed 10 m inshore and approximately parallel to the first two. Transects were located in coral reef habitats, at depths ranging from 5–15 m with a mean of 11.6 m. All species of fish and selected invertebrates were recorded along the transects.

A pair of divers each surveyed a 25 m × 2 m belt on either side of the transect line with counts combined to generate the full transect width. The divers first performed a “high” swim, swimming down the transect and recording larger mobile species, mid-water species and uncommon species. On return, the divers performed a “low” swim quantifying the remaining fishes, including cryptic species and recruits. During the “low” swim, sea urchins were also counted but only included “free” urchin species (Tripneustes gratilla, Chondrocidaris gigantea, Diadema paucispinum, Echinothrix calamaris, and Heterocentrotus mammillatus) and not rock boring species. All fishes were classified into functional groups, with herbivores further subdivided into scrapers, grazers and browsers, based on trophic classifications and diet information from Donovan et al. (2018). Only herbivorous fish species were considered in the analysis.

Co-located information on macroalgal, algal turf, hard coral, and CCA cover was available for 6 years of the DAR surveys; 2003, 2007, 2011, 2014, 2016, and 2017. Surveys were conducted following a standard operating procedure used by the Hawai‘i Department of Land and Natural Resources’ DAR (Williams et al., 2006), where photographs were used to survey coral reef benthos, taken at 1 m intervals along the length of fixed monitoring transects. DAR staff analyzed images using Coral Point Count with Extensions (Kohler and Gill, 2006) where benthic habitat cover was averaged across the four transects per site.

The United States National Oceanic and Atmosphere Administration (NOAA) Coral Reef Watch program (CRW) generates daily, satellite-derived global 5 km resolution SST dating back to 1985 (NOAA Coral Reef Watch, 2017). From the SST observations, CRW also generates sea surface temperature anomaly (SSTA) by subtracting the long-term (7 years) SST mean from the SST at that given location at that same time of year, with a positive anomaly indicating the SST is warmer than average. SSTA data were extracted for each of the survey sites from 1999 to 2019 with yearly SSTA means derived from daily SSTA measurements (Supplementary Figure 2). CRW temperature data were accessed at https://coralreefwatch.noaa.gov/product/5km/index.php#data_access.

To determine how benthic cover and herbivore densities were affected by temperature, survey data were matched with the SSTA conditions of the previous year, i.e., we expected changes in benthic cover and densities would lag behind changes in the environment. Previous research has shown that a 1-year lag in seawater temperature data (rather than current or 2-year lagged data) best explains changes in fish and benthic structures (Halford and Caley, 2009; Edmunds and Lasker, 2016; Ribeiro et al., 2018; Petraitis and Dudgeon, 2020), where 1-year lags in physical conditions are often used to assess changes in coral reef composition (Cooper et al., 2015; Donovan et al., 2021). Here we have adopted a similar approach with our analysis. If the previous yearly SSTA mean at a given survey site was positive, this was considered a “warm” year and if the SSTA was negative, this was considered a “cool” year for the analysis.

All statistical procedures were performed in R software (v 4.0.2). To assess whether benthic cover differed between positive and negative SSTA, we used bootstrapping with 10,000 iterations of site-level observations for the 6 years of benthic data to generate 95% confidence intervals (CI) for percentage benthic cover across warm and cool years. Non-overlapping 95% CIs was an indication that means were significantly different between depths (Cumming and Finch, 2005). Since herbivore counts were performed annually, to assess the impact of temperature on herbivore density, the mean annual density was compared with mean annual SSTA using Spearman’s correlations.

We used generalized additive models (GAM) to identify habitat-abundance associations and account for non-linear relationships observed between benthic cover and herbivore groups using the mgcv package in R. For these models, only fish and urchin data for which co-located benthic data were available were used (n = 148). Eight models were developed where we included five variables in each model; sea urchins, browsers, grazers, scrapers, and site. Conventional random effects can be incorporated into GAM models (Wood, 2006) where “site” was specified as a random effect in our models to account for possible spatial autocorrelation of sites and population sample bias (Zuur et al., 2009). A separate model was conducted for each benthic type – algal turf, macroalgae, CCA, and coral. To compare reefs exposed to positive and negative SSTA, we categorized all surveys into “warm” (n = 73) and “cool” (n = 75) years, performing separate analyses depending on temperature. Collinearity among variables was assessed prior to all analyses via regression analysis and using variance inflation factors (VIF) where the highest VIF for any driver was less than 2, lower than our a priori cut-off value of 3 (Zuur, 2012).

All models were performed using the following procedure. Models were fitted with a beta distribution, suitable for proportional data which cannot be modeled as binomial. Here, the percentage of benthic cover was represented as a proportion between 0 and 1. Model validation was performed by assessing the heterogeneity of the error distribution and histograms of the residuals. Given our small number of response values, we limited the number of knots to five to prevent overfitting of our models.

We used model selection and averaging procedures from the MuMIn R package (Barton, 2013) calculating Akaike’s information criterion (AIC) corrected for small sample size (AICc) (Anderson, 2008) and the AICc-based relative importance weights (w). Here, all possible models through various combinations of variables were run using the dredge procedure. P-values derived from GAM models are usually approximates and typically low, therefore, we used a weight of evidence approach to test for driver significance (Burnham and Anderson, 2002). The variable importance output gives the total weight of each driver across all possible models, that is, the summed weights across all models. As models sum to 1, values close to 1 indicate drivers that occur in large portions of the models and the higher the probability that that driver is important (Burnham and Anderson, 2002).

To display the relationships between benthic cover and each herbivore group, we used the predict procedure from the mgvc R package, to generate visualizations of smoothers of the top ranked models. The top ranked models included all models with a weight > 0.05. Model averaging generated coefficients that we then used for prediction. The predict function takes the averaged, fitted GAM model and produces predictions given model covariate values, where here the predictor variable of interest had values equally spaced between the variable’s minimum and maximum values, and all other drivers were set to their means. We used the resulting output to generate smoothers with the response of benthic cover shown against untransformed variables to ease interpretation of the effect of each herbivore group on benthic cover.

Mean coral cover was ∼37% from 2003 to 2014. In 2016, there was a large decrease where coral cover was 18.9 ± 1.4% (mean ± standard error). In 2017, coral cover dropped to a mean of 17.8 ± 1.4% (Figure 1A). Macroalgal (MA) cover was low, with values ranging from 0.2 to 1.0% of benthic cover between 2003 and 2017. Algal turf cover represented the dominant benthic class during the survey period, showing a general increase in mean cover across the permanent monitoring sites from 42.3 ± 2.3% in 2003 to 55.3 ± 2.4% in 2017. CCA cover was 8.3 ± 0.8% in 2003 and increased to 17.8 ± 2.4 in 2017 (Figure 1A).

From 1999 to 2019, there was a steady increase in sea urchin density across the permanent plots in West Hawai‘i, with a mean density increase of 0.15 in 1999 to 1.09 sea urchins per m² in 2019 (Figure 1B). Scrapers showed more variable mean densities across the survey period with a mean of 0.017 ± 0.003 in 1999 and 0.024 ± 0.004 in 2019 (Figure 1C). Browsers had a mean density of 0.013 ± 0.002 in 1999 and 0.016 ± 0.002 in 2019 (Figure 1C). Grazers showed an increasing trend across the survey period with a mean density of 0.28 ± 0.02 in 1999 and 0.56 ± 0.06 in 2019 (Figure 1D).

Bootstrapped estimates of CI for benthic cover show that there was significantly higher coral cover in cool years (36%) than in warm years (24%). Comparing this result to Figure 1, the majority of this coral loss occurred between 2014 and 2016. On the other hand, there was greater CCA cover in warm years (14%) than in cool years (6%). Both algal turf and MA cover were similar between warm and cool years (Figure 2).

When comparing mean herbivore density and mean SSTA, there were significant, positive correlations for three of the herbivore groups: sea urchins, grazers, and browsers, where there were greater densities with increases in SSTA (Figure 3).

When considering the relationship between algal turf cover and each herbivore group in cool years, the top models contained all factors at least once (Supplementary Table1 and Figure 4). Positive relationships were seen between algal turf cover and sea urchin and scraper densities, where algal turf cover increased with increasing herbivore density (Figure 4). A negative relationship between algal turf cover and grazer density was seen with a flatter relationship between algal turf cover and browser density. Site had the highest variable importance weight (w = 1) indicating significantly different site differences influencing the relationship between algal turf cover and herbivore group. Grazers had the second highest variable importance weight across all models (w = 0.973; Supplementary Table 1 and Figure 5) indicating that this variable was more strongly coupled to algal turf cover than other herbivore groups, which had much lower weights. Additionally, the top model containing grazers and site as the only two variables explained the greatest amount of variance (adj-R² = 0.9).

For MA, the top models contained all factors except site (Supplementary Table 1 and Figure 4). A negative relationship was seen between macroalgal cover and grazer density with a slightly less steep negative relationship observed between MA and sea urchin density. A slightly positive relationship of MA cover was seen with browser density until a density of ∼0.025 individuals/m² where cover began to decrease. A flat relationship between MA cover and scraper density existed until a density of ∼0.03 individuals/m² where cover began to increase steeply (Figure 4). Scrapers and browsers had the highest weights across all models (w = 0.990 and 0.862 respectively; Figure 5) with the top model explaining a moderate proportion of variance (adj-R² = 0.42).

For CCA, the top models contained all factors except grazers (Supplementary Table 1 and Figure 4). Positive relationships were seen between CCA cover and both sea urchin and browser densities, with a slightly negative relationship between CCA cover and scraper density (Figure 4). Sea urchins and browsers had the highest weights across all models (w = 0.999 and 0.986 respectively; Figure 5) with the top model explaining a high amount of variance (adj-R² = 0.65).

For coral, the top models contained all factors except browsers (Supplementary Table 1 and Figure 4). A positive relationship between coral cover and grazer density was seen with a slightly negative relationship between coral cover and scraper density. There was a steeper negative relationship between coral cover and sea urchin density (Figure 4). Similar to algal turf cover, site had the highest variable importance weight (w = 1), followed closely by sea urchins (w = 0.992) and grazers (0.938; Figure 5) with the top model explaining a very high amount of variance in coral cover (adj-R² = 0.93).

When considering the relationship between algal turf cover and herbivore group across warm years, the top models contained all factors except site (Supplementary Table 1 and Figure 4). Positive relationships were seen between algal turf cover and sea urchin, grazer and browser densities, where algal turf cover increased with increasing herbivore density (Figure 4). A slight positive relationship was seen between algal turf cover and scrapers. Sea urchins had the highest weight across all models (w = 0.947; Figure 5) however only a low amount of variance was explained by the top model (adj-R² = 0.23).

For MA, the top models contained all factors except site (Supplementary Table 1 and Figure 4). Positive relationships between MA and both grazer and browser densities were seen with a slightly positive relationship between MA and sea urchin density. On the other hand, a slightly negative relationship between MA and scraper density was seen (Figure 4). Browsers had the highest weight across all models (w = 0.986; Figure 5). A low amount of variance was explained by the top model (adj-R² = 0.29).

For CCA, the top models contained all factors except site (Supplementary Table 1 and Figure 4). Positive relationships between CCA cover and both sea urchin and grazer densities were seen, with flatter relationships between CCA cover and both browser and scraper density (Figure 4). Browsers had the highest weight across all models (w = 0.999; Figure 5) where, similar to the models with algal turf and MA, only a low amount of variance was explained by the top model (adj-R² = 0.22).

For coral, the top models contained all factors except browsers (Supplementary Table 1 and Figure 4). Negative relationships between coral cover and both sea urchin and grazer density were seen with a positive relationship between coral cover and scraper density (Figure 4). Sea urchins had the highest weight across all models (w = 1) closely followed by grazers and site (w = 0.999 for both; Figure 5) with a high amount of variance explained by the top model (adj-R² = 0.84).

Over the past 20 years in West Hawai‘i, hard coral cover has decreased, with the steepest decline in cover after 2015. Decrease in coral cover are associated with positive SSTA, illustrating the impact of the severe marine heatwave in 2015 where nearly 50% of coral cover was lost to mortality (Kramer et al., 2016). Macroalgal cover represents a minimal benthic component of coral reefs surveyed in West Hawai‘i with no significant differences in cover between warm and cool SSTA. CCA cover has shown a consistent increase since 2011, likely driven by positive SSTA where there is significantly higher CCA cover in warm years in comparison to cool years. Positive impacts of warming on CCA cover have been noted in other studies (O’Leary and McClanahan, 2010). Algal turf cover has also increased over time, representing a mean of 55% of cover in 2017. For algal turf, there are similar amounts of cover across warm and cool years indicating that increased temperatures are not driving changes in cover, although cover is not necessarily reflective of algal productivity, which is known to increase with increased temperature (Connell and Russell, 2010; Bender et al., 2014). Algal turf cover is highly negatively correlated with coral cover (Supplementary Figure 3), where the increases in algal turf cover are likely related to the increased space available for colonization after coral death (Kramer et al., 2016), with opposite trajectories seen in algal turf cover and coral cover from 2014 onward.

For all herbivore groups, there was a higher density in 2019 in comparison to 1999, where sea urchins and grazers show the clearest upward trends over the past 20 years. For sea urchins, grazers and browsers, these increases in density are influenced by increasing temperature, with significant positive relationships between density and SSTA seen for these three groups. For scrapers, a slight positive relationship between density and SSTA was observed but this was not strong or significant. The increases in scraper density after 2015 could have been influenced by factors other than temperature. For example, a ban on the use of SCUBA gear while spearfishing in West Hawai‘i was enacted in 2013. This may have helped boost scraper densities as the scraper group is mainly composed of parrotfishes, a family often targeted by spearfishers. The increased densities of sea urchins, grazers and browsers observed with positive SSTA indicate that thermal tolerance limits have not yet been reached for these herbivores on Hawaiian reefs, although it is important to note that biomass and fish size were not assessed here. With ongoing ocean warming, it will be important to assess whether densities remain on an upward trajectory, especially if thermal tolerances are breached.

Our analyses revealed herbivore-specific relationships with benthic types, where the direction of the relationship differed depending on temperature. In cool years, algal turf cover decreased with increasing grazer density. The very high importance of grazers in the model suggests that this group of herbivores are most tightly linked to algal turf cover. These results suggest that during cool years, grazers are either inhabiting areas with lower algal turf cover, or that grazers are contributing to reductions in algal turf cover. This relationship is reversed in warm years where increases in algal turf cover become coupled with increases in grazer density, a pattern observed across all herbivore groups.

Ectotherms are physiologically plastic, likely tolerating increases in temperature through herbivory (Pörtner, 2002; Seebacher et al., 2015). Studies have shown aquatic animals increase herbivory in response to high temperatures, including zooplankton (Boersma et al., 2016), freshwater fish (Prejs, 1984), and parrotfish (Smith, 2008). Positive correlations between algal turf cover and herbivore densities have also been observed in other studies (Russ et al., 2015; Morais et al., 2020). This study shows that in warm years, herbivores are responding to spatial patterns in algal distributions where grazing surface area is a strong predictor of the abundance of herbivorous fishes (Tootell and Steele, 2016; Oakley-Cogan et al., 2020).

Similar relationships were seen with MA, where increases in MA cover were coupled with increases in grazer and browser densities in warm years. In cool years, we see a negative relationship between MA cover and grazer and browser densities, a role that has been identified across many reef systems (Williams et al., 2019), especially for browsers whose dominant food source is macroalgae (Hoey and Bellwood, 2008). With macroalgae being the dominant food source, the strong positive relationship with MA cover and highest model weight in warm years suggests that during warmer years, browsers are resource-limited and are tracking nutritional resources closely. Similar, but much weaker relationships between MA and sea urchin density were also seen. In these surveys, the dominant species of sea urchin is T. gratilla, a species which shows a greater preference for brown algae such as Turbinaria spp. (Lewis et al., 2018), which is rare on the reefs surveyed here (M. Lamson, pers. comm.). MA cover is very low in West Hawaii, which could explain the weak correlations observed here between MA and sea urchin density.

For CCA, the positive relationship between CCA cover and sea urchins across warm and cool years are steeper than any other herbivore, with sea urchins having the highest importance for explaining CCA cover. This could indicate that sea urchins are clearing some space for CCA to grow, a response also documented on other reefs (Sammarco et al., 1974; Johansson et al., 2010). Similar results in the Caribbean have been noted with CCA cover increasing with increasing grazer density due to removal of competitive fleshy algae (O’Leary and McClanahan, 2010).

For coral, we see decreases in cover with increasing sea urchin density across both warm and cold years. Whether sea urchins are causing coral cover decline or merely associating with low cover coral reefs in Hawai’i will be important to determine, especially as sea urchins show a strong upward population increase. It is important to note that of all sea urchins surveyed here, only T. gratilla shows similar abundances from day to night, while all other species show an increase in night-time surveys (Walsh et al., 2004), indicating that the reported urchin densities are likely underestimated.

On the other hand, coral cover increases with grazer density in cool years, but decreases with grazer density in warm years. A positive relationship between coral cover and scraper density emerges in warm years. Previous research has shown that highly abundant sea urchins reduce the density of coral recruits and increase the mortality of small coral fragments (Dang et al., 2020) where sea urchins prey on living coral tissue of adult colonies (Ruiz-Ramos et al., 2011). There is, however, the confounding strong correlation between algal turf and coral cover which could be influencing the negative relationships observed here, where these negative relationships could be driven by herbivores favoring nutritional resources and areas with greater algal turf rather than coral cover.

Additionally, in cool years, site variability greatly influenced algal turf cover and CCA cover, but barely in warm years, suggesting a homogenizing impact of increased temperature. For MA, site had low weight across both warm and cool models, likely due to low MA cover throughout all West Hawai’i. On the other hand, site had high weight across both warm and cool models for coral cover, suggesting important roles of local adaptation and local conditions influencing cover.

Despite the large increases in herbivore densities over the past two decades, algal turf continues to increase in cover. Herbivore abundance seems to be responding accordingly to increases in food availability where increased algal cover is driving energetic shifts towards greater herbivore populations in warm years. While algal turf cover would increase at a much greater rate without herbivores, the results suggest that current grazing levels may not be enough to cause reductions in algal turf cover, where algal production is still exceeding grazing capacity. The increased metabolism and consumption of herbivorous fish driven by warming are not outpacing algal turf productivity but is still an important mechanism contributing to control of algal turf in a warming ocean. The relationships between herbivore groups and benthic cover are conceptualized in Figure 6. A mismatch between the thermal limits of herbivorous fish and that of algal turf species could drive a greater competitive advantage for algae during warm years.

Other factors are likely contributing to increases in algal turf cover, where eutrophication and sedimentation can be just as important as a reduction in herbivory in causing phase shifts (Gilby et al., 2015). Increases in nutrients and sediment loads can stimulate algal growth with bottom-up processes precipitating coral-algal phase shifts on Mesoamerican reefs despite healthy herbivore populations (Arias-Gonzalez et al., 2017). As nitrogen pollution strongly impacts Hawaiian resource fish populations, consisting of mainly herbivores (Foo et al., 2021), addressing water quality in Hawai’i should be a priority in conjunction with herbivore management. Additionally, across algal turf, MA, and CCA benthic types, the amount of variance explained by the models greatly decreased when comparing warm to cool years with a shift in the relative importance of different herbivore groups depending on temperature, indicating that during warm years, factors other than herbivory are important in driving cover. Warming reduces the predictability of relationships between benthic types and herbivore groups.

Further, there are clear differences between the impacts of herbivorous fish and sea urchins on Hawaiian coral reefs, with critical implications for herbivore management. The groups are not functionally equivalent where grazer densities are more strongly associated with decreased algal turf cover and greater sea urchin densities with increasing CCA cover. A balance of sea urchin abundance relative to herbivorous fishes will be important to ensure a positive influence of sea urchins on reef health and resilience, and prevention of high substrate erosion (McClanahan, 1995).

It is important to highlight the limitations of this study as well as critical areas of future research. Here, we generalize the impacts of warm and cool years on the benthic-herbivore relationship (Figure 6), while it is likely that both the extent of heat stress experienced and previous heatwave events would influence the trajectory of change (Maynard et al., 2008; Hughes et al., 2018b). Future studies, which monitor both benthic cover and herbivore abundance annually could utilize rates of change alongside various temperature metrics, to understand how severe heatwaves and ecological memory potentially modify the benthic-herbivore relationship.

Investigating algae-fish relationships has important consequences for local management. Our analyses show that increasing temperature disrupts the normal habitat associations of herbivore densities on coral reefs. Temperature changes should, therefore, be integrated into herbivore management strategies. The future health of coral reefs and recovery from bleaching will largely be determined by local herbivore population effectiveness in algal control (Williams et al., 2019), particularly in West Hawai‘i where algal turf already dominates reefs (Jouffray et al., 2015). Our results show that herbivore populations, despite increasing in density over the past decade, are likely impacting algal turf but not enough to cause a reduction in algal turf cover over time. Increasing temperatures may be establishing a higher threshold of herbivore densities required to cause a reduction in algal turf cover on Hawaiian reefs. In conjunction with bottom-up drivers of algal turf cover, i.e., mitigation of land-use impacts to improve water quality, algal-dominated reefs will benefit from protection measures to increase herbivore populations and to encourage the increased resilience and expansion of corals. Reducing coastal eutrophication will also be critical in mitigating future algal turf growth in a warming ocean, where herbivores are likely to be pushed outside of their optimal thermal windows.