El Niño-driven phase shift to algal dominance on Isla del Caño’s coral reefs: implications for urgent restoration

Sea surface temperature

Daily SST data for Isla del Caño (8°42′59″N 83°53′06″W) were obtained from Coral Reef Watch 5 km Heat Stress Product (1985–2025; CoralTemp V3.1; Heron et al., 2015; Skirving et al., 2020). The Maximum of Monthly Means of climatological SST (MMM) was computed for 1985–1990 and 1993, excluding 1991–1992 due to anomalies caused by the Mt. Pinatubo eruption (Heron et al., 2016). Degree Heating Weeks (DHW), indicative of thermal stress accumulation, were calculated as days exceeding MMM + 0.5 °C, using a 12-week rolling sum and dividing this by seven to express thermal stress in °C-weeks (Liu, Strong & Skirving, 2003; Liu et al., 2018). MMM + 0.5 °C was used instead of MMM + 1 °C to increase the sensitivity for capturing localized bleaching events (Lachs et al., 2023). The bleaching threshold was defined as MMM + 0.5 °C and the cooling threshold was defined as the climatological monthly mean (1985–1990 and 1993) minus 0.5 °C. This value represents a lower thermal bound used to identify days that may contribute to cold-water exposure or thermal relief, potentially moderating heat stress impacts. Changes in the frequency of days breaching this cooling threshold were analyzed to assess whether such thermal buffering has declined over time. A linear regression model was used to assess the long-term SST trends, with date as the predictor, and the slope estimate was used to calculate annual and decadal warming. The influence of ENSO on SST variability was also explored using linear regression, using date and Oceanic Niño Index (ONI) as the predictor variables. Diagnostic residual vs. time and normal Q-Q plots indicated that heteroscedasticity and outliers had minimal impact on the model. The change in annual frequency of days breaching bleaching and cooling thresholds between 1985 and 2025 were explored using linear regression models.

Study area, survey methods and environmental parameters

Between 2019 and 2025, nine coral community sites were surveyed within the Reserva Biológica Isla del Caño (8°42′45.70″N 83°53′23.20″W) and on the northern coast of the Osa Peninsula within the Osa Conservation Area (ACOSA), Costa Rica (Table S1; Fig. 1). Of the investigated coral communities, which ranged in depth from <2 m to 16 m (Table S1), eight sites were surveyed for benthic cover and coral health between February 2024 and February 2025, including Ancla, Barco Somero, Barco Profundo, Cueva de Tiburones “Cueva”, Esquina, Este Intermedio, San Josecito, and Tina (Fig. 1). A ninth site, Chorro, was additionally surveyed for benthic cover in 2019, 2020 and both Chorro and Cueva in 2021 (Table S2). Landslides were observed at Isla del Caño in November 2025, with large sediment plumes still present in January 2025 (Fig. S2).

Benthic cover was recorded using different, but comparable, methods, including the categories: live coral (coral), sand, rubble, rock, CCA, recently dead coral (dead coral), turf algae, Caulerpa, cyanobacteria, and macroalgae. Recently dead coral was defined as an absence of tissue with recognizable skeletal features and absent or minimal overgrowth. Caulerpa was included as a distinct category because of its documented proliferation and impact on Costa Rican reefs (Fernández, 2007). For surveys in 2019, 2020 and 2021, three haphazardly-laid 10 m point intercept transects were used, recording the benthic cover category every 20 cm Sistema (Nacional de Áreas de Conservación (SINAC), 2021). Surveys conducted in February and March 2024 used 1 × 1 m quadrats to record the proportion of each benthic category over the full length of three 30 m (San Josecito) and three 10 m transects (Barco Somero, and Cueva; Table S2). For the remainder of the surveys, point and line intercept transects were conducted, recording substrate categories every 10 cm (in April 2024), and then the continuous 10 m line for three or six transects at each site (Tn = 270 transects; Table S2). Coral diversity, abundance and health were recorded from April 2024 using three or six 1 ×10 m belt transects, separated by 5 m. Corals were identified to genus, or species when possible, and categorized as healthy, pale, partially bleached (patches), mostly bleached (>50% of the colony) and bleached (whole colony). An individual coral colony was defined as having distinct area of tissue, meaning that multiple colonies could arise from a coral undergoing partial mortality. Benthic cover data were converted into proportions, enabling comparison of data collected using different survey methods. Total coral abundance, the abundance of each coral taxa, total bleaching prevalence and bleaching prevalence per species were calculated for each belt transect. Temperature loggers (ONSET HOBO MX2202) were programmed to record every 30 min and deployed at Ancla 15 m, Barco Profundo 16 m, Barco Somero 11 m and Tina (8 m), though remained temporarily and variably functional between April and August 2024.

Fieldwork was conducted under permit numbers: SINAC-ACOSA-DT-PI-R-006-2019, SINAC-ACOSA-DASP-PI-R-067-2021, SINAC-ACOSA-D-PI-R-013-2023, SINAC-ACOSA-DR-PI-R-024-2024.

Data analysis

All analyses were conducted in R (R Core Team, 2021) and are summarized in the supplementary section (Table S4). Bayesian models were implemented in brms (Bürkner, 2017) using Stan (Carpenter et al., 2017) Hamiltonian Monte Carlo (HMC) sampling via the No-U-Turn Sampler (NUTS). Unless otherwise stated, models were fit with four chains of 4,000 iterations each (2,000 warm-up), using an adapt_delta of 0.99 for convergence and four cores for parallel computation. Model convergence was confirmed using R-hat values (≈1.00) and effective sample sizes (Bulk_ESS, Tail_ESS). Posterior predictive checks and leave-one-out cross-validation (LOO) were used to assess model fit and predictive accuracy. Proportion data (e.g., coral cover, bleaching prevalence) were bounded between 0.0001 and 0.9999 to avoid boundary issues. Continuous predictors such as date were standardized using z-scores (e.g., z_date) to improve numerical stability and allow comparability of effect sizes. Models involving proportional data were fitted using beta or zero-inflated beta (ZIB) distributions using a logit link for the mean and an identity link for zero inflation, where applicable. All p-values for non-Bayesian tests were Bonferroni-adjusted where appropriate.

Coral bleaching prevalence and taxa-specific models

