Behavioral responses buffer seasonal variation more strongly than endocrine and chemical responses in poison frogs

Shirley J. Serrano-Rojas*^a,b, Andrius Pašukonis^c, Mabel Gonzalez^a, Camilo Rodriguez^a, Rodrigo F. Calvo Usto^b, Andrea Carazas Muñiz^b, Cristel Sandoval García^d, Jean Pier Zolorzano^b, Luisa F. Arcila-Pérez^e, Sergio Boluarte-Salinas^f, Esau D. Baldarrago^b, Alfredo Sosa-Salazar^g, Lauren A. O’Connell^a

^a Department of Biology, Stanford University, Stanford, CA 94305, USA
^b Department of Biology, Universidad Nacional de San Antonio Abad del Cusco, Cusco 08003, Perú
^c Institute of Biosciences, Vilnius University Life Sciences Center, Saulėtekio al. 7, LT-10257 Vilnius, Lithuania
^d Department of Biology, Universidad Nacional Mayor de San Marcos, Lima 15081, Perú
^e Department of Biological Sciences, Universidad de los Andes, Bogotá 111711, Colombia
^f Department of Forestry, Universidad Nacional Agraria La Molina, Lima 15024, Perú
^g Department of Biology, Universidad Nacional Federico Villarreal, Lima 15007, Perú

*Corresponding author: shirley.serrano25@gmail.com

Abstract

Seasonal rainfall shapes biological responses in tropical ecosystems, yet how tropical organisms integrate behavioral and physiological responses to cope with seasonality remains poorly understood. We assessed how four poison frog species with contrasting reproductive strategies respond to dry and wet season environmental conditions. We quantified spatial behavior, microhabitat use, hormone release rates, and chemical defenses in two seasonal breeders (Allobates femoralis and Ameerega trivittata) and two year-round breeders (Ameerega macero and Ameerega shihuemoy). Seasonal breeders exhibited pronounced sex-specific shifts in space use, where males expanded their space use during the wet season, likely to track reproductive opportunities, while A. femoralis females increased their spatial use during the dry season, likely responding to foraging demands when prey resources are sparse. Year-round breeders maintained similar space use across seasons, likely reflecting their ability to access key resources within the same space to reproduce year-round. Microhabitat use was flexible, as seasonal breeders shifted toward humid refugia during the dry season and reproduction-associated microhabitats during the wet season, whereas year-round breeders selected microhabitats that facilitate continuous reproduction across seasons. Despite these behavioral responses, water-borne corticosterone and testosterone, as well as chemical defenses showed no consistent seasonal variation, suggesting that behavioral responses may be partially decoupled from shifts in endocrine and chemical defenses. These results support a role for behavioral buffering in mediating responses to seasonal environmental variation. However, given the increasing unpredictability in rainfall timing and intensity as a result of climate change, how these coping strategies will function in the long term is uncertain.

Introduction

Climate change drives biodiversity loss (Pfenning-Butterworth et al. 2024), partially by altering the timing, intensity, and predictability of environmental factors like temperature and rainfall (Webb et al. 2005). For some species, seasonal rainfall strongly shapes resource availability, reproductive timing, and environmental stress (Feng et al. 2013). Tropical ecosystems exhibit pronounced rainfall seasonality, with distinct dry and wet seasons that shape biological and ecological processes (Brown and Shine 2006; Mandl et al. 2018). Thus, seasonal habitats provide a natural context to understand how organisms coordinate behavioral and physiological responses to climatic variability (Varpe 2017; Gotthard et al. 2025). Understanding these mechanisms is critical for predicting resilience or vulnerability to climate change.

Contrasting life-history traits provide a comparative framework for understanding how species respond to seasonal variability (Hjernquist et al. 2012; Hall and MacPherson 2025). Reproductive timing, for example, ranges from reproducing during specific periods of favorable conditions in a season to reproducing continuously throughout the year (Ndithia et al. 2017; Heldstab et al. 2021). Seasonal breeders minimize exposure to unfavorable conditions by limiting their reproduction to periods of high resource availability (Hall et al. 2018; Sun et al. 2020; Soanes et al. 2021). In contrast, year-round breeders exhibit behaviors that may reduce exposure to seasonal constraints while maintaining reproductive activity across fluctuating environmental conditions (Prado et al. 2005; Heinermann et al. 2015). Thus far, relatively few tropical studies have investigated seasonal responses through a framework that integrates behavior, spatial ecology, and physiology across species with contrasting life-history strategies (Lopez et al. 2023). Behavioral responses to seasonal variability may also occur independently of broad endocrine shifts, highlighting the importance of assessing both behavioral and physiological responses simultaneously. Glucocorticoids and gonadal steroids are often used to assess physiological responses to environmental variation because of their roles on energy allocation, reproduction, and behavior (Crespi et al. 2013; Husak et al. 2021). In amphibians, corticosterone is one of the main glucocorticoids involved in metabolic regulation and responses to environmental change, whereas testosterone influences reproductive behaviors, spatial activity, and aggression (Wingfield et al. 1998; Moore and Jessop 2003; Denver 2009; Eikenaar et al. 2012; Husak et al. 2021). Assessing seasonal variation in these hormones alongside behavioral responses can therefore reveal whether seasonal environmental change is accompanied by endocrine adjustments. Filling this critical gap is important in the context of global climate change because climate projections for tropical regions consistently predict alterations in the timing, intensity, and predictability of rainfall (Chadwick et al. 2016; Pendergrass et al. 2017; Song et al. 2023).

Across frog species, seasonal rainfall shapes life-history traits (Watling and Donnelly 2002; Ficetola and Maiorano 2016), alters microhabitat selection, movement, and diet (Born et al. 2010), and drives plasticity in parental care (Schulte and Lötters 2013). Poison frogs (Dendrobatidae and Aromobatidae) offer contrasting life-history traits for assessing seasonal responses. Some species breed seasonally and rely on temporary ponds or phytotelmata for tadpole development, while others breed year-round using similar habitats during the wet season and pools along stream margins during the dry season (Wells 2010). Despite these differences in reproduction timing, poison frogs broadly share parental care behaviors and rely on fine-scale spatial behavior and specific microhabitats for reproduction and thermal refuges (Brown et al. 2008). Poison frogs also sequester alkaloids from arthropod prey as chemical defenses (Saporito et al. 2003; Daly et al. 2005), and the composition of these defenses can vary across space and time, reflecting changes in prey communities across seasons (Saporito et al. 2006). Seasonal changes in rainfall may influence chemical defenses by altering prey availability and foraging opportunities. Thus, assessing chemical defenses can provide further insight into how seasonal environmental variation affects poison frogs via trophic interactions and into the potential consequences of predicted rainfall changes in the tropics.

In this study, we investigated how seasonal rainfall affects poison frogs with contrasting reproductive timing. We focused on two seasonal breeders (Allobates femoralis and Ameerega trivittata), which reproduce mostly during the wet season, and two year-round breeders (Ameerega macero and Ameerega shihuemoy), which reproduce throughout the year. Here, we define behavioral buffering as behavioral adjustments that reduce exposure to unfavorable environmental conditions, including shifts in space use area, or microhabitat preferences across seasons. We hypothesize that seasonal breeders will show stronger behavioral, physiological, and chemical defense responses to rainfall shifts than year-round breeders. Specifically, we predict that seasonal breeders will reduce space use and seek refuge in humid microhabitats during the dry season. Seasonal environmental changes will be accompanied by shifts in corticosterone and testosterone release rates, thus reflecting physiological responses to changing environmental and reproductive conditions. These environmental changes, together with reduced space use, may also alter access to alkaloid-rich prey and lead to seasonal shifts in chemical defenses. In contrast, year-round breeders are expected to exhibit behavioral buffering, maintaining similar movement patterns, space use, hormone release rates, and chemical defenses across seasons despite variation in rainfall.

seasons.

Figure 1: Fig. 1. Study sites, environmental conditions, and focal species. (A) Map of southeastern Peru showing the two study sites: the Manu Learning Centre (purple), where Allobates femoralis, Ameerega macero, and Ameerega shihuemoy were studied, and the Los Amigos Biological Station (blue), where Ameerega trivittata was studied. (B) Mean temperature and (C) daily precipitation (mm/day) during the dry (yellow) and wet (blue) seasons at the Manu Learning Centre. (D–G) Field photographs of tagged individuals of the focal species: (D) A. femoralis (brilliant-thighed poison frog), (E) A. trivittata (three-striped poison frog), (F) A. macero (Manu poison frog), and (G) A. shihuemoy (Amarakaeri poison frog). Comparable environmental data for Los Amigos Biological Station are provided in Supplementary Fig. S1.

Materials and methods

Study site and species

We conducted this study in two sites within the Madre de Dios Region, southeastern Peru (Fig. 1A). Across both sites, seasonal variation is primarily driven by rainfall rather than temperature, reflecting the strong precipitation seasonality typical of tropical forest systems.

Three species (Allobates femoralis, Ameerega macero, and Ameerega shihuemoy) were studied in a regenerating forest within the buffer zone of the Manu Biosphere Reserve (The Manu Learning Centre; 450 – 700 m asl). Based on weather data collected on site during the sampling periods (dry season: August–early October 2022; wet season: January–March 2023), mean daily temperature differed between seasons (Wilcoxon rank-sum: W=1288, p<0.001; Fig. 1B), with the wet season slightly warmer (22.82 °C) than the dry season (22.04 °C). Rainfall probability and intensity were also higher in the wet season (3.75× and ~3.3× greater, respectively; p < 0.001), , and daily rainfall patterns showed pronounced seasonal differences (Fig. 1C).

Ameerega trivittata was studied at Los Amigos Biological Station (300 m asl), 142 km from Manu. Based on weather data collected on site during the sampling periods (dry season: November 2022; wet season: January–February 2024), this floodplain forest showed no seasonal difference in mean daily temperature (Welch t-test: t=0.10, p=0.92; dry = 25.14 °C, wet = 25.11 °C), but rainfall probability was 6.72× higher in the wet season ( p < 0.001; Supplementary Methods S1 and Supplementary Fig. S1).

All four species share similar mating and parental behaviors (Rodríguez and Myers 1993; Roithmair 1994; Kaefer et al. 2012; Serrano-Rojas et al. 2017); Fig. 1D–G). Based on breeding duration, we classified A. femoralis and A. trivittata as seasonal breeders and A. macero and A. shihuemoy as year-round breeders (see Supplementary Methods S2 and Fig. S2). The four focal species belong to the superfamily Dendrobatoidea but differ in their phylogenetic relatedness. Three species (A. macero, A. shihuemoy, and A. trivittata) are members of the family Dendrobatidae and the genus Ameerega (Guillory et al. 2020). In contrast, A. femoralis belongs to the family Aromobatidae, which diverged from Dendrobatidae approximately 40 million years ago (Santos et al. 2009). While we do not account for phylogenetic relatedness in our analyses, this context provides a comparative framework for understanding species differences.

Telemetry, space-use quantification, and microhabitat use

We selected telemetry methods based on body size. Smaller species (A. femoralis, A. macero, and A. shihuemoy) were tracked using harmonic direction finding (HDF) with passive transponders, whereas the larger A. trivittata was tracked using radio-telemetry with very high frequency (VHF) transmitters following (Pašukonis et al. 2014, 2018, 2022). See Supplementary Methods S3 for details and a step-by-step protocol available on protocols.io (Serrano-Rojas et al. 2026).

