Spatial and temporal variability of Honduran benthic communities and reef functionality throughout the fourth global coral bleaching event

Mass coral bleaching events are intensifying under climate change, yet their ecological consequences remain highly context dependent. We quantified shifts in benthic community structure, coral diversity and reef functionality during the 2023–2025 Fourth Global Mass Coral Bleaching Event (GCBE) in Tela Bay and Roatán, Honduras, two traditionally high-cover reef systems within the Mesoamerican Reef (MAR). Following three years of elevated heat stress (maximum DHW of 18 °C-weeks and 12 °C-weeks) both locations exhibited marked declines up to 60% in live coral cover, accompanied by increases of up to 58% in fleshy macroalgae cover, indicating shifts toward more algal-dominated benthic states. These shifts were accompanied by declines in coral diversity and reef functional index (RFI), suggesting reduced contribution of structurally important coral assemblages to key ecological processes. However, benthic community responses to bleaching were not uniform across sites. Increased multivariate dispersion following the event revealed substantial site-level heterogeneity, particularly in Tela Bay, where some sites retained greater coral cover and higher post-disturbance diversity indices, suggesting localized differences in resistance or early recovery potential. Coral community shifts were also strongly associated with declines in historically abundant taxa such as Agaricia tenuifolia, indicating that prolonged thermal stress can reduce the functional contribution of common and previously persistent coral taxa, not only structurally complex reef-building species. Our results show that even reef locations with high coral cover remain vulnerable to bleaching-driven reorganization of benthic communities and associated functional declines. Nevertheless, strong spatial variability among sites suggests that local environmental conditions and early-bleaching assemblage structure may mediate severity and recovery potential. These findings highlight the need to prioritize site-specific monitoring and management strategies in the MAR that identify and protect reefs with greater persistence potential under accelerating climate stress.

Coral bleaching represents one of the most immediate and visible symptoms of climate-induced stress in marine ecosystems (Pratchett et al., 2008). It is primarily driven by prolonged periods of anomalously high ocean temperatures, often associated with marine heatwaves (MHWs) (Hobday et al., 2016). Record high water temperatures and associated global bleaching events are becoming increasingly more frequent, leaving limited time to adequately document and assess their ecological impacts (Frölicher et al., 2018; Hughes et al., 2018; Muñiz-Castillo et al., 2019). The synergistic effects of global and local pressures have intensified bleaching events since 2003, culminating in the most recent 2023–2025 Fourth Global Coral Bleaching Event (GCBE) (Spady et al., 2026). Influenced by the 2023 El Niño event, approximately 77.2% of the world’s reef locations experienced bleaching-level thermal stress (≥4 Degree Heating Weeks, DHW), nearly four times the extent documented during the first global bleaching event in 1998 (Spady et al., 2026). Given that coral reefs are among the most productive and biodiverse ecosystems on Earth, supporting extensive biodiversity and critical ecological processes, these trends raise major concerns regarding the future condition and persistence of coral ecosystems globally (Stuart-Smith et al., 2018).

When MHWs persist for extended periods, corals may lose their symbiotic photosynthetic microalgae (Symbiodiniaceae), leading to the loss of pigment (Glynn, 1996). Bleached corals exhibit reduced energy reserves, growth, calcification rates, and reproductive potential (Leinbach et al., 2021) compromising reef functionality, including the ability of corals to maintain reef productivity through structural complexity and to support biodiversity (Graham and Nash, 2013; Bellwood et al., 2019; Brandl et al., 2019). Sustained thermal stress drives coral bleaching and mortality, reducing net calcium carbonate production and vertical accretion, and promoting community shifts toward opportunistic, non-framework-building corals. Thereby reef 3D structure is simplified, weakening key ecosystem services such as coastal protection, fisheries and ecotourism (Reaser et al., 2000; Magel et al., 2019; Eddy et al., 2021; Álvarez-Filip et al., 2026).

Furthermore, these structural losses can drive phase shifts from coral to macroalgae-dominated states (McManus and Polsenberg, 2004). Reduced structural complexity limits shelter for juvenile fish and alters grazing dynamics, weakening top-down control of algal-growth (Agudo-Adriani et al., 2016; Richardson et al., 2020). Additionally, increased macroalgal cover and limited substrate availability can alter competitive interactions among benthic communities, reducing coral larval settlement and recovery potential and reinforcing feedback loops that can stabilize algal-dominated states (van de Leemput et al., 2016; Ford et al., 2024). While measures such as live coral cover and diversity indices can describe changes in community structure, they do not fully capture the consequences of those changes for ecological functioning. Because coral taxa differ in their contribution to reef accretion, structural complexity, and habitat provision, evaluating reef functionality can provide a more ecologically meaningful assessment of disturbance impacts (González-Barrios and Álvarez-Filip, 2018; Estrada-Saldívar et al., 2019).

The Mesoamerican Reef (MAR), the largest barrier reef system in the Atlantic Ocean, extends along the coastlines of Mexico, Belize, Guatemala and Honduras, and supports exceptionally high biodiversity and ecosystem connectivity (Gress et al., 2019). This study focused on two ecologically important reef systems within Honduras, Capiro Bank in Tela Bay and Cordelia Banks, Roatán. These reefs are considered among the healthiest within the MAR, with reported live coral cover amongst the highest in the region (McField et al., 2022). Beyond high coral cover, these locations have also been reported as areas of high reef functionality, supporting ecological functioning and ecosystem services (González-Barrios and Álvarez-Filip, 2018). Honduran reefs function as important reservoirs of biodiversity, supporting over 500 fish species, 73 coral species, and diverse sponge and hydrozoan assemblages (SERNA, 2019). These reefs support ecological and commercially important species, provide coral larvae and donor sites for restoration initiatives, and sustain critical ecosystem services through fisheries, tourism, recreation, and coastal protection (Medina, 2011; Canty et al., 2021; Trau Rodas, 2024). Tourism and recreation generate an estimated USD 162.63 million annually, while fisheries contribute approximately USD 99.55 million annually, representing the highest economic values within the MAR (Ruiz De Gauna et al., 2021). However, despite their ecological and economic significance, the MAR is increasingly threatened by recurrent marine heatwaves, coral disease outbreaks, declining water quality, and local anthropogenic stressors that collectively reduce reef resilience and recovery capacity (Helm, 2016; Tuholske, 2016; Muñiz-Castillo et al., 2019; McField et al., 2022; Truc et al., 2023).

While several studies in the Caribbean have quantified effects from the Fourth GCBE (Goreau and Hayes, 2024; Birkart and Álvarez-Filip, 2025; Doherty et al., 2025; Manzello et al., 2025; Green et al., 2026), there remains limited understanding of how bleaching impacts vary spatially across reef systems within the region. Addressing this gap is essential for identifying patterns of reef persistence, evaluating functional consequences of bleaching and establishing baselines for assessing long-term reef recovery. In this study, we used benthic and coral percent cover data to investigate whether the 2023–2025 mass bleaching event resulted in spatially variable benthic community and phase-shift responses across two reef systems in Honduras: Tela Bay and Cordelia Banks, Roatán. Our main objectives were to 1) assess temporal variation in benthic assemblages following bleaching, 2) determine if each location responded differently to bleaching, and 3) evaluate changes in coral community structure, diversity, and reef functionality.

Study area

Surveys were conducted within two marine protected areas (MPAs) located along the Honduran coastline of the Mesoamerican Reef System (MAR). We focused on two locations to assess spatial and temporal variations in response to prolonged thermal stress. Roatán (R), part of the Bay Islands and located approximately 55 km offshore from northern Honduras, encompasses the Cordelia Banks reef complex (17 km²), located 1.6 km southwest of the island (Medina, 2011) (Figure 1). Cordelia Banks lies within the Bay Islands Marine National Park and is recognized as an Important Wildlife Site due to its ecological significance for coral reef biodiversity. The monitored sites were located within reef conservation zones where extractive activities are prohibited (ICF, 2025). Tela Bay (T) is situated in the northern mainland coast within the Tela Bay Marine Wildlife Refuge, a marine protected area established to conserve coral reefs, seagrass beds, and fisheries resources. All surveyed sites in Tela Bay were located within the Capiro Reef, which are designated fisheries recovery zones. This zone was established to protect critical marine habitats and species and restrict most extractive activities, allowing only limited hook-and-line fishing under specific management regulations (Rivera, 2018; Anariba et al., 2021) (Figure 1).

To provide ecological context for the study locations, we summarized long-term benthic monitoring data from the AGRRA benthic monitoring database, collected by Healthy Reefs for Healthy People and partners between 2006 and 2023. These summaries were used descriptively to characterize historical reef condition and coral community composition in Tela Bay and Cordelia Banks prior to the 2023–2025 bleaching event (McField et al., 2024; AGRRA, 2025). According to this dataset, the two regions have historically exhibited contrasting benthic assemblages. Reefs in the Cordelia Banks, Roatán, maintained consistently high live coral cover (46–73%), relatively low to moderate fleshy macroalgal cover (3–18%), and stable crustose coralline algae cover (7–12%). In contrast, Tela Bay reefs showed greater temporal variability, with live coral cover declining from 69% in 2011 to 17% in 2023, accompanied by an increase in fleshy macroalgae from 1% to 29%, while crustose coralline algae cover decreased from 17% to approximately 7% over the same period. These long-term patterns suggest that, prior to the bleaching event, Cordelia Banks remained largely coral-dominated and relatively stable, whereas Tela Bay reefs had already shown signs of increasing macroalgal dominance. The two study regions also differed in species composition. Coral communities at Cordelia Banks were primarily characterized by Agaricia tenuifolia, Porites porites, and Acropora cervicornis, whereas coral communities in Tela Bay were characterized by Agaricia tenuifolia, Porites astreoides, and Agaricia agaricites (McField et al., 2024; AGRRA, 2025). These taxa represent dominant coral species recorded within the study areas prior to the recent bleaching events and provide important context for interpreting subsequent changes in coral community structure and reef functionality.

SST and DHW data were obtained from CoralTemp (NOAA Coral Reef Watch, 2019) via ERDDAP (Liu et al., 2017). Average values in each of the regions considered in this study revealed that over the course of the Fourth GCBE, three distinct heat stress events occurred in the region during 2023, 2024, and 2025 (Figure 2). For the purposes of this study, these events are defined as three discrete analytical periods (T1, T2, and T3) which affected the two study locations during different months. A total of 9 sites from Tela Bay were surveyed in May 2023, September 2023, and May 2025 while 3 sites from Cordelia Banks were surveyed in September 2023, February 2024, and September 2025 (Supplementary Table 1). Survey timing was aligned with the most relevant thermal stress periods for each location, based on satellite-derived temperature anomalies, to account for spatial differences in the timing and intensity of heat stress exposure.

Field surveys

All surveys followed standardized long-term monitoring protocols by Atlantic Gulf Rapid Reef Assessment (AGRRA) (Roth et al., 2024) in collaboration with local comanagers, and Healthy Reef Initiative partners (HRI). Benthic cover was quantified by SCUBA divers using point-intercept sampling along a 10 m transect with observations recorded at 10 cm intervals for a total of 100 points per transect. At each point, trained divers recorded the benthic category directly beneath the transect line using AGRRA Benthos Underwater Codes and CORAL Codes-V5.4 (Roth et al., 2024). When multiple benthic components were present at a point, the two most common were recorded; if no living cover was present, the underlying substratum type was noted. Although surveys were conducted over multiple periods of the bleaching event, bleached corals were not recorded within this survey methodology. Between 18 and 20 transects were conducted at each site across the three sampling periods, each transect was separated by more than 5 m, resulting in a total of 224 transects across 12 sites, at a mean depth of 9.8 m (Supplementary Table 1).

For each transect, we calculated percent cover from point-intercept counts of major benthic categories, including total live coral, turf algae, fleshy macroalgae, and crustose coralline algae. Each benthic category was calculated as a proportion by summing the number of points for all its codes and dividing by the total number of available points in the transect length; this proportion was then multiplied by 100 to obtain the percent cover per transect. Transect-level values were then averaged at the site level to obtain mean percent cover for each benthic category per period. Similarly, coral species-level percent cover was calculated per transect and subsequently averaged across transects within each site and sampling period.

Reef functionality index

Reef functionality was assessed before and after the bleaching event across each location using the Reef Functional Index (RFI), which integrates species-specific coral percent cover with morpho-functional traits following González-Barrios and Álvarez-Filip (2018). The RFI incorporates three primary components for each coral species: (1) live coral cover, (2) calcification rate, and (3) structural complexity. Calcification rate (kg CaCO3 m −2 yr −1), rugosity index and height (cm) were obtained from González-Barrios and Álvarez-Filip (2018) to create a species or genus-specific functional coefficient (Fc) (Supplementary Table 2). Live coral cover was directly obtained from benthic surveys and expressed as species or genus-level percent cover for each site and sampling periodWhen percent cover was identified only to the genus level, functional coefficients were averaged across the species within that genus (Supplementary Table 2). The RFI was calculated using: ∑ (L c c i 100) (F c i), where Lcc i is the percentage of live coral cover of species i and Fc i represents the functional coefficient of each species.Overall, RFI values were averaged and analyzed at the site level with higher values indicating dominance of structurally complex, fast-calcifying species. RFI values were fourth root transformed to facilitate visualization within figures (Supplementary Table 3).

Data analysis

Univariate analysis

Univariate analyses were used to assess spatial and temporal variation in benthic variables, including percent live coral cover, turf algae, fleshy macroalgae, and crustose coralline algae (CCA) in addition to coral diversity metrics: species richness, Pielou’s evenness and Shannon diversity (Supplementary Table 4), and reef functional index (RFI). To meet the normality of residuals assumption, turf algae were log10-transformed and Shannon diversity index was transformed using the orderNorm method which transforms the data based off of a rank mapping to the normal distribution. Transformations were chosen as the most appropriate following simulations from the bestNormalize package (Peterson, 2021). Two sites had a species richness of 1 and were therefore excluded from Pielou’s evenness calculations. We used Spearman rank correlation test to examine the relationship between benthic variables (Supplementary Figure 1). Variables were then analyzed using linear mixed-effect models, with location (Tela Bay and Roatán) and period (T1, T2, T3) included as fixed effects and site was included as a random effect. For RFI, coral species richness and Shannon diversity were included as a fixed effect to account for its influence on reef functional structure. Model assumptions were assessed using residual diagnostics. When needed transformations were applied to improve model assumptions. Type III ANOVA was used to evaluate the significance of the fixed effects and their interaction. Differences in response variables were calculated for each Location × Period combination using estimated marginal means (EMMs) derived from the fitted models, followed by pairwise comparisons among sampling periods.

Multivariate analysis

Multivariate analyses were used to assess spatial and temporal variation in benthic community assemblages, coral species composition, and coral functionality. Benthic assemblages were defined using percent cover of major benthic categories, including total live coral cover, fleshy macroalgae, and crustose coralline algae. Non-metric multidimensional scaling (nMDS) based on Bray-Curtis dissimilarities were used to visualize patterns in both benthic and coral community structure. Percent cover data were square root transformed prior to analysis to reduce the influence of dominant taxa and zero inflation (Anderson, 2001). Biological variables and coral species were overlaid on ordinations using the ‘Envfit’ function from the package vegan (Oksanen, 2013).

Differences in spatial and temporal benthic assemblages, coral species compositions, and community weighted coral functionality traits were tested using a three-way permutational multivariate analysis of variance (PERMANOVA) with 9999 permutations using adonis2 in the vegan package (Anderson, 2001). Location and sampling period (T1, T2, T3) were included as fixed effects. The site was treated as nested within location to account for local spatial structure (Green, 1993). Where appropriate, permutations were constrained by site to account for repeated sampling structure using the strata argument in the vegan package. Differences in community composition among sampling periods were assessed using pairwise PERMANOVA based on Bray–Curtis dissimilarities calculated from square-root-transformed community data. Pairwise comparisons were conducted for T1 versus T2, T1 versus T3, and T2 versus T3, with permutations constrained within site to account for the repeated-measures sampling design. Differences in multivariate dispersion among sampling periods were assessed using the function betadisper (PERMDISP) based on the same Bray–Curtis dissimilarities (Anderson et al., 2006). This approach evaluates the homogeneity of group variances by measuring the average distance of each sample to its group centroid in multivariate space. Statistical significance was tested for using permutation tests with 999 permutations (permutest) (Anderson et al., 2006), allowing us to assess whether dispersion differed between periods. Mean distances to group centroids were calculated for each period to quantify the direction and magnitude of changes within location heterogeneity.

SIMPER analyses were conducted to identify the benthic categories and coral species contributing most to dissimilarities between groups (Clarke, 1993) with permutation tests providing an indication of the consistency of each species’ contribution across samples. Species with higher average contributions and low permutation p-values were considered the most reliable contributors to observed compositional differences. All statistical analysis was conducted in R (version 4.4).

Thermal exposure context

Both Tela Bay and Roatán experienced elevated sea surface temperatures and accumulated Degree Heating Weeks during the 2023–2025 marine heatwave and associated bleaching event, with survey periods aligned to early, mid-, and prolonged thermal stress phases (Figures 2A, B). Although the magnitude of exposure differed between locations, the dominant prolonged thermal stress trend was similar across both reef systems. In 2023, DHW at Tela Bay began to increase in early June, reaching over 12 °C-weeks by early September and peaking at 18 °C-weeks in December. Conversely, while heat stress also started in June in Roatán, it remained below 4 °C-weeks for ~3 months until September, before increasing rapidly to 12 °C-weeks in October. In 2024, both regions experienced an even more severe event, reaching 19 °C-weeks in Tela Bay and 14 °C-weeks in Roatán. The 2025 event was milder but still substantial, with Tela Bay peaking at 16 °C-weeks and Roatán at 8 °C-weeks. Together, these three consecutive years of heat stress represent the most significant on record for both regions.

Changes in benthic reef assemblages

Benthic reef assemblages changed markedly across the bleaching event, with temporal variation explaining most of the observed community shift.

Across both Tela Bay and Roatán, live coral cover declined substantially while fleshy macroalgae increased throughout the sampled periods (coral: F 2, 20 = 62.51, p < 0.001, Figure 3, Supplementary Table 5; fleshy macroalgae: F 2, 20 = 65.50, p < 0.001, Figure 3; Supplementary Table 5). Coral cover response also differed by location (F 1, 10 = 6.33, p = 0.031, Supplementary Table 5). Although both locations followed similar trends where coral cover significantly declined from T1 to T2 (Roatán: T1–T2: estimate = 0.411, t 20 = 4.60, p = 0.0005; Tela Bay: T1–T2: estimate = 0.528, t 20 = 10.25, p < 0.0001, Supplementary Table 7), neither location showed a significant decline in coral cover from T2–T3 (Roatán: estimate = 0.033, t 20 = 0.374, p = 0.93; Tela Bay: T1–T2: estimate = 0.077, t 20 = 1.49, p = 0.32, Supplementary Table 7). Fleshy macroalgae at Roatán significantly increased from T1 to T2 but not from T2–T3 (T1–T2: estimate = −0.218, t 20 = −2.98, p = 0.020; T2–T3: estimate = −0.163, t 20 = −2.231, p = 0.09, Supplementary Table 7), whereas at Tela Bay, fleshy macroalgae increased significantly across all periods (T1–T2: estimate = −0.218, t 20 = −5.17, p < 0.001; T1–T3: estimate = −0.585, t 20 = −13.84, p < 0.001; T2–T3: estimate = −0.367, t 20 = −8.67, p < 0.001, Supplementary Table 7). Demonstrated by a significant Location × Period interaction (F 2, 20 = 3.86, p = 0.038, Supplementary 5).

Turf algae also exhibited significant temporal effects (turf algae (log10 transformed): adjusted F 2, 20 = 33.08, p < 0.001, Figure 3, Supplementary 5);. Across both locations, turf algae initially peaked at T2 and then declined again at T3, Roatán indicated significant differences among all periods (T1–T2: estimate = −1.545, t 20 = −3.28, p = 0.010; T1–T3: estimate = 1.278, t 20 = 2.71, p = 0.034; T2–T3: estimate = 2.823, t 20 = 5.99, p < 0.0001, Supplementary 7), whereas at Tela Bay, T1 differed from T2 (estimate = −1.355, t 20 = −4.98, p = 0.0002) and T2 differed from T3 (estimate = 1.522, t 20 = 5.60, p = 0.0001), but T1 and T3 did not differ (p = 0.815, Supplementary Table 7). Contrastingly, crustose coralline algae did not differ significantly among study periods (F 2, 20 = 1.45, p = 0.259) or between locations (F 1, 10 = 0.29, p = 0.603), but there was a significant Location × Period interaction (F 2, 20 = 7.39, p = 0.004, Figure 3, Supplementary Table 5). CCA significantly decreased from T1 to T2 at Roatán (estimate = 0.062, t 20 = 3.09, p = 0.015, Supplementary Table 7), whereas no other within-location temporal comparisons were significant.

Multivariate analyses showed significant differences in benthic assemblage composition across sampling periods and between locations, although period explained substantially more variation than location (PERMANOVA; period: R² = 0.713, p < 0.001; location: R² = 0.023, p < 0.001; Supplementary Table 8). These temporal changes followed a similar pattern across locations, as the interaction between sampling period and location was not significant (R² = 0.025, p = 0.284). Pairwise PERMANOVA comparisons indicated significant differences between all sampling periods, including T1 versus T2 (F 1, 22 = 40.39, R ² = 0.647, p = 0.0005), T1 versus T3 (F 1, 22 = 89.90, R ² = 0.803, p = 0.0005), and T2 versus T3 (F 1, 22 = 13.84, R ² = 0.386, p = 0.0005, Supplementary Table 9). Multivariate dispersion differed significantly among sampling periods for the benthic assemblage (PERMDISP; F 2, 33 = 4.88, p = 0.020), with mean dispersion increasing from 0.068 at T1 to 0.138 at T2 before declining to 0.105 at T3. Thus, among-site heterogeneity increased following the disturbance but partially declined by T3 while remaining elevated relative to T1.

The nMDS ordination revealed clear differences among surveyed periods throughout the bleaching event. However, three sites in Tela Bay surveyed during T2 clustered more closely with T1 assemblages than with other T2 or T3 samples, suggesting some similarity to pre-disturbance community structure (Figure 4). SIMPER analysis identified live coral cover as the main contribution to dissimilarity among communities between T1 to T2, accounting for 30.3% of total variation, followed by fleshy macroalgae (12.5%) and turf algae (11.2%). Between T2 and T3, shifts were largely driven by continued increases in fleshy macroalgae (20.4%) and declines in turf algae (14.2%) (Figure 5, Supplementary Table 10).,,,

Coral community structure and diversity

Pielou’s evenness and Shannon Diversity varied across sampling periods, independent of location (Pielou’s Evenness: F 2, 19. 6 = 18.14, p < 0.001, Figure 6; Shannon Diversity Index (orderNorm transformed): F 2, 20 = 16.02, p < 0.001, Figure 7; Supplementary Table 5). Both diversity metrics exhibited similar patterns across both study locations where they measured significantly higher at T1 than at T2 and T3, while T2 and T3 did not differ at either location (Supplementary Table 7). In contrast, species richness differed significantly among sampling periods (F 2, 20 = 16.32, p < 0.001, Figure 8), but not between locations (F 1, 10 = 1.55, p = 0.242), and the Location × Period interaction was not significant (F 2, 20 = 2.74, p = 0.089, Supplementary Table 5). Pairwise comparisons showed that richness declined significantly from T1 to T2 and T3 at Tela Bay (T1–T2: estimate = 5.78, t 20 = 6.02, p < 0.0001; T1–T3: estimate = 6.78, t 20 = 7.07, p < 0.0001), whereas no pairwise temporal comparisons were significant at Roatán.

Coral species community composition differed significantly across periods and locations. PERMANOVA indicated that period explained more variation in coral community composition than location (period: R² = 0.238, p< 0.001; location: R² = 0.130, p< 0.001, Supplementary Table 8). These temporal changes followed a similar pattern across locations, as the interaction between time and location was not significant (R² = 0.055, p = 0.197). Pairwise PERMANOVA comparisons indicated significant differences between T1 and T2 (F 1, 22 = 6.00, R ² = 0.214, p < 0.001) and between T1 and T3 (F 1, 22 = 10.24, R ² = 0.318, p < 0.001), whereas coral species composition did not differ significantly between T2 and T3 (F 1, 22 = 0.88, R ² = 0.039, p = 0.257; p > 0.05, Supplementary Table 9). Multivariate dispersion also differed significantly among periods (PERMDISP; F 2, 33 = 18.00, p = 0.001), with mean dispersion increasing from 0.255 at T1 to 0.538 at T2 and remaining similarly elevated at T3 (0.536). Thus, the shift in coral assemblage structure was accompanied by a substantial increase in among-site heterogeneity following the disturbance, which persisted through T3.

The nMDS ordination of coral species composition showed separation between T1 assemblages to those observed in periods T2 and T3, primarily driven by site differences (Figure 9). Between T1 and T2, the largest contributors to community dissimilarity were Agaricia tenuifolia (ATEN) and Agaricia agaricites (AGAR), accounting for 38.9% and 21.2% of total dissimilarity. In contrast, differences between T2 and T3 were primarily driven by Agaricia tenuifolia (ATEN) (22.6%), Porites astreoides (PAST) (19.7%), and Acropora cervicornis (ACER) (9.57%) (Figure 10, Supplementary Table 11).

Reef functional index

Reef Functional Index (RFI) differed significantly among sampling periods, independent of location (F 2, 20 = 76.49, p < 0.001, Figure 11; Supplementary Table 5). RFI decreased significantly from T1 to T2 and T3 at both Roatán (T1–T2: estimate = 0.00103, t 20 = 4.93, p < 0.001; T1–T3: estimate = 0.00121, t 20 = 5.75, p < 0.001) and Tela Bay (T1–T2: estimate = 0.00132, t 20 = 10.92, p < 0.001; T1–T3: estimate = 0.00158, t 20 = 13.00, p < 0.001, Supplementary Table 7), while T2 and T3 did not differ significantly at either location.

A positive relationship was observed between species richness and RFI in Figure 12, but the relationship was a marginal association (F 1, 19. 8 = 3.86, p = 0.064, Supplementary Table 5). However, Shannon diversity index was significantly associated with RFI (F 1, 26. 3 = 40.20, p < 0.001, Supplementary Table 5).

The 2023–2025 Fourth Global Coral Bleaching Event (GCBE) produced substantial changes in benthic community structure across both Tela Bay and Roatán, primarily driven by temporal variation. Both locations showed substantial declines in live coral cover (T: -60.3%, R: -44.0%), diversity (species richness: R: -4.0, T: -6.7; Shannon diversity: R: -0.3, T: -0.2), and Reef Functional Index (R: -0.00121, T: -0.00157), accompanied by increases in fleshy macroalgae (T: +58.6%, R: +38.2%) and an initial increase in turf algae (T: +16.4%, R: +21.5%), which then declined after prolonged exposure to thermal stress (T: -17.2%%, R: -30.2%). These patterns are consistent with recent observations across multiple regions throughout the Caribbean and tropical eastern Pacific where prolonged thermal stress during the Fourth GCBE caused substantial coral mortality ranging from 40–93%, alongside major shifts toward more algal-dominated reef states (López-Pérez et al., 2024; Alvarado et al., 2025; Doherty et al., 2025; Goodbody-Gringley and Chequer, 2025; Liñán‐Cabello et al., 2025; Palmer et al., 2025). The severity of thermal exposure across these regions was also reflected in Degree Heating Week (DHW) accumulation, with maximum values of approximately 12-20 °C-weeks in Panama, ~15 °C-weeks in Veracruz, and > 17 °C-weeks in Little Cayman (Doherty et al., 2025; Glass et al., 2026; Horta-Puga et al., 2026). These values greatly exceeded the 8 °C-weeks threshold associated with severe coral bleaching and mortality, while remaining within the range of 8-19 °C-weeks reported across our study region during the Fourth GCBE, indicating that thermal stress experiences at our study sites were comparable to other regions. However, areas with long-term increases in DHW experience greater bleaching severity, suggesting that reefs increasingly exposed to elevated thermal stress may reduce their capacity to withstand subsequent marine heatwaves (Muñiz-Castillo et al., 2024). Importantly, these changes were not limited to temporary declines in coral cover alone. The concurrent declines in coral diversity and RFI indicate that bleaching also altered coral assemblage structure and reduced the contribution of coral communities to key reef functions. Our findings suggest that the changes observed in Tela Bay and Roatán represent broader consequences of consecutive thermal stress events, not only a single bleaching event.

The severe ecological and structural losses observed in the historically robust reefs of Tela Bay and Roatán are consistent with a broader Western Atlantic trajectory of functional decline and carbonate loss. This shift is characterized by not only a decline in coral cover, but also by reorganization of reef communities. For example, in the Southwestern Atlantic of Brazil, long-term monitoring has revealed a 45% decline in dominant, framework building corals, which are being replaced by a >150% expansion of opportunistic corals (Carlos-Júnior et al., 2026). Furthermore, both Florida and Mexico experienced large functional declines due to >97% mortality of Acropora cervicornis and Acropora palmata (Birkart and Álvarez-Filip, 2025; Manzello et al., 2025). This is particularly concerning where prolonged coral loss and functionality results in reduced carbonate production and a transition toward net-erosional reef states. Region-wide carbonate models across the western Atlantic predict more than 70% of reefs to transition to a net-erosional state by 2040 (Perry et al., 2025). This prediction is consistent with recent assessments in the Caribbean showing that bleaching and disease associated with the 2023–2025 marine heatwave pushed more than 70% of reef sites into net-erosional states, where carbonate loss exceeds reef accretion (Álvarez-Filip et al., 2026). Even historically healthy areas like Bonaire exhibit average coral cover under 10% and nearly 30% of its reef sites are already net eroding (Meesters et al., 2025). Mexico’s Mesoamerican reef displays a similar decline in accretion rate with 38.4% of sites net eroding (Perry et al., 2025). Therefore, based on similar regional trends, declines in coral cover and functionality will likely lead towards decreased carbonate production, indicating a lower ability of natural reef recovery may result.

Although temporal variation accounted for the majority of the observed community change, spatial variation between locations was also statistically significant for live coral cover and RFI, while significant period × location interactions were observed for fleshy macroalgae and CCA throughout the disturbance period. This suggests that while bleaching drove an overall phase-shift in benthic communities, local environmental conditions, such as wave exposure, turbidity, and nutrient input, may have influenced the magnitude and trajectory of these changes and should be investigated in future studies (Wiedenmann et al., 2013; Freitas et al., 2019; Cunio and Mercado-Molina, 2026). As seen in Figure 10, coral assemblages showed increasing spatial variability over time, indicating that sites followed divergent trajectories in response to repeated thermal stress. This pattern suggests that bleaching impacts may be shaped by site-specific local environmental conditions and ecological memory. Ecological memory describes how coral life-history traits, physiological adaptations, and behavioral strategies shaped by past environmental conditions or disturbance histories influence the capacity of ecosystems to respond to future stressors such as thermal stress (Ogle et al., 2015; Hughes et al., 2019; Schweiger et al., 2019). This concept may help explain the contrast between our two focal sites. Tela Bay and Roatán have completely different historical patterns to exposure to heat stress. Corals in Tela Bay had been more exposed to different heat stress events (Muñiz-Castillo et al., 2019; McField et al., 2024), yet Tela Bay showed higher coral cover and richness at the start of our study, opposite to what a purely cumulative damage concept predicts. This fits the ecological memory hypothesis: earlier repeated exposure may increase the resistance to the consequent heat stress events (Hughes et al., 2019; Muñiz-Castillo et al., 2024; González-Barrios and Álvarez-Filip, 2018, 2023). However, loss can also arrive late or because of a compounded effect of multiple stressors, in the Caribbean, much mortality follows disease outbreaks in combination with coral bleaching events (Álvarez-Filip et al., 2026). Disentangling long-term memory from local stressors will require finer-scale thermal reconstructions and the inclusion of the thermal exposure context or history as priority for future works.

Changes in coral species composition were linked to decline in reef functionality. Community dissimilarity was associated with changes in Agaricia tenuifolia, Agaricia agaricites, Porites astreoides and Acropora cervicornis, while RFI declined across sampling periods. Structurally complex, fast-calcifying corals such as Acropora are widely recognized as major contributors to reef accretion and three-dimensional habitat structure (González-Barrios and Álvarez-Filip, 2018). However, our results also show that abundant taxa with lower per-colony structural contributions can strongly influence community-level function when they dominate local assemblages. Importantly, RFI incorporates live coral percent cover into its calculation; meaning declines in RFI likely reflect not only shifts in the functional composition of the coral community but also overall reductions in live coral cover. Consequently, some of the observed decline in reef functionality is attributable to coral loss itself, in addition to changes in the relative abundance of functionally important taxa. Declines in Agaricia spp. were important because these taxa were historically abundant in the study sites and contributed substantially to overall coral cover and RFI (McField et al., 2022). Although Agaricia and Porites are generally considered dominant, fast-growing and prolific recruiting taxa within Caribbean assemblages (Doherty et al., 2025), their decline during this event suggests that even common and previously persistent corals can lose their buffering role under extreme disturbance.

The relationship between coral diversity and reef response to bleaching remains complex and may vary across ecological scales. Previous regional-scale studies have suggested that more diverse assemblages may exhibit increased bleaching severity because they include a greater range of coral taxa, including thermally sensitive species (Muñiz-Castillo et al., 2024). In contrast, our findings suggest that diversity may play a different role at local scales under prolonged stress. Specifically, the positive relationship between species richness and RFI suggests that more species-rich assemblages may partially buffer the functional consequences of bleaching, even if greater diversity does not necessarily reduce bleaching susceptibility itself. As found by Baum et al. (2023), coral community losses following a heatwave depended on pre-heatwave community structure, highlighting how differences in the composition of coral assemblages can shape responses to thermal disturbance. In this context, species-rich assemblages may provide greater functional redundancy, increasing the likelihood that some taxa can persist and continue contributing to key ecosystem processes when other taxa are disproportionately affected by bleaching (Carturan et al., 2022). Because coral thermal tolerance varies with symbiont composition, previous thermal exposure and environmental conditions, coral taxa can differ substantially in their susceptibility to bleaching (Grottoli et al., 2014; Thomas et al., 2018; Denis et al., 2024). For example, in Palau, specific genotypes of Porites cf. lobata show enhanced thermal tolerance across different habitats suggesting that variation in thermal tolerance may be influenced by acclimatization and genetic differences (Rivera et al., 2022). Similarly, changes in coral host-symbiont associations have been linked to increased thermal tolerance in Eastern Pacific corals (Turnham et al., 2023). Together, these findings suggest that acclimatization, genetic variation, or host–microbiome interactions may enhance coral thermal tolerance and contribute to future recovery. Thus, greater species richness may help mediate the effects of bleaching by increasing the likelihood that functionally important taxa persist and maintain key ecosystem processes such as calcium carbonate secretion and reef three-dimensionality despite changes in community composition (Briggs et al., 2020). Consistent with this interpretation, Carturan et al. (2022) found that functional richness was positively associated with coral resilience through its effects on coral cover and rugosity 10 years after disturbance. Similar patterns have also been documented in terrestrial ecosystems, where greater species richness can increase resistance to drought by reducing the functional consequences of species-specific sensitivity to environmental stress (Bai and Tang, 2024). However, the absence of comparable species-specific bleaching and mortality data limits our ability to determine whether the observed relationship reflects differences in bleaching susceptibility, post-bleaching mortality, functional redundancy or a combination of these mechanisms. Consequently, although diversity may influence bleaching susceptibility at regional-scales, local assemblages characterized by higher species richness prior to disturbance may enhance the likelihood that coral communities retain their capacity to provide necessary ecological functions post disturbance by retaining functional redundancy.

The increase in fleshy macroalgae and the variable responses of turf algae and crustose coralline algae (CCA) observed in both regions will have important implications for recovery trajectories. Coral loss can open space for rapid algal colonization, particularly where herbivory is insufficient or nutrient inputs favor algal growth (McManus and Polsenberg, 2004; Graham et al., 2013). Once macroalgae become established, they can suppress coral recruitment, intensify competition for space, and reinforce feedbacks that slows or prevents recovery toward coral-dominated states (Mumby et al., 2007; Mumby, 2009; van de Leemput et al., 2016). Similarly, the increasing prevalence of cyanobacterial mats across reef benthic communities indicates the need to consider this group in future studies, particularly their potential role in influencing reef recovery following disturbances such as bleaching (de Bakker et al., 2017; Ribeiro et al., 2018). At the same time, not all benthic groups have the same implications for recovery. CCA and sparse turf algae contribute to carbonate production or substrate stabilization and recruitment success (Vermeij et al., 2011; Lang et al., 2016; Jorissen et al., 2021), whereas dense turf and fleshy macroalgae can inhibit coral recruitment and favor algal persistence (Kuffner et al., 2006; Arnold et al., 2010). Furthermore, we observed significant period x location variation in CCA cover, with an initial increase during T2 at Tela Bay while declining in Roatán, before both locations returned to pre-disturbance levels under prolonged thermal stress. These results suggest that variations in site responses to thermal stress highlight the importance of long-term monitoring efforts following bleaching impacts to better understand site-specific recovery trajectories.

Overall, our findings indicate that even high-cover reef systems in the Honduran Caribbean are vulnerable to bleaching-driven reorganization of benthic communities and associated functional decline under prolonged thermal stress. However, the increase in mean dispersion among-sites for both benthic and coral assemblages also suggest that some reefs may have experienced lower impacts than others, highlighting their differences in response to prolonged bleaching impacts. Management should therefore move beyond broad regional assessments and prioritize site-specific monitoring that identifies reefs retaining coral cover, diversity, and functional contribution after a disturbance. Reducing local stressors, particularly poor water quality and nutrient inputs, remains essential because these pressures can compound the effects thermal stress and promote macroalgal expansion (Fung et al., 2011; Graham et al., 2013, 2015; Mumby et al., 2007; Sura et al., 2025). Effective implementation of herbivore protections is also critical to maintaining grazing processes that limit algal dominance and support coral recovery (Hempson et al., 2018; Richardson et al., 2018; Elma et al., 2023). However, the functional role of parrotfishes involves important trade-offs, as grazing that reduces algal competition may occur alongside reef-framework bioerosion (Molina-Hernández and Álvarez-Filip, 2024). The relative contribution of these functions depends on species composition, body size structure, and reef condition, highlighting the importance of considering herbivore functional roles rather than biomass alone when evaluating reef recovery and resilience. Conservation and restoration planning should consider both structurally complex reef-building taxa and abundant coral taxa that currently support local reef function. As marine heatwaves intensify, protecting reefs with greater persistence potential while reducing local stressors will be central to maintaining reef structure, fisheries habitat, coastal protection, and tourism benefits across the MAR.

This study evaluated the responses of benthic and coral communities in Tela Bay and Roatán to the 2023–2025 Fourth GCBE through analyses of coral cover, community diversity, and reef functional index. Our findings show that three consecutive years of elevated thermal stress were associated with substantial declines in live coral cover, reductions in coral diversity, promoted shifts toward more algal-dominated states, and reduced overall reef functionality. These changes indicate that bleaching impacts extended beyond visible coral paling or mortality and affected broader aspects of reef community structure and ecological function.

Importantly, reef responses varied considerably among sites, as indicated by increased multivariate dispersion following the disturbance. This spatial variability suggests that disturbance trajectories were likely influenced by a combination of pre-bleaching community structure mediating effects, site-specific ecological memory and adaptations, and local environmental conditions.

The spatial variability observed in community response highlights the need to move beyond generalized assessments of bleaching impact towards a deeper understanding of what drives resilience and recovery at the local level. Identifying and protecting reefs with greater persistence potential must become a central priority for reef conservation and management, requiring coordinated efforts to integrate site-specific monitoring, reduce local stressors, and direct resources toward the reefs with greatest persistence potential.