Published January 23, 2026 | Version v1
Dataset Open

Beyond bleaching: Collapse of net coral reef carbonate budgets in the Eastern Tropical Pacific after the 4th global coral bleaching event.

Authors/Creators

Description

Description of the data and file structure

Title:** Beyond bleaching: Collapse of net coral reef carbonate budgets in the Eastern Tropical Pacific after the 4th global coral bleaching event.**

Francisco Medellín-Maldonado1,2*, Rebeca Granja-Fernández3, Tania González-Mendoza4, Rafael Andrés Cabral-Tena5, Lorenzo Álvarez-Filip1, Andrés López-Pérez2

Abstract

Marine heatwaves' increasing frequency and severity drive unprecedented ecological changes in coral reefs. Although extensive research has documented mass coral bleaching events and subsequent declines in live coral cover worldwide, the impacts on structural functionality and the disruption of key mechanisms governing reef accretion remain poorly understood. Here, we investigate the effects of the fourth global bleaching event on carbonate production and degradation processes in coral reefs of the southern Eastern Tropical Pacific. Using a combination of field surveys, carbonate production models, and thermal stress analyses, we show that prolonged SST anomalies peaked at 2.24 °C in August 2023, with degree heating weeks reaching a record-breaking 16, far surpassing established bleaching and mortality thresholds. This extreme thermal event led to regional coral mortality exceeding 76%, with some reefs experiencing complete loss of live coral cover. Consequently, these reefs shifted from highly productive carbonate states (4.25 ± 2.13 kg·CaCO₃·m⁻²·yr⁻¹) to net erosional states (-2.67 ± 2.82 kg·CaCO₃·m⁻²·yr⁻¹). Furthermore, we observed significant disruptions in the mechanisms underlying both carbonate production and erosion. Our findings underscore the vulnerability of reef structural functionality to climate change-induced thermal stress and provide crucial insights into the potential status of coral reefs affected by 4th global coral bleaching event.

Repository contains the following:

R script, and raw data

  1. Figures 
    • data.txt
    • code.txt
  2. Statics code
    • code.txt
    • data.txt
  3. Files
    • DATA_RSPB-2025-0482.xlsx
    • CODE_RSPB-2025-0482_reviwers_2.R
    • data_cover_SST_3_c.txt
    • cobertura_huatus_23vs24.txt
    • SST_Huatulco_promedio_2016.txt
    • data_grupos.txt
    • H23vs24.txt

Files and variables

File: DATA_RSPB-2025-0482.xlsx

Description: Data of SST and DHW from 1985 to 2023, Coral cover from 2016 to 2024

**Files: **data_cover_SST_3_c.txt to H23vs24.txt is the data necesar to run de code in R

Variables
  • Data base of surveys and historical monthly sea surface temperature (SST, DHW), records for the study area were obtained from the NOAA Coral Reef Watch Coral Bleaching Thermal Stress dataset (NOAA/CRW), which is based on satellite observations with a 5-km resolution and covers the period from January 1, 1985 to December 31, 2023 (accessed on January 12, 2024).

Code/software

################ CODE RSPB-2025-0482 #############################

##### #####################

############################################33

# FIT MODELS#

library(glmmTMB)

library(ggeffects)

library(DHARMa)

library(ggplot2)

###########################################################

data_cover_SST_3_c

names(data_cover_SST_3_c)

view(data_cover_SST_3_c)

data_cover_SST_3_c$Year<- as.factor(data_cover_SST_3_c$Year)

data_cover_SST_3_c$Reef<- as.factor(data_cover_SST_3_c$Reef)

data_cover_SST_3_c$Seson<- as.factor(data_cover_SST_3_c$Seson)

data_cover_SST_3_c<-within(data_cover_SST_3_c, anidado <- factor(Reef:Year))## C O N C A T E N A R

data_cover_SST_3_c$anidado<- as.factor(data_cover_SST_3_c$anidado)

hist(data_cover_SST_3_c$coral.cover, breaks = 20, main = "Distribución de cobertura coralina")

boxplot(coral.cover ~ Seson, data = data_cover_SST_3_c)

shapiro.test(data_cover_SST_3_c$coral.cover)

boxplot(Max_Av._DHW ~ Reef, data = data_cover_SST_3_c)

levels(data_cover_SST_3_c$Seson)

data_cover_SST_3_c$Seson<- factor(data_cover_SST_3_c$Seson, levels = c ("Before Niño", "After Niño" ))

levels(data_cover_SST_3_c$Reef)

data_cover_SST_3_c$Reef <- factor(data_cover_SST_3_c$Reef, levels = c ("San Agistín", "Riscalillo" ,"Jicaral" , "Dos Hermanas", "La India","Cacaluta", "Maguey","Organo", "Violín" , "La Entrega", "Chahué" ))

##########

#########################################33

data_cover_SST_3_c$prop_cover <- (data_cover_SST_3_c$coral.cover + 0.001) / 100.002

mod_beta <- glmmTMB(prop_cover ~ Reef + Year,

                family = beta_family(link = "logit"), 

                data = data_cover_SST_3_c)

summary(mod_beta)

simulateResiduals(mod_beta, n=250, plot= T) #

Anova(mod_beta)

Reef_unit<- ggpredict(mod_beta, terms = c("Reef", "Year"), type = "fixed") %>%

plot() +

theme_classic(base_size = 25)+

labs(x = "", y= "% Coral cover")+

guides()+ggtitle("") +

theme(

axis.line = element_line(color = "black"),  # Asegurar que los ejes se vean

panel.border = element_rect(color = "black", fill = NA, size = 0.5)  # Añadir borde alrededor

)+

scale_color_manual(values = c("#F38D30", "#B2D8EE","#87B6D9","#5F90BB","#3B6895", "navy", "blue","red"))

Reef_unit #1200*550

#============================================

# Predicciones

pred_beta <- ggpredict(mod_beta, terms = c("Reef", "Year"), type = "fixed")

# Pasar a data.frame y escalar a %

pred_df <- as.data.frame(pred_beta)

pred_df$predicted <- pred_df$predicted * 100

pred_df$conf.low <- pred_df$conf.low * 100

pred_df$conf.high <- pred_df$conf.high * 100

# Gráfico categórico con error bars

Reef_unit2 <- ggplot(pred_df, aes(x = x, y = predicted, color = group)) +

geom_pointrange(aes(ymin = conf.low, ymax = conf.high),

              position = position_dodge(width = 0.6),

              size = 1) +

geom_point(data = data_cover_SST_3_c,

         aes(x = Reef, y = prop_cover * 100, color = Year),

         position = position_jitter(width = 0.2, height = 0),

         size = 2, alpha = 0.6) +

theme_classic(base_size = 20) +

labs(x = "", y = "Coral cover (%)", color = "year") +

theme(

axis.line = element_line(color = "black"),

panel.border = element_rect(color = "black", fill = NA, size = 0.5),

axis.text.x = element_text(angle = 40, hjust = 1),

legend.position = "top",

legend.direction = "horizontal",   # leyenda en horizontal

legend.justification = "center"    # centrada

) +

#guides(colour = "none")+

scale_color_manual(values = c("#F38D30", "#B2D8EE","#87B6D9","#5F90BB","#3B6895", "navy", "blue","red"))

Reef_unit2

#############################################################################################################

Year_plot<- ggpredict(mod_beta, terms = c("Year"), type = "fixed") %>%

plot() +

theme_classic(base_size = 25)+

labs(x = "Year", y= "% Coral cover")+

guides()+ggtitle(" ") +

geom_line(color = "#8B0000", size = 1) +

annotate("rect",xmin=2023,xmax=2024,ymin=0,ymax=Inf, alpha=0.4, fill='#FF4500', )+

annotate("rect",xmin=2015,xmax=2016,ymin=0,ymax=Inf, alpha=0.4, fill='#FF4500', )+

scale_x_continuous(breaks=c(2015,2016,2017,2018,2019,2020,2021,2022,2023,2024))+

scale_y_continuous(breaks=c(0,0.20,0.40,0.60,0.80, 1.10))+

theme(axis.line = element_line(color = "black"), # Asegurar que los ejes se vean

    panel.border = element_rect(color = "black", fill = NA, size = 0.5))

Year_plot #1200*550

#============================================

#############################################################################################################

view(data_cover_SST_3_c)

mod_beta_DHW <- glmmTMB(prop_cover ~ Max_Es._DHW +

                     (1|Reef),

             family = beta_family(link = "logit"),

             data = data_cover_SST_3_c)

mod_beta_DHW1 <- glmmTMB(prop_cover ~ Max_Es._DHW +

                       (1|Reef)+(1|Year),

                     family = beta_family(link = "logit"),

                     data = data_cover_SST_3_c)

mod_beta_DHW3 <- glmmTMB(prop_cover ~ Max_Es._DHW +Reef+

                       (1|Year),

                     family = beta_family(link = "logit"),

                     data = data_cover_SST_3_c)

mod_beta_DHW2 <- glmmTMB(prop_cover ~ Max_Av._DHW +

                    (1|Reef),

              family = beta_family(link = "logit"),

              data = data_cover_SST_3_c)

AIC(mod_beta_DHW, mod_beta_DHW1, mod_beta_DHW2, mod_beta_DHW3)

summary(mod_beta_DHW3)

Anova(mod_beta_DHW1)

simulateResiduals(mod_beta_DHW, n=250, plot= T) #

simulateResiduals(mod_beta_DHW1, n=250, plot= T) # the best

simulateResiduals(mod_beta_DHW3, n=250, plot= T) #

ggpredict(mod_beta_DHW1, terms = c("Max_Es._DHW"), type = "fixed") %>%

plot(data = T, colors = "#8B0000") +

theme_classic(base_size = 25)+

labs(x = "Max annual DHW", y= "% Coral cover")+

guides()+ggtitle(" ") +

theme(

axis.line = element_line(color = "black"),  # Asegurar que los ejes se vean

panel.border = element_rect(color = "black", fill = NA, size = 0.5))

ggpredict(mod_beta_DHW3, terms = c("Max_Es._DHW","Reef"), type = "fixed") %>%

plot() +

theme_classic(base_size = 25)+

labs(x = "Max annual DHW", y= "% Coral cover")+

guides()+ggtitle(" ") +

scale_colour_manual(values = c("#8B2500", "#CD6839", "#FF0000","#FF9932",

                             "#FFCC65", "#FFEE99", "#54FF9F", "#CCFFFF",

                             "#99EDFF","#79CDCD", "#528B8B","#0054FF"),

                  name = "Reef")+

theme(

axis.line = element_line(color = "black"),  # Asegurar que los ejes se vean

panel.border = element_rect(color = "black", fill = NA, size = 0.5))

ggpredict(mod_beta_DHW2, terms = c("Max_Av._DHW"), type = "fixed") %>%

plot(data = T, colors = "#8B0000") +

theme_classic(base_size = 25)+

labs(x = "Max annual DHW", y= "% Coral cover")+

guides()+ggtitle(" ") +

theme(

axis.line = element_line(color = "black"),  # Asegurar que los ejes se vean

panel.border = element_rect(color = "black", fill = NA, size = 0.5))

pred_DHW <- ggpredict(mod_beta_DHW1, terms = "Max_Es._DHW")

# --- 4. Gráfico final ---

ggplot() +

# Datos observados (prop_cover en porcentaje)

geom_point(data = data_cover_SST_3_c,

         aes(x = Max_Es._DHW, y = prop_cover*100, color = Reef),

         alpha = 0.7, size = 3) +

# Curva del modelo

geom_line(data = pred_DHW,

        aes(x = x, y = predicted*100),

        color = "#8B0000", size = 1.2) +

# Banda de confianza

geom_ribbon(data = pred_DHW,

          aes(x = x, ymin = conf.low*100, ymax = conf.high*100),

          fill = "#8B0000", alpha = 0.4) +

geom_hline(yintercept =50, color= "red",linetype="dashed", lwd=1)+ ###media

geom_vline(xintercept =5, color= "black",linetype="dashed", lwd=1)+ ###media

scale_x_continuous(breaks=c(0,2,4,6,8,10,12,14,16,18))+

# Estilo

#guides(colour = "none")+

theme_classic(base_size = 30) +

theme(axis.line = element_line(color = "black"),

panel.border = element_rect(color = "black", fill = NA, size = 1))+

labs(x = "Annual maximum DHW (C°-weeks)",

   y = "% Coral cover") +

scale_colour_manual(values = c("#8B2500", "#CD6839", "#FF0000","#FF9932",

                           "#FFCC65", "#FFEE99", "#54FF9F", "#CCFFFF",

                           "#99EDFF","#79CDCD", "#528B8B","#0054FF"),

                name = "Reef")

###########################################################################################

###################################################################################################################################

data_cober2<-cobertura_huatus_23vs24

names(data_cober2)

view(data_cober2)

data_cober2$Loc<- as.factor(data_cober2$Loc)

data_cober2$Año<- as.factor(data_cober2$Año)

data_cober2$Mes<- as.factor(data_cober2$Mes)

data_cober2$date<-as.Date(data_cober2$date)

levels(data_cober2$Loc)

data_cober2$Loc <- factor(data_cober2$Loc, levels = c ("San Agustín", "Riscalillo","Jicaral", "Dos Hermanas", "La India","Cacaluta", "Maguey","Organo", "Violín", "La Entrega", "Chahué" ))

Cober_hut<- ggplot(data_cober2, aes(x = Año, y = Vivo, group =Loc, color=Loc)) +

geom_point(size=7, pch = 21, stroke = 3) +

geom_line(size=1.2, linewidth = 1.5,position = position_dodge(width = 0.01)) +

scale_colour_manual(values = c("#8B2500", "#CD6839", "#FF0000","#FF9932", "#FFCC65", "#FFEE99", "#54FF9F", "#CCFFFF", "#99EDFF","#79CDCD", "#528B8B","#0054FF"),name = "Reef") +

scale_fill_manual(values =c("#8B2500", "#CD6839", "#FF0000","#FF9932", "#FFCC65", "#FFEE99", "#54FF9F", "#CCFFFF", "#99EDFF","#79CDCD", "#528B8B","#0054FF")) +

scale_linetype_manual(values =c("#8B2500", "#CD6839", "#FF0000","#FF9932", "#FFCC65", "#FFEE99", "#54FF9F", "#CCFFFF", "#99EDFF","#79CDCD", "#528B8B","#0054FF")) +

xlab('')+

scale_y_continuous(breaks=c(0,20,40,60,80,100))+

ylab('% Coral cover')+

theme_classic(base_size = 25) +

theme(

axis.line = element_line(color = "black"),  # Asegura que los ejes se vean

panel.border = element_rect(color = "black", fill = NA, size = 0.5)  # Añade un borde cerrado

) +

labs(title = "")

Cober_hut ## 1150*1650

Cober_hut2<- ggplot(data_cober2, aes(x = Año, y = Vivo, group =Loc, fill=Loc, color=Loc)) +

geom_point(size=7, pch = 21, stroke = 3) +

geom_smooth(method = "lm")+

theme_classic(base_size = 20) +

scale_colour_manual(values = c("#8B2500", "#CD6839", "#FF0000","#FF9932", "#FFCC65", "#FFEE99", "#54FF9F", "#CCFFFF", "#99EDFF","#79CDCD", "#528B8B","#0054FF"),name = "Arrecife") +

scale_fill_manual(values =c("#8B2500", "#CD6839", "#FF0000","#FF9932", "#FFCC65", "#FFEE99", "#54FF9F", "#CCFFFF", "#99EDFF","#79CDCD", "#528B8B","#0054FF")) +

scale_linetype_manual(values =c("#8B2500", "#CD6839", "#FF0000","#FF9932", "#FFCC65", "#FFEE99", "#54FF9F", "#CCFFFF", "#99EDFF","#79CDCD", "#528B8B","#0054FF")) +

xlab('')+

ylab('Cobertura de coral (% diferencia)')

###############################################################################################

mea_dat_SST<-SST_Huatulco_promedio_2016

names(mea_dat_SST)

summary(mea_dat_SST)

str(mea_dat_SST)

mea_dat_SST$day <-as.factor(mea_dat_SST$day)

mea_dat_SST$Month <-as.factor(mea_dat_SST$Month)

mea_dat_SST$Year <-as.factor(mea_dat_SST$Year)

mea_dat_SST$dates<-as.Date(mea_dat_SST$dates)

mea_dat_SST<-within(mea_dat_SST, Time <- factor(day:Month))## C O N C A T E N A R

mea_dat_SST$Month<-as.ts (mea_dat_SST$Month)

SST_huatus_mean <- ggplot(mea_dat_SST, aes(x = Month, y = SST, colour=Year)) +

#facet_wrap(~year, nrow=2) + # a panel for each mountain range

geom_point(alpha = 0.3) +

theme_classic(base_size = 25) +

theme(

axis.line = element_line(color = "black"),  # Asegura que los ejes se vean

panel.border = element_rect(color = "black", fill = NA, size = 0.5)  # Añade un borde cerrado

) +

geom_smooth(method = "gam", aes(color = Year))+

scale_x_continuous(breaks=c(0,2,4,6,8,10,12))+

scale_y_continuous(breaks=c(22,24,26,28,30,32,34))+

xlab("Month of each year") +

ylab("SST")+

labs(title = " ")+

annotate("text", x = 6.5, y = 31.7, label = "max SST 2023", size = 5, color = "#22292F") +

scale_color_manual(values = c("#1E90FF","#1E90FF","#1C86EE","#87CEFA", "#104E8B","#00BFFF","#00BFFF","red","#009ACD","#00688B","#87CEFA","#B0E2FF","#A4D3EE",

                            "#CCD6EB","#D7E8FF","#CEE2FF","#C6DDFD","#BDD7FB","#AACAF5", "#F4F4FD","#F4F3FF","#E9EBF5","#F0EEFF","#EBE9FF","#E6E4FE","#E1DEFC",

                            "#DBD7FB","#FA4E5A","#F53F4D","#F02B3F","#E42237","#D61F33","#C81B2E","#BA182A","#AC1425","#9E1021","#910B1C","#760312","red"))+

scale_fill_manual(values = c("#1E90FF","#1E90FF","#1C86EE","#87CEFA", "#104E8B","#00BFFF","#00BFFF","red","#009ACD","#00688B","#87CEFA","#B0E2FF","#A4D3EE",

                           "#CCD6EB","#D7E8FF","#CEE2FF","#C6DDFD","#BDD7FB","#AACAF5", "#F4F4FD","#F4F3FF","#E9EBF5","#F0EEFF","#EBE9FF","#E6E4FE","#E1DEFC",

                           "#DBD7FB","#FA4E5A","#F53F4D","#F02B3F","#E42237","#D61F33","#C81B2E","#BA182A","#AC1425","#9E1021","#910B1C","#760312","red"))+

geom_hline(yintercept =31.2 , color= "black",linetype="dashed", lwd=1.2)+ ###

guides(colour = "none")

SST_huatus_mean ### 550*1000

DHW_huatus_mean <- ggplot(mea_dat_SST, aes(x = Month, y = DHW, colour=Year)) +

#facet_wrap(~year, nrow=2) +

geom_line(alpha = 2.5) +

theme_classic(base_size = 25) +

theme(

axis.line = element_line(color = "black"),  # Asegura que los ejes se vean

panel.border = element_rect(color = "black", fill = NA, size = 0.5)  # Añade un borde cerrado

) +

#geom_smooth(method = "gam")+

scale_x_continuous(breaks=c(0,2,4,6,8,10,12))+

scale_y_continuous(breaks=c(0,2,4,6,8,10,12,14,16))+

xlab("Month of each year") +

ylab("DHW")+

labs(title = "")+

guides(colour = "none")+ guides(fill = "none")+

geom_hline(yintercept =16.2 , color= "black",linetype="dashed", lwd=1.2)+

#geom_hline(yintercept =4.9 , color= "#760312",linetype="dashed", lwd=1.2)+

#annotate("text", x = 3, y = 5.5, label = "avarege DHW 2023", size = 4, color = "#22292F") +

annotate("text", x = 8, y = 17.2, label = "max DHW 2023", size = 5, color = "#22292F") +

scale_color_manual(values = c("#1E90FF","#1E90FF","#1C86EE","#87CEFA", "#104E8B","#00BFFF","#00BFFF","red","#009ACD","#00688B","#87CEFA","#B0E2FF","#A4D3EE",

                            "#CCD6EB","#D7E8FF","#CEE2FF","#C6DDFD","#BDD7FB","#AACAF5", "#F4F4FD","#F4F3FF","#E9EBF5","#F0EEFF","#EBE9FF","#E6E4FE","#E1DEFC",

                            "#DBD7FB","#FA4E5A","#F53F4D","#F02B3F","#E42237","#D61F33","#C81B2E","#BA182A","#AC1425","#9E1021","#910B1C","#760312","red"))+

scale_fill_manual(values = c("#1E90FF","#1E90FF","#1C86EE","#87CEFA", "#104E8B","#00BFFF","#00BFFF","red","#009ACD","#00688B","#87CEFA","#B0E2FF","#A4D3EE",

                           "#CCD6EB","#D7E8FF","#CEE2FF","#C6DDFD","#BDD7FB","#AACAF5", "#F4F4FD","#F4F3FF","#E9EBF5","#F0EEFF","#EBE9FF","#E6E4FE","#E1DEFC",

                           "#DBD7FB","#FA4E5A","#F53F4D","#F02B3F","#E42237","#D61F33","#C81B2E","#BA182A","#AC1425","#9E1021","#910B1C","#760312","red"))

DHW_huatus_mean ### 550*1000

panel1bc<-SST_huatus_mean/DHW_huatus_mean

panel1bc<-panel1bc+plot_layout(guides = "collect")

panel1bc

str()

###################3333

names (data_grupos)

data_grupos <- data_grupos %>%

mutate(Year = as.factor(Year))

Hut_grupos<-ggplot(data_grupos, aes(width = 0.9, fill=Group, y=as.numeric(Cover), x=Year ))+

geom_bar(position="stack", stat="identity" )+

#scale_x_continuous(breaks = c(2023, 2024)) +

labs(x = "", y = "% Cover groups") +

labs(title = "")+

scale_fill_manual(values = c("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")) +

theme_classic(base_size = 22)+theme(

axis.line = element_line(color = "black"),  # Asegura que los ejes se vean

panel.border = element_rect(color = "black", fill = NA, size = 0.5)  # Añade un borde cerrado

)

zplor<-Hut_grupos+coord_polar(theta = "x",start=0) + theme_bw() +

theme(legend.key.size = unit(1.8, "cm")) + ylim(c(-100,100))+

theme(legend.position = c(0.1, 0.1),axis.text.y = element_text(size = 30,colour="black"),axis.text.x=element_text(size = 22,colour="black"))+

geom_hline(yintercept =0 , color= "black",linetype="dashed", lwd=1)+

geom_hline(yintercept =100 , color= "black",linetype="dashed", lwd=1)+

theme(axis.text.y=element_blank(), axis.ticks.y=element_blank())+

theme(plot.title=element_text(size=22),

    axis.title=element_text(size=22,face="bold"),

    legend.text = element_text(size = 18, face="bold"))

Cober_hut<- ggplot(data_cober2, aes(x = Año, y = Vivo, group =Loc, color=Loc)) +

geom_point(size=7, pch = 21, stroke = 3) +

geom_line(size=1.2, linewidth = 1.5,position = position_dodge(width = 0.01)) +

theme_classic(base_size = 22) +

scale_colour_manual(values = c("#8B2500", "#CD6839", "#FF0000","#FF9932", "#FFCC65", "#FFEE99", "#54FF9F", "#CCFFFF", "#99EDFF","#79CDCD", "#528B8B","#0054FF"),name = "Arrecife") +

scale_fill_manual(values =c("#8B2500", "#CD6839", "#FF0000","#FF9932", "#FFCC65", "#FFEE99", "#54FF9F", "#CCFFFF", "#99EDFF","#79CDCD", "#528B8B","#0054FF")) +

scale_linetype_manual(values =c("#8B2500", "#CD6839", "#FF0000","#FF9932", "#FFCC65", "#FFEE99", "#54FF9F", "#CCFFFF", "#99EDFF","#79CDCD", "#528B8B","#0054FF")) +

xlab('')+

ylab('% Coral cover')+

labs(title = "b)")

Hut_grupos/Cober_hut ## 1150*1650

###############################################################################################

#########################

datH<-H23vs24

names(datH)

View(datH)

datH$Reef<- as.factor(datH$Reef)

datH$Year<- as.factor(datH$Year)

datH$Status<- as.factor(datH$Status)

datH<-within(datH, Status_Year <- factor(Status:Year))

datH<-within(datH, Status_Reef <- factor(Status:Reef))

levels(datH$Reef)

datH$Reef <- factor(datH$Reef, levels = c ("Sn Agustin", "Riscalillo","Jicaral", "Dos Hermanas", "La india " ,"Cacaluta", "Maguey","Organo", "Violin" , "La Entrega", "Chahué" ))

levels(datH$Status)

datH$Status <- factor(datH$Status, levels = c ("Before ENSO", "After ENSO"))

levels(datH$Status_Year)

datH$Status_Reef <- factor(datH$Status_Reef, levels = c ( "Before ENSO:2023", "After ENSO:2024" ))

######################### Net G #######################

mod.HG<- lmer(Net.G~Year+(1|Reef), data=datH, REML = FALSE)

Anova(mod.HG)

lsmeans(mod.HG,pairwise ~ Year)

Gnet<-ggpredict(mod.HG, terms = c("Year" ), type = "fixed") %>%

plot() +

theme_classic(base_size = 25)+

xlab("") +

labs(y=expression(GkgCaCO[3]~m{-2}~year{-1})) +

guides()+ggtitle(" ") +

geom_hline(yintercept =0, color= "red",linetype="dashed", lwd=1)+ ###media

annotate("rect",xmin=2023,xmax=2024,ymin=0,ymax=-Inf, alpha=0.3, fill="#FF4500", )+

theme(

axis.line = element_line(color = "black"),  # Asegurar que los ejes se vean

panel.border = element_rect(color = "black", fill = NA, size = 0.5)  # Añadir borde alrededor

)+

scale_colour_manual(values = c("#8B2500", "#CD6839", "#FF0000","#FF9932", "#FFCC65", "#FFEE99", "#54FF9F", "#CCFFFF", "#99EDFF","#79CDCD", "#528B8B","#0054FF"),name = "Reef") +

scale_fill_manual(values =c("#8B2500", "#CD6839", "#FF0000","#FF9932", "#FFCC65", "#FFEE99", "#54FF9F", "#CCFFFF", "#99EDFF","#79CDCD", "#528B8B","#0054FF")) +

scale_linetype_manual(values =c("#8B2500", "#CD6839", "#FF0000","#FF9932", "#FFCC65", "#FFEE99", "#54FF9F", "#CCFFFF", "#99EDFF","#79CDCD", "#528B8B","#0054FF"))

Gnet

############################################################################

######################### CORAL #######################

mod.HGC<- lmer(Coral~Status+(1|Reef), data=datH, REML = FALSE)

Anova(mod.HG)

lsmeans(mod.HG,pairwise ~ Status)

simulateResiduals(mod.HGC, n=250, plot= T) #

Coral_NCC<-ggpredict(mod.HGC, terms = c("Status" ), type = "fixed") %>%

plot() +

theme_classic(base_size = 25)+

xlab("") +

labs(y=expression(CoralkgCaCO[3]~m{-2}~year{-1})) +

guides()+ggtitle(" ") +

geom_hline(yintercept =0, color= "red",linetype="dashed", lwd=1)+ ###media

theme(

axis.line = element_line(color = "black"),  # Asegurar que los ejes se vean

panel.border = element_rect(color = "black", fill = NA, size = 0.5)  # Añadir borde alrededor

)+

scale_y_continuous(breaks=c(0,5,10))+

scale_colour_manual(values = c("#8B2500", "#CD6839", "#FF0000","#FF9932", "#FFCC65", "#FFEE99", "#54FF9F", "#CCFFFF", "#99EDFF","#79CDCD", "#528B8B","#0054FF"),name = "Reef") +

scale_fill_manual(values =c("#8B2500", "#CD6839", "#FF0000","#FF9932", "#FFCC65", "#FFEE99", "#54FF9F", "#CCFFFF", "#99EDFF","#79CDCD", "#528B8B","#0054FF")) +

scale_linetype_manual(values =c("#8B2500", "#CD6839", "#FF0000","#FF9932", "#FFCC65", "#FFEE99", "#54FF9F", "#CCFFFF", "#99EDFF","#79CDCD", "#528B8B","#0054FF"))

Coral_NCC

# Predicciones

gp <- ggpredict(mod.HGC, terms = "Status")

Coral_NCC <- ggplot(gp, aes(x = x, y = predicted)) +

# 1. Intervalos de confianza como líneas verticales sin "caps"

geom_linerange(aes(ymin = conf.low, ymax = conf.high),

             color = "#99EDFF", linewidth = 1) +

# 2. Conectar las predicciones con línea

geom_line(aes(group = 1), color = "#99EDFF", linewidth = 1) +

# 3. Puntos del modelo

geom_point(color = "#99EDFF", size = 3, alpha = 0.5) +

# 4. Agregar los datos crudos (add.data = TRUE)

geom_jitter(data = datH, aes(x = Status, y = Coral),

          width = 0.1, alpha = 0.3, color = "gray30") +

# Línea horizontal = 0

geom_hline(yintercept = 0, color = "red", linetype = "dashed", linewidth = 1) +

theme_classic(base_size = 25) +

#scale_y_continuous(breaks = c(0,2,6,8,10)) +

labs(x = "",

   y = expression(Coral~kg~CaCO[3]~m^{-2}~year^{-1})) +

theme(axis.line = element_line(color = "black"),

    panel.border = element_rect(color = "black", fill = NA, size = 0.5))

Coral_NCC

############################################################################

######################### CCA #######################

str(datH)

mod.HGCCA<- lmer(CCA~Status+(1|Reef), data=datH, REML = FALSE)

simulateResiduals(mod.HGCCA, n=250, plot= T) #

Anova(mod.HGCCA)

lsmeans(mod.HGCCA,pairwise ~ Status)

mod_beta <- glmmTMB(CCA ~ Status + (1|Reef),

                family = beta_family(link = "logit"), 

                data = datH)

summary(mod_beta)

simulateResiduals(mod_beta, n=250, plot= T) #

CCA_NCC<-ggpredict(mod_beta, terms = c("Status"), type = "fixed") %>%

plot() +

theme_classic(base_size = 25)+

xlab("") +

#scale_y_continuous(breaks=c(0,5,10))+

labs(y=expression(CCAandotherscalcifierskg~CaCO[3]~m{-2}~year{-1})) +

guides()+ggtitle(" ") +

geom_hline(yintercept =0, color= "red",linetype="dashed", lwd=1)+ ###media

annotate("rect",xmin=0.85,xmax=Inf,ymin=0,ymax=-Inf, alpha=0.3, fill='#7a5d7e', )+

theme(

axis.line = element_line(color = "black"),  # Asegurar que los ejes se vean

panel.border = element_rect(color = "black", fill = NA, size = 0.5)  # Añadir borde alrededor

)+

scale_colour_manual(values = c("#8B2500", "#CD6839", "#FF0000","#FF9932", "#FFCC65", "#FFEE99", "#54FF9F", "#CCFFFF", "#99EDFF","#79CDCD", "#528B8B","#0054FF"),name = "Reef") +

scale_fill_manual(values =c("#8B2500", "#CD6839", "#FF0000","#FF9932", "#FFCC65", "#FFEE99", "#54FF9F", "#CCFFFF", "#99EDFF","#79CDCD", "#528B8B","#0054FF")) +

scale_linetype_manual(values =c("#8B2500", "#CD6839", "#FF0000","#FF9932", "#FFCC65", "#FFEE99", "#54FF9F", "#CCFFFF", "#99EDFF","#79CDCD", "#528B8B","#0054FF"))

CCA_NCC

#============================================

#############################################################################################################

##### ######## ###

library(ggeffects)

library(ggplot2)

library(dplyr)

# Predicciones

gp <- ggpredict(mod_beta, terms = "Status")

CCA_NCC <- ggplot(gp, aes(x = x, y = predicted)) +

# 1. Intervalos de confianza como líneas verticales sin "caps"

geom_linerange(aes(ymin = conf.low, ymax = conf.high),

             color = "#FFFF00", linewidth = 1) +

# 2. Conectar las predicciones con línea

geom_line(aes(group = 1), color = "#FFFF00", linewidth = 1) +

# 3. Puntos del modelo

geom_point(color = "#FFFF00", size = 3, alpha = 0.5) +

# 4. Agregar los datos crudos (add.data = TRUE)

geom_jitter(data = datH, aes(x = Status, y = CCA),

          width = 0.1, alpha = 0.3, color = "gray30") +

# Línea horizontal = 0

geom_hline(yintercept = 0, color = "red", linetype = "dashed", linewidth = 1) +

# Sombreado

annotate("rect", xmin = 0.85, xmax = Inf, ymin = 0, ymax = -Inf,

       alpha = 0.3, fill = "#7a5d7e") +

theme_classic(base_size = 25) +

#scale_y_continuous(breaks = c(0,2,6,8,10)) +

labs(

x = "",

y = expression(CCA~and~others~calcifiers~kg~CaCO[3]~m^{-2}~year^{-1})

) +

theme(

axis.line = element_line(color = "black"),

panel.border = element_rect(color = "black", fill = NA, size = 0.5)

)

CCA_NCC

############################################################################

######################### Micro #######################

mod.HGCMicro<- lmer(Micro~Status+(1|Reef), data=datH, REML = FALSE)

Anova(mod.HGCMicro)

lsmeans(mod.HG,pairwise ~ Status)

Micro_NCC<-ggpredict(mod.HGCMicro, terms = c("Status" ), type = "fixed") %>%

plot()+

xlab("") +

labs(y=expression(CoralkgCaCO[3]~m{-2}~año{-1})) +

guides()+ggtitle(" ") +

geom_hline(yintercept =0, color= "red",linetype="dashed", lwd=1)+ ###media

annotate("rect",xmin=0.85,xmax=Inf,ymin=0,ymax=-Inf, alpha=0.3, fill='#7a5d7e', )+

theme(

axis.line = element_line(color = "black"),  # Asegurar que los ejes se vean

panel.border = element_rect(color = "black", fill = NA, size = 0.5)  # Añadir borde alrededor

)+

scale_colour_manual(values = c("#8B2500", "#CD6839", "#FF0000","#FF9932", "#FFCC65", "#FFEE99", "#54FF9F", "#CCFFFF", "#99EDFF","#79CDCD", "#528B8B","#0054FF"),name = "Reef") +

scale_fill_manual(values =c("#8B2500", "#CD6839", "#FF0000","#FF9932", "#FFCC65", "#FFEE99", "#54FF9F", "#CCFFFF", "#99EDFF","#79CDCD", "#528B8B","#0054FF")) +

scale_linetype_manual(values =c("#8B2500", "#CD6839", "#FF0000","#FF9932", "#FFCC65", "#FFEE99", "#54FF9F", "#CCFFFF", "#99EDFF","#79CDCD", "#528B8B","#0054FF"))

Micro_NCC

############################################################################

######################### Macro #######################

mod.HGCMacro<- lmer(Macro~Status+(1|Reef), data=datH, REML = FALSE)

Anova(mod.HGCMacro)

Macro_NCC<-ggpredict(mod.HGCMacro, terms = c("Status" ), type = "fixed") %>%

plot()+

xlab("") +

labs(y=expression(CoralkgCaCO[3]~m{-2}~año{-1})) +

guides()+ggtitle(" ") +

geom_hline(yintercept =0, color= "red",linetype="dashed", lwd=1)+ ###media

annotate("rect",xmin=0.85,xmax=Inf,ymin=0,ymax=-Inf, alpha=0.3, fill='#7a5d7e', )+

theme(

axis.line = element_line(color = "black"),  # Asegurar que los ejes se vean

panel.border = element_rect(color = "black", fill = NA, size = 0.5)  # Añadir borde alrededor

)+

scale_colour_manual(values = c("#8B2500", "#CD6839", "#FF0000","#FF9932", "#FFCC65", "#FFEE99", "#54FF9F", "#CCFFFF", "#99EDFF","#79CDCD", "#528B8B","#0054FF"),name = "Reef") +

scale_fill_manual(values =c("#8B2500", "#CD6839", "#FF0000","#FF9932", "#FFCC65", "#FFEE99", "#54FF9F", "#CCFFFF", "#99EDFF","#79CDCD", "#528B8B","#0054FF")) +

scale_linetype_manual(values =c("#8B2500", "#CD6839", "#FF0000","#FF9932", "#FFCC65", "#FFEE99", "#54FF9F", "#CCFFFF", "#99EDFF","#79CDCD", "#528B8B","#0054FF"))

Macro_NCC

############################################################################

######################### Parrotfish #######################

names(datH)

View(datH)

mod.HGCParrotfish<- lmer(Parrotfish~Status+(1|Reef), data=datH)

summary(mod.HGCParrotfish2)

Anova(mod.HGCParrotfish2)

Parrotfish_NCC<-ggpredict(mod.HGCParrotfish, terms = c("Status" ), type = "fe") %>%

plot(connect.lines = T, add.data = T, line.size=1, dot.size=3, dot.alpha = 0.5, colors = '#993404') +

theme_classic(base_size = 25)+

xlab("") +

labs(y=expression(CoralkgCaCO[3]~m{-2}~año{-1})) +

guides()+ggtitle(" ") +

geom_hline(yintercept =0, color= "red",linetype="dashed", lwd=1)+ ###media

annotate("rect",xmin=0.85,xmax=Inf,ymin=0,ymax=-Inf, alpha=0.3, fill='#7a5d7e', )+

theme(

axis.line = element_line(color = "black"),  # Asegurar que los ejes se vean

panel.border = element_rect(color = "black", fill = NA, size = 0.5)  # Añadir borde alrededor

)+

scale_colour_manual(values = c("#8B2500", "#CD6839", "#FF0000","#FF9932", "#FFCC65", "#FFEE99", "#54FF9F", "#CCFFFF", "#99EDFF","#79CDCD", "#528B8B","#0054FF"),name = "Reef") +

scale_fill_manual(values =c("#8B2500", "#CD6839", "#FF0000","#FF9932", "#FFCC65", "#FFEE99", "#54FF9F", "#CCFFFF", "#99EDFF","#79CDCD", "#528B8B","#0054FF")) +

scale_linetype_manual(values =c("#8B2500", "#CD6839", "#FF0000","#FF9932", "#FFCC65", "#FFEE99", "#54FF9F", "#CCFFFF", "#99EDFF","#79CDCD", "#528B8B","#0054FF"))

Parrotfish_NCC

############################################################################

######################### Urchin #######################

mod.HGCUrchin<- lmer(Urchin~Status+(1|Reef), data=datH)

Anova(mod.HGCUrchin)

Urchin_NCC<-ggpredict(mod.HGCUrchin, terms = c("Status" ), type = "fixed") %>%

plot()+

xlab("") +

labs(y=expression(CoralkgCaCO[3]~m{-2}~año{-1})) +

guides()+ggtitle(" ") +

geom_hline(yintercept =0, color= "red",linetype="dashed", lwd=1)+ ###media

annotate("rect",xmin=0.85,xmax=Inf,ymin=0,ymax=-Inf, alpha=0.3, fill='#7a5d7e', )+

theme(

axis.line = element_line(color = "black"),  # Asegurar que los ejes se vean

panel.border = element_rect(color = "black", fill = NA, size = 0.5)  # Añadir borde alrededor

)+

scale_colour_manual(values = c("#8B2500", "#CD6839", "#FF0000","#FF9932", "#FFCC65", "#FFEE99", "#54FF9F", "#CCFFFF", "#99EDFF","#79CDCD", "#528B8B","#0054FF"),name = "Reef") +

scale_fill_manual(values =c("#8B2500", "#CD6839", "#FF0000","#FF9932", "#FFCC65", "#FFEE99", "#54FF9F", "#CCFFFF", "#99EDFF","#79CDCD", "#528B8B","#0054FF")) +

scale_linetype_manual(values =c("#8B2500", "#CD6839", "#FF0000","#FF9932", "#FFCC65", "#FFEE99", "#54FF9F", "#CCFFFF", "#99EDFF","#79CDCD", "#528B8B","#0054FF"))

Urchin_NCC

#############################################################

Files

cobertura_huatus_23vs24.txt

Files (9.1 MB)

Name Size Download all
md5:a51015cb3bd0b58710a032e5820ef192
3.1 kB Preview Download
md5:492273585377de23cee328c8207f747f
31.4 kB Download
md5:3d585a8008e19a5e47a656e9cb695dc9
7.5 kB Preview Download
md5:fbcb099359fe1f06aef9fdba53d04c0c
301 Bytes Preview Download
md5:ccbb3dfe180ec5d02d1e3c027d896db0
7.9 MB Download
md5:ba5967b7a093de5e98849ea672f47ba4
4.1 kB Preview Download
md5:137e649c28a0dc3698db2720b674d729
34.3 kB Preview Download
md5:a12efe0f885918dae350c7450c664086
198.0 kB Preview Download
md5:99afaec1877438a7848ee3b6ebcb9fed
879.0 kB Preview Download