We tracked 212 frogs across species. Individuals were located 5–7 times daily, recording position and perching substrate. When direct observation was not possible due to dense vegetation, we estimated locations within ~1 m. The tracking area was mapped using compasses and laser distance meters to establish reference points (Ringler et al. 2016). Distances and bearings were recorded using the Epicollect5 app (Centre for Genomic Pathogen Surveillance 2025). We excluded all locations associated with tadpole transport because these movements are highly directional, faster, and cover longer distances than routine movements, which would inflate space-use estimates. Additionally, we excluded tracking days with fewer than four locations and individuals tracked for fewer than five full days, because sampling rate and tracking duration influence spatial parameters. 

We defined microhabitat as the structural substrate and spatial position occupied at each relocation. We categorized substrates into 13 categories. Exposed substrates (e.g., leaf litter, leaves, branches, logs, palm roots, rocks) were frequently used for calling, perching, and moving. Potential clutch sites were those previously observed to be used for egg deposition, such as palm bracts and rolled leaves that form sheltered cavities. Refuge or concealed microhabitats (e.g., root cavities, tree fall areas, rock or soil cavities, and spaces under logs or leaf litter) were primarily used during the hottest periods of the day, although recent treefalls also served as calling sites.

Water-borne hormones

We measured water-borne corticosterone and testosterone release rates in individuals tracked for ≥5 days following an indirect non-invasive water-borne hormone sampling method originally developed for amphibians by Gabor et al. (2013) and subsequently adapted for other frog species, including poison frogs (Baugh et al. 2018; Rodríguez et al. 2022). Briefly, we captured frogs between 15:00 and 18:00 h and placed them individually in glass containers containing 40 mL of distilled water for 1 h, allowing hormones released through the skin and urine to accumulate in the water. Water samples were then filtered through C18 cartridges and hormones were eluted with ethanol for subsequent laboratory analysis. See Supplementary Methods S4 for extraction and processing details and Table S1 for sample sizes.

We quantified corticosterone and testosterone using commercial enzyme immunoassay kits (corticosterone: ADI-900-097; testosterone: ADI-900-065, Enzo Life Sciences, Farmingdale, NY, USA). We converted the hormone concentrations to total hormone released and expressed as water-borne hormone release rates (pg/h; see Supplementary Methods S4). To minimize batch effects, we randomized the order of sample assays, ensuring each batch included individuals from all species and seasons within the same geographical location. Due to differences in sampling timing, A. trivittata samples were processed separately. Since our analysis focuses on within-species seasonal comparisons, this approach does not affect our main inferences. Samples with a coefficient of variation (CV) > 20% were excluded (corticosterone: 38/145; testosterone: 29/145). Intra-assay CVs were 8.27% and 7.38%, and inter-assay CVs were 6.77% and 10.17%, respectively.

Alkaloid collection, extraction, and annotation

We used non-lethal metabolite profiling of skin secretions following Gonzalez (2021) to extract and quantify alkaloids in frogs tracked for ≥5 days (see Supplementary Methods S5 for details). Alkaloids were analyzed using gas chromatography/mass spectrometry (GC/MS) following Gonzalez et al. (2021), and annotated using the Global Natural Products Social Molecular Networking (GNPS) pipeline (Supplementary Methods S6; Supplementary Table S2).

Statistical analyses

All statistics were conducted in R (version 2025.09.0). While statistical methods are briefly described here, a complete description can be found in the Supplementary Materials. Code is available on the GitHub repository (https://github.com/laurenoconnelllab/Behavioral-responses-buffer-seasonal-variation-more-strongly-than-endocrine-and-chemical-responses).

#rm(list = ls())
library(broom)
library(car)
library(circlize)
library(conflicted)
library(cowplot)
library(performance)
library(DHARMa)
library(emmeans)
library(glmmTMB)
library(ggeffects)
library(ggtext)
library(ggridges)
library(gt)
library(ggfortify)
library(ggh4x)
library(ggpubr)
library(here)
library(kableExtra)
library(lme4)
library(magick)
library(patchwork)
library(purrr)
library(RColorBrewer)
library(rstatix)
library(scales)
library(see)
library(stringr)
library(sf)
library(tidyverse)
library(vegan)
library(wesanderson)
conflicts_prefer(dplyr::select)# Will prefer dplyr::select over any other package. 
conflicts_prefer(dplyr::filter)# Will prefer dplyr::filter over any other package. 
conflicts_prefer(dplyr::lag)# Will prefer dplyr::lag over any other package.
conflicts_prefer(dplyr::recode)# Will prefer dplyr::recode over any other package.
source(here("Functions", "Functions.r"))

# Space use data
tracking_all <- read_csv(here("Data/all_tracking_spaceuse.csv"))
daily_movement_envir <- read_csv(here("Data/daily_movement_envir.csv"))

# Microhabitat data
microhabitat_data <- read_csv(here("Data/clean_microhabitat_data.csv"))

# Hormone data
testosterone_data <- read_csv(here("Data/clean_testosterone_data.csv"))
corticosterone_data <- read_csv(here("Data/clean_corticosterone_data.csv"))
spaceuse_hormones<- read_csv(here("Data/clean_spaceuse_hormone_data.csv"))

# Alkaloid data
alk_data <- read_csv(here("Data/alkaloids_filtered_26012026.csv"))
select_alk <- read_csv(here("Data/revised_alkaloid_annotations.csv"))
metadata <- read_csv(here("Data/metadata_alkaloids.csv"))

Space use and daily movement

We estimated individual space use (m²) as the 95% utilization distribution (UD) derived from kernel density estimation (KDE), using the kernelUD() function from the adehabitatHR R package (Calenge 2006) and a conservative plug-in bandwidth selection method with the hpi() function from the ks R package (Chacón and Duong 2018). Most individuals were sampled in only one season, although fourteen individuals were recaptured across seasons. Therefore, seasonal comparisons primarily represent population-level changes in space use rather than within-individual seasonal changes. For each species, we tested the effects of season and sex, and their interaction, on space use area using linear models, followed by post hoc tests (Supplementary Methods S7). Models included snout–vent length (SVL) and tracking duration as covariates to account for body-size differences that may influence movement capacity and variation in sampling effort among individuals, respectively. To explore seasonal shifts in location, we quantified seasonal space use overlap and centroid shifts for recaptured individuals (see detailed methods in Supplementary Methods S8).

As tadpole transport events were recorded unevenly across individuals, including them in space use area estimates would have introduced sampling bias and reduced comparability among individuals. However, we provide complementary analysis including tadpole transport in Supplementary Methods S9.

Additionally, we estimated individual daily movement and assessed how rainfall and temperature influence daily movement using linear mixed-effect models for each species (Supplementary Methods S10).

# Space use models

#Clean data
tracking_all_noNA <- tracking_all %>%
  filter(!is.na(HPI95)) %>% #NA are those frogs with too few tracking days
  mutate(
    season = factor(season, levels = c("dry", "wet")),
    sex = factor(sex, levels = c("f", "m")),
    species = factor(species),
    log_HPI95 = log(HPI95))

#Subset the data
tracking_af <- tracking_all_noNA %>% filter(species == "Allobates femoralis")
tracking_am <- tracking_all_noNA %>% filter(species == "Ameerega macero")
tracking_at <- tracking_all_noNA %>% filter(species == "Ameerega trivittata")
tracking_as <- tracking_all_noNA %>% filter(species == "Ameerega shihuemoy")

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

#Allobates femoralis
#Normality
#hist(tracking_af$log_HPI95)

#Rename data
df <- tracking_af

#Selected linear model
lm1 <- lm(log_HPI95 ~ season * sex + svl + days_tracked, data = df)
af_model_output_df <- summary(lm1)$coefficients %>% #this info is for TableS3
  as.data.frame() %>%
  tibble::rownames_to_column("Predictor")
car::Anova(lm1, type = 3)

#Post-hoc comparisons

#Pair comparisons
# estimated marginal means of season by sex
emm <- emmeans(lm1, ~ season | sex) #marginal effect of season within each level of sex
emm_season <- emmeans(lm1, ~ season) # just for the plot

emm_df_af <- as.data.frame(emm)
emm_df_af_season <- as.data.frame(emm_season)

#Combine results
emm_df_af_all <- bind_rows(emm_df_af, 
                           emm_df_af_season %>%mutate(sex = "All"))
#Pairwise contrasts
emm_contrasts_af <- contrast(emm, method = "pairwise", by= "sex" ,adjust = "fdr") %>%
  summary(infer = TRUE) %>% #this is to add additional info such as conf interv, p valies, sign, inference
  as.data.frame()

emm_contrasts_season_af <- contrast(emm_season, method = "pairwise", adjust = "fdr") %>%
  summary(infer = TRUE) %>% 
  as.data.frame()

#Combine results
emm_contrast_af_all <- bind_rows(emm_contrasts_af, 
                                 emm_contrasts_season_af %>%mutate(sex = "All"))
#Add species name
emm_df_af_all$species <- "Allobates femoralis"
emm_contrast_af_all$species <- "Allobates femoralis"
af_model_output_df$species <- "Allobates femoralis"

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

#Ameerega trivittata

#Normality
#hist(tracking_at$log_HPI95)

#Rename data
df <- tracking_at

#Selected linear model
lm5 <- lm(log_HPI95 ~ season+sex, data = df) #selected
at_model_output_df <- summary(lm5)$coefficients %>% 
  as.data.frame() %>%
  tibble::rownames_to_column("Predictor")
lm3 <- lm(log_HPI95 ~ season * sex  + days_tracked, data = df) #just for the plot
summary(lm5)
car::Anova(lm5, type = 2)

#Post-hoc comparisons
#Pair comparisons
#Estimated marginal means by season
emm_season <- emmeans(lm5, ~ season) #this is the important pairwise
emm <- emmeans(lm3, ~ season | sex) #just for the plot

emm_df_at <- as.data.frame(emm)
emm_df_at_season <- as.data.frame(emm_season)

#Combine results
emm_df_at_all <- bind_rows(emm_df_at, 
                           emm_df_at_season %>%mutate(sex = "All"))
#Pairwise contrasts
emm_contrasts_season_at <- contrast(emm_season, method = "pairwise", adjust = "fdr") %>%
  summary(infer = TRUE) %>%
  as.data.frame()

emm_contrasts_at <- contrast(emm, method = "pairwise", by= "sex" ,adjust = "fdr") %>%
  summary(infer = TRUE) %>% #this is to add additional info such as conf interv, p valies, sign, inference
  as.data.frame()

#Combine results
emm_contrast_at_all <- bind_rows(emm_contrasts_at, 
                                 emm_contrasts_season_at %>%mutate(sex = "All"))

#Add species name
emm_df_at_all$species <- "Ameerega trivittata"
emm_contrast_at_all$species <- "Ameerega trivittata"
at_model_output_df$species <- "Ameerega trivittata"

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

#Ameerega macero

#Normality
#hist(tracking_am$log_HPI95)

#Rename data
df <- tracking_am

#Selected linear model
lm3 <- lm(log_HPI95 ~ season * sex  + svl, data = df)
am_model_output_df <- summary(lm3)$coefficients %>% 
  as.data.frame() %>%
  tibble::rownames_to_column("Predictor")
summary(lm3)
car::Anova(lm3, type = 3)

#Post-hoc comparisons
#Pair comparisons
#Estimated marginal means of season by sex
emm <- emmeans(lm3, ~ season | sex)
emm_season <- emmeans(lm3, ~ season) # just for the plot

emm_df_am <- as.data.frame(emm)
emm_df_am_season <- as.data.frame(emm_season)

#Combine results
emm_df_am_all <- bind_rows(emm_df_am, 
                           emm_df_am_season %>%mutate(sex = "All"))
#Pairwise contrasts
emm_contrasts_am <- contrast(emm, method = "pairwise", by= "sex", adjust = "fdr") %>%
  summary(infer = TRUE) %>%
  as.data.frame()

emm_contrasts_season_am <- contrast(emm_season, method = "pairwise", adjust = "fdr") %>%
  summary(infer = TRUE) %>% 
  as.data.frame()

#Combine results
emm_contrast_am_all <- bind_rows(emm_contrasts_am, 
                                 emm_contrasts_season_am %>%mutate(sex = "All"))

#Add species name
emm_df_am_all$species <- "Ameerega macero"
emm_contrast_am_all$species <- "Ameerega macero"
am_model_output_df$species <- "Ameerega macero"

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

#Ameerega shihuemoy

#Normality
#hist(tracking_as$log_HPI95)

#Rename data
df <- tracking_as

#Selected linear model
lm3 <- lm(log_HPI95 ~ season + svl, data = df) #selected model
as_model_output_df <- summary(lm3)$coefficients %>% 
  as.data.frame() %>%
  tibble::rownames_to_column("Predictor")
lm6 <- lm(log_HPI95 ~ season*sex, data = df) #just for the plot
summary(lm3)
car::Anova(lm3, type = 2)

#No posthoc because not significant difference
emm_season <- emmeans(lm3, ~ season) #just for the plot
emm <- emmeans(lm6, ~ season | sex) #just for the plot

emm_df_as <- as.data.frame(emm)
emm_df_as_season <- as.data.frame(emm_season)

#Combine results
emm_df_as_all <- bind_rows(emm_df_as, 
                           emm_df_as_season %>%mutate(sex = "All"))
#Pairwise contrasts
emm_contrasts_season_as <- contrast(emm_season, method = "pairwise") %>%
  summary(infer = TRUE) %>%
  as.data.frame()

emm_contrasts_as <- contrast(emm, method = "pairwise", by= "sex" ,adjust = "fdr") %>%
  summary(infer = TRUE) %>%
  as.data.frame()

#Combine results
emm_contrast_as_all <- bind_rows(emm_contrasts_as, 
                                 emm_contrasts_season_as %>%mutate(sex = "All"))

#Add species name to both
emm_df_as_all$species <- "Ameerega shihuemoy"
emm_contrast_as_all$species <- "Ameerega shihuemoy"
as_model_output_df$species <- "Ameerega shihuemoy"

#Merge data for plot
emm_df <- bind_rows(
  emm_df_af_all,
  emm_df_at_all,
  emm_df_am_all,
  emm_df_as_all)

emm_contrasts <- bind_rows(
  emm_contrast_af_all,
  emm_contrast_at_all,
  emm_contrast_am_all,
  emm_contrast_as_all)

#============================================================================================================
# Figure 2
# Species-specific seasonal shifts in space use area reflect contrasting breeding strategies in poison frogs 
#============================================================================================================

# Plot settings
species_levels <- c(
  "Allobates femoralis",
  "Ameerega trivittata",
  "Ameerega macero",
  "Ameerega shihuemoy")

species_colors <- c(
  "Allobates femoralis" = "#CBC3E3",
  "Ameerega trivittata" = "#CBC3E3",
  "Ameerega macero"     = "#b5d1a1",
  "Ameerega shihuemoy"  = "#b5d1a1")

season_cols <- c(
  dry = "#E7B800",
  wet = "#0072B2")

species_labels <- c(
  "Allobates femoralis" = "Allobates\nfemoralis",
  "Ameerega trivittata" = "Ameerega\ntrivittata",
  "Ameerega macero"     = "Ameerega\nmacero",
  "Ameerega shihuemoy"  = "Ameerega\nshihuemoy")

# Factor ordering
for (dat in list(emm_df, df)) {
  dat$season  <- factor(dat$season, levels = c("dry", "wet"))
  dat$sex     <- factor(dat$sex, levels = c("f", "m", "All"))
}
emm_df$species        <- factor(emm_df$species, levels = species_levels)
emm_contrasts$species <- factor(emm_contrasts$species, levels = species_levels)

#Extract p-values per group and make a label
emm_labels <- emm_contrasts %>%
  as.data.frame() %>%
  mutate(
    p_label = case_when(
      p.value < 0.001 ~ "***",
      p.value < 0.01  ~ "**",
      p.value < 0.05  ~ "*",
      TRUE ~ " " )) %>%
  select(species, sex, p_label)

# Remove singleton observation
emm_df_plot <- emm_df %>%
  filter(!(species == "Ameerega shihuemoy" &
             sex == "m" &
             season == "dry"))

# Plot
plot_figure2 <- ggplot() +
  geom_linerange(data = emm_df_plot,
                 aes(season, ymin = lower.CL, ymax = upper.CL, color = season),
                 linewidth = 5) +
  geom_point(data = emm_df_plot, aes(season, emmean), shape = 21, 
             fill = "black", color = "black", size = 3) +
  scale_color_manual(values = season_cols) +
  scale_fill_manual(values = season_cols) +
  ggh4x::facet_grid2(sex ~ species,
                     labeller = labeller(species = species_labels,
                                         sex = c(f = "Female", m = "Male")),
                     scales = "free",
                     strip = ggh4x::strip_themed(
                       background_x = elem_list_rect(fill = species_colors, alpha = 0.5),
                       text_x = elem_list_text(face = "italic", size = 14))) +
  geom_text(
    data = emm_labels,
    aes(x = 1.5, y = 6.2, label = p_label),
    inherit.aes = FALSE,
    size = 14) +
  scale_x_discrete(
    labels = c(dry = "Dry", wet = "Wet"),
    expand = expansion(add = 1)) +
  scale_y_continuous(limits = c(0.5, 7)) +
  labs(y = expression(log ~ Space ~ Use)) +
  theme_bw(base_size = 20) +
  theme(
    legend.position = "none",
    axis.title.x = element_blank(),
    axis.text.x = element_text(color = "black", size = 12),
    aspect.ratio = 2)


#==========================================================================
# Supplementary Table S3
# Predictors of space use (HPI95) 
#==========================================================================
#Merge data for Table S3
model_output <- bind_rows(
  af_model_output_df,
  at_model_output_df,
  am_model_output_df,
  as_model_output_df)

table_results_s3 <- model_output %>%
  filter(Predictor != "(Intercept)")%>%
  transmute(
    Species = paste0("*", species, "*"), #this is for italics
    Predictor = recode(
      Predictor,
      seasonwet = "Season (wet)",
      sexm = "Sex (male)",
      svl = "SVL",
      days_tracked = "Days tracked",
      `seasonwet:sexm` = "Season (wet) × Sex (male)"),
    β = Estimate,
    SE = round(`Std. Error`, 3),
    t = round(`t value`, 3),
    p = round(`Pr(>|t|)`, 3),
    Sig = case_when(
      `Pr(>|t|)` < 0.001 ~ "***",
      `Pr(>|t|)` < 0.01 ~ "**",
      `Pr(>|t|)` < 0.05 ~ "*",
      TRUE ~ ""))

table_results_s3 %>%
  gt() %>%
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_column_labels()) %>%
  fmt_markdown(columns = Species) %>% #this is for italics
  fmt_number(
    columns = c(β, SE, t, p),
    decimals = 3) %>%
  cols_label(
    Sig = "")

#==========================================================================
# Supplementary Table S4
# Pairwise contrasts of seasonal differences in space use(log-transformed) 
#==========================================================================

table_results_s4 <- emm_contrasts %>%
  filter(
    (species == "Allobates femoralis" & sex %in% c("f", "m")) |
    (species == "Ameerega trivittata" & sex == "All") |
    (species == "Ameerega macero" & sex %in% c("f", "m")) |
    (species == "Ameerega shihuemoy" & sex == "All")) %>%
  transmute(
    Species = paste0("*", species, "*"), #this is for italics
    Sex = sex,
    Contrast = gsub("dry - wet", "Dry − Wet", contrast),
    Estimate = round(estimate, 3),
    SE = round(SE, 3),
    `t-value` = round(t.ratio, 2),
    `p-value` = round(p.value, 3),
    Sig = case_when(
      p.value < 0.001 ~ "***",
      p.value < 0.01 ~ "**",
      p.value < 0.05 ~ "*",
      TRUE ~ ""))

table_results_s4 %>%
  gt() %>%
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_column_labels()) %>%
  fmt_markdown(columns = Species) %>% #this is for italics
  fmt_number(
    columns = c(Estimate, SE, `t-value`, `p-value`),
    decimals = 3) %>%
  cols_label(
    Sig = "")

# Space use overlap and centroid shifts

# Allobates femoralis

#recapture1
af_03f_dry <- st_read(here("Data/Individual_shp_recaptures/tracking2022/Allobates_femoralis/HPI_95_femotrack03f.shp"))
af_34f_wet <- st_read(here("Data/Individual_shp_recaptures/tracking2023/Allobates_femoralis/HPI_95_femotrack34f.shp"))

#recapture2
af_30f_dry <- st_read(here("Data/Individual_shp_recaptures/tracking2022/Allobates_femoralis/HPI_95_femotrack30f.shp"))
af_33f_wet <- st_read(here("Data/Individual_shp_recaptures/tracking2023/Allobates_femoralis/HPI_95_femotrack33f.shp"))

#recapture3
af_06m_dry <- st_read(here("Data/Individual_shp_recaptures/tracking2022/Allobates_femoralis/HPI_95_femotrack06m.shp"))
af_05m_wet <- st_read(here("Data/Individual_shp_recaptures/tracking2023/Allobates_femoralis/HPI_95_femotrack05m.shp"))

#recapture4
af_13m_dry <- st_read(here("Data/Individual_shp_recaptures/tracking2022/Allobates_femoralis/HPI_95_femotrack13m.shp"))
af_02m_wet <- st_read(here("Data/Individual_shp_recaptures/tracking2023/Allobates_femoralis/HPI_95_femotrack02m.shp"))

# Ameerega trivittata

#recapture 5
at_08f_dry <- st_read(here("Data/Individual_shp_recaptures/tracking2022/Ameerega_trivittata/HPI_95_trivitrack08.shp"))
at_06f_wet <- st_read(here("Data/Individual_shp_recaptures/tracking2024/Ameerega_trivittata/HPI_95_trivitrack06.shp"))

#recapture 6
at_06m_dry <- st_read(here("Data/Individual_shp_recaptures/tracking2022/Ameerega_trivittata/HPI_95_trivitrack06.shp"))
at_07m_wet <- st_read(here("Data/Individual_shp_recaptures/tracking2024/Ameerega_trivittata/HPI_95_trivitrack07.shp"))

# Ameerega macero

#recapture7
am_03f_dry <- st_read(here("Data/Individual_shp_recaptures/tracking2022/Ameerega_macero/HPI_95_macetrack03f.shp"))
am_32f_wet <- st_read(here("Data/Individual_shp_recaptures/tracking2023/Ameerega_macero/HPI_95_macetrack32f.shp"))

#recapture8
am_09f_dry <- st_read(here("Data/Individual_shp_recaptures/tracking2022/Ameerega_macero/HPI_95_macetrack09f.shp"))
am_03f_wet <- st_read(here("Data/Individual_shp_recaptures/tracking2023/Ameerega_macero/HPI_95_macetrack03f.shp"))

#recapture9
am_29f_dry <- st_read(here("Data/Individual_shp_recaptures/tracking2022/Ameerega_macero/HPI_95_macetrack29f.shp"))
am_25f_wet <- st_read(here("Data/Individual_shp_recaptures/tracking2023/Ameerega_macero/HPI_95_macetrack25f.shp"))

#recapture10
am_01m_dry <- st_read(here("Data/Individual_shp_recaptures/tracking2022/Ameerega_macero/HPI_95_macetrack01m.shp"))
am_22m_wet <- st_read(here("Data/Individual_shp_recaptures/tracking2023/Ameerega_macero/HPI_95_macetrack22m.shp"))

#recapture11
am_08m_dry <- st_read(here("Data/Individual_shp_recaptures/tracking2022/Ameerega_macero/HPI_95_macetrack08m.shp"))
am_01m_wet <- st_read(here("Data/Individual_shp_recaptures/tracking2023/Ameerega_macero/HPI_95_macetrack01m.shp"))

#recapture12
am_12m_dry <- st_read(here("Data/Individual_shp_recaptures/tracking2022/Ameerega_macero/HPI_95_macetrack12m.shp"))
am_06m_wet <- st_read(here("Data/Individual_shp_recaptures/tracking2023/Ameerega_macero/HPI_95_macetrack06m.shp"))

#recapture13
am_34m_dry <- st_read(here("Data/Individual_shp_recaptures/tracking2022/Ameerega_macero/HPI_95_macetrack34m.shp"))
am_27m_wet <- st_read(here("Data/Individual_shp_recaptures/tracking2023/Ameerega_macero/HPI_95_macetrack27m.shp"))

# Ameerega shihuemoy

#recapture14
as_12_dry <- st_read(here("Data/Individual_shp_recaptures/tracking2022/Ameerega_shihuemoy/HPI_95_shihuetrack12f.shp"))
as_17_wet <- st_read(here("Data/Individual_shp_recaptures/tracking2023/Ameerega_shihuemoy/HPI_95_shihuetrack17f.shp"))

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

# Define the recapture pairs
recaptures <- list(
  list(id="af_f1", species="Allobates_femoralis", sex="F", dry=af_03f_dry, wet=af_34f_wet),
  list(id="af_f2", species="Allobates_femoralis", sex="F", dry=af_30f_dry, wet=af_33f_wet),
  list(id="af_m1", species="Allobates_femoralis", sex="M", dry=af_06m_dry, wet=af_05m_wet),
  list(id="af_m2", species="Allobates_femoralis", sex="M", dry=af_13m_dry, wet=af_02m_wet),

  list(id="at_f1", species="Ameerega_trivittata", sex="F", dry=at_08f_dry, wet=at_06f_wet),
  list(id="at_m1", species="Ameerega_trivittata", sex="M", dry=at_06m_dry, wet=at_07m_wet),

  list(id="am_f1", species="Ameerega_macero", sex="F", dry=am_03f_dry, wet=am_32f_wet),
  list(id="am_f2", species="Ameerega_macero", sex="F", dry=am_09f_dry, wet=am_03f_wet),
  list(id="am_f3", species="Ameerega_macero", sex="F", dry=am_29f_dry, wet=am_25f_wet),

  list(id="am_m1", species="Ameerega_macero", sex="M", dry=am_01m_dry, wet=am_22m_wet),
  list(id="am_m2", species="Ameerega_macero", sex="M", dry=am_08m_dry, wet=am_01m_wet),
  list(id="am_m3", species="Ameerega_macero", sex="M", dry=am_12m_dry, wet=am_06m_wet),
  list(id="am_m4", species="Ameerega_macero", sex="M", dry=am_34m_dry, wet=am_27m_wet),

  list(id="as_f1", species="Ameerega_shihuemoy", sex="F", dry=as_12_dry, wet=as_17_wet))

#===================================================================================================
# Supplementary Table S5
# Individual-level space use overlap and centroid shifts between recaptures for each species and sex
#===================================================================================================

# Compute overlap for all pairs
results <- lapply(recaptures, function(recapture) {
  overlap_results <- compute_overlap_poly(
    dry_polygon = recapture$dry, 
    wet_polygon = recapture$wet)
  data.frame(
    id = recapture$id,
    species = recapture$species,
    sex = recapture$sex,
    overlap = overlap_results$overlap,
    centroid_dist_m = overlap_results$dist)
})
results <- do.call(rbind, results) #merge results for each species

table_results_s5 <- results%>%
  transmute(
    Species = paste0("*", gsub("_", " ", species), "*"), #to correct sp format
    Sex = recode(
      sex,
      `F` = "female",
      `M` = "male"),
    `Recapture ID` = id,
    `Space Use Overlap` = round(overlap,3),
    `Centroid shifts (m)` = round(centroid_dist_m,3))

table_results_s5 %>%
  gt() %>%
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_column_labels()) %>%
  fmt_markdown(columns = Species) %>% #this is for italics
  fmt_number(
    columns = c(`Space Use Overlap`, `Centroid shifts (m)`),
    decimals = 3)

#=============================================================================================
# Supplementary Table S6
# Mean ± SD of space use overlap and centroid shifts (m) calculated for each species and sex
#=============================================================================================
# Summarize the results for Table S6
summary_results <- results %>% 
  group_by(species, sex) %>%
  summarise(
    n = n(),
    mean_overlap = mean(overlap),
    sd_overlap = sd(overlap),
    mean_dist = mean(centroid_dist_m),
    sd_dist = sd(centroid_dist_m),
    .groups = "drop")

summary_results <- summary_results %>%
  mutate(
    species = factor(species, levels = unique(results$species))) %>%
  arrange(species)

table_results_s6 <- summary_results %>%
  transmute(
    Species = paste0("*", gsub("_", " ", species), "*"),
    Sex = recode(
      sex,
      `F` = "female",
      `M` = "male"),
    n = n,
    `Space Use Overlap mean` = round(mean_overlap,3),
    `Space Use Overlap SD` = round(sd_overlap,3),
    `Centroid shifts (m) mean` = round(mean_dist,3),
    `Centroid shifts (m) SD` = round(sd_dist,3))

table_results_s6 %>%
  gt() %>%
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_column_labels()) %>%
  fmt_markdown(columns = Species) %>% #this is for italics
  fmt_number(
    columns = c(`Space Use Overlap mean`,
                `Space Use Overlap SD`,
                `Centroid shifts (m) mean`,
                `Centroid shifts (m) SD`),
    decimals = 3)

# Daily movement models

# Subset data
daily_movement_envir_af <- daily_movement_envir%>%filter(species=="Allobates femoralis")
daily_movement_envir_at <- daily_movement_envir%>%filter(species=="Ameerega trivittata")
daily_movement_envir_am <- daily_movement_envir%>%filter(species=="Ameerega macero")
daily_movement_envir_as <- daily_movement_envir%>%filter(species=="Ameerega shihuemoy")

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

# Allobates femoralis
model_af <- lmer(log(daily_dist) ~ daily_rainfall+mean_temp_daily+(1|ID_season), data = daily_movement_envir_af) #selected model: model4 = model_af
summary(model_af) # model coeficients
car::Anova(model_af, type = 2) # fixed effects

#Extract coefficients and Wald statistics for summary tables
af_results <- extract_model_results(model_af, "Allobates femoralis")
#Model diagnosis
#plot(simulateResiduals(model_af)) # all good

# Plot for Supplementary Fig. S3 Effect of rainfall and temperature on daily movement across poison frog species
plot_af <- plot_effect_dailymovement(model_af,
                                     predictor = "daily_rainfall", #significant variable
                                     x_label = "Daily rainfall",
                                     color = "steelblue")

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

# Ameerega trivittata
model_at <- lmer(log(daily_dist) ~ daily_rainfall+mean_temp_daily+(1|ID_season), data = daily_movement_envir_at)#selected model: model4 = model_at
summary(model_at) # model coeficients
car::Anova(model_at, type = 2) # fixed effects

#Extract coefficients and Wald statistics for summary tables
at_results <- extract_model_results(model_at, "Ameerega trivittata")
#Model diagnosis
#plot(simulateResiduals(model_at)) # best fit

# Plot for Supplementary Fig. S3 Effect of rainfall and temperature on daily movement across poison frog species
plot_at <- plot_effect_dailymovement(model_at,
                                     predictor = "mean_temp_daily", #significant variable
                                     x_label = "Mean daily temperature",
                                     color = "#E7B800")

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

# Ameerega macero
model_am <- lmer(log(daily_dist) ~ daily_rainfall+mean_temp_daily+sex+svl+(1|ID_season), data = daily_movement_envir_am) #selected model: model2 = model_am
summary(model_am) # model coeficients
car::Anova(model_am, type = 2) # fixed effects

#Extract coefficients and Wald statistics for summary tables
am_results <- extract_model_results(model_am, "Ameerega macero")
#Model diagnosis
#plot(simulateResiduals(model_am)) # all good

# Plot for Supplementary Fig. S3 Effect of rainfall and temperature on daily movement across poison frog species.
plot_am <- plot_effect_dailymovement(model_am,
                                     predictor = "mean_temp_daily", #significant variable
                                     x_label = "Mean daily temperature",
                                     color = "#E7B800")
#================================================================================================================

# Ameerega shihuemoy
model_as <- lmer(log(daily_dist) ~ daily_rainfall+mean_temp_daily+(1|ID_season), data = daily_movement_envir_as)#selected model: model4 = model_as
summary(model_as) # model coeficients
car::Anova(model_as, type = 2) # fixed effects

#Extract coefficients and Wald statistics for summary tables
as_results <- extract_model_results(model_as, "Ameerega shihuemoy")
#Model diagnosis
#plot(simulateResiduals(model_as)) # all good

# Plot for Supplementary Fig. S3 Effect of rainfall and temperature on daily movement across poison frog species.
plot_as <- plot_effect_dailymovement(model_as,
                                     predictor = "mean_temp_daily", #non - significant
                                     x_label = "Mean daily temperature",
                                     color = "#E7B800")


#=====================================================
# Supplementary Table S8
# Environmental predictors of daily movement distance
#=====================================================
# Merge model outputs for Table S8
model_results <- bind_rows(af_results, at_results, am_results, as_results)

table_results_s8 <- model_results %>%
  filter(Predictor != "(Intercept)") %>%
  transmute(
    Species = paste0("*", gsub("_", " ", species), "*"),
    Predictor = recode(
      Predictor,
      daily_rainfall = "Daily rainfall",
      mean_temp_daily = "Mean daily temperature",
      sex = "Sex (male)",
      svl = "Body size (SVL)"),
    `Estimate (β)` = round(Estimate, 3),
    SE = round(`Std. Error`, 3),
    `Wald χ²` = round(Chisq, 3),
    `p-value` = round(`Pr(>Chisq)`, 3),
    Sig = case_when(
      `Pr(>Chisq)` < 0.001 ~ "***",
      `Pr(>Chisq)` < 0.01 ~ "**",
      `Pr(>Chisq)` < 0.05 ~ "*",
      TRUE ~ ""))

table_results_s8 %>%
  gt() %>%
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_column_labels()) %>%
  fmt_markdown(columns = Species)

#=====================================================================================================
# Supplementary Figure S3
# Predicted effects of rainfall and temperature on daily movement distance across poison frog species
#=====================================================================================================
supp_fig_s3 <- (plot_af | plot_at) /(plot_am | plot_as)
# ggsave(
#   filename = "Figures/Supplementary_Figure_S3.tiff",
#   plot = supp_fig_s3,
#   width = 260,
#   height = 180,
#   units = "mm",
#   dpi = 600,
#   compression = "lzw"
# )

Microhabitat use

To assess seasonal differences in microhabitat use, we fitted generalized linear mixed-effects models (GLMMs) with negative binomial error distribution (nbinom2) using the glmmTMB() function from the glmmTMB R package (Brooks et al. 2026) for each species. Full model details are in Supplementary Methods S11.

# Microhabitat use models

# Subset data
af_microh_data <- microhabitat_data %>% filter(species == "Allobates femoralis")
at_microh_data <- microhabitat_data %>% filter(species == "Ameerega trivittata")
am_microh_data <- microhabitat_data %>% filter(species == "Ameerega macero")
as_microh_data <- microhabitat_data %>% filter(species == "Ameerega shihuemoy")

#=====================================================================================================
# Allobates femoralis

# Counts and proportions
af_counts_prop <- prepare_microhabitat_data(af_microh_data)

# Remove microhabitat variables not present in all combinations
af_filtered <- remove_problem_levels(af_counts_prop, 
                                     ~ season + sex + catg_subtrates)

# Selected negative binomial model 
model_af <- glmmTMB(
  count ~ season * sex * catg_subtrates + (1 | ID_new),
  family = nbinom2(),
  offset = log(total_obs),
  data = af_filtered)

# Main effects
main_effects_af <- glmmTMB:::Anova.glmmTMB(model_af, type = "III")

# PostHoc
# Important substrates to keep (only the ones with sig in at least one of the species)
keep_subs <- c("On_ll", "Treefall_refuge", "Under_ll", "Under_log",
               "Ins_roots", "Ll_refuge", "On_log", "Soil_rock_refuge")

contr_af <- extract_posthoc(model_af, by_vars = c("sex", "catg_subtrates"),
  keep_subs = keep_subs)

contr_breeding_af_female <- contr_af %>% filter(sex == "f")

contr_breeding_af_male <- contr_af %>% filter(sex == "m")

#=====================================================================================================
# Ameerega trivittata

# Counts and proportions
at_counts_prop <- prepare_microhabitat_data(at_microh_data)

# Remove microhabitat variables not present in all combinations
at_filtered <- remove_problem_levels(at_counts_prop, 
                                     ~ season + sex + catg_subtrates)

# Selected negative binomial model 
model_at <- glmmTMB(
  count ~ season * sex * catg_subtrates + (1 | ID_new),
  family = nbinom2(),
  offset = log(total_obs),
  data = at_filtered)

# Main effects
main_effects_at <- glmmTMB:::Anova.glmmTMB(model_at, type = "III")

# PostHoc
# Important substrates to keep (only the ones with sig in at least one of the species)
contr_at <- extract_posthoc(model_at, by_vars = c("catg_subtrates"),
                            keep_subs = keep_subs)

#=====================================================================================================
# Ameerega macero

# Counts and proportions
am_counts_prop <- prepare_microhabitat_data(am_microh_data)

# Remove microhabitat variables not present in all combinations
am_filtered <- remove_problem_levels(am_counts_prop, 
                                     ~ season + sex + catg_subtrates)

# Selected negative binomial model 
model_am <- glmmTMB(
  count ~ season * sex * catg_subtrates + (1 | ID_new),
  family = nbinom2(),
  offset = log(total_obs),
  data = am_filtered)

# Main effects
main_effects_am <- glmmTMB:::Anova.glmmTMB(model_am, type = "III")

# PostHoc
# Important substrates to keep (only the ones with sig in at least one of the species)
contr_am <- extract_posthoc(model_am, by_vars = c("catg_subtrates"),
                            keep_subs = keep_subs)

#=============================================================================================
# Ameerega shihuemoy
# Counts and proportions
as_counts_prop <- prepare_microhabitat_data(as_microh_data)

# Remove microhabitat variables not present in all combinations
as_filtered <- remove_problem_levels(as_counts_prop, 
                                     ~ season + sex + catg_subtrates)

# Selected negative binomial model 
model_as <- glmmTMB(
  count ~ season * sex * catg_subtrates + (1 | ID_new),
  family = nbinom2(),
  offset = log(total_obs),
  data = as_filtered)

# Main effects
main_effects_as <- glmmTMB:::Anova.glmmTMB(model_as, type = "III")

# PostHoc
# Important substrates to keep (only the ones with sig in at least one of the species)
contr_as <- extract_posthoc(model_as, by_vars = c("catg_subtrates"),
                            keep_subs = keep_subs)

#=============================================================================================
# Combine Posthoc results
contr_af_female <- contr_breeding_af_female %>%
  mutate(sp = "Allobates femoralis ♀") %>% select(-sex)

contr_af_male <- contr_breeding_af_male %>%
  mutate(sp = "Allobates femoralis ♂") %>% select(-sex)

contr_at <- contr_at %>%
  mutate(sp = "Ameerega trivittata")

contr_am <- contr_am %>%
  mutate(sp = "Ameerega macero")

contr_as <- contr_as %>%
  mutate(sp = "Ameerega shihuemoy")


combined_df <- bind_rows(
  contr_af_female, contr_af_male, contr_at, contr_am, contr_as)

#=====================================================================================================
# Plot Figure 3
# Seasonal shifts in microhabitat use across four poison frog species
#=====================================================================================================

# Formatting labels by storing them in objects
substrate_labels <- c(
  "On_ll" = "On leaf litter",
  "Under_ll" = "Under leaf litter",
  "Treefall_refuge" = "Treefall refuge",
  "Under_log" = "Under log",
  "Ins_roots" = "Inside roots",
  "Ll_refuge" = "Potential clutch site",
  "On_log" = "On log",
  "Soil_rock_refuge" = "Rock refuge")

species_order <- c(
  "Allobates femoralis ♀",
  "Allobates femoralis ♂",
  "Ameerega trivittata",
  "Ameerega macero",
  "Ameerega shihuemoy")

facet_labels <- c(
  "Allobates femoralis ♀" = "Allobates\nfemoralis ♀",
  "Allobates femoralis ♂" = "Allobates\nfemoralis ♂",
  "Ameerega trivittata"   = "Ameerega\ntrivittata",
  "Ameerega macero"       = "Ameerega\nmacero",
  "Ameerega shihuemoy"    = "Ameerega\nshihuemoy")

# Change labels
combined_df <- combined_df %>% mutate(
  catg_subtrates = recode(catg_subtrates, !!!substrate_labels), #!!!inserts each named element of substrate_labels as its own argument
  sp = factor(sp, levels = species_order))

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

# Facet Strip Styling
strip_bg <- strip_themed(
  background_x = elem_list_rect(
    fill = c("#CBC3E3", "#CBC3E3", "#CBC3E3",
             "#b5d1a1", "#b5d1a1")),
  text_x = elem_list_text(face = "italic", size = 12))

#Plot
plot <- combined_df %>%
  ggplot(aes(y = catg_subtrates,
             x = log_ratio,
             color = signif,
             shape = nonEst)) +
  annotate("rect",
           xmin = -Inf, xmax = 0,
           ymin = -Inf, ymax = Inf,
           fill = "#0072B2", alpha = 0.6) +
  annotate("rect",
           xmin = 0, xmax = Inf,
           ymin = -Inf, ymax = Inf,
           fill = "#E7B800", alpha = 0.6) +
  geom_point(size = 3) +
  geom_errorbar(aes(xmin = lower, xmax = upper),
                 height = 0.2,
                 linewidth = 0.7) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray17") +
  scale_shape_manual(values = c(`FALSE` = 16, `TRUE` = 1)) +
  scale_color_manual(values = c(`FALSE` = "gray46", `TRUE` = "black"))+
  labs(
    x = "Wet season ← 0 → Dry season", 
    y = ""  ) +
  ggh4x::facet_grid2(. ~ sp, scales = "free_y", 
                     labeller = labeller(sp = facet_labels),
                     strip = strip_bg) +
  theme_classic()+
  theme(legend.position = "none",
        axis.text.y = element_text(size = 12),   # increase y-axis label size
        axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14),
        panel.spacing = unit(1, "lines"),
        strip.text = element_text(face = "bold.italic", size = 12))

#=============================================================================================
# Table S9
# Type III Wald χ² tests from generalized linear mixed models for microhabitat 
# (negative binomial, log offset)
#=============================================================================================

af_df <- extract_typeIII_main_effects(main_effects_af, "Allobates femoralis")
at_df <- extract_typeIII_main_effects(main_effects_at, "Ameerega trivittata")
am_df <- extract_typeIII_main_effects(main_effects_am, "Ameerega macero")
as_df <- extract_typeIII_main_effects(main_effects_as, "Ameerega shihuemoy")

final_table <- bind_rows(af_df, at_df, am_df, as_df)

#Format table names
final_table <- final_table %>%
  mutate(
    Effect = recode(Effect,
                    "(Intercept)" = "(Intercept)",
                    "season" = "season",
                    "sex" = "sex",
                    "catg_subtrates" = "substrate categories",
                    "season:sex" = "season*sex",
                    "season:catg_subtrates" = "season*substrate categories",
                    "sex:catg_subtrates" = "sex*substrate categories",
                    "season:sex:catg_subtrates" = "season*sex*substrate categories"),
    Species = paste0("*", Species, "*"))  # italics

final_table %>%
  gt(groupname_col = "Species") %>%
  fmt_markdown(columns = Species) %>%
  cols_label(
    Effect = "Effect",
    Chi2 = "χ²",
    df = "df",
    p_value = "p-value",
    Sig = "") %>%
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_column_labels())%>%
  tab_source_note(
    source_note = "*** p < 0.001, ** p < 0.01, * p < 0.05, . p < 0.1")


#=======================================================================================
# Table S10
# Post-Hoc pairwise comparison of substrate use between dry and wet 
# seasons for species with   # significant interactions
#=======================================================================================

sig_table <- combined_df %>% filter(signif == TRUE)

#Separate species and sex
sig_table <- sig_table %>%
  mutate(
    Sex = case_when(
      str_detect(sp, "♀") ~ "female",
      str_detect(sp, "♂") ~ "male",
      TRUE ~ "—"), 
    Species = str_remove(sp, " ♀| ♂"))

# Format p-values and significance codes
sig_table <- sig_table %>%
  mutate(
    p_label = case_when(
      p.value < 0.001 ~ "<0.001",
      TRUE ~ sprintf("%.4f", round(p.value, 4))
    ),
    Sig = case_when(
      p.value < 0.001 ~ "***",
      p.value < 0.01  ~ "**",
      p.value < 0.05  ~ "*",
      p.value < 0.1   ~ ".",
      TRUE ~ ""))

# Create formatted table
final_sig_table <- sig_table %>%
  transmute(
    Species = paste0("*", Species, "*"),
    Substrate = catg_subtrates,
    Sex = Sex,
    Contrast = contrast,
    Ratio = round(ratio, 3),
    SE = round(SE, 3),
    z = round(z.ratio, 3),
    `p-value` = p_label,
    Sig = Sig)

# Table
final_sig_table %>%
  gt(groupname_col = "Species") %>%
  fmt_markdown(columns = Species) %>%
  cols_label(
    Sig = "") %>%
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_column_labels()) %>%
  tab_source_note(
    source_note = "*** p < 0.001, ** p < 0.01, * p < 0.05, . p < 0.1")

Water-borne hormone release rates

We assessed seasonal differences in hormonal release rates using species-specific models. Testosterone and corticosterone release rates (pg/h) were log-transformed and analyzed with linear models, except for A. macero (testosterone) and A. trivittata (corticosterone), which required Gamma GLMs with log links. For species with data for both sexes, we assessed whether seasonal effects depended on sex while accounting for body size (log(hormone) ~ season*sex + SVL); otherwise, we only included season and SVL. Significant effects were evaluated using estimated marginal means with FDR-corrected comparisons. Results were back-transformed for visualization. Additionally, we assessed the influence of testosterone release rates, sex and their interaction on space use area by using linear models for each species.

#Testosterone 

#Subset data per species
testo_data_af <- testosterone_data %>% filter(species =="Allobates femoralis")
testo_data_at <- testosterone_data %>% filter(species =="Ameerega trivittata")
testo_data_am <- testosterone_data %>% filter(species =="Ameerega macero")
testo_data_as <- testosterone_data %>% filter(species =="Ameerega shihuemoy")

#=======================================================================================
# Allobates femoralis

#hist(log(testo_data_af$mean_value))
# Selected linear model
lm1 <- lm(log(mean_value) ~ season * sex + svl, data = testo_data_af)

af_model_output_df <- summary(lm1)$coefficients %>% 
  as.data.frame() %>%
  tibble::rownames_to_column("Predictor") %>%
  dplyr::mutate(species = "Allobates femoralis")

# Main effects
car::Anova(lm1, type = "III") 

# Post-hoc - see Functions.R for code
af_results <- get_species_emm(lm1, "Allobates femoralis")
af_results_df <- as.data.frame(af_results)

#=======================================================================================
# Ameerega trivittata

#hist(log(testo_data_at$mean_value))
# Selected linear model
lm1 <- lm(log(mean_value) ~ season * sex + svl, data = testo_data_at)

at_model_output_df <- summary(lm1)$coefficients %>% 
  as.data.frame() %>%
  tibble::rownames_to_column("Predictor") %>%
  dplyr::mutate(species = "Ameerega trivittata")

# Main effects
car::Anova(lm1, type = "III") 

# Post-hoc - see Functions.R for code
at_results <- get_species_emm(lm1, "Ameerega trivittata")
at_results_df <- as.data.frame(at_results)

#=======================================================================================
# Ameerega macero

#hist(log(testo_data_am$mean_value)) #doesnt follow normality
# Selected gamma generalized linear model with log link, which is appropriate for positive continuous data
glm_am <- glm(
  mean_value ~ season*sex + svl,
  family = Gamma(link = "log"),
  data = testo_data_am)

am_model_output_df <- summary(glm_am)$coefficients %>% 
  as.data.frame() %>%
  tibble::rownames_to_column("Predictor") %>%
  dplyr::mutate(species = "Ameerega macero")

# Main effects
car::Anova(glm_am, type = "III") #non-significant effects

#Post-Hoc - This is just for the plot
am_results <- get_species_emm(glm_am, "Ameerega macero")
am_results_df <- as.data.frame(am_results)

#=======================================================================================
# Ameerega shihuemoy

#hist(log(testo_data_as$mean_value))
# Selected linear model
lm1 <- lm(log(mean_value) ~ season + svl, data = testo_data_as)
as_model_output_df <- summary(lm1)$coefficients %>% 
  as.data.frame() %>%
  tibble::rownames_to_column("Predictor") %>%
  dplyr::mutate(species = "Ameerega shihuemoy")

# Main effects
car::Anova(lm1, type = "II")

#PostHoc - This is just for the plot
lm2 <- lm(log(mean_value) ~ season*sex + svl, data = testo_data_as)
as_results <- get_species_emm(lm2, "Ameerega shihuemoy")
as_results_df <- as.data.frame(as_results)

#=======================================================================================
# Figure 4A
# Seasonal variation in testosterone release rates across four poison frog species
#=======================================================================================

# Merge all testosterone data for the plot
results_list <- list(
  af_results, at_results, am_results, as_results)
all_emmeans_df <- bind_rows(map(results_list, "emmeans"))
all_sig_df     <- bind_rows(map(results_list, "sig"))

# Plot
plot_testosterone <- plot_hormone_emmeans(
  all_emmeans_df,
  all_sig_df,
  ylab = expression("Testosterone levels (pg H"^{-1}*")"))

#=======================================================================================
# Table S11
# Species-specific ANOVA/GLM results for testosterone concentrations
#=======================================================================================
#Merge data for Table S11
model_output <- bind_rows(
  af_model_output_df,
  at_model_output_df,
  am_model_output_df,
  as_model_output_df)

table_results_s11 <- model_output %>%
  filter(Predictor != "(Intercept)")%>%
  transmute(
    Species = paste0("*", species, "*"), #this is for italics
    Predictor = recode(
      Predictor,
      seasonwet = "Season (wet)",
      sexMale = "Sex (male)",
      svl = "SVL",
      `seasonwet:sexMale` = "Season (wet) × Sex (male)"),
    β = Estimate,
    SE = round(`Std. Error`, 3),
    t = round(`t value`, 3),
    p = round(`Pr(>|t|)`, 3),
    Sig = case_when(
      `Pr(>|t|)` < 0.001 ~ "***",
      `Pr(>|t|)` < 0.01 ~ "**",
      `Pr(>|t|)` < 0.05 ~ "*",
      TRUE ~ ""))

table_results_s11 %>%
  gt() %>%
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_column_labels()) %>%
  fmt_markdown(columns = Species) %>% #this is for italics
  fmt_number(
    columns = c(β, SE, t, p),
    decimals = 3) %>%
  cols_label(
    Sig = "")

# Models assessing the effect of testosterone on space use
  
#Subset data
space_hormones_af <- spaceuse_hormones%>%filter(species == "Allobates femoralis")
space_hormones_at <- spaceuse_hormones%>%filter(species == "Ameerega trivittata")
space_hormones_am <- spaceuse_hormones%>%filter(species == "Ameerega macero")
space_hormones_as <- spaceuse_hormones%>%filter(species == "Ameerega shihuemoy")

#=======================================================================================
# Allobates femoralis
# Selected model
af_model <- lm(log(HPI95) ~ t_mean_value * sex, data = space_hormones_af)
af_model_output_df <- summary(af_model)$coefficients %>% 
  as.data.frame() %>%
  tibble::rownames_to_column("Predictor") %>%
  dplyr::mutate(species = "Allobates femoralis")
car::Anova(af_model, type = 3)

#=======================================================================================
# Ameerega trivittata
# Selected model
at_model <- lm(log(HPI95) ~ t_mean_value * sex, data = space_hormones_at)
at_model_output_df <- summary(at_model)$coefficients %>% 
  as.data.frame() %>%
  tibble::rownames_to_column("Predictor") %>%
  dplyr::mutate(species = "Ameerega trivittata")
car::Anova(at_model, type = 3)

#=======================================================================================
# Ameerega macero
# Selected model
am_model <- lm(log(HPI95) ~ t_mean_value * sex, data = space_hormones_am)
am_model_output_df <- summary(am_model)$coefficients %>% 
  as.data.frame() %>%
  tibble::rownames_to_column("Predictor") %>%
  dplyr::mutate(species = "Ameerega macero")
car::Anova(am_model, type = 3)


#=======================================================================================
# Ameerega shihuemoy
# Selected model
as_model <- lm(log(HPI95) ~ t_mean_value + sex, data = space_hormones_as)
as_model_output_df <- summary(as_model)$coefficients %>% 
  as.data.frame() %>%
  tibble::rownames_to_column("Predictor") %>%
  dplyr::mutate(species = "Ameerega shihuemoy")
car::Anova(as_model, type = 2)

#==============================================================================================
# Table S12
# Predictors of space use (HPI95) from testosterone levels. 
# Linear models testing effects of testosterone and sex on log-transformed HPI95 across species
#==============================================================================================
#Merge data for Table S12
model_output <- bind_rows(
  af_model_output_df,
  at_model_output_df,
  am_model_output_df,
  as_model_output_df)

table_results_s12 <- model_output %>%
  filter(Predictor != "(Intercept)")%>%
  transmute(
    Species = paste0("*", species, "*"), #this is for italics
    Predictor = recode(
      Predictor,
      t_mean_value = "Testosterone",
      sexm = "Sex (male)",
      `t_mean_value:sexm` = "Testosterone x Sex (male)"),
    β = Estimate,
    SE = round(`Std. Error`, 3),
    t = round(`t value`, 3),
    p = round(`Pr(>|t|)`, 3),
    Sig = case_when(
      `Pr(>|t|)` < 0.001 ~ "***",
      `Pr(>|t|)` < 0.01 ~ "**",
      `Pr(>|t|)` < 0.05 ~ "*",
      TRUE ~ ""))

table_results_s12 %>%
  gt() %>%
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_column_labels()) %>%
  fmt_markdown(columns = Species) %>% #this is for italics
  fmt_number(
    columns = c(β, SE, t, p),
    decimals = 3) %>%
  cols_label(
    Sig = "")

# Corticosterone

#Subset data per species
cort_data_af <- corticosterone_data %>% filter(species =="Allobates femoralis")
cort_data_at <- corticosterone_data %>% filter(species =="Ameerega trivittata")
cort_data_am <- corticosterone_data %>% filter(species =="Ameerega macero")
cort_data_as <- corticosterone_data %>% filter(species =="Ameerega shihuemoy")

#=======================================================================================
# Allobates femoralis

#hist(log(hormone_data_af$mean_value))
# Selected linear model
lm1 <- lm(log(mean_value) ~ season * sex + svl, data = cort_data_af)
af_model_output_df <- summary(lm1)$coefficients %>% 
  as.data.frame() %>%
  tibble::rownames_to_column("Predictor") %>%
  dplyr::mutate(species = "Allobates femoralis")

# Main effects
car::Anova(lm1, type = "III") 

# Post-hoc - see Functions.R for code
af_results <- get_species_emm(lm1, "Allobates femoralis")
af_results_df <- as.data.frame(af_results)

#=======================================================================================
# Ameerega trivittata

#hist(log(cort_data_at$mean_value)) #doesnt follow normality
# Selected gamma generalized linear model with log link, which is appropriate for positive continuous data
glm_at <- glm(
  mean_value ~ season*sex + svl,
  family = Gamma(link = "log"),
  data = cort_data_at)

at_model_output_df <- summary(glm_at)$coefficients %>% 
  as.data.frame() %>%
  tibble::rownames_to_column("Predictor") %>%
  dplyr::mutate(species = "Ameerega trivittata")

# Main effects
car::Anova(glm_at, type = "III") #non-significant effects

# Post-hoc - see Functions.R for code
at_results <- get_species_emm(glm_at, "Ameerega trivittata")
at_results_df <- as.data.frame(at_results)


#=======================================================================================
# Ameerega macero

#hist(log(cort_data_am$mean_value)) #doesnt follow normality
# Selected linear model
lm1 <- lm(log(mean_value) ~ season * sex + svl, data = cort_data_am)
am_model_output_df <- summary(lm1)$coefficients %>% 
  as.data.frame() %>%
  tibble::rownames_to_column("Predictor") %>%
  dplyr::mutate(species = "Ameerega macero")

# Main effects
car::Anova(lm1, type = "III")

#PostHoc - This is just for the plot
am_results <- get_species_emm(lm1, "Ameerega macero")
am_results_df <- as.data.frame(am_results)

#=======================================================================================
# Ameerega shihuemoy

#hist(log(testo_data_as$mean_value))
# Selected linear model
lm1 <- lm(log(mean_value) ~ season + svl, data = cort_data_as)
as_model_output_df <- summary(lm1)$coefficients %>% 
  as.data.frame() %>%
  tibble::rownames_to_column("Predictor") %>%
  dplyr::mutate(species = "Ameerega shihuemoy")

# Main effects
car::Anova(lm1, type = "II")

#PostHoc - This is just for the plot
lm2 <- lm(log(mean_value) ~ season*sex + svl, data = cort_data_as)
as_results <- get_species_emm(lm2, "Ameerega shihuemoy")
as_results_df <- as.data.frame(as_results)

#=======================================================================================
# Figure 4B
# Seasonal variation in Corticosterone release rates across four poison frog species
#=======================================================================================

# Merge all testosterone data for the plot
results_list <- list(
  af_results, at_results, am_results, as_results)
all_emmeans_df <- bind_rows(map(results_list, "emmeans"))
all_sig_df     <- bind_rows(map(results_list, "sig"))

# Plot
plot_corticosterone <- plot_hormone_emmeans(
  all_emmeans_df,
  all_sig_df,
  ylab = expression("Corticosterone levels (pg H"^{-1}*")"))

#=======================================================================================
# Table S13
# Species-specific ANOVA/GLM results for corticosterone concentrations
#=======================================================================================
#Merge data for Table S3
model_output <- bind_rows(
  af_model_output_df,
  at_model_output_df,
  am_model_output_df,
  as_model_output_df
)

table_results_s13 <- model_output %>%
  filter(Predictor != "(Intercept)")%>%
  transmute(
    Species = paste0("*", species, "*"), #this is for italics
    Predictor = recode(
      Predictor,
      seasonwet = "Season (wet)",
      sexMale = "Sex (male)",
      svl = "SVL",
      `seasonwet:sexMale` = "Season (wet) × Sex (male)"),
    β = Estimate,
    SE = round(`Std. Error`, 3),
    t = round(`t value`, 3),
    p = round(`Pr(>|t|)`, 3),
    Sig = case_when(
      `Pr(>|t|)` < 0.001 ~ "***",
      `Pr(>|t|)` < 0.01 ~ "**",
      `Pr(>|t|)` < 0.05 ~ "*",
      TRUE ~ ""))

table_results_s13 %>%
  gt() %>%
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_column_labels()) %>%
  fmt_markdown(columns = Species) %>% #this is for italics
  fmt_number(
    columns = c(β, SE, t, p),
    decimals = 3) %>%
  cols_label(
    Sig = "")

Alkaloids

We normalized alkaloid peak intensities using nicotine internal standard (IS) (Grant et al. 2012; Caty et al. 2025). We then grouped alkaloids into structural families based on revised annotations (Supplementary Table S2). We quantified alkaloid family-level abundance (ng), alkaloid richness (number of detected alkaloids), and total alkaloid abundance for each individual frog. We also quantified detection frequency per alkaloid family within species and seasons. To evaluate seasonal differences in overall alkaloid composition, we constructed a Bray–Curtis dissimilarity matrices using alkaloid family-level abundances and performed species-specific permutational multivariate analyses of variance (PERMANOVA) using the adonis2() function in the vegan R package (Oksanen et al. 2025) with 999 permutations. Homogeneity of multivariate dispersion was evaluated using betadisper(). To assess seasonal differences in individual alkaloid families within species, we used Wilcoxon rank-sum tests and adjusted p-values for multiple comparisons using the Benjamini–Hochberg procedure. For visualization, we summarized alkaloid abundance using mean abundance and 95% confidence intervals and included individual-level observations for frogs with detectable alkaloids.

#Merge alkaloid data with alkaloid metadata
alk_data <- alk_data %>% left_join(metadata %>% select(Sample, ID), by = "Sample")

#Normalize using alkaloid standard
alk_data_normalized <- alk_data %>%
  mutate(across(12:38, ~ ifelse(.x < 150000, 0, .x)))%>% # convert to 0 all values less than 150000
  mutate(across(12:38, ~ .x / `476`))%>%
  mutate(across(13:38, ~ .x * 10)) #convert to nanograms

#Count the number of features per sample
counts <- alk_data_normalized %>%
  mutate(n_features = rowSums(across(13:38, ~ . > 0), na.rm = TRUE))

#Dataset including only samples with alkaloids
samples_alkaloids <- counts %>% filter(n_features >0)


#Quantities - abundance and richness per sample
frog_summary <- samples_alkaloids %>%
  mutate(
    total_ng = rowSums(select(., 13:38), na.rm = TRUE),
    alkaloid_richness = rowSums(select(., 13:38) > 0, na.rm = TRUE))

#Count number of features per species and season
alk_cols <- c("552", "619", "621", "655", "689", "709", "788", "811", "812", "822",
              "838", "845", "886", "895", "899", "907", "921", "939", "944", "982",
              "998", "999", "1003", "1011", "1012", "1015")

# Number of alkaloids/features per species and season
population_richness <- frog_summary %>%
  group_by(ATTRIBUTE_season_species) %>%
  summarise(
    n_alkaloids = sum(colSums(across(all_of(alk_cols)) > 0) > 0),
    .groups = "drop")

# Sample size per species and season
sample_size_species_season <-frog_summary%>%
  group_by(ATTRIBUTE_season_species)%>%
  summarise(n())

# Number of alkaloids per species
population_richness_sp <- frog_summary %>%
  group_by(ATTRIBUTE_species) %>%
  summarise(
    n_unique_alkaloids = sum(colSums(across(all_of(alk_cols)) > 0) > 0),
    .groups = "drop")

# Mean chemical load per species and season
species_season_summary <- frog_summary %>%
  group_by(ATTRIBUTE_species, ATTRIBUTE_season) %>%
  summarise(
    mean_total_ng = mean(total_ng, na.rm = TRUE),
    sd_total_ng   = sd(total_ng, na.rm = TRUE),
    mean_richness = mean(alkaloid_richness, na.rm = TRUE),
    sd_richness   = sd(alkaloid_richness, na.rm = TRUE),
    n_frogs       = n(),
    .groups = "drop")

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

# Seasonal comparison between alkaloid families
season_abund_feature <- samples_alkaloids %>% #samples_alkaloids includes only samples with alkaloids
  select(-c(3,5:12,40)) %>% 
  pivot_longer(
    cols = where(is.numeric),
    names_to = "alkaloid",
    values_to = "abundance") %>%
  mutate(alkaloid = as.numeric(alkaloid)) %>%
  left_join(
    select_alk %>% select(alkaloid, Mabel_family),
    by = "alkaloid") %>%
  mutate(
    season  = ATTRIBUTE_season,
    species = ATTRIBUTE_species)

season_abund_feature

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

# Prepare data for analysis

# 1) Collapse alkaloid families within each frog

# Some frogs may contain multiple alkaloids belonging to the same structural family.
# Here, we sum them WITHIN each frog so that each frog has one abundance value per alkaloid family.
profile_summary_ng_collapsed <- season_abund_feature %>%
  group_by(Sample,
           ATTRIBUTE_season,
           ATTRIBUTE_species,
           Mabel_family) %>%
  summarise(
    abundance_fam = sum(abundance, na.rm = TRUE),
    .groups = "drop")

# 2) Create a wide frog x family matrix
# Each row = one frog
# Each column = one alkaloid family
# Values = abundance (ng)
frog_wide <- profile_summary_ng_collapsed %>%
  pivot_wider(names_from  = Mabel_family,
              values_from = abundance_fam,
              values_fill = 0)

# 3) Extract metadata
metadata <- frog_wide %>%
  select(Sample, ATTRIBUTE_species, ATTRIBUTE_season) %>%
  mutate(ATTRIBUTE_season = factor(ATTRIBUTE_season,levels = c("dry", "wet")),
         ATTRIBUTE_species = factor(ATTRIBUTE_species),
         sp_season = paste(ATTRIBUTE_species, ATTRIBUTE_season, sep = "_"))

# 4) Abundance matrix
# Remove metadata columns and keep only alkaloid abundances
alk_matrix <- frog_wide %>% select(-Sample, -ATTRIBUTE_species, -ATTRIBUTE_season) %>%
  as.data.frame()
rownames(alk_matrix) <- frog_wide$Sample

# 5) Remove frogs with no alkaloids
# Frogs with all zeros cannot be used in Bray-Curtis analyses
keep <- rowSums(alk_matrix) > 0
alk_matrix <- alk_matrix[keep, ]
metadata <- metadata[keep, ]

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

# Statistical test

#=======================================================================================
# Table S14
# Species-specific PERMANOVA testing for seasonal differences in alkaloid composition 
# using Bray–Curtis dissimilarities of alkaloid family abundances (ng)
#=======================================================================================

# Permanova (within species), all my analysis are done by species
species_list <- unique(metadata$ATTRIBUTE_species)
permanova_results <- lapply(species_list, function(sp) {
  samples_sp <- metadata %>%
    filter(ATTRIBUTE_species == sp) %>%
    pull(Sample)
  adonis2(
    alk_matrix[samples_sp, ] ~ ATTRIBUTE_season,
    data = metadata %>% filter(Sample %in% samples_sp),
    method = "bray",
    permutations = 999)
})

names(permanova_results) <- species_list

# Table results
table_results_s14 <- bind_rows(
  lapply(names(permanova_results), function(sp) {
    res <- permanova_results[[sp]]
    tibble(
      Species = paste0("*", sp, "*"),  # markdown italics
      Df = res["Model", "Df"],
      F = res["Model", "F"],
      R2 = res["Model", "R2"],
      P_value = res["Model", "Pr(>F)"]
    )}))

table_results_s14 <- table_results_s14 %>%
  mutate(Species = factor(Species,
                          levels = c(
                            "*Ameerega trivittata*",
                            "*Ameerega macero*",
                            "*Ameerega shihuemoy*"))) %>% arrange(Species)

table_results_s14 %>%
  gt() %>%
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_column_labels()) %>%
  fmt_markdown(columns = Species) %>%
  fmt_number(
    columns = c(F, R2, P_value),
    decimals = 3) %>%
  cols_label(
    Species = "Species",
    Df = "Df",
    F = "F",
    R2 = html("R<sup>2</sup>"),
    P_value = "P-value")

#=======================================================================================
# Table S15
# Seasonal comparisons of alkaloid family abundance within each species. Differences between 
# dry and wet seasons were assessed using pairwise Wilcoxon rank-sum tests. Reported values 
# include sample sizes (n), test statistic (W), raw and Benjamini–Hochberg-adjusted p-values
#=======================================================================================

# Univariate analysis (families)
# Wilcoxon tests
family_tests <- profile_summary_ng_collapsed %>%
  group_by(ATTRIBUTE_species, Mabel_family) %>%
  filter(sum(abundance_fam) > 0) %>%   # remove all-zero groups
  wilcox_test(
    abundance_fam ~ ATTRIBUTE_season) %>%
  adjust_pvalue(method = "BH") %>%
  add_significance("p.adj")

# Table results
table_results_s15 <- family_tests %>%
  transmute(
    Species = paste0("*", ATTRIBUTE_species, "*"),
    `Alkaloid family` = Mabel_family,
    Dry = n1,
    Wet = n2,
    W = statistic,
    `P-value` = round(p, 3),
    `P-adjusted` = round(p.adj, 3))

table_results_s15 <- table_results_s15 %>%
  mutate(Species = factor(Species,
                          levels = c(
                            "*Ameerega trivittata*",
                            "*Ameerega macero*",
                            "*Ameerega shihuemoy*"))) %>% arrange(Species)
table_results_s15 %>%
  gt() %>%
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_column_labels()) %>%
  fmt_markdown(columns = Species) #this is for italics


#=======================================================================================
# Figure 5B
# Mean abundance (ng) of alkaloid families within each species–season combination (sample sizes indicated as samples with detectable alkaloids/total sampled). 
#=======================================================================================

# Plot settings

species_levels <- c("Ameerega trivittata",
                    "Ameerega macero",
                    "Ameerega shihuemoy")

family_levels <- c("HTX", "DHQ", "N-methyl-DHQ",
                   "5,8-I", "5,6,8-I", "3,5-P",
                   "Pyr", "Tricyclic", "Unclass", "Unknown")

species_colors <- c("Ameerega trivittata" = "#CBC3E3",
                    "Ameerega macero"     = "#b5d1a1",
                    "Ameerega shihuemoy"  = "#b5d1a1")

season_colors <- c(dry = "#E7B800", wet = "#0072B2")
season_alpha  <- c(dry = 0.4, wet = 0.9)

# Family palette
family_palette <- wes_palette("Zissou1", length(family_levels), type = "continuous")

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

# Family-level summary statistics (mean abundance and 95% confidence intervals)

profile_summary_ng <- profile_summary_ng_collapsed %>%
  group_by(ATTRIBUTE_season, ATTRIBUTE_species, Mabel_family) %>%
  summarise(n = n_distinct(Sample), 
            mean_ng = mean(abundance_fam, na.rm = TRUE),
            sd_ng = sd(abundance_fam,na.rm = TRUE),
            se_ng = sd_ng / sqrt(n), 
            ci_ng = qt(0.975, df = n - 1) * sd_ng / sqrt(n),
            .groups = "drop")
#=================================================================

# Factor ordering

profile_summary_ng <- profile_summary_ng %>%
  mutate(Mabel_family = factor(Mabel_family, levels = family_levels),
         ATTRIBUTE_species = factor(ATTRIBUTE_species, levels = species_levels))

profile_summary_ng_collapsed <- profile_summary_ng_collapsed %>%
  mutate(ATTRIBUTE_species = factor(ATTRIBUTE_species, levels = species_levels))

# Individual observations shown as jittered points
raw_points <- profile_summary_ng_collapsed %>% filter(abundance_fam > 0)

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

# Plot mean abundance by alkaloid family

plot_families <- ggplot(
  profile_summary_ng,
  aes(x = Mabel_family, y = mean_ng, fill = Mabel_family)) +
  # mean abundance bars
  geom_col(aes(alpha = ATTRIBUTE_season), position = position_dodge(width = 0.9), 
           color = "black") +
  # family colors
  scale_fill_manual(values = family_palette, guide = "none") +
  ggnewscale::new_scale_fill() +
  # season aesthetics
  scale_fill_manual(values = season_colors, name = "Season")+
  scale_alpha_manual(values = season_alpha, name = "Season") +
  # 95% confidence intervals
  geom_errorbar(aes(group = ATTRIBUTE_season, 
                    ymin = pmax(0, mean_ng - ci_ng),
                    ymax = mean_ng + ci_ng),
                position = position_dodge(width = 0.9), 
                width = 0.25, color = "black") +
  # Individual sample abundances
  geom_jitter(data = raw_points,
              aes(x = Mabel_family, y = abundance_fam, fill = ATTRIBUTE_season), 
              position = position_jitterdodge(jitter.width = 0.15, 
                                              dodge.width = 0.9),
              shape = 21, color = "black", stroke = 0.4, alpha = 0.5, size = 2.2,
              inherit.aes = FALSE) +
  # Species panels
  ggh4x::facet_wrap2(~ ATTRIBUTE_species, ncol = 1,scales = "fixed",
                     strip = ggh4x::strip_themed(
                       background_x = elem_list_rect(fill = species_colors,
                                                     alpha = 0.5), 
                       text_x = elem_list_text(face = "italic", size = 14))) +
  # Axes
  scale_y_sqrt(limits = c(0, 70), labels = comma) +
  labs(
    x = "Alkaloid Family",
    y = "Mean abundance (ng)") +
  theme_bw() +
  theme(
    axis.title = element_text(size = 16),
    axis.text.x = element_text(size = 12, angle = 45, hjust = 1, vjust = 1),
    axis.text.y = element_text(size = 12),
    strip.text = element_text(face = "italic"),
    legend.title = element_text(size = 14),
    legend.text = element_text(size = 12),
    panel.grid.minor = element_blank())

#=======================================================================================
# Figure S4
# Heatmap showing the prevalence (%) of alkaloid families
#=======================================================================================

# Calculate prevalence among ALL sampled frogs

# Species-season labels (positive samples / total samples)
axis_labels <- c(
  "dry_Ameerega trivittata" = "<span style='color:#E6B800;'><i>Ameerega trivittata</i><br>n=12/13</span>",
  "wet_Ameerega trivittata" = "<span style='color:#0072B2;'><i>Ameerega trivittata</i><br>n=22/26</span>",
  "dry_Ameerega macero"     = "<span style='color:#E6B800;'><i>Ameerega macero</i><br>n=3/19</span>",
  "wet_Ameerega macero"     = "<span style='color:#0072B2;'><i>Ameerega macero</i><br>n=8/26</span>",
  "dry_Ameerega shihuemoy"  = "<span style='color:#E6B800;'><i>Ameerega shihuemoy</i><br>n=4/6</span>",
  "wet_Ameerega shihuemoy"  = "<span style='color:#0072B2;'><i>Ameerega shihuemoy</i><br>n=7/16</span>"
)

# Order of species-season combinations and alkaloid families
species_levels <- rev(c(
  "dry_Ameerega trivittata",
  "wet_Ameerega trivittata",
  "dry_Ameerega macero",
  "wet_Ameerega macero",
  "dry_Ameerega shihuemoy",
  "wet_Ameerega shihuemoy"
))

family_levels <- c(
  "HTX", "DHQ", "N-methyl-DHQ", "5,8-I", "5,6,8-I",
  "3,5-P", "Pyr", "Tricyclic", "Unclass", "Unknown"
)

# Calculate prevalence (% frogs containing each alkaloid family)
prevalence_summary <- counts %>%
  select(-c(3, 5:12, 40)) %>%
  pivot_longer(
    cols = where(is.numeric),
    names_to = "alkaloid",
    values_to = "abundance") %>%
  mutate(alkaloid = as.numeric(alkaloid)) %>%
  left_join(
    select(select_alk, alkaloid, Mabel_family),
    by = "alkaloid") %>%
  group_by(
    Sample,
    ATTRIBUTE_species,
    ATTRIBUTE_season,
    Mabel_family) %>%
  summarise(
    abundance_fam = sum(abundance),
    .groups = "drop") %>%
  filter(ATTRIBUTE_species != "Allobates femoralis") %>%
  group_by(
    Sample,
    ATTRIBUTE_species,
    ATTRIBUTE_season,
    Mabel_family) %>%
  summarise(
    present = any(abundance_fam > 0),
    .groups = "drop") %>%
  group_by(
    ATTRIBUTE_species,
    ATTRIBUTE_season,
    Mabel_family) %>%
  summarise(
    prevalence = mean(present) * 100,
    n_positive = sum(present),
    n_total = n(),
    .groups = "drop") %>%
  mutate(
    species_season = factor(
      paste(ATTRIBUTE_season, ATTRIBUTE_species, sep = "_"),
      levels = species_levels),
    Mabel_family = factor(
      Mabel_family,
      levels = family_levels))

# Heatmap
plot_prevalence <- ggplot(prevalence_summary, aes(x = Mabel_family,
                                                  y = species_season,
                                                  fill = prevalence)) +
  geom_tile(color = "white") +
  geom_text(aes(label = round(prevalence, 1)), size = 3) +
  scale_y_discrete(labels = axis_labels) +
  scale_fill_gradient(
    low = "white",
    high = "#0072B2",
    limits = c(0, 100),
    name = "% frogs\nwith alkaloid") +
  labs(
    x = "Alkaloid Family",
    y = NULL) +
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    axis.text.y = ggtext::element_markdown(),
    panel.grid = element_blank(),
    axis.title = element_text(size = 14),
    legend.title = element_text(size = 12),
    legend.text = element_text(size = 10))

Results

Seasonal and sex differences in space use reflect breeding strategies

We evaluated whether season influences space use area in species with contrasting breeding strategies (Fig. 2A). When considering all individuals (i.e., both males and females; Fig. 2A, upper panel), A. trivittata and A. macero showed slightly larger space use area in the wet season, whereas A. femoralis and A. shihuemoy showed slightly larger space use area in the dry season. However, statistically significant seasonal differences were only detected in A. trivittata. In the seasonal breeder A. femoralis, space use area showed sex-specific seasonal patterns (LM: season*sex: F1,50=12.99, p < 0.001). Tracking duration, included as a covariate in the model, also showed a significant positive effect on space use area (p < 0.05). Females showed larger space use area in the dry season (emmeans estimate = 1.28 ± 0.50, padj = 0.014), while males exhibited smaller space use area in the dry season (estimate = –0.84 ± 0.35, padj = 0.021). In the seasonal breeder A. trivittata, space use area differed between seasons regardless of sex (LM: season: F1,33=6.02, p < 0.020), with individuals exhibiting smaller space use area in the dry season (estimate = –0.75 ± 0.30, padj = 0.017). In the year-round breeder A. macero, space use area showed a season*sex interaction (LM: F1,49=4.71, p = 0.035); however, we detected no significant seasonal shift for either sex. The other year-round breeder, A. shihuemoy, showed no seasonal variation in space use area. Full model outputs and post hoc comparisons are provided in Supplementary Tables S3 and S4. Seasonal space use overlap and centroid shifts varied among species and sexes (Supplementary Results S1; Supplementary Tables S5 and S6). However, because the number of recaptures was limited, these results should be interpreted as exploratory only. The inclusion of tadpole transport movements in space use area estimates did not alter the overall conclusions (Supplementary Results S2; Supplementary Table S7).

Environmental variables also influenced daily movement differently among species: rainfall was associated with movement in A. femoralis, temperature influenced movement in A. trivittata and A. macero, and no effects were detected for A. shihuemoy (Supplementary Table S8 and Fig. S3).

Figure 2: Species-specific seasonal shifts in space use reflect contrasting breeding strategies in poison frogs. (A) Seasonal differences in space use area between sexes and across species. Points represent estimated marginal means, and bars represent the 95% confidence intervals from the selected linear models. Asterisks indicate significant differences ( * < 0.05). (B–E) Kernel density estimates illustrating the spatial distribution and overlap of seasonal space use in recaptured individual: (B) male Allobates femoralis, (C) male Ameerega trivittata, (D) male Ameerega macero, and (E) female Ameerega shihuemoy (no male recaptures were obtained for this species). Yellow indicates the dry season and blue indicates the wet season. Purple panels represent seasonal breeders and green panels represent year-round breeders. Sample sizes for panel A were: A. femoralis dry season (12 females, 14 males) and wet season (7 females, 23 males); A. trivittata dry (7 females, 9 males) and wet (8 females, 12 males); A. macero dry (8 females, 13 males) and wet (11 females, 22 males); and A. shihuemoy dry (6 females, 1 male) and wet (8 females, 9 males).

Species-specific seasonal microhabitat shifts

In the seasonal breeder A. femoralis, seasonal changes in microhabitat differed between sexes for specific substrate types (GLMM: season*sex*microhabitat: χ²(8)= 22.37, p = 0.004; Supplementary Table S9). Post hoc contrasts (Supplementary Table S10) revealed that females increased use of exposed structures (e.g., logs) during the wet season (dry/wet ratio = 0.56, p_adj = 0.033). Males used treefall refuges more frequently in the wet season (dry/wet ratio = 0.18, p_adj = 0.038) but increased their use of hidden microhabitats during the dry season, including under leaf litter (dry/wet ratio =1.37, p_adj = 0.036) and under logs (dry/wet ratio =1.79, p_adj = 0.041). A. trivittata, also showed seasonal shifts in microhabitat use independent of sex (GLMM: season*microhabitat: χ²(8)=31.92, p < 0.001; Fig. 3A, Table S9), with a greater use of hidden microhabitats during the dry season, including under leaf litter (dry/wet ratio =2.01, p_adj < 0.001), inside tree fall refuges (dry/wet ratio =7.40, p_adj < 0.001), and inside root cavities (dry/wet ratio =1.59, p_adj = 0.034). In contrast, during the wet season, they showed greater use of exposed microhabitats, such as on top of the leaf litter (dry/wet ratio = 0.33, p_adj < 0.001), and within leaf litter cavities (dry/wet ratio = 0.53, p_adj = 0.007), classified here as potential clutch sites.

Among year-round breeders, shifts were species-specific. Ameerega macero showed significant seasonal shifts in microhabitat use (GLMM: season*microhabitat: χ²(9)=19.10, p = 0.024; Fig. 3A, Table S9), with greater use of log and rock-associated microhabitat during the dry season, including on top of logs (dry/wet ratio = 1.74, p_adj = 0.046), under logs (dry/wet ratio =2.27, p_adj < 0.001), and inside rock refuges (dry/wet ratio = 4.23, p_adj = 0.035). In contrast, A. shihuemoy showed no seasonal variation in microhabitat use (GLMM: season*microhabitat: χ²(7)=6.35, p = 0.499).