Data and code associated with Goodbody-Gringley et al 2026 Coral Reefs

Reproducible R Code for Goodbody-Gringley et al. (2025) Brett D. Jameson November 09, 2025 Contents 1 Overview 2 2 Setup 2 3 Data Import and Cleaning 2 3.1 Photophysiology ........................................... 2 3.2 MortalityandSurvivalData..................................... 4 3.3 CoralHealthMetrics......................................... 6 4 Statistical Analyses 7 4.1 Photophysiology ........................................... 7 4.2 MortalityandSurvival........................................ 17 4.3 HostHealthMetrics ......................................... 22 5 Data Visualizations 26 5.1 ManuscriptFigures.......................................... 26 5.2 MortalityandSurvival........................................ 30 6 Supplementary Tables 37 6.1 Table S1 - GLMM Summary for Endosymbiont Density . . . . . . . . . . . . . . . . . . . . . 37 6.2 Table S2 - Beta GLMM Summary for EQY . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 6.3 TableS3-AFTModelSummary.................................. 40 6.4 Table S4 - Mortality GLMM Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 6.5 Table S5 - Predicted Probabilties Table . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 7 Session Info 44 1 1 Overview This document reproduces all analyses and visualizations associated with: Goodbody-Gringley et al. (2025) —Fate tracking of Acropora cervicornis colonies postrestorative out-planting reveals genotype-specific responses to marine heatwaves. All analyses were performed in R (v4.5.1) and compiled into this single R Markdown archive for open-access distribution on Zenodo. The file contains fully reproducible code for data import, cleaning, modeling, statistical testing, and figure generation used in the manuscript. All raw data required to run these analyses are included in the /data directory. Derived outputs (model summaries, figures, tables) are exported to /figures and /tables. 2 Setup This section initializes the analysis environment for reproducibility. It sets global chunk options, installs and loads all required packages, defines consistent plotting themes, and seeds the random number generator for replicable results. All file paths are relative to the project root. 3 Data Import and Cleaning This section imports, inspects, and formats all datasets used in the analysis. Each subsection standardizes key variables (e.g., genotype codes, sampling dates) and performs minimal cleaning to ensure data integrity before modeling. Data are stored in the “/data” directory. 3.1 Photophysiology 3.1.1 Endosymbiont Densities These data quantify symbiont cell densities before and after the 2023 marine heatwave. Replicate measurements per fragment are summarized to means and standard deviations, with genotypes standardized to uppercase codes and domes extracted from sample IDs. # Load data zoox <- read_csv("data/resembid.zoox.density.master.csv") # Extract dome ID from sample_id (A–Z code) and standardize genotype naming zoox <- zoox %>% mutate( dome = str_extract(sample_id, "[A-Za-z]"), genotype = str_to_upper(genotype), time = case_when( sample_date == "2022-06-01" ~"Before", sample_date == "2023-09-04" ~"After", 2 TRUE ~NA_character_ ), time = forcats::fct_relevel(factor(time), "Before","After") ) # Summarize replicate zoox counts into a single mean per fragment zoox_summary <- zoox %>% group_by(sample_date, time, genotype, sample_id, dome) %>% summarise( mean_zoox = mean(zoox_cm2, na.rm = TRUE), sd_zoox = sd(zoox_cm2, na.rm = TRUE), n = n(), .groups = "drop" ) # Prepare the main analysis dataframe zoox_df <- zoox_summary %>% rename(date = sample_date) %>% mutate( genotype = factor(genotype, levels = c("OB","B","G","K","KW","LG","S","Y")), dome = factor(dome)) # Sanity checks stopifnot(all(zoox_df$mean_zoox >0)) # required for Gamma/log models table(zoox_df$time, zoox_df$genotype) # check balance ## ## OB B G K KW LG S Y ## Before 3 3 3 3 3 3 3 3 ## After 9 9 8 8 4 6 8 5 glimpse(zoox_df) ## Rows: 81 ## Columns: 8 ## $ date <date> 2022-06-01, 2022-06-01, 2022-06-01, 2022-06-01, 2022-06-01,~ ## $ time <fct> Before, Before, Before, Before, Before, Before, Before, Befo~ ## $ genotype <fct> B, B, B, G, G, G, K, K, K, KW, KW, KW, LG, LG, LG, OB, OB, O~ ## $ sample_id <chr> "103_b", "106_d", "202_b", "103_a", "106_e", "202_a", "103_c~ ## $ dome <fct> b, d, b, a, e, a, c, c, h, c, e, i, b, a, b, a, b, c, e, a, ~ ## $ mean_zoox <dbl> 840913.6, 969236.5, 628105.1, 443574.8, 630967.3, 655365.1, ~ ## $ sd_zoox <dbl> 128957.57, 50972.71, 52889.10, 58288.23, 40705.66, 50731.55,~ ## $ n <int> 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, ~ 3.1.2 Effective Quantum Yield (FV/Fm) The effective quantum yield (EQY; Fv/Fm) dataset tracks photosynthetic efficiency across genotypes and timepoints. Measurements are averaged per fragment per date, converted to proportions within (0, 1), and formatted for Beta regression analysis. 3 fire <- read_csv("data/resembid.fire.csv")%>% mutate( date = paste(month, year, sep = "_"), date = factor( date, levels = c("May_2022","April_2023","September_2023"), labels = c("May-22","Apr-23","Sep-23") ), dome = as.character(dome), tag = as.character(tag), genotype = factor(genotype), fragment_id = paste(dome, tag, sep = "_") )%>% select(date, dome, tag, fragment_id, genotype, fv.fm) %>% filter(genotype != "C") # Make genotype factor with OB reference fire$genotype <- factor(fire$genotype, levels = c("OB",setdiff(unique(fire$genotype), "OB"))) # Average duplicate fragment measurements per timepoint fire_avg <- fire %>% group_by(fragment_id, date, genotype, dome) %>% summarise(fv.fm = mean(fv.fm), .groups = "drop") # Ensure beta regression values strictly in (0,1) n<- nrow(fire_avg) epsilon <- 1e-4 fire_avg$fv.fm <- (fire_avg$fv.fm *(n -1)+0.5)/n fire_avg$fv.fm <- pmin(pmax(fire_avg$fv.fm, epsilon), 1-epsilon) glimpse(fire_avg) ## Rows: 165 ## Columns: 5 ## $ fragment_id <chr> "A_1", "A_1", "A_1", "A_10", "A_10", "A_10", "A_11", "A_11~ ## $ date <fct> May-22, Apr-23, Sep-23, May-22, Apr-23, Sep-23, Apr-23, Se~ ## $ genotype <fct> G, G, G, KW, KW, KW, OB, OB, G, G, G, B, B, B, K, K, Y, Y,~ ## $ dome <chr> "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A"~ ## $ fv.fm <dbl> 0.1486424, 0.2525091, 0.1580848, 0.1441697, 0.2530061, 0.2~ 3.2 Mortality and Survival Data This dataset integrates fragment survival, mortality percentages, and time-to-event outcomes. survival_df <- read_csv("data/resembid.tle.master.csv")%>% mutate( Month = factor( Month, levels = c("May_2022","April_2023","September_2023"), labels = c("May-22","Apr-23","Sep-23") ), month_date = as.Date(parse_date_time(Month, orders = "my")), health = factor( 4 health, levels = c("healthy","pale","bleached","disease","predation","dead","missing") ), dome = as.character(dome), tag = as.character(tag), genotype = factor(str_to_upper(genotype)), perc_mortality = parse_number(perc_mortality), tle_alive = as.numeric(tle_alive), tle_dead = as.numeric(tle_dead), max_tle = as.numeric(max_tle), coral_id = paste(dome, tag, sep = "_") ) # Restrict mortality data btw (0, 1) mort_df <- survival_df %>% mutate(mort_prop = (perc_mortality /100)*0.99 +0.005)%>% select(-c(tle_alive, tle_dead, max_tle, health)) # Subset dataset for survival time/event analysis time_event_df <- mort_df %>% arrange(coral_id, month_date) %>% group_by(coral_id, genotype, dome) %>% summarise( date_death = if (any(mort_prop == 1,na.rm = TRUE)) min(month_date[mort_prop == 1]) else as.Date(NA), last_obs = max(month_date), .groups = "drop" )%>% mutate( event = ifelse(!is.na(date_death), 1,0), event_date = ifelse(event == 1, date_death, last_obs), time_months = as.numeric( difftime(as.Date(event_date), as.Date("2022-05-01"), units = "days") )/30 )%>% filter(!is.na(time_months) &time_months >= 0) # Quick check table(time_event_df$event) ## ## 0 ## 72 table(time_event_df$genotype) ## ## B G K KW LG OB S Y ##99999999 glimpse(survival_df) 5 ## Rows: 216 ## Columns: 11 ## $ Month <fct> Sep-23, Sep-23, Sep-23, Sep-23, Sep-23, Sep-23, Sep-23,~ ## $ dome <chr> "a", "a", "a", "a", "a", "a", "a", "a", "a", "a", "a", ~ ## $ tag <chr> "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11"~ ## $ genotype <fct> G, B, KW, KW, OB, LG, LG, OB, LG, KW, OB, G, B, K, Y, S~ ## $ tle_alive <dbl> 990, 330, 0, 480, 970, 510, 0, 380, 800, 660, 610, 345,~ ## $ tle_dead <dbl> 35, 0, 160, 120, 0, 90, 260, 0, 70, 109, 710, 280, 350,~ ## $ max_tle <dbl> 1025, 330, 160, 600, 970, 600, 260, 380, 870, 769, 1320~ ## $ perc_mortality <dbl> 3, 0, 100, 20, 0, 15, 100, 0, 8, 14, 54, 45, 49, 100, 1~ ## $ health <fct> bleached, bleached, dead, bleached, pale, bleached, dea~ ## $ month_date <date> 2023-09-01, 2023-09-01, 2023-09-01, 2023-09-01, 2023-0~ ## $ coral_id <chr> "a_1", "a_2", "a_3", "a_4", "a_5", "a_6", "a_7", "a_8",~ glimpse(time_event_df) ## Rows: 72 ## Columns: 8 ## $ coral_id <chr> "a_1", "a_10", "a_11", "a_12", "a_13", "a_14", "a_15", "a_~ ## $ genotype <fct> G, KW, OB, G, B, K, Y, S, S, Y, K, B, B, G, S, Y, K, KW, K~ ## $ dome <chr> "a", "a", "a", "a", "a", "a", "a", "a", "a", "a", "a", "a"~ ## $ date_death <date> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N~ ## $ last_obs <date> 2023-09-01, 2023-09-01, 2023-09-01, 2023-09-01, 2023-09-0~ ## $ event <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0~ ## $ event_date <dbl> 19601, 19601, 19601, 19601, 19601, 19601, 19601, 19601, 19~ ## $ time_months <dbl> 16.26667, 16.26667, 16.26667, 16.26667, 16.26667, 16.26667~ 3.3 Coral Health Metrics Categorical health observations are prepared for multinomial modeling in this section. Genotype and health categories are releveled for interpretability, with counts tabulated by month and genotype. health <- survival_df %>% mutate(fragment_id = paste(dome, tag, sep = "_"), health = relevel(health, ref = "healthy"), genotype = relevel(genotype, ref = "OB")) %>% select(Month, month_date, genotype, health, fragment_id, dome) %>% drop_na(health) counts_tbl <- health %>% count(Month, genotype, health) glimpse(counts_tbl) ## Rows: 50 ## Columns: 4 ## $ Month <fct> May-22, May-22, May-22, May-22, May-22, May-22, May-22, May-2~ ## $ genotype <fct> OB, B, G, K, KW, LG, S, Y, OB, OB, B, G, G, G, K, K, KW, KW, ~ ## $ health <fct> healthy, healthy, healthy, healthy, healthy, healthy, healthy~ ## $ n <int> 9, 9, 9, 9, 9, 9, 9, 9, 8, 1, 9, 7, 1, 1, 8, 1, 4, 5, 3, 1, 2~ 6 4 Statistical Analyses This section describes the statistical modeling workflows used to evaluate treatment effects on coral physiology, survival, and health. Each subsection includes model specifications, diagnostics, and post-hoc comparisons corresponding to analyses presented in the main text and figures. 4.1 Photophysiology 4.1.1 Endosymbiont Density (Gamma GLMM) Endosymbiont densities were modeled using a Gamma generalized linear mixed model (GLMM) with a log link. Fixed effects included sampling period (Before vs. After MHW) and genotype, with a random intercept for dome. # Fit Gamma GLMM (log link) zoox_df$genotype <- factor(zoox_df$genotype, levels = c("OB",setdiff(unique(zoox_df$genotype), "OB"))) m_zoox <- glmer( mean_zoox ~time +genotype +(1|dome), data = zoox_df, family = Gamma(link = "log"), control = glmerControl(optimizer = "bobyqa",optCtrl = list(maxfun = 1e5)) ) summary(m_zoox) ## Generalized linear mixed model fit by maximum likelihood (Laplace ## Approximation) [glmerMod] ## Family: Gamma ( log ) ## Formula: mean_zoox ~ time + genotype + (1 | dome) ## Data: zoox_df ## Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05)) ## ## AIC BIC logLik -2*log(L) df.resid ## 2176.0 2202.3 -1077.0 2154.0 70 ## ## Scaled residuals: ## Min 1Q Median 3Q Max ## -2.25032 -0.82108 -0.09437 0.51559 2.76482 ## ## Random effects: ## Groups Name Variance Std.Dev. ## dome (Intercept) 0.009581 0.09788 ## Residual 0.126987 0.35635 ## Number of obs: 81, groups: dome, 9 ## ## Fixed effects: ## Estimate Std. Error t value Pr(>|z|) ## (Intercept) 13.63717 0.14392 94.757 < 2e-16 *** ## timeAfter -0.94352 0.09611 -9.817 < 2e-16 *** ## genotypeB 0.14123 0.14786 0.955 0.33950 ## genotypeG -0.26979 0.15093 -1.788 0.07385 . ## genotypeK -0.42043 0.15042 -2.795 0.00519 ** 7 ## genotypeKW 0.15513 0.17361 0.894 0.37158 ## genotypeLG -0.28802 0.15816 -1.821 0.06860 . ## genotypeS -0.20601 0.15192 -1.356 0.17508 ## genotypeY 0.18526 0.16641 1.113 0.26558 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Correlation of Fixed Effects: ## (Intr) tmAftr gntypB gntypG gntypK gntyKW gntyLG gntypS ## timeAfter -0.249 ## genotypeB -0.519 -0.083 ## genotypeG -0.507 -0.070 0.490 ## genotypeK -0.511 0.017 0.477 0.468 ## genotypeKW -0.507 0.001 0.423 0.425 0.410 ## genotypeLG -0.471 -0.003 0.458 0.450 0.450 0.391 ## genotypeS -0.534 -0.064 0.486 0.485 0.468 0.429 0.447 ## genotypeY -0.522 -0.001 0.447 0.433 0.431 0.386 0.408 0.437 # Optional: test interaction (not retained) m_zoox_int <- update(m_zoox, . ~.+time:genotype) anova(m_zoox, m_zoox_int, test = "Chisq") npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq) m_zoox 11 2175.989 2202.328 -1076.994 2153.989 NA NA NA m_zoox_int 18 2185.432 2228.532 -1074.716 2149.432 4.557194 7 0.7138197 # Model Diagnostics sim_zoox <- simulateResiduals(m_zoox, n = 2000) plot(sim_zoox) 8 0.0 0.4 0.8 0.0 0.4 0.8 QQ plot residuals Expected Observed KS test: p= 0.38133 Deviation n.s. Outlier test: p= 1 Deviation n.s. Dispersion test: p= 0.122 Deviation n.s. Model predictions (rank transformed) DHARMa residual 0.0 0.4 0.8 0.00 0.50 1.00 DHARMa residual vs. predicted Quantile deviations detected (red curves) Combined adjusted quantile test n.s. DHARMa residual testDispersion(sim_zoox) DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated Simulated values, red line = fitted model. p−value (two.sided) = 0.122 Frequency 0.0 0.5 1.0 1.5 2.0 2.5 3.0 0 10 30 50 ## ## DHARMa nonparametric dispersion test via sd of residuals fitted vs. 9 ## K / S 1.043 0.0999 Inf 1 0.437 0.9999 ## Y / S 1.038 0.1130 Inf 1 0.341 1.0000 ## ## date = Sep-23: ## contrast odds.ratio SE df null z.ratio p.value ## OB / G 1.853 0.2160 Inf 1 5.295 <.0001 ## OB / B 1.790 0.2170 Inf 1 4.806 <.0001 ## OB / KW 0.984 0.1280 Inf 1 -0.124 1.0000 ## OB / LG 0.750 0.0876 Inf 1 -2.464 0.2110 ## OB / K 1.364 0.1520 Inf 1 2.791 0.0972 ## OB / Y 1.430 0.1800 Inf 1 2.840 0.0854 ## OB / S 1.207 0.1380 Inf 1 1.645 0.7228 ## G / B 0.966 0.1250 Inf 1 -0.267 1.0000 ## G / KW 0.531 0.0737 Inf 1 -4.560 0.0001 ## G / LG 0.405 0.0508 Inf 1 -7.202 <.0001 ## G / K 0.736 0.0883 Inf 1 -2.554 0.1731 ## G / Y 0.772 0.1030 Inf 1 -1.938 0.5245 ## G / S 0.652 0.0806 Inf 1 -3.463 0.0125 ## B / KW 0.550 0.0785 Inf 1 -4.192 0.0007 ## B / LG 0.419 0.0542 Inf 1 -6.728 <.0001 ## B / K 0.762 0.0957 Inf 1 -2.165 0.3732 ## B / Y 0.799 0.1100 Inf 1 -1.627 0.7340 ## B / S 0.675 0.0868 Inf 1 -3.060 0.0459 ## KW / LG 0.762 0.1060 Inf 1 -1.958 0.5105 ## KW / K 1.386 0.1860 Inf 1 2.431 0.2260 ## KW / Y 1.453 0.2130 Inf 1 2.546 0.1761 ## KW / S 1.227 0.1680 Inf 1 1.497 0.8092 ## LG / K 1.819 0.2210 Inf 1 4.931 <.0001 ## LG / Y 1.907 0.2560 Inf 1 4.817 <.0001 ## LG / S 1.610 0.2000 Inf 1 3.834 0.0032 ## K / Y 1.049 0.1360 Inf 1 0.366 1.0000 ## K / S 0.885 0.1050 Inf 1 -1.029 0.9702 ## Y / S 0.844 0.1120 Inf 1 -1.274 0.9087 ## ## P value adjustment: tukey method for comparing a family of 8 estimates ## Tests are performed on the log odds ratio scale # Temporal contrasts within each genotype time_contrasts <- contrast(emm_fire, method = "pairwise",by = "genotype",adjust = "tukey") time_contrasts ## genotype = OB: ## contrast odds.ratio SE df null z.ratio p.value ## (May-22) / (Apr-23) 0.690 0.0657 Inf 1 -3.898 0.0003 ## (May-22) / (Sep-23) 0.640 0.0689 Inf 1 -4.146 0.0001 ## (Apr-23) / (Sep-23) 0.928 0.0822 Inf 1 -0.846 0.6741 ## ## genotype = G: ## contrast odds.ratio SE df null z.ratio p.value ## (May-22) / (Apr-23) 0.542 0.0484 Inf 1 -6.860 <.0001 ## (May-22) / (Sep-23) 0.929 0.1010 Inf 1 -0.677 0.7768 ## (Apr-23) / (Sep-23) 1.716 0.1750 Inf 1 5.304 <.0001 ## ## genotype = B: 16 ## contrast odds.ratio SE df null z.ratio p.value ## (May-22) / (Apr-23) 0.541 0.0536 Inf 1 -6.193 <.0001 ## (May-22) / (Sep-23) 0.951 0.1080 Inf 1 -0.445 0.8968 ## (Apr-23) / (Sep-23) 1.756 0.2050 Inf 1 4.832 <.0001 ## ## genotype = KW: ## contrast odds.ratio SE df null z.ratio p.value ## (May-22) / (Apr-23) 0.586 0.0672 Inf 1 -4.659 <.0001 ## (May-22) / (Sep-23) 0.500 0.0617 Inf 1 -5.614 <.0001 ## (Apr-23) / (Sep-23) 0.853 0.1130 Inf 1 -1.202 0.4519 ## ## genotype = LG: ## contrast odds.ratio SE df null z.ratio p.value ## (May-22) / (Apr-23) 0.567 0.0559 Inf 1 -5.747 <.0001 ## (May-22) / (Sep-23) 0.401 0.0435 Inf 1 -8.429 <.0001 ## (Apr-23) / (Sep-23) 0.706 0.0766 Inf 1 -3.211 0.0038 ## ## genotype = K: ## contrast odds.ratio SE df null z.ratio p.value ## (May-22) / (Apr-23) 0.604 0.0549 Inf 1 -5.548 <.0001 ## (May-22) / (Sep-23) 0.750 0.0762 Inf 1 -2.837 0.0127 ## (Apr-23) / (Sep-23) 1.240 0.1210 Inf 1 2.208 0.0698 ## ## genotype = Y: ## contrast odds.ratio SE df null z.ratio p.value ## (May-22) / (Apr-23) 0.580 0.0611 Inf 1 -5.168 <.0001 ## (May-22) / (Sep-23) 0.751 0.0887 Inf 1 -2.421 0.0409 ## (Apr-23) / (Sep-23) 1.294 0.1580 Inf 1 2.121 0.0857 ## ## genotype = S: ## contrast odds.ratio SE df null z.ratio p.value ## (May-22) / (Apr-23) 0.598 0.0551 Inf 1 -5.577 <.0001 ## (May-22) / (Sep-23) 0.630 0.0669 Inf 1 -4.350 <.0001 ## (Apr-23) / (Sep-23) 1.053 0.1080 Inf 1 0.506 0.8684 ## ## P value adjustment: tukey method for comparing a family of 3 estimates ## Tests are performed on the log odds ratio scale 4.2 Mortality and Survival 4.2.1 Accelerated Failure Time (AFT) Model To compare survival duration among genotypes, we fitted parametric Accelerated Failure Time (AFT) models using both Weibull and log-normal distributions. The model with the lowest AIC was retained (log-normal in this case). Model coefficients were exponentiated to obtain time ratios, where values < 1 indicate shorter survival relative to the reference genotype (OB). # Ensure reference genotype is OB survival_df$genotype <- relevel(factor(survival_df$genotype), ref = "OB") aft_weib <- survreg(Surv(time_months, event) ~genotype, data = time_event_df, dist = "weibull") 17 aft_logn <- survreg(Surv(time_months, event) ~genotype, data = time_event_df, dist = "lognormal") print(AIC(aft_weib, aft_logn)) # choose lowest AIC (expected: lognormal) ## df AIC ## aft_weib 9 18 ## aft_logn 9 18 # Extract time-ratios ± 95 % CI get_time_ratios <- function(model) { s<- summary(model) est <- s$table[, "Value"]# coefficient estimates (log time scale) se <- s$table[, "Std. Error"]# standard errors TR <- exp(est) lower <- exp(est -1.96 *se) upper <- exp(est +1.96 *se) data.frame( term = rownames(s$table), estimate = est, time_ratio = TR, lower = lower, upper = upper, p = s$table[, "p"] ) } aft_logn <- survreg(Surv(time_months, event) ~genotype, data = time_event_df, dist = "lognormal") tr_logn <- get_time_ratios(aft_logn) 4.2.2 Mortality (Beta GLMM) Mortality proportions were modeled using a Beta GLMM with a logit link, including genotype, month, and their interaction as fixed effects, and a random intercept for dome. Model diagnostics assessed dispersion, uniformity, and multicollinearity. # Make genotype a factor with OB as the reference mort_df$genotype <- factor(mort_df$genotype, levels = c("OB",setdiff(unique(mort_df$genotype), "OB"))) mod_zoib <- glmmTMB( mort_prop ~genotype *Month +(1|dome), family = beta_family(link = "logit"), data = mort_df ) summary(mod_zoib) ## Family: beta ( logit ) ## Formula: mort_prop ~ genotype * Month + (1 | dome) ## Data: mort_df 18 ## ## AIC BIC logLik -2*log(L) df.resid ## -686.2 -598.5 369.1 -738.2 190 ## ## Random effects: ## ## Conditional model: ## Groups Name Variance Std.Dev. ## dome (Intercept) 0.09049 0.3008 ## Number of obs: 216, groups: dome, 3 ## ## Dispersion parameter for beta family (): 1.12 ## ## Conditional model: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) -1.695e+00 4.274e-01 -3.966 7.32e-05 *** ## genotypeG -2.449e-05 5.398e-01 0.000 0.99996 ## genotypeB -6.237e-06 5.398e-01 0.000 0.99999 ## genotypeKW -1.257e-05 5.398e-01 0.000 0.99998 ## genotypeLG -2.727e-06 5.398e-01 0.000 1.00000 ## genotypeK -2.207e-05 5.398e-01 0.000 0.99997 ## genotypeY -1.122e-05 5.398e-01 0.000 0.99998 ## genotypeS -1.684e-05 5.398e-01 0.000 0.99998 ## MonthApr-23 5.259e-01 5.483e-01 0.959 0.33747 ## MonthSep-23 1.763e+00 5.550e-01 3.177 0.00149 ** ## genotypeG:MonthApr-23 3.743e-01 7.760e-01 0.482 0.62957 ## genotypeB:MonthApr-23 -2.087e-01 7.720e-01 -0.270 0.78691 ## genotypeKW:MonthApr-23 1.741e+00 7.790e-01 2.235 0.02539 * ## genotypeLG:MonthApr-23 1.045e+00 7.849e-01 1.332 0.18302 ## genotypeK:MonthApr-23 -2.362e-01 7.720e-01 -0.306 0.75966 ## genotypeY:MonthApr-23 1.549e+00 7.834e-01 1.977 0.04800 * ## genotypeS:MonthApr-23 -3.644e-01 7.706e-01 -0.473 0.63632 ## genotypeG:MonthSep-23 5.973e-01 7.836e-01 0.762 0.44594 ## genotypeB:MonthSep-23 3.454e-01 7.735e-01 0.446 0.65526 ## genotypeKW:MonthSep-23 1.206e+00 7.734e-01 1.560 0.11887 ## genotypeLG:MonthSep-23 9.284e-01 7.785e-01 1.192 0.23308 ## genotypeK:MonthSep-23 3.316e-01 7.843e-01 0.423 0.67244 ## genotypeY:MonthSep-23 1.302e+00 7.769e-01 1.676 0.09365 . ## genotypeS:MonthSep-23 1.217e+00 7.778e-01 1.564 0.11770 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 # Simulate residuals and assess model assumptions res_mort <- DHARMa::simulateResiduals(mod_zoib, n = 2000) # Residual plots and tests plot(res_mort) # residual vs predicted, QQ-plot 19 0.0 0.4 0.8 0.0 0.4 0.8 QQ plot residuals Expected Observed KS test: p= 0 Deviation significant Outlier test: p= 1 Deviation n.s. Dispersion test: p= 0.695 Deviation n.s. Model predictions (rank transformed) DHARMa residual 0.0 0.4 0.8 0.00 0.50 1.00 DHARMa residual vs. predicted Quantile deviations detected (red curves) Combined adjusted quantile test significant DHARMa residual testDispersion(res_mort) # dispersion � 1 indicates good fit DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated Simulated values, red line = fitted model. p−value (two.sided) = 0.695 Frequency 0.4 0.5 0.6 0.7 0.8 0 5 10 15 20 ## ## DHARMa nonparametric dispersion test via sd of residuals fitted vs. 20 ## simulated ## ## data: simulationOutput ## dispersion = 0.96112, p-value = 0.695 ## alternative hypothesis: two.sided testUniformity(res_mort) # KS test for residual uniformity 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 QQ plot residuals Expected Observed KS test: p= 0 Deviation significant Outlier test: p= 1 Deviation n.s. Dispersion test: p= 0.695 Deviation n.s. ## ## Asymptotic one-sample Kolmogorov-Smirnov test ## ## data: simulationOutput$scaledResiduals ## D = 0.23315, p-value = 1.267e-10 ## alternative hypothesis: two-sided performance::check_collinearity(mod_zoib) Term VIF VIF_CI_lowVIF_CI_highSE_factor Tolerance Tolerance_CI_lowTolerance_CI_high genotype 1533.89801 1215.19535 1936.25392 1.688717 0.0006519 0.0005165 0.0008229 Month 47.05976 37.39414 59.29239 2.619162 0.0212496 0.0168656 0.0267422 genotype:Month31115.46597 24648.28213 39279.56985 176.395765 0.0000321 0.0000255 0.0000406 performance::check_singularity(mod_zoib) ## [1] FALSE 21 4.3 Host Health Metrics 4.3.1 Biased-Reduced Multinomial Logistic Regression Coral health categories (e.g., healthy, pale, bleached, dead) were modeled using a bias-reduced multinomial logistic regression to address potential small-sample bias in certain states. The model tested effects of genotype and sampling month on the probability of each health outcome. Predicted probabilities are also extracted for visualization and multivariate analysis. # Fit bias-reduced multinomial logistic regression m_health <- brmultinom(health ~Month +genotype, data = health) summary(m_health) ## Call: ## brmultinom(formula = health ~ Month + genotype, data = health) ## ## Coefficients: ## (Intercept) MonthApr-23 MonthSep-23 genotypeB genotypeG genotypeK ## pale -4.018123 0.3180008 6.542381 -1.8026280 0.43563647 -1.6422430 ## bleached -9.310781 2.1986065 9.446074 2.9810181 4.24134735 0.7450345 ## disease -4.419854 2.8081640 5.367596 -1.4602715 -0.05872965 -0.4600601 ## predation -4.609449 0.3422697 4.621647 0.4026286 0.94396376 0.3224372 ## dead -7.950194 3.9482140 10.014915 0.8156513 1.97755513 -1.1904781 ## missing -6.998419 1.7299950 7.125984 0.5592091 1.62528269 -0.8548391 ## genotypeKW genotypeLG genotypeS genotypeY ## pale 0.3628353 0.5104544 -0.6628913 0.3080688 ## bleached 4.8160572 5.8317173 4.3771492 5.0528957 ## disease -0.6143925 1.1759528 -1.2548931 -0.6591577 ## predation 1.3033918 1.3922107 0.8939818 1.2920207 ## dead 4.1607542 4.0119910 1.9372006 3.8098326 ## missing 2.5939865 2.7899803 1.5781350 4.0892621 ## ## Std. Errors: ## (Intercept) MonthApr-23 MonthSep-23 genotypeB genotypeG genotypeK ## pale 1.404824 1.420628 1.524391 2.007666 1.452929 1.430918 ## bleached 2.377241 1.500643 1.803428 2.165269 2.091230 2.129267 ## disease 1.426438 1.277581 1.682546 1.653342 1.218865 1.142875 ## predation 1.773167 1.181972 1.559956 2.088649 2.083588 1.865254 ## dead 1.881982 1.467087 1.859907 1.532113 1.422001 1.496591 ## missing 2.263077 1.404775 1.809589 2.453424 2.360847 2.486337 ## genotypeKW genotypeLG genotypeS genotypeY ## pale 1.768945 1.781814 1.844093 1.776284 ## bleached 2.143853 2.095230 2.115199 2.113890 ## disease 1.744500 1.148720 1.712733 1.751565 ## predation 2.157161 2.179013 2.082408 2.157319 ## dead 1.378891 1.395146 1.456078 1.381037 ## missing 2.365349 2.369226 2.371687 2.046216 ## ## Residual Deviance: 254.3719 ## Log-likelihood: -127.186 ## AIC: 374.3719 ## ## Type of estimator: AS_mixed (mixed bias-reducing adjusted score equations) 22 ## Number of Fisher Scoring iterations: 45 # Model fit indices cat("\n--- Model fit ---\n") ## ## --- Model fit --- cat("Residual deviance:",deviance(m_health), "\n") ## Residual deviance: 254.3719 cat("LogLik:",as.numeric(logLik(m_health)), "\n") ## LogLik: -127.186 cat("AIC:",AIC(m_health), "\n") ## AIC: 374.3719 # Predicted probabilities by genotype × month newdat <- expand.grid( Month = unique(health$Month), genotype = unique(health$genotype) ) pred_df <- cbind(newdat, predict(m_health, newdata = newdat, type = "probs")) # Tidy for plotting plot_df <- pred_df %>% pivot_longer( cols = -c(Month, genotype), names_to = "health_status", values_to = "probability" ) # Preview head(plot_df) Month genotype health_status probability Sep-23 G healthy 0.0059656 Sep-23 G pale 0.1151121 Sep-23 G bleached 0.4746851 Sep-23 G disease 0.0145127 Sep-23 G predation 0.0155207 Sep-23 G dead 0.3397761 23 # Plot predicted probabilities ggplot(plot_df, aes(x = Month, y = probability, color = genotype, group = genotype)) + geom_line(size = 1)+ geom_point(size = 2)+ facet_wrap(~health_status, scales = "free_y")+ labs( x = "Month", y = "Predicted probability", color = "Genotype", title = "Bias-reduced multinomial regression of coral health status" )+ theme_minimal(base_size = 14) predation healthy missing pale bleached dead disease May−22Apr−23Sep−23 May−22Apr−23Sep−23 May−22Apr−23Sep−23 0.00 0.05 0.10 0.15 0.20 0.0 0.1 0.2 0.3 0.4 0.0 0.2 0.4 0.6 0.00 0.04 0.08 0.12 0.0 0.2 0.4 0.00 0.25 0.50 0.75 1.00 0.000 0.025 0.050 0.075 0.100 Month Predicted probability Genotype OB B G K KW LG S Y Bias−reduced multinomial regression of coral health status 4.3.2 NMDS ordination of predicted probabilities (exploratory) To visualize similarities in predicted coral health composition among genotypes and timepoints, a non-metric multidimensional scaling (NMDS) ordination was applied to the matrix of predicted probabilities. Health-state vectors were fitted to the ordination using envfit to identify directions of strongest association. nmds_mat <- pred_df %>% select(healthy, pale, bleached, disease, predation, dead, missing) %>% as.matrix() rownames(nmds_mat) <- paste(pred_df$genotype, pred_df$Month, sep = "_") nmds_res <- metaMDS(nmds_mat, distance = "bray",k = 2,trymax = 100) 24 ## Run 0 stress 0.03252542 ## Run 1 stress 0.05633943 ## Run 2 stress 0.05773263 ## Run 3 stress 0.0424488 ## Run 4 stress 0.04471079 ## Run 5 stress 0.04481317 ## Run 6 stress 0.05641732 ## Run 7 stress 0.04244883 ## Run 8 stress 0.03672854 ## Run 9 stress 0.05850263 ## Run 10 stress 0.04294394 ## Run 11 stress 0.04481317 ## Run 12 stress 0.05633943 ## Run 13 stress 0.05741067 ## Run 14 stress 0.05633943 ## Run 15 stress 0.05641733 ## Run 16 stress 0.03252702 ## ... Procrustes: rmse 0.0008197671 max resid 0.003527119 ## ... Similar to previous best ## Run 17 stress 0.05784905 ## Run 18 stress 0.03252542 ## ... Procrustes: rmse 1.239708e-05 max resid 3.0579e-05 ## ... Similar to previous best ## Run 19 stress 0.05735467 ## Run 20 stress 0.03252542 ## ... Procrustes: rmse 6.156161e-06 max resid 1.754569e-05 ## ... Similar to previous best ## *** Best solution repeated 3 times # Extract stress value and ordination coordinates nmds_stress <- round(nmds_res$stress, 3) nmds_coords <- as.data.frame(scores(nmds_res, display = "sites")) nmds_coords <- nmds_coords %>% mutate( genotype = sapply(strsplit(rownames(nmds_coords), "_"), `[`,1), Month = sapply(strsplit(rownames(nmds_coords), "_"), `[`,2) ) cat("\nNMDS 2D stress =", nmds_stress, "\n") ## ## NMDS 2D stress = 0.033 # Fit vectors for all health categories health_vars <- pred_df %>% select(healthy, pale, bleached, disease, predation, dead, missing) efit <- envfit(nmds_res, health_vars, perm = 999) vectors <- as.data.frame(efit$vectors$arrows *sqrt(efit$vectors$r)) %>% rownames_to_column("health_status") 25 "disease" ="orange", "predation" ="purple", "dead" ="red", "missing" ="grey" ) alluvial_plot <- ggplot(data_filtered, aes(x = Month, stratum = health, alluvium = coral_id, fill = health)) + geom_flow(stat = "alluvium",lode.guidance = "frontback",alpha = 0.7)+ geom_stratum(alpha = 0.7)+ scale_fill_manual(values = status_colors) + scale_y_continuous(limits = c(0,10), breaks = seq(0,10,2)) + facet_wrap(~genotype, ncol = 4)+ labs(y = "Fragment count (n)",x = NULL,fill = "Health Status")+ plot_theme + theme(axis.text.x = element_text(angle = 45,hjust = 1)) alluvial_plot LG OB S Y B G K KW May−22 Apr−23 Sep−23 May−22 Apr−23 Sep−23 May−22 Apr−23 Sep−23 May−22 Apr−23 Sep−23 0 2 4 6 8 10 0 2 4 6 8 10 Fragment count (n) Health Status healthy pale bleached disease predation dead missing # Panel B - Ranked probabilities rank_df <- pred_df %>% select(Month, genotype, healthy, pale, bleached, disease, dead) %>% pivot_longer( cols = c(healthy, pale, bleached, disease, dead), names_to = "health_status", values_to = "probability" )%>% group_by(genotype, health_status) %>% summarize(mean_prob = mean(probability), .groups = "drop") 32 # Get mean healthy probabilities per genotype healthy_order <- rank_df %>% filter(health_status == "healthy")%>% arrange(desc(mean_prob)) %>% pull(genotype) # Reorder genotype factor by healthy probability rank_df$genotype <- factor(rank_df$genotype, levels = healthy_order) # Ensure health status order rank_df$health_status <- factor( rank_df$health_status, levels = c("healthy","pale","bleached","disease","dead") ) # Ensure genotypes are in desired order gen_order <- c("ob",sort(setdiff(unique(rank_df$genotype), "ob"))) rank_df$genotype <- factor(rank_df$genotype, levels = gen_order) # Custom colors to match alluvial status_colors_bar <- c( "healthy" ="darkgreen", "pale" ="lightgreen", "bleached" ="yellow", "disease" ="orange", "dead" ="red" ) # Plot rank_plot <- ggplot(rank_df, aes(x = genotype, y = mean_prob, fill = health_status)) + geom_col(position = "stack",width = 0.8,color = "black",linewidth = 0.2,alpha = 0.7)+ scale_fill_manual(values = status_colors_bar) + labs(x = "Genotype",y = "Mean predicted prob.",fill = "Health Status")+ plot_theme + theme(axis.text.x = element_text(angle = 45,hjust = 1), legend.position = "none") rank_plot 33 0.00 0.25 0.50 0.75 1.00 B G K KW LG OB S Y Genotype Mean predicted prob. # Panel C - NMDS ## Keep only selected health categories for vectors vectors_filtered <- vectors %>% filter(!health_status %in% c("predation","missing")) coords_filtered <- nmds_coords %>% filter(!Month %in% c("May-22")) # Extract stress stress_val <- round(nmds_res$stress, 3) # Reorder Month factor nmds_coords$Month <- factor(nmds_coords$Month, levels = c("May-22","Apr-23","Sep-23")) # Now plot nmds_plot <- ggplot(coords_filtered, aes(x = NMDS1, y = NMDS2, color = genotype, shape = Month)) + geom_point(size = 4,stroke = 0.3,alpha = 0.6)+ geom_segment( data = vectors_filtered, aes(x = 0,y = 0,xend = NMDS1, yend = NMDS2), arrow = arrow(length = unit(0.3,"cm")), color = "black", inherit.aes = FALSE )+ geom_text( data = vectors_filtered, aes(x = NMDS1, y = NMDS2, label = health_status), color = "black", 34 nudge_y = 0.075, nudge_x = -0.05, size = 3, inherit.aes = FALSE )+ annotate("text",x = -0.15,y = 0.6, label = paste0("2D stress = ",round(nmds_res$stress, 2)), hjust = 0,vjust = 1,size = 3,fontface = "italic")+ scale_color_brewer(palette = "Dark2",name = "Genotype")+ labs(x = "NMDS1",y = "NMDS2")+ plot_theme + theme(legend.position = "none", axis.text.x = element_text(angle = 45,hjust = 1),) nmds_plot healthy pale bleached disease dead 2D stress = 0.03 −0.5 0.0 0.5 −1.0 −0.5 0.0 0.5 1.0 NMDS1 NMDS2 # Figure 7 - Composite leg_health <- get_legend( alluvial_plot + theme( legend.position = "right", legend.title = element_text(size = 12,face = "plain"), legend.text = element_text(size = 10), legend.key.size = unit(4.5,"mm"), # slightly larger boxes legend.box.spacing = unit(1,"mm") ) ) leg_nmds <- get_legend( 35 nmds_plot + theme( legend.position = "right", legend.title = element_text(size = 10), legend.text = element_text(size = 9), legend.key.size = unit(3.5,"mm"), legend.margin = margin(t = -4,b = 0,unit = "mm")# � upward shift ) ) # Remove legends from plots A<- alluvial_plot +theme(legend.position = "none") B<- rank_plot +theme(legend.position = "none") C<- nmds_plot +theme(legend.position = "none") # Adjust margins directly A_adj <- A+theme(plot.margin = margin(3,5,3,5,"mm")) B_adj <- B+theme(plot.margin = margin(5,3,3,7,"mm")) C_adj <- C+theme(plot.margin = margin(7,3,3,3,"mm")) # Assemble composite figure ## Top row (A + legend) row_top <- plot_grid( A_adj, leg_health, nrow = 1, rel_widths = c(1,0.19), align = "v",axis = "tb", labels = "A", label_size = 14,label_fontface = "plain", label_x = 0.02,label_y = 1.02 ) ## Bottom row (B + C side by side) row_bottom_panels <- plot_grid( B_adj, C_adj, nrow = 1, rel_widths = c(1,1), align = "h",axis = "tb", labels = c("B","C"), label_size = 14,label_fontface = "plain", label_x = c(0.02,0.02), label_y = c(1.03,1.03) ) ## Add NMDS legend on far right row_bottom <- plot_grid( row_bottom_panels, leg_nmds, nrow = 1, rel_widths = c(1,0.19), align = "v",axis = "tb" ) ## Final stack 36 final_plot <- plot_grid( row_top, row_bottom, ncol = 1, rel_heights = c(1.7,1.4), align = "hv",axis = "tblr" ) final_plot LG OB S Y B G K KW May−22 Apr−23 Sep−23 May−22 Apr−23 Sep−23 May−22 Apr−23 Sep−23 May−22 Apr−23 Sep−23 0 2 4 6 8 10 0 2 4 6 8 10 Fragment count (n) A Health Status healthy pale bleached disease predation dead missing 0.00 0.25 0.50 0.75 1.00 BGK KW LG OBSY Genotype Mean predicted prob. B healthy pale bleached disease dead 2D stress = 0.03 −0.5 0.0 0.5 −1.0 −0.5 0.0 0.5 1.0 NMDS1 NMDS2 CMonth Apr−23 Sep−23 Genotype B G K KW LG OB S Y # Save ggsave( filename = "./figures/Figure_7_health_composite.png", plot = final_plot, width = 7,height = 6.5, units = "in",dpi = 600,bg = "white" ) 6 Supplementary Tables This section generates formatted tables (Tables S1–S5) corresponding to model outputs described in the manuscript. Tables are produced using the gt package and exported as Word .docx files for easy integration into supplementary materials. 6.1 Table S1 - GLMM Summary for Endosymbiont Density Fixed-effect estimates and 95 % CIs from the Gamma GLMM modeling mean endosymbiont density with a random intercept for dome. 37 tidy_mod <- broom.mixed::tidy(m_zoox, effects = "fixed",conf.int = TRUE,conf.level = 0.95) tbl_zoox <- tidy_mod %>% mutate( term = dplyr::recode( term, "(Intercept)" ="Intercept (OB, Before)", "timeAfter" ="Time (After vs Before)", "genotypeB" ="Genotype B", "genotypeG" ="Genotype G", "genotypeK" ="Genotype K", "genotypeKW" ="Genotype KW", "genotypeLG" ="Genotype LG", "genotypeS" ="Genotype S", "genotypeY" ="Genotype Y" ), `95% CI`=paste0( "[",formatC(conf.low, digits = 2,format = "f"), ", ", formatC(conf.high, digits = 2,format = "f"), "]" ), p.value.display = ifelse(p.value <0.001,"<0.001",formatC(p.value, digits = 3,format = "f")), p.value.numeric = ifelse(p.value <0.001,0.001, p.value) )%>% select(term, estimate, std.error, statistic, p.value.display, `95% CI`, p.value.numeric) %>% gt() %>% fmt_number(columns = c(estimate, std.error, statistic), decimals = 2)%>% cols_label( term = "Term", estimate = "Estimate", std.error = "SE", statistic = "t-value", p.value.display = "p-value", `95% CI`="95% CI" )%>% cols_hide(columns = p.value.numeric) %>% cols_align(align = "left",columns = term) %>% cols_align( align = "right", columns = c(estimate, std.error, statistic, p.value.display, `95% CI`) )%>% tab_options(table.font.size = 12,data_row.padding = px(4), table.align = "center")%>% tab_style( style = cell_text(weight = "bold"), locations = cells_body(columns = p.value.display, rows = p.value.numeric <0.05) )%>% tab_source_note( source_note = "Model: Gamma GLMM (log link); random intercept for dome." ) tbl_zoox 38 Term Estimate SE t-value p-value 95% CI Intercept (OB, Before) 13.64 0.14 94.76 <0.001 [13.36, 13.92] Time (After vs Before) -0.94 0.10 -9.82 <0.001 [-1.13, -0.76] Genotype B 0.14 0.15 0.96 0.339 [-0.15, 0.43] Genotype G -0.27 0.15 -1.79 0.074 [-0.57, 0.03] Genotype K -0.42 0.15 -2.79 0.005 [-0.72, -0.13] Genotype KW 0.16 0.17 0.89 0.372 [-0.19, 0.50] Genotype LG -0.29 0.16 -1.82 0.069 [-0.60, 0.02] Genotype S -0.21 0.15 -1.36 0.175 [-0.50, 0.09] Genotype Y 0.19 0.17 1.11 0.266 [-0.14, 0.51] Model: Gamma GLMM (log link); random intercept for dome. # Optional: Save to word doc. # gtsave(tbl_zoox, filename = "tables/Table_S1_Zoox_GammaGLMM.docx") 6.2 Table S2 - Beta GLMM Summary for EQY Parameter estimates, odds ratios, and 95 % CIs from the Beta GLMM (logit link) on Fv/Fm, including random intercepts for dome and fragment ID. tidy_fire <- broom.mixed::tidy(m_fire, effects = "fixed",conf.int = TRUE,conf.level = 0.95) tbl_fire <- tidy_fire %>% mutate( term = case_when( term == "(Intercept)" ~"Intercept (baseline: May-22, OB)", term == "dateApr-23" ~"Date: Apr-23", term == "dateSep-23" ~"Date: Sep-23", TRUE ~term ), OR = round(exp(estimate), 2), `95% CI`=paste0("[",round(conf.low, 2), ", ",round(conf.high, 2), "]"), p.display = ifelse(p.value <0.001,"<0.001",formatC(p.value, digits = 3,format = "f")), p.numeric = ifelse(p.value <0.001,0.001, p.value) )%>% select(term, estimate, OR, std.error, statistic, p.display, `95% CI`, p.numeric) %>% gt() %>% fmt_number(columns = c(estimate, OR, std.error, statistic), decimals = 2)%>% cols_label( term = "Term", estimate = "Estimate", OR = "Odds Ratio (OR)", std.error = "SE", statistic = "z-value", p.display = "p-value", `95% CI`="95% CI" )%>% cols_hide(columns = p.numeric) %>% cols_align(align = "left",columns = term) %>% cols_align(align = "right",columns = c(estimate, OR, std.error, statistic, p.display, `95% CI`)) %>% tab_options(table.font.size = 12,data_row.padding = px(6), table.align = "center")%>% 39 Term Estimate Odds Ratio (OR) SE z-value p-value 95% CI Intercept (baseline: May-22, OB) -1.46 0.23 0.08 -17.41 <0.001 [-1.62, -1.3] Date: Apr-23 0.37 1.45 0.10 3.90 <0.001 [0.18, 0.56] Date: Sep-23 0.45 1.56 0.11 4.15 <0.001 [0.24, 0.66] genotypeG -0.24 0.78 0.11 -2.18 0.029 [-0.46, -0.03] genotypeB -0.19 0.83 0.11 -1.69 0.092 [-0.4, 0.03] genotypeKW -0.23 0.79 0.11 -2.07 0.038 [-0.45, -0.01] genotypeLG -0.18 0.83 0.11 -1.63 0.102 [-0.4, 0.04] genotypeK -0.15 0.86 0.11 -1.38 0.167 [-0.37, 0.06] genotypeY -0.20 0.82 0.11 -1.78 0.075 [-0.41, 0.02] genotypeS -0.20 0.82 0.11 -1.84 0.066 [-0.42, 0.01] dateApr-23:genotypeG 0.24 1.27 0.13 1.85 0.064 [-0.01, 0.5] dateSep-23:genotypeG -0.37 0.69 0.15 -2.44 0.015 [-0.67, -0.07] dateApr-23:genotypeB 0.24 1.27 0.14 1.75 0.080 [-0.03, 0.51] dateSep-23:genotypeB -0.40 0.67 0.16 -2.53 0.011 [-0.7, -0.09] dateApr-23:genotypeKW 0.16 1.18 0.15 1.09 0.275 [-0.13, 0.45] dateSep-23:genotypeKW 0.25 1.28 0.16 1.51 0.132 [-0.07, 0.57] dateApr-23:genotypeLG 0.20 1.22 0.14 1.43 0.154 [-0.07, 0.46] dateSep-23:genotypeLG 0.47 1.60 0.15 3.07 0.002 [0.17, 0.77] dateApr-23:genotypeK 0.13 1.14 0.13 1.01 0.314 [-0.13, 0.39] dateSep-23:genotypeK -0.16 0.85 0.15 -1.07 0.286 [-0.45, 0.13] dateApr-23:genotypeY 0.17 1.19 0.14 1.21 0.225 [-0.11, 0.45] dateSep-23:genotypeY -0.16 0.85 0.16 -1.01 0.315 [-0.47, 0.15] dateApr-23:genotypeS 0.14 1.15 0.13 1.07 0.283 [-0.12, 0.4] dateSep-23:genotypeS 0.02 1.02 0.15 0.10 0.919 [-0.28, 0.31] Model: Beta GLMM (logit link); random intercepts for fragment_id and dome. tab_style( style = cell_text(weight = "bold"), locations = cells_body(columns = p.display, rows = p.numeric <0.05) )%>% tab_source_note(source_note = "Model: Beta GLMM (logit link); random intercepts for fragment_id and dome.") tbl_fire gtsave(tbl_fire, filename = "tables/Table_S2_EQY_BetaGLMM.docx") 6.3 Table S3 - AFT Model Summary Time-ratio estimates and confidence intervals from the log-normal AFT model. Ratios < 1 indicate shorter survival relative to reference genotype OB. # Identify reference genotype reference <- "genotypeOB" # Create a row for reference genotype ref_row <- data.frame( term = reference, estimate = NA, time_ratio = 1, lower = NA, 40 upper = NA, p = NA ) # Combine reference with other genotypes tr_logn_full <- bind_rows(ref_row, tr_logn %>% filter(term != "(Intercept)")) # Build polished table tbl_aft <- tr_logn_full %>% mutate( term = case_when( term == "genotypeob" ~"Genotype OB", term == "genotypeg" ~"Genotype G", term == "genotypeb" ~"Genotype B", term == "genotypekw" ~"Genotype KW", term == "genotypelg" ~"Genotype LG", term == "genotypek" ~"Genotype K", term == "genotypey" ~"Genotype Y", term == "genotypes" ~"Genotype S", TRUE ~term ), estimate = ifelse(is.na(estimate), NA,round(estimate, 2)), `Time Ratio`=round(time_ratio, 2), `95% CI`=ifelse( !is.na(lower), paste0("[",round(lower, 2), ", ",round(upper, 2), "]"), "-" ), p.display = ifelse( is.na(p), "-", ifelse(p <0.001,"<0.001",formatC(p, digits = 3,format = "f")) ), p.numeric = ifelse(is.na(p), NA,ifelse(p <0.001,0.001, p)) )%>% select(term, estimate, `Time Ratio`, p.display, `95% CI`, p.numeric) %>% gt() %>% fmt_number(columns = c(estimate, `Time Ratio`), decimals = 2)%>% cols_label( term = "Genotype", estimate = "Estimate", `Time Ratio`="Time Ratio", p.display = "p-value", `95% CI`="95% CI" )%>% cols_hide(columns = p.numeric) %>% cols_align(align = "left",columns = term) %>% cols_align(align = "right",columns = c(estimate, `Time Ratio`, p.display, `95% CI`)) %>% tab_options( table.font.size = 12, data_row.padding = px(6), table.align = "center" )%>% # Bold significant p-values (< 0.05) 41