Site-level variation in bleaching prevalence was modelled using a Bayesian beta regression with logit link and a random intercept for site. Logit-scale intercept was used to estimate a baseline bleaching prevalence through inverse logit transformation. A multivariate ZIB regression model was used to analyze bleaching prevalence among different coral genera, incorporating structural zeros to account for sites where certain coral species were absent. The response variables included bleaching prevalence for coral taxon categories of Pocillopora, Porites, Pavona, Psammocora, and “Other” corals, with each modelled as a function of their respective structural zeros and those of other genera.

Coral diversity

The Shannon Diversity Index was computed using total abundance of 13 coral taxa, including the Other category with the vegan package. Rows with missing species count values were removed. Due to significant deviations from normality and variance homogeneity (Shapiro–Wilk and Levene’s p < 0.05), Kruskal–Wallis tests were conducted to determine differences in taxa abundance among sites. Dunn’s post-hoc tests (Bonferroni-adjusted) were performed to identify pairwise differences between sites.

Benthic composition shifts

Shifts in benthic composition were assessed using similarity percentage (SIMPER) analysis (vegan) conducted on “During Bleaching” (March to end of September 2024) and “After Disturbance” (December 2024 to February 2025) time periods. Ecologically relevant benthic variables were included; coral cover, macroalgae, turf, CCA, cyanobacteria, bleached coral and dead coral. SIMPER was performed at two levels: (1) across all sites combined to evaluate broad-scale compositional changes and (2) separately for each site to highlight site-specific patterns. Both analyses used 100 permutations to ensure robust estimates. Variables were ranked based on their contribution to compositional differences, and their mean before and after values were extracted to evaluate temporal trends. P-values were derived from permutation tests comparing within-group and between-group Bray–Curtis dissimilarities. These metrics identify the variables contributing most to observed shifts in community composition. To avoid artificial inflation of change, missing values were replaced with the mean for the time period at that site, or the global mean if the variable was entirely absent. Percentage change and p-values were used to identify the strongest drivers of composition shifts.

Temporal trends in coral cover and benthic covariates

To assess temporal changes in coral cover, we used Bayesian ZIB regression models implemented in brms, with a logit link for the beta mean and identity link for the zero-inflation component. All models included site as a random intercept. Candidate models differed in how time was modelled (e.g., linear vs. polynomial), and we compared their predictive performance using LOO, based on expected log predictive density (ELPD). The linear time model provided the highest predictive accuracy and was thus selected for inference. Based on this, we evaluated site-specific coral cover trends using simple linear regressions over scaled date (z_date), offering interpretable measures of change at each site. While more complex models could capture non-linear or heterogeneous dynamics, the linear model proved most appropriate given the temporal resolution and consistency across sites.

Temporal trends in benthic cover for algal groups (macroalgae, Caulerpa, turf algae, cyanobacteria, and CCA) were analyzed using Generalized additive models (GAMs) with beta error distribution and logit link. GAMs allowed flexibility in capturing potential nonlinear trends while controlling for site-level random effects. Penalized smoothing in mgcv allowed interpretability while avoiding overfitting. Models were fitted using REML, and diagnostics (residual checks, smooth term significance) were conducted to ensure model validity. A linear model was used to assess the relationship between coral cover and turf algae while accounting for site and temporal effects. Coral cover was modeled as a function of turf algae cover, site, and standardized time (z_date). The model was fitted using ordinary least squares regression, with model assumptions checked via residual diagnostics.

Coral cover and local temperature

HOBO logger data from four locations: Barco Profundo, Barco Somero, Tina, and Ancla was used to evaluate differences in local temperature conditions across sites. Raw temperature data were processed to calculate monthly mean, maximum, and minimum temperatures per site. All values were aggregated at monthly resolution after cleaning and aligning timestamps to ensure maximum overlap in sampling dates among sites. A one-way ANOVA was conducted separately for mean, max, and min monthly temperature data to test for differences among sites. Tukey’s Honest Significant Difference (HSD) post hoc tests were performed to identify pairwise differences. All comparisons were adjusted for family-wise error rate, and significance was determined at α = 0.05. A beta regression model was used to analyze the effects of temperature and site on coral cover. The initial model included mean, minimum, and maximum temperature as predictors but was refined based on statistical significance and model performance. Minimum and maximum temperature were removed due to their non-significant effects (p = 0.585 and p = 0.429, respectively), resulting in a final model with mean temperature and site as predictors. The beta regression framework, with a logit link function for the mean and an identity link for the precision parameter (φ), was chosen to account for heteroscedasticity and the bounded distribution of coral cover. Maximum likelihood estimation (BFGS optimization) was applied, and model performance was assessed using pseudo R².

Temporal trends in coral taxa

Temporal trends in the proportional cover of Pocillopora, Porites, Pavona, Psammocora, and Other coral were modelled using a multivariate beta regression model with random slopes. This allowed proportional data to be modelled while accounting for variation in coral cover (random intercept) and divergent temporal trends among sites (random slope for z_date). Alternative models were tested to determine the best structure. Model comparison showed that the random slopes model provided a better balance between complexity and interpretability than more parameter-rich interaction models, which increased uncertainty without improving fit. Models were implemented in brms, and convergence was confirmed through standard diagnostics.

Multivariate benthic analysis

Principal component analyses (PCAs) were performed to explore multivariate gradients of key benthic parameters. An initial, full model PCA was conducted including biological and benthic variables, as well as total coral cover and the abundance each of Pocillopora spp., Porites spp., Pavona spp., Psammocora spp., and Other coral, and environmental parameters (Supplementary File). Variables with minimal influence on components were then removed from the analysis, including pH, salinity, phosphate, and oxygen. In the first two PCA analyses, Caulerpa and dead coral contributed to PC1; however, they were not observed in the most recent surveys and thus were absent from the final PCA. This absence reflects an ecological shift rather than a methodological change, indicating that these components were no longer present at detectable levels during the final sampling period.

PCAs were run on three datasets: (i) During Bleaching (all survey data up to July), (ii) After Bleaching (August to December 2024), and (iii) the Latest dataset (January or February 2025) for all sites, to capture potential shifts post-La Niña rains and landslide events. Each PCA used standardized variables (mean = 0, SD = 1), and prcomp in R was employed with centring and scaling. Zero-variance variables were removed, and incomplete cases were filtered out. PCA loadings were used to identify which benthic or environmental variables contributed most strongly to disparities among sites. To assess temporal trends in benthic community composition, PCA loadings were extracted and standardized to facilitate direct comparison and site scores were compared across the three time periods. Key benthic variables structuring PC1 were examined to identify benthic community composition shifts. All variables were transformed to z-scores (mean = 0, SD = 1) before PCA to prevent biases from differences in magnitude or measurement units. Only variables present in all three PCA results were included, excluding those filtered out due to low variance or collinearity. PC1 orientations were aligned across time periods to ensure consistent interpretation, with eigenvectors adjusted by multiplying by -1 when necessary.

Ecological recovery feasibility index

Composite, multi-metric indices are well established in restoration science, e.g., expert-driven priority scores for forests (Orsi, Geneletti & Newton, 2011) and mangroves (Su et al., 2022), and leading coral-reef restoration guides recommend site ranking, suggesting a call for analogous tools (Goergen et al., 2025; Shaver et al., 2020). However, no existing index integrates site-specific benthic data with data-driven weights and continuous scoring of both positive (coral cover, CCA, diversity; Nyström et al., 2008; Tebben et al., 2015) and negative (turf algae, cyanobacteria, bleaching; Ford et al., 2018) drivers. We designed the Ecological Recovery Feasibility Index (ERFI) as an initial, generalized screening tool to rank reefs by their baseline benthic-community conditions and identify those most likely to support natural recovery. The ERFI therefore produces a transparent, replicable restoration feasibility score.

To obtain PCA-derived weights we performed a PCA on six benthic variables (live coral cover, CCA cover, turf-algal cover, cyanobacterial cover, bleaching prevalence, Shannon diversity) using the Latest Data (January–February 2025). For each variable j, we extracted its loadings lj ₁, lj ₂, lj ₃ on principal components 1–3 and summed their absolute values to obtain: ${\displaystyle L_j=\left|l_{\left\{j1\right\}}\right|+\left|l_{\left\{j2\right\}}\right|+\left|l_{\left\{j3\right\}}\right|.}$

Expert direction weights, Dj, were assigned to each variable based on ecological effect: ${\displaystyle D_j=\left\{\begin{array}{c}\hfill +4ifj=coral;\hfill \\ \hfill +2ifj=CCAordiversity;\hfill \\ \hfill -2ifj=cyanobacteria;\hfill \\ \hfill -2ifj=bleaching;\hfill \\ \hfill -3ifj=turfalgae\hfill \end{array}\right\}.}$

For each site i and variable j, we computed the z-score: ${\displaystyle z_{\left\{ij\right\}}=\frac{\left(x_{\left\{ij\right\}}-{\overline{x}}_j\right)}{s_j}}$

where xij is the observed value, and x⁻j and sj are the mean and standard deviation of variable j across all sites.

To calculate composite weights and raw ERFI, we multiplied the PCA weight and the direction weight to obtain: ${\displaystyle C_j=L_j\times D_j}$

and then summed across variables to compute each site’s raw score: ${\displaystyle ERF{I}_i=\sum _j\left(C_j\times z_{\left\{ij\right\}}\right).}$

Finally, we completed a normalization step [0, 1] to facilitate interpretation. We rescaled the raw scores by: ${\displaystyle ERF{I}_i\ast =\frac{\left(ERF{I}_i-min\left(ERF{I}_i\right)\right)}{\left(max\left(ERF{I}_i\right)-min\left(ERF{I}_i\right)\right)}.}$

A worked numeric example of the above steps is provided in Table S3.

SST at Isla del Caño

SST increased significantly over the 40-year period (1985–2025) at Isla del Caño (R² = 0.095; F_(1,14642) = 1544, p = 2.2e−16; Fig. 2) with an MMM of 29.08 °C based on the highest monthly mean SST (1985–1990 and 1993). The estimated rate of SST increase was approximately 0.23 °C per decade, equating to a 0.92 °C warming over the 40-year period (Fig. 2). Incorporating the ONI into the regression model significantly improved the fit (R² = 0.28; F_(2,14606) = 2872, p = 2.2e−16), confirming that ENSO as a strong predictor of SST and modulator of SST anomalies. El Niño periods are visible in the Isla del Caño SST record, with peaks or clustered peaks in SST and DHW most prominent in 1997–98, 2010, 2014–17, 2020 and 2023–24 (Fig. 2). The triple peaks of DHW during the 2023–24 El Niño reached 1.6-fold higher than the peak of the 2016 event, and 2-fold greater than the peak of 1998 (Fig. 2). The highest SST was also recorded during the 2023–24 El Niño, reaching 31.2 °C, with peaks in 1998 and 2016 remaining below 31 °C (Fig. 2).

Linear regression models indicate that the annual frequency of days breaching the bleaching threshold (MMM + 0.5 °C) has significantly increased through time (Fig. S1), with an estimated additional 2.07 bleaching days per year (F_(1,39) = 12.94 ; p = 0.0009; R² = 0.23). Conversely, that the frequency of cool days has significantly decreased, by approximately 0.72 days per year (F_(1,39) = 17.62 ; p = 0.0002; R² = 0.29). Years prior to 2012 predominantly experienced frequent cool days, including in El Niño years when warm (bleaching) days were also frequent, including 1998 and 2010 (Fig. S1). After 2012, the frequency of cool days is markedly less, with no substantial cool periods occurring since 2021. Warm days were most frequent in 2023, amounting to two thirds of the year with no cool days (Fig. S1).

Local impacts, bleaching prevalence and coral susceptibility

Severe, mortality-inducing coral bleaching was observed at Isla del Caño from 2023 and throughout 2024 (Fig. 3), prior to the study period (Table S1). Numerous corals bleached and were observed to progress to a state of partial or complete mortality (Figs. 3A–3E), suggesting that the 2024 to 2025 survey data span the tail-end of an extreme bleaching event (Fig. 2). Some corals were photographed severely bleached and subsequently recovered (Figs. 3F–3H). Additionally, in late November 2024 during the rainy season, a series of landslides occurred on Isla del Caño (Figs. S2A and S2B), with a significant amount of sediment washed onto reef sites (particularly Cueva and Ancla; Fig. S2C). During this time the Biological Reserve was closed to the public due to unsafe conditions, and it is unclear exactly when the landslides occurred and the immediate extent of the sediment plumes. Crown-of-thorns starfish were also observed sporadically at several sites, including predating on massive coral species (Fig. S3).

There was a baseline bleaching prevalence of approximately 23% across all surveys, as calculated through logit transformation (−1.17; 95% CI [−1.54 to −0.79]), with significant variation in bleaching prevalence among sites (beta regression with logit link; SD 0.48; 95% CI [0.22–0.91]) and moderately low within site variability (phi 1.94; 95% CI [1.61–2.32]). Ancla, Barco Profundo, and Chorro demonstrated increasing bleaching prevalence through time (Fig. 4), whilst decreasing prevalence was observed at San Josecito. The remaining sites, Barco Somero, Tina, Cueva, Este Intermedio and Esquina had approximately constant bleaching prevalence through time (Fig. 4).

Baseline bleaching prevalence varied among coral genera. Model (ZIB) estimates indicate that Pocillopora spp. exhibited the highest bleaching prevalence at 33.9%, followed by Porites lobata (24.0%), Psammocora spp. (15.6%), and Pavona spp. (10.1%), while Other corals exhibited minimal bleaching at just 0.25% (Table 1). The estimated standard deviation of site-level variation in bleaching prevalence was highest for Porites lobata (SD = 0.19), suggesting greater among-site variability than other taxa, and was lowest for Other coral (SD = 0.08). Pocillopora spp., Psammocora spp. and Pavona spp. exhibited intermediate variability. The dispersion parameter (phi) estimates suggest moderate consistency in bleaching prevalence for most genera, while Pavona spp. showed moderate dispersion and Other coral had high dispersion, indicating close to zero bleaching in this group across all sites (Table 1). Zero-inflation was negligible across all coral taxa (Table 1), indicating that observed zeros were consistent with model expectations and that zero-inflated components were not necessary.

Coral diversity

A total of 12 coral taxa were observed and identified at Isla del Caño, including Pocillopora spp., Porites lobata, Pavona clavus, Pavona gigantea, Pavona varians, Pavona maldiviensis, Pavona chiriquiensis, Pavona frondifera, Psammacora stellata, Psammacora profundacella, Tubastraea coccinea and Gardineroseris planulata. Pocillopora spp. and Porites lobata were present at all sites, and Pavona clavus and Pavona gigantea were present at all sites except Este Intermedio and San Josecito, respectively. Gardineroseris planulata was the rarest coral taxon, observed only at Ancla, Barco Somero and Cueva (Table S5).

Shannon Diversity Index varied significantly among sites (Kruskal-Wallis test: χ² = 92.623, df = 8, p < 2.2e − 16; Fig. 5), with Cueva exhibiting the highest mean Shannon Index (1.28 ± 0.43), followed by Ancla (1.14 ± 0.39) and Barco Somero (0.90 ± 0.44). San Josecito had the lowest diversity (0.12 ± 0.16), with Este Intermedio also showing relatively low values (0.43 ± 0.31). Post-hoc Dunn’s test confirmed significant differences among sites (Table S6). For instance, Cueva had significantly higher diversity than San Josecito (p adj. = 6.30E−12), Este Intermedio (p adj. = 0.000002), and Esquina (p adj. = 0.02). Ancla also had significantly higher diversity than San Josecito (p adj. = 2.94E−12) and Este Intermedio (p adj. = 2.99E−06) (Table S6).

Community composition change

The SIMPER analysis identified key contributors to benthic composition changes across sites between the “During Bleaching” (March to August 2024) and “After Disturbance” (December 2024 to February 2025) periods (Tables 2 and 3). Turf algae and dead coral contributed most to overall community dissimilarity, with turf increasing by 70.62% (p = 0.010) and dead coral decreasing by 81.30% (p = 0.010; Table 2). Coral cover declined significantly by 40.42% (p = 0.010), from a mean of 19.33% to 11.46%. CCA increased modestly by 13.0%, though not significantly (p = 0.347). Macroalgae and cyanobacteria declined by 34.96% (p = 0.069) and 74.37% (p = 0.010), respectively. Bleaching also contributed to dissimilarity, increasing by 35.30%, but was not significant (p = 0.881).

Site-level patterns revealed strong declines in coral cover at several locations (Table 3). Tina experienced a 52.44% decline (p = 0.010), while Ancla and Esquina showed declines of 46.0% (p = 0.040) and 50.56% (p = 0.010), respectively. Coral cover also declined at Barco Somero (−8.4%) and Este Intermedio (−45.34%), though these were not significant (Table 3). Turf algae increased at Ancla, Barco Profundo, Chorro, Cueva, and Esquina, with changes ranging from 43.2% to 336.7%, but was only significant at Cueva (p = 0.010), Barco Profundo (p = 0.010), and San Josecito (p = 0.030; Table 3; Fig. S3). Macroalgae declined at most sites and was a notable contributor to benthic change at Barco Somero, Tina, Barco Profundo, San Josecito, and Este Intermedio (Table 3; Fig. S4). Recently dead coral at San Josecito declined by 99.75% (p = 0.010) and similarly at Cueva (−97.1%, p = 0.010), Chorro (−29.8%, p = 0.277), and Este Intermedio (−97.64%, p = 0.010). The mean percentage coral cover after disturbance was lowest at Barco Profundo (3.52%) and highest at Chorro (16.86%) (Table S7). Four sites had less than 10% coral cover after disturbance (Cueva, Barco Profundo, Ancla, and San Josecito), while Esquina had 10.59% and Este Intermedio 11.12%.

Temporal changes in coral cover

The initial proportion of coral cover varied among sites, with Tina exhibiting the highest mean cover at 55% (SD ± 0.09, n = 3) in mid-April 2024, and Barco Profundo the lowest, with less than 1% (SD ± 0.115, n = 3; Fig. 6). Bayesian ZIB regression model revealed a general decline in coral cover over time. The best-supported model included a fixed linear effect of time and a random intercept for site. This model estimated a significant negative linear effect (Estimate = −0.19, 95% CI [−0.29 to −0.08]), indicating a consistent decline in coral cover across the study period. The overall baseline coral cover was estimated at approximately 11.4% (Intercept = −2.05, 95% CI [−2.45 to −1.66]), with moderate variation among sites (SD = 0.53, 95% CI [0.28–0.98]). There was no evidence of substantial zero-inflation (zi = 0.00, 95% CI [0.00–0.01]), suggesting that complete coral absence was rare. Linear regressions exploring site-level variability in coral trajectories highlighted the strongest significant decline in coral cover at Cueva (slope = −0.059, SE = 0.014, p < 0.001; Fig. 6), followed by Chorro (−0.022, SE = 0.008, p = 0.009) and Tina (−0.199, SE = 0.085, p = 0.024). Coral cover declines at Este Intermedio (−0.137, p = 0.096), Esquina (−0.131, p = 0.204), Ancla (−0.079, p = 0.204), Barco Profundo (−0.046, p = 0.277), San Josecito (−0.050, p = 0.287) and Barco Somero (−0.012, p = 0.858) were not statistically significant.

Temporal trends in algal cover

GAMs used to assess temporal trends in algal cover demonstrated a significant non-linear temporal trend in macroalgae (edf = 5.012, p = 0.0055) with site-level variability (edf = 5.299, p = 0.0013), with an adjusted R² of 0.0139. Caulerpa exhibited a weaker temporal pattern (edf = 1.746, p = 0.0113), though site effects were negligible (p = 0.8621). Turf algae displayed a strong temporal trend (edf = 1.665, p < 0.0001) without a site effect (p = 0.88), explaining 33% of deviance. Cyanobacteria and CCA showed no significant temporal trends (p = 0.187 and p = 0.48, respectively), although CCA exhibited significant site-level variation (edf = 7.120, p < 0.0001).

A linear model exploring the relationship between coral cover and turf algae demonstrated a significant negative association (β = −0.131, SE = 0.039, t = −3.40, p = 0.00081), indicating that higher turf algae cover corresponded with lower coral cover. Site effects were significant at Barco Profundo (β = −0.087, SE = 0.034, t = −2.53, p = 0.012) and Tina (β = 0.059, SE = 0.030, t = 1.98, p = 0.049), highlighting site-specific differences in coral cover. Temporal trends were not significant (β = −0.035, SE = 0.033, t = −1.08, p = 0.282; Fig. S5).

Site temperature and coral cover

Mean monthly temperatures differed significantly among the four sites (ANOVA: F_(3,10) = 8.39, p = 0.0044; Fig. S6), whereas no significant differences were detected for maximum (F_(3,10) = 0.281, p = 0.838) or minimum temperature (F_(3,10) = 2.95, p = 0.085). Post hoc comparisons revealed that Tina had significantly higher mean temperatures than both Ancla (difference = +1.12 °C, p = 0.023) and Barco Profundo (+1.63 °C, p = 0.0035; Table S8). No other pairwise comparisons were statistically significant. Coral cover was significantly influenced by mean temperature and varied across sites. Mean temperature had a positive effect (Estimate = 0.877, SE = 0.275, z = 3.187, p = 0.0014), indicating higher coral cover at warmer sites. The model explained 32.91% of the variation in coral cover (Pseudo R² = 0.3291), suggesting that while temperature influenced coral cover, site-specific factors played a substantial role in shaping spatial patterns. The precision parameter (φ = 5.542, SE = 1.058, p < 0.0001) indicated low dispersion.

Temporal trends in the proportion cover of each coral taxa

The multivariate beta regression model with random slopes indicated variability in coral cover trends among sites. Posterior estimates showed no strong overall trends in the proportion cover of each coral taxa through time, though site-level differences may be driving variation in coral composition over time, as indicated by site-level variation in coral proportion changes (Table S9; Fig. S7).

Shift from coral to algal dominance

PCA were used to assess changes in reef structure across three sequential time points (Tables S10 and S11): During bleaching (March to July 2024), After bleaching (August to December 2024), and the Latest period (January to February 2025). Each PCA was interpreted independently and captured a high proportion of total variance in reef composition and collectively indicate shifts in community structure from coral dominance toward turf algae dominance at most sites (Fig. 7). While some sites have stabilized with remaining coral cover, others have continued to degrade, with San Josecito experiencing a near-total collapse. Pocillopora spp. was more predominant at shallow sites and was associated with high bleaching prevalence (Fig. 7). Pavona spp. and Psammocora spp. were initially associated with deeper reefs (Cueva and Ancla; Fig. 7A), though lost prominence as turf algae and bleaching increased (Figs. 7A and 7C). Similarly, Porites spp., while less strongly associated with deep sites, also reduced in dominance as turf and bleaching increased (Fig. 7). The absence of dead coral by the Latest Data PCA (Fig. 7C) suggests that sites initially characterized by recent mortality transitioned to a different state, likely dominated by algal overgrowth.

During bleaching (March to late July 2024)

The first four principal components (PCs) of survey data from March to late July 2024 explained 80.6% of the total variance, with PC1 (30.9%) and PC2 (21.6%) as the dominant axes. PC1 primarily reflected a depth-related gradient, with deeper sites associated with Pavona spp., Psammacora spp., and Other coral (positive loadings), contrasting with sites experiencing high bleaching and dead coral. Ancla, Cueva, and Barco Profundo aligned with deeper reefs, while San Josecito exhibited extreme bleaching and coral mortality. For PC2, Pocillopora spp. and dead coral were positively associated, and Porites spp., turf, macroalgae, and CCA negatively associated. San Josecito had the highest positive score, while Barco Somero had the most negative, followed by Este Intermedio and Tina, indicating high turf cover at the latter sites. PC3 (16.9%) explained site variance based Caulerpa (positive), contrasting with coral cover and CCA.

After bleaching (August to December 2024)

The PCA for August to December 2024 distinguished sites with high bleaching prevalence, turf algae, and coral loss from deeper sites, explaining 86% of the total variance among the first four PCs (PC1 32.6%, PC2 24.2%, PC3 15.3%, and PC4 13.9%). PC1 separated deeper sites (0.42), associated with higher coral cover, from sites with high bleaching prevalence (−0.32) and declining coral cover (−0.23). Ancla and Cueva clustered with deeper reefs, while San Josecito and Tina exhibited high bleaching and loss of coral cover. PC2 differentiated sites based on coral cover (−0.36) and CCA (−0.24), contrasting with turf algae (+0.40) and dead coral (−0.35). PC3 (15.3%) captured site differences in coral cover (−0.37) versus increasing turf algae dominance. PC4 (13.9%) reflected variation in macroalgae (−0.48) and total Other coral (−0.49), distinguishing sites with increased macroalgal presence from those with remaining coral dominance.

Latest (January and February 2025)

The PCA for the latest data demonstrated a strong gradient distinguishing coral-dominated and algal-dominated reef sites. The first four PCs explained 86.74% of the total variance, with PC1 (34.5%) and PC2 (22.1%) as the dominant axes. PC1 primarily reflected a coral-to-algae gradient, with coral cover (−0.426) negatively associated and turf algae (0.347), CCA (0.279), and depth (0.366) positively associated. Barco Somero and Esquina aligned with coral-dominated reefs, while Barco Profundo and Cueva exhibited higher algal cover. PC2 (22.1%) captured variation in coral cover (0.279) and depth (0.304), further distinguishing site-level differences. Ancla and San Josecito showed strong divergence, between axes, but bleaching was not a major structuring variable in this period. San Josecito exhibited low coral cover (−0.43), high turf algae (0.35), consistent with post-bleaching impacts. Barco Somero and Esquina had relatively higher coral cover and Barco Profundo and Cueva were positioned towards algal dominance.

Comparing PCA outcomes across time periods, there was a shift in the dominant variables associated with structuring community composition. For example, coral-associated variables contributed strongly to PC1 during the bleaching period but were less influential in the most recent PCA, where turf algae had higher loadings. Turf algae became more strongly aligned with sites in ordination space post-bleaching, reflecting a continued shift in the principal structuring variables. Bleached coral initially played a key role in explaining community variation but was less strongly associated with PCA axes in later analyses, suggesting either coral mortality or stabilization of post-disturbance conditions. Pavona spp. showed more stable contributions to PC1 across all time periods, suggesting a consistent relationship with community composition, which may reflect resilience to disturbance and bleaching events.

Ecological recovery feasibility index

The Ecological Recovery Feasibility Index (ERFI), which prioritizes coral cover and ecologically relevant PCA loadings (Fig. 8), revealed significant variation in reef restoration potential across sites. Chorro (21.18) and Esquina (12.89) exhibited the highest feasibility scores, indicating favorable conditions for coral recovery. Este Intermedio (5.49) and Tina (4.56) showed moderate restoration potential, suggesting a balance between coral persistence and algal encroachment. Cueva (3.15) and Barco Somero (3.11) also ranked positively but with lower scores, likely influenced by higher macroalgae and cyanobacteria cover. Conversely, Ancla (−10.41), San Josecito (−14.26), and Barco Profundo (−25.72) ranked lowest, suggesting degraded reef conditions that may hinder self-recovery. As expected given the index construction, sites with higher feasibility scores also exhibited greater coral cover and lower turf and cyanobacteria prevalence, whereas low-ranking sites were characterized by persistent algal dominance.

The 2023–24 El Niño surpassed historical SST at Isla del Caño

The 2023–24 El Niño marked the most extreme thermal event at Isla del Caño for 40 years, with SSTs reaching unprecedented levels, a triple-peaked DHW anomaly, and the lowest frequency of cool days. This surpasses prior El Niño events documented to affect Costa Rica’s ETP coral reefs since the 1980s (Alvarado et al., 2020; Cortés et al., 2010; Cortés & Jiménez, 2003; Guzmán et al., 1987; Jiménez et al., 2001) and signals a loss of climate variability that may have previously buffered reefs from thermal stress (Pearce-Kelly et al., 2025; Trew & Maclean, 2021). While extreme cooling events can cause coral die-off, such as in Costa Rica’s north Pacific (Palmer et al., 2022), cool periods reduce bleaching impact and aid reef recovery (Green et al., 2019). With the exception of the 1983 El Niño, Isla del Caño’s coral communities have previously escaped devastation from thermal anomalies (Alvarado et al., 2020; Jiménez et al., 2001), leading to the suggestion of local resilience (Alvarado et al., 2020; Friedlander et al., 2022; Guzmán & Cortés, 2007; Romero-Torres et al., 2020). However, retrospective analysis shows that locally the thermal events were not prolonged or extreme and were accompanied by cool days. Over the 40-year period from 1985 to 2025, SSTs at Isla del Caño have risen by 0.92 °C, slightly above the global average of 0.7 °C (Samset et al., 2023), and are consistent with global climate change and ENSO-induced variability (Intergovernmental Panel on Climate Change (IPCC), 2021). This change in thermal regime, with a warming trend decoupled from cool periods, suggests that thermal stress is likely to exceed the pace of natural coral reef recovery at Isla del Caño (Guzmán & Cortés, 2007), increasing the risk of a persistent shift to algal dominated reefs (Nyström et al., 2012; Pearce-Kelly et al., 2025).

Coral to algal dominance

The 2023–24 El Niño triggered a rapid and severe shift from coral to algal dominance at Isla del Caño and the Osa Peninsula, signaling widespread reef degradation (Armstrong McKay et al., 2022; Bland et al., 2018; Darling et al., 2019; Nyström et al., 2012), a likely breach of a tipping point threshold (Pearce-Kelly et al., 2025) and the potential onset of an alternative stable state (Bruno et al., 2009; Crisp, Tebbett & Bellwood, 2022; Folke et al., 2004). Coral cover declined by over 40% and algal cover increased by more than 70%, though as reef sites were surveyed from the tail-end of the thermal anomaly (Reimer et al., 2024), it is likely that the coral loss documented is a gross underestimate. Furthermore, six of the nine sites had approximately 10% coral cover or less in the most recent benthic surveys (January and February 2025) commensurate with levels reported for ecosystem collapse, which is defined as the endpoint of ecosystem degradation (Bland et al., 2018). While the relatively short-term data does not confirm the stability of an alternative state, the observed ecological patterns are consistent with early indicators or such transitions (Fung, Seymour & Johnson, 2011). For example, the magnitude of change and potential for positive feedback loops, whereby persistent algal dominance suppresses coral growth and recruitment, suggest that recovery trajectories may be impaired (Eddy et al., 2021; Norström et al., 2009). The benthic community composition shift at Isla del Caño suggests a significant digression from its previously reported status as a resilient biodiversity hotspot (Cortés et al., 2010; Friedlander et al., 2022; Glynn et al., 2017a).

Spatial and taxonomic trends of reef degradation

Variability in spatial and taxonomic trends highlights the influence of local conditions in buffering or promoting reef degradation in response to the thermal event (Winslow et al., 2024). The most substantial site-level changes were observed at San Josecito, Ancla, Barco Profundo, and Cueva, where significant coral loss, turf proliferation, and reductions in macroalgae and cyanobacteria occurred. Shallow reefs such as San Josecito and Tina, where Pocillopora spp. dominate (Cortés & Jiménez, 2003), exhibited early signs of degradation, with high bleaching prevalence and coral mortality, consistent with the genus’s thermal sensitivity (Marshall & Baird, 2000; Palmer, Bythell & Willis, 2010). In contrast, deeper reefs such as Ancla, Cueva, and Barco Profundo initially appeared more resistant, likely benefiting from cooler water suggestive of a depth-related bleaching refugia (Muir et al., 2017) and a higher prevalence of more thermally tolerant genera Pavona and Psammocora (Marshall & Baird, 2000; Palmer, 2018b; Palmer, Bythell & Willis, 2010). However, these deeper sites were not immune to degradation. Pavona spp., in particular, exhibited partial mortality, demonstrating breaches of damage thresholds, where accumulated stress exceed repair capacity (Palmer, 2018a). The resultant fragmentation inflated colony counts, creating a misleading perception of resistance. It is likely that the timings of the surveys, conducted after the peak bleaching, missed much of the acute mortality phase, particularly for Pavona clavus, which local observation indicate it had already severely declined. This highlights the need for continuous assessment to accurately track community changes and mortality timing (Obura et al., 2019).

In addition to thermal stress, site-specific degradation was compounded by synergistic local stressors. Cueva and Ancla were likely exposed to intense sedimentation, given their proximity to the November 2024 landslides, which likely contributed to coral declines by blocking light, altering microbial communities, and increasing the immune and energetic demands for sediment clearance (Sheridan et al., 2014). Additionally, nutrient influx from fine sediments likely promoted turf algae proliferation (Fourney & Figueiredo, 2017; Vermeij et al., 2010). These impacts were reflected in the SIMPER analysis, which identified turf algae as the primary contributor to benthic community change at Ancla, Cueva, and Barco Profundo.

Several sites, including Esquina, Este Intermedio, and Chorro, exhibited more gradual transitions, initially maintaining higher coral cover before progressively shifting toward turf dominance. Tina and Barco Somero also showed signs of decline, though degradation at Barco Somero was less severe, with turf algae not reaching the dominance levels observed at San Josecito or Barco Profundo. These patterns suggest that site-specific factors, such as initial community composition or local conditions, may have modulated the pace and extent of degradation (Norström et al., 2009).

Across all sites, only 12 coral taxa were recorded in coral reef surveys (2024 to 2025), a reduction from the 22 species historically reported at Isla del Caño (Cortés & Guzmán, 1998; Guzmán & Cortés, 2001). Porites panamensis was notably absent and may be locally extirpated, while Gardineroseris planulata was rare. Whilst our surveys grouped all Pocillopora spp. together, due to difficulties with identification to species level in the field (Johnston, Forsman & Toonen, 2018), their overall decline points to a loss of diversity that threatens long-term resilience (Folke et al., 2004; Nyström, Folke & Moberg, 2000). However, the persistence of some Pocillopora spp. colonies raises the possibility that these individuals may possess enhanced thermal tolerance, potentially through endosymbiont shuffling or the selective survival of more resilient species or genotypes (Palacio-Castro et al., 2023; Romero-Torres et al., 2020). Further impacting the system, crown-of-thorns starfish (Acanthaster planci) were observed preying on surviving coral, including Porites lobata, reducing reef recovery potential.

The cumulative effects of bleaching, sedimentation, and predation appear to have driven an ecological simplification at Isla del Caño and the north Osa Peninsula. Although bleaching was visually prevalent in 2023, its contribution to post-disturbance dissimilarity was low, and dead coral decreased at most sites, reflecting overgrowth of coral skeletons (Romanó De Orte et al., 2021). This suggests that the most significant coral loss occurred prior to the 2024 and 2025 surveys and that the true extent of mortality, and among-species variation, may be underestimated. The widespread decline in live coral cover reinforces this conclusion and indicates a substantial ecological shift, with many reefs converging towards a simplified, turf-dominated state that appears to be beginning to stabilize.

Ecological recovery feasibility

Composite indices, which combine multiple ecological indicators into a single decision-support tool, are widely used in terrestrial and aquatic restoration ecology to synthesize complex ecological data into decision-support tools for prioritizing interventions and evaluating site-specific recovery potential (Orsi, Geneletti & Newton, 2011; Riedler et al., 2015). These frameworks provide transparent, replicable means to integrate multiple ecological indicators and have been instrumental in informing triage-based management in degraded ecosystems (Smith et al., 2023b). In coral reef restoration, the need for analogous tools has been emphasized to align restoration efforts with ecological baselines and resilience-based objectives (Hein et al., 2021; Shaver et al., 2020; Shaver et al., 2022). Responding to this gap, the ERFI introduced here offers a preliminary framework for evaluating the relative recovery potential of reef communities in Costa Rica. By integrating science-informed ecological trajectories, current benthic composition, and PCA-derived ecological gradients, the index can be used to facilitate triage-based restoration decisions at both local and regional scales (Kleypas et al., 2021; Obura et al., 2019). As with other composite indices, the resolution and accuracy of the ERFI will benefit from ongoing refinement and validation, particularly when applied to other reef systems (Smith et al., 2023b).

The application of the ERFI to Isla del Caño and northern Osa Peninsula reefs reveals substantial site-level variation in recovery potential. Critically, the index highlights sites such as San Josecito, Ancla, and Barco Profundo as having limited feasibility for natural recovery, underscoring the need for reactive restoration at these locations (Hein et al., 2021; Kleypas et al., 2021). In contrast, sites like Chorro and Esquina demonstrated relatively higher recovery scores yet remain potentially vulnerable to accelerating thermal anomalies.

These findings affirm that passive recovery alone is unlikely to restore ecological function at the most degraded reef sites. Proactive, site-specific interventions that are tailored to the ecological conditions and disturbance history of each reef are increasingly essential to sustaining biodiversity and ecosystem services in the face of compounding stressors (Alvarado et al., 2025; Pearce-Kelly et al., 2025).

Management and restoration implications

The ERFI highlights that several reef sites in Isla del Caño and the northern Osa Peninsula, especially San Josecito, Ancla, and Barco Profundo, have crossed ecological thresholds, marked by steep coral loss and turf algae proliferation. These changes exceed permissible variability defined by SINAC’s ecological monitoring protocol, PRONAMEC, which calls for comprehensive impact mitigation but lacks clear ecological response strategies (Sistema Nacional de Áreas de Conservación (SINAC), 2016; Nacional de Áreas de Conservación (SINAC), 2021). The transition to turf dominance, particularly at high-use tourism sites such as Ancla, Cueva, and Barco Profundo, aligns with previous calls for more conservation efforts and sustainable use of the Osa Conservation Area (Alvarado et al., 2015), and underscores the need for urgent site-specific restoration and enhanced management (BIOMARCC-SINAC-GIZ, 2016; Naranjo-Arriola, 2021).

Coral gardening remains the most widely implemented restoration method globally, showing moderate success in boosting local coral cover, especially for fast-growing genera like Pocillopora and Acropora (Boström-Einarsson et al., 2020; Shaver et al., 2022; Suggett & Van Oppen, 2022). In Costa Rica, low-cost methods such as floating and benthic nurseries have enabled local-scale outplanting of both branching and massive corals, with promising growth (∼6–9 cm/year) and survival (60–90%) rates (Alvarado et al., 2025). However, these gains are fragile without addressing underlying stressors like turf algae overgrowth and thermal stress.

To stabilize benthic trajectories, ecological management must address both biotic and abiotic stressors. While manual macroalgae removal has proven effective in some contexts (Smith et al., 2023a), Isla del Caño’s primary issue is turf algae. Small-scale turf removal around coral colonies can improve coral health and offers an entry point for community-based stewardship (Cetz-Navarro et al., 2013). Meanwhile, strengthening top-down control through herbivore protection, especially of key grazing fish and invertebrates, can suppress algal dominance and support recovery (Alvarado, Cortés & Reyes-Bonilla, 2015; Mumby, Hastings & Edwards, 2007) though the negative impacts of bioerosion should be considered (Alvarado et al., 2016; Alvarado, Cortés & Reyes-Bonilla, 2015). Reinforcing the no-take status of Isla del Caño (ACOSA-TNC-UCI-ELAP, 2008) and expanding herbivore biomass monitoring should be prioritized. In parallel, addressing land-based sediment and nutrient inputs remains critical (Alvarado, Fernández & Cortés, 2009), especially amid ongoing coastal development (Hall, 2024).

In Costa Rica, where restoration is led by NGOs and community-based organizations within the SINAC framework (Alvarado et al., 2025; SINAC-GIZ, 2020), low-tech, scalable methods remain the most feasible. However, long-term reef resilience depends not only on managing present threats but also on maintaining key ecosystem processes, such as coral reproduction and larval dispersal (Kleypas et al., 2021; Romero-Torres et al., 2018; Shaver et al., 2022). Further investigation into these fundamental processes is needed within the region (Glynn et al., 2017b; Lequeux et al., 2018; Santiago-Valentín et al., 2019) to enable a strategic long-term restoration plan to be developed in line with Costa Rica’s national restoration protocol (SINAC-GIZ, 2020). Experimental strategies, such as larval enhancement, cryopreservation, microbiome manipulation, and assisted evolution, hold promise (Barshis et al., 2013; Bay et al., 2023; Suggett et al., 2024), though are currently limited to pilot stages in high-income settings due to cost and infrastructure demands (Bayraktarov et al., 2016; McLeod et al., 2022; Schmidt-Roach et al., 2025), with limited field-scale validation (Pearce-Kelly et al., 2025). Immune profiling, as proposed under the damage threshold hypothesis of coral susceptibility (Palmer, 2018a), offers a practical and comparatively low-cost strategy for understanding coral health dynamics and informing selective propagation. Rather than aiming to manipulate immunity, this approach focuses on identifying coral phenotypes with inherently higher baseline immune responses, such as antioxidants and phenoloxidase activity (Palmer, 2018b; Palmer, Mydlarz & Willis, 2008; Palmer, Bythell & Willis, 2010), which can be prioritized for nursery rearing and outplanting. This aligns with resilience-based restoration goals (Shaver et al., 2022; SINAC-GIZ, 2020) by enhancing the likelihood of survival under recurrent stress events, and may be particularly valuable in guiding cost-effective restoration in resource-limited settings (Klepac et al., 2024).

Reef recovery depends on coordinated governance and reducing local anthropogenic pressures, including from tourism (Alvarado et al., 2015). Isla del Caño, though a biological reserve (Sistema Nacional de Áreas de Conservación (SINAC), 2009), supports intensive marine tourism that benefits local economies and conservation (BIOMARCC-SINAC-GIZ, 2016; Naranjo-Arriola, 2021). Although reduced visitor numbers during COVID-19 did not influence the abundance of predatory fish (Valverde, Cambra & Espinoza, 2025), tourism visitation has led to coral breakage and potentially other negative impacts, such as from sunscreen-derived toxins (He et al., 2019; Lamb et al., 2014), at sites including Ancla, Barco Profundo, Cueva, and Esquina (Naranjo-Arriola, 2021). Management responses, such as rotating site access, redistributing visitor pressure, setting carrying capacities, and providing targeted guide and tourist training, should be urgently reviewed (ACOSA-TNC-UCI-ELAP, 2008; Alvarado et al., 2015). With reef health underpinning regional tourism and conservation funding, ongoing collaboration between SINAC, NGOs, and local tourism operators is essential to align reef use with long-term ecosystem resilience (Alvarado et al., 2025; Kleypas et al., 2021).