Phenology and environmental determinants of explosive breeding in gliding treefrogs: diel timing of rainfall matters

Abstract

The influence of abiotic and biotic factors on the temporal pattern of calling and breeding for many temperate anurans is well understood. However, few studies have documented patterns of reproduction in explosive-breeding tropical frogs or incorporated multiple environmental factors in their analyses, especially across multiple breeding seasons. We combine long-term natural history observations and automated data collection methods with boosted regression tree (BRT) analysis to determine the phenology and determinants of explosive breeding in the gliding treefrog, Agalychnis spurrelli. We monitored breeding for a total of 418 days across three breeding seasons and determined the relative importance of several environmental factors on the probability of calling and breeding activity. Our study population of A. spurrelli on Costa Rica’s Osa Peninsula forms breeding aggregations up to 11 times per year during 1–2-day long explosive-breeding events, from late May to mid-September. Calling and breeding activity are strongly and positively related to accumulated rainfall during the previous 24 and 48–24 h before, particularly rainfall during the afternoon and evening. Day-of-year, days since breeding occurred, and lunar phase also influence reproductive activity. This study provides the first description and analysis of the phenology and factors that predict explosive breeding in A. spurrelli and illustrates the value of using automated data collection paired with BRTs for the analysis of complex ecological data.

La influencia de factores abióticos y bióticos en el patrón temporal de cantos y reproducción de muchos anuros de zonas templadas es bien conocida. Sin embargo, pocos estudios han documentado patrones de reproducción en ranas tropicales, que se reproducen de manera explosiva, o han incorporado múltiples factores ambientales en sus análisis, especialmente a lo largo de múltiples temporadas de reproducción. Combinamos observaciones de historia natural a largo plazo y métodos automatizados de recolectar datos con análisis de árboles de regresión y técnicas de remuestreo (‘boosted regression trees’, BRT) para determinar la fenología y los determinantes de la reproducción explosiva en la rana planeadora, Agalychnis spurrelli. Monitoreamos la reproducción durante un total de 418 días a través de tres temporadas de reproducción y determinamos la importancia relativa de varios factores ambientales en la probabilidad de cantos y reproducción. Nuestra población de estudio de A. spurrelli en la Península de Osa en Costa Rica forma agregaciones reproductivas hasta 11 veces por año, los eventos de reproducción explosiva tienen 1 a 2 días de duración, desde finales de mayo hasta mediados de septiembre. Los cantos y la reproducción están fuerte y positivamente relacionados con la acumulación de lluvia en las 24 y 48 a 24 horas anteriores, particularmente durante la tarde y la noche. El día del año, los días desde la última reproducción, y la fase lunar también influyen la reproducción. Este estudio es la primera descripción y análisis de la fenología y los factores que predicen la reproducción explosiva en A. spurrelli e ilustra el valor de utilizar la recolección de datos automatizada junto con BRT para el análisis de datos ecológicos complejos.

INTRODUCTION

The environment plays a key role in determining the temporal patterns of anuran reproductive activity, including calling and breeding. Studies have documented the influences of multiple abiotic factors, including temperature, rainfall, humidity, and light intensity (Brooke et al. 2000; Oseen and Wassersug 2002; Saenz et al. 2006; Hartel et al. 2007; Wells 2007; Llusia et al. 2013; Schalk and Saenz 2016). Recent work also suggests that photoperiod and lunar cycles are important for the timing of migrations to breeding sites and for calling activity (Both et al. 2008; Canavero and Arim 2009; Grant et al. 2009, 2012; Schalk and Saenz 2016; Underhill and Höbel 2018; Jarvis et al. 2021). However, how environmental factors influence breeding patterns vary greatly between and, sometimes, within species, particularly across latitudes and between species with different reproductive modes (Wells 1977, 2007; Oseen and Wassersug 2002; Gottsberger and Gruber 2004; Saenz et al. 2006; Canavero and Arim 2009; Schalk and Saenz 2016).

In temperate zones, temperature and rainfall are considered the prevailing determinants of calling and breeding activity (Hartel et al. 2007; Wells 2007; Llusia et al. 2013). These factors are particularly important for species that breed in temporary ponds that fill or thaw during warmer months with greater amounts of rainfall. Breeding periods are generally longer in the wet tropics than in temperate regions (Duellman and Trueb 1986; Wells 2007), but while some tropical species show strong responses to rainfall, photoperiod, and other abiotic factors, others do not (Inger and Bacon Jr 1968; Crump 1974; Brooke et al. 2000; Gottsberger and Gruber 2004; Medeiros et al. 2016; Schalk and Saenz 2016; Ulloa et al. 2019).

Anuran breeding patterns are generally categorized as either explosive or prolonged, ranging along a spectrum from a single night of breeding annually to year-round breeding (Wells 1977, 2007). “Prolonged” typically refers to breeding periods of several months, while “explosive” covers breeding periods from a single night in some species to several days or weeks of breeding in others (Wells 1977, 2007). These distinctions are most applicable to “breeding periods” of populations, rather than “breeding seasons” of species, since populations may vary in their breeding periods within a defined breeding season. For instance, northern populations of the edible frog (Pelophylax esculentus complex) have a single short breeding period each year (Forselius 1962), while those in warmer regions show multiple peaks of reproductive activity over a breeding period of 2–3 months (Wells 1977). The effects of environmental factors on the temporal pattern of breeding in temperate-zone explosive breeders have been well studied (reviewed in Miwa 2007). However, relatively few studies have documented general breeding patterns in explosive-breeding tropical anurans (Aichinger 1987; Duellman 1995; Bevier 1997; Prado et al. 2005), and even fewer have assessed the importance of multiple environmental factors on reproductive timing at finer scales and/or across multiple breeding seasons (but see Gottsberger and Gruber 2004; Schalk and Saenz 2016; Ulloa et al. 2019).

In this study, we combine natural history observations and automated acoustic monitoring with machine-learning statistical analyses to assess the phenology and environmental determinants of reproduction in a large, explosive-breeding population of gliding treefrogs, Agalychnis spurrelli, on Costa Rica’s Osa Peninsula. Agalychnis spurrelli is a largely nocturnal species that inhabits the canopy of humid lowland forests from southern Costa Rica to northwestern Ecuador (Duellman 1970; Savage 2002). Adults typically reproduce during explosive-breeding aggregations, where hundreds to thousands of individuals aggregate on vegetation overhanging forest ponds; males engage in intense physical competition while females collectively lay thousands to millions of eggs in communal masses (Scott and Starrett 1974; Savage 2002; Thompson et al. 2016; Güell et al. 2019). Typical breeding aggregations at our study site consist of ca. 2000 individuals (Güell et al. 2019). However, aggregations can be smaller or larger, ranging from a few hundred to multiple thousands of adults (BAG personal observation). Information on the reproductive ecology of this species is sparse and often anecdotal; breeding is thought to occur only a few times per year in response to heavy rains early in the rainy season (Scott and Starrett 1974; Savage 2002; Ortega-Andrade et al. 2011). However, long-term systematic monitoring and formal analyses of the timing and frequency of reproduction, or the factors that influence breeding, have not been conducted. Thus, the specific objectives of our study were to 1) describe the phenology of explosive breeding in A. spurrelli and 2) determine what factors trigger calling and explosive breeding activity.

MATERIALS AND METHODS

Study population and survey methods

We investigated the breeding pattern of A. spurrelli at Shampoo Pond (8°24ʹ55″N, 83°20ʹ45″W) on Costa Rica’s Osa Peninsula across three breeding seasons. Our surveys included 418 days (105, 113, and 200 consecutive days from 14 May to 26 August 2018, 25 May–14 September 2019, and 19 May–4 December 2021). Breeding aggregations of A. spurrelli are concentrated on one side of the pond and detectable both visually and acoustically from any location at the pond (Thompson et al. 2016; Güell et al. 2019). We conducted daily visual surveys at dawn (~06:00 h) by walking around the pond’s perimeter and its surrounding swamp areas and recording the presence or absence of explosive-breeding aggregations of A. spurrelli. To assess the breeding phenology and environmental factors that trigger breeding for most individuals within our study population, we focused our analysis on explosive-breeding aggregations—that is, when hundreds to thousands of frogs were present in a dense aggregation and females laid large communal egg masses (Güell et al. 2019). We occasionally observed a few dispersed pairs breeding non-communally, without an aggregation, particularly on days before and after the larger explosive-breeding aggregations. We estimated that the frogs participating in these non-communal breeding events represented <1% of all breeding activity, based on point count estimates of typical aggregations sizes (Güell et al. 2019). We, therefore, did not consider these cases representative of the breeding patterns of the population as a whole and did not include them in our analyses. In 2021, we deployed a single automated recording device at Shampoo Pond to monitor calling activity and infer breeding beyond the period of our visual surveys; we identified explosive-breeding aggregations using audio recordings only in 2021.

Acoustic recordings and analyses

Traditional in-person, survey-based methods for collecting reproductive phenology data are time-consuming, expensive, and often invasive (Digby et al. 2013). An alternative is to use automated recording devices to enhance data collection over space and time (Llusia et al. 2013; Schalk and Saenz 2016; Ulloa et al. 2019; Larsen et al. 2021). To complement our surveys and assess A. spurrelli calling and breeding for a longer period in 2021, we used an automated device (SM MINI, Wildlife Acoustics, Maynard, MA, USA) to record 5 min of 44.1 kHz audio at the start of each hour (120 min/day) from 1 June–16 December 2021. We attached the device to a large tree trunk, 3 m above the water surface and within 3 m from where A. spurrelli breeding aggregations are densest, with its microphone oriented toward the pond center (Thompson et al. 2016; Güell et al. 2019).

We listened to every recording for A. spurrelli vocalizations to assess calling and breeding activity. At our study site, A. spurrelli vocalizations are similar to those described by Scott and Starrett (1974) as a low-frequency, weak call (also see “aggressive call B” in Cossio and Medina-Barcenas 2020). Explosive-breeding choruses are audible in person and in our audio recordings as a continuous low-frequency murmur (BAG personal observation). Our analysis of the predictors of breeding activity includes only days with explosive-breeding choruses; our analysis of the predictors of calling activity includes all days when A. spurrelli called, in aggregation or individually. On days with explosive-breeding choruses, we considered the hourly recording with the highest chorus amplitude to represent the peak of breeding activity and estimated activity duration to span the first to last hour when A. spurrelli choruses were audible. Heavy rain between 20:00 and 00:00 h was common during our survey period in 2021 and sometimes masked frog vocalizations in our recordings. Because of this, our audio analysis sometimes underestimates breeding activity durations.

Environmental data

We obtained rainfall, air temperature, and relative humidity measurements from weather stations at Osa Conservation’s Piro Biological Station (Davis, Vantage Pro 2, 8°24ʹ14″N, 83°20ʹ11″W, 1.4 m altitude) for all dates in 2019 and 21 August–4 December 2021, and from the Instituto Meteorológico Nacional de Costa Rica’s weather station at the Greg Gund Conservation Center (https://www.imn.ac.cr/estaciones-automaticas, 8°24ʹ40″N, 83°19ʹ13″W, 235 m altitude) for all dates in 2018 and 19 May–20 August 2021. We obtained dates of full moons from the US Naval Observatory Astronomical Applications Department (https://www.willbell.com/almanacs/almanac_mica.htm). We used the number of days since the last full moon (DFM; 0–29) as a predictor in our models. We did not measure ambient light intensity or record whether the moon was visible.

Statistical/machine-learning analyses

We used boosted regression trees (BRTs; “dismo” package; Hijmans et al. 2017) to examine the major predictors of A. spurrelli calling and breeding (Friedman et al. 2000; Schapire 2003; Elith et al. 2006, 2008; Leathwick et al. 2006; De’Ath 2007; Buston and Elith 2011). BRTs are a powerful, assumption-free modeling technique that combines the benefits of standard regression and boosting approaches. They accommodate any type of variable (continuous, categorical, missing and non-independent data) and can include correlated sets of independent variables (Elith et al. 2008). In comparison to standard regression, they offer increased predictive power by basing predictions on multiple model fits via boosting, rather than on one simple model. BRTs determine the relative importance of predictor variables and allow for complex functions and their interactions (Elith et al. 2008). Here, we used BRTs to investigate the major predictors of calling and breeding activity in A. spurrelli. These analyses allowed us to identify and compare the relative influence of our abiotic and biotic factors on the probability of calling and breeding activity in A. spurrelli with high accuracy.

Two parameters must be chosen for BRTs: the “learning rate” (lr) and “tree complexity” (tc). The lr is used to adjust the contribution of each tree as it is added to the model; decreasing (or slowing) the lr increases the number of trees required. The tc is used to adjust the number of nodes in a tree, or the number of interactions within the model; fitting more complex trees decreases the number of trees required for minimum error. In general, fitting a model with a smaller lr and a higher tc is preferable and typically fits models with > 1000 trees (see Elith et al. 2008 for details and guidelines for selecting optimal parameter values). Here, we used lr’s between 0.001 and 0.0005 and tc’s of 5 to fit models that were slow enough to estimate responses reliably yet able to include complex interactions.

We compared different combinations of lr and tc and used a 10-fold cross-validation (CV) method for model development and evaluation, following recommendations by Elith et al. (2008). This approach progressively tests the model on withheld portions of the data and estimates the optimal number of trees needed to achieve minimum predictive error. We evaluated the success of models using different parameters of lr and tc by comparing the change in deviance explained and chose the model with at least 1000 trees that achieved the smallest predictive error, as suggested by Elith et al. (2008). We also assessed model performance by calculating the discrimination ability, measured as the area under the receiver operating characteristic curve (AUC). AUC values are commonly used in machine learning for model comparison; they combine sensitivity and specificity to indicate the discrimination ability of the model (Ferrier and Watson 1997). BRTs do not provide P values, but instead provide estimates for the relative importance of each predictor, based on how often it is selected and whether it improves the model when selected. The relative influence of predictors is scaled to sum to 100% and higher numbers indicate a stronger influence on the response variable. BRTs also automatically model interactions between predictors and provide a measure of relative strength for each interaction (see Elith et al. 2008 for an explanation of interaction strength quantification). To focus on the major factors that influence breeding and calling in A. spurrelli, we limited our interpretations of predictor variables, and their interactions, to those with a relative influence >5% (as in Buston and Elith 2011).

We visualized the relative influence of all our predictors using partial dependence plots, which show the effect of each predictor on the response variable while controlling for the average effect of all other variables in the model. We used three-dimensional partial dependence plots to visualize the most important (i.e., strongest) interaction between predictors in our models to predict explosive breeding and calling. Environmental conditions were considered predictor variables and A. spurrelli breeding or calling was the binary response variable. Since individuals begin arriving at breeding sites before midnight (as early at 20:00 h) and peak calling and breeding activity occurs primarily near dawn (i.e., 06:00 h) (Scott and Starrett 1974; Thompson et al. 2016; Güell et al. 2019), we used the total rainfall and the mean temperature and relative humidity during the previous 24 h as predictors in our analyses (i.e., from midnight to midnight). We also included rainfall accumulated during the periods 48–24 and 72–48 h before, day-of-year (DOY), days since the last full moon (DFM), and days since breeding occurred (DBO) as predictors in our models. Then, because rainfall during the previous 24 h was the strongest predictor of breeding (see “Results” section), we binned rainfall data into four 6-h periods to further investigate the relative importance of the exact timing of rainfall on the probability of breeding. Considering that choruses are loudest at, and most eggs are laid by, ca. 06:00 h, we analyzed the importance of five periods of rainfall preceding 06:00 h on the probability of breeding (i.e., 6–0 h before: midnight to dawn; 12–6 h before: dusk to midnight; 18–12 h before: noon to dusk; 24–18 h before: dawn to noon; and 30–24 h before: midnight to dawn of the previous night). Finally, to assess support for specific interpretations of our results, we discuss the performance of our full models (above) versus reduced models excluding rainfall during the previous 48 h, DOY, or DBO as predictor variables. All statistical tests were performed in the R statistical environment (R Core Development Team 2021) using RStudio (version 1.0.143).

RESULTS

Field surveys

We recorded 35 days of explosive breeding during our survey period (N = 14, 12, and 9 in 2018, 2019, and 2021, respectively). The earliest breeding occurred on 22 May in 2018 and 2021, and on 28 May in 2019; the latest breeding occurred on 12 August in 2018, 13 September in 2019, and 16 August in 2021. The mean period from first to last day of breeding, across seasons, was 93 days (83, 109, and 87 days in 2018, 2019, and 2021). The explosive-breeding choruses of A. spurrelli recorded by our automated device in 2021 (N = 6) began as early as 20:00 h, but on average began at 23:10 ± 1:56 h (mean ± SD in hours:minutes). Peak calling consistently occurred at ca. 06:00 h (5:49 ± 0:24 h) and calling lasted until 9:19 ± 0:49 h. The mean duration of choruses was 9.83 ± 2.04 h.

Most explosive breeding events lasted one day. However, on six occasions (N = 3, 1, and 2 in 2018, 2019, and 2021), explosive breeding took place on two consecutive days (Figure 1a), with hundreds to thousands of frogs chorusing and breeding communally. These 2-day-long events represented 34% of all explosive breeding during our survey period. We observed paired and unpaired females, some of whom were visibly gravid, sleeping near the breeding site after the first day of breeding activity, presumably so they could resume breeding the next day. The relative size of aggregations—based on simple visual estimates—and total duration of reproductive activity on the second day of breeding events were often less than on the first day (BAG personal observation).

Environmental variables

The total rainfall recorded across the entire survey period was 6754.7 mm and the highest daily rainfall was 300 mm on 14 July 2019 (Figure 1a). It rained more, and was warmer, during our 2019 survey period than either 2018 or 2021 (total rainfall accumulations: 1143.0, 3247.0, and 2364.7 mm; mean temperatures: 24.8 ± 0.6, 26.3 ± 0.7, and 24.9 ± 1.3 °C in 2018, 2019, and 2021, respectively; mean ± SD here and throughout the text). There was on average more rain during the afternoon and evening on days before breeding occurred than on other days (Figure 1b). Across our survey period, the mean rainfall during the 24 h before breeding was 50.0 ± 66.2 mm (N = 35), compared to 11.1 ± 33.5 mm (N = 383) before days without breeding.

Probability of explosive breeding

The probability of explosive breeding on any day during our survey period was 0.084. The final model to predict it explained 58% of the deviance using training data and 48% of the deviance using 10-fold CV. These estimates tell us how well the model explains the observed data and predicts withheld data, respectively. The AUC was 0.97 in training and 0.80 in CV; this metric indicates the discrimination ability of the model.

Seven variables were important predictors of explosive breeding: rainfall during the previous 24 h, rainfall 48–24 h before, DOY, DFM, humidity during the previous 24 h, rainfall 72–48 h before, and DBO (Table 1). The relative influence of temperature during the previous 24 h and of breeding season (year) was ≤5%, so we will not interpret them further.

The partial response plots for these seven variables (Figure 2) indicate that explosive breeding is most likely to occur after days with very high rainfall (over ~95 mm), higher mean humidity (more than ~97% RH), and early in the rainy season (before 24 June or DOY 175). The rainfall 48–24 and 72–48 h earlier also influenced the likelihood of explosive breeding (Figure 2b and f); more rain 48–24 h earlier and little-to-no rain 72–48 h earlier increased it. The rug plots for rainfall during the previous 24 h and the period 48–24 h earlier indicate that 90–91% of days had insufficient rain to elicit breeding (Figure 2a and b). Explosive breeding was more likely to occur 16–29 days after the full moon—that is, between waxing crescent and waxing gibbous lunar phases—and least likely during full and waning gibbous phases (Figures 2d and 3a). Finally, the probability of explosive breeding varied strongly with days since breeding occurred (DBO). It was more likely one day after breeding, less likely 2–10 days after, about equally likely 10–17 days after, and subsequently more likely (Figure 2g).

There were four interactions between predictors of explosive breeding, the most important being between DOY and rainfall during the previous 24 h (Table 2). The effect of DOY (Figure 2c) was only evident when rainfall was high, as breeding was always unlikely with little or no rain (Figure 4a). The effects of humidity and of rainfall 48–24 h before were similarly dependent on rain during the previous 24 h; when the latter was low (ca. <25 mm) breeding was unlikely and there was little or no effect of humidity or rainfall 48–24 h before. The relative influence of breeding season was <5% so we do not interpret its interaction with DFM.

Probability of explosive breeding: diel timing of rainfall

Rainfall amounts during all five 6-h periods during the 30 h before peak breeding time (i.e., 06:00 h) were important predictors of breeding (Figure 5). The model explained 59% of the deviance using training data and 54% of the deviance using CV. The AUC was 0.88 in training and 0.73 in CV. The partial response plots for these variables (Figure 5) indicate that more rainfall during all these periods increased the likelihood of breeding. The strongest predictors of breeding were rain 18–12 and 12–6 h before peak breeding time—that is, noon to dusk and dusk to midnight (Figure 5f). The rug plots for rainfall during these time periods show that 88–90% of days had insufficient rainfall to elicit breeding (Figure 5a and b).

Probability of A. spurrelli calling

The probability of A. spurrelli calling on any day during our automated acoustic recording period in 2021 was 0.19. We recorded their vocalizations on 36 days, including 29 without an explosive-breeding chorus. Rainfall during the 24 h before days with calling was 34 ± 35.2 mm versus 7 ± 13.3 mm before days without calling (N = 152).

Six variables were important predictors of A. spurrelli calling: rainfall during the previous 24 h, mean humidity, DFM, DBO, DOY, and rainfall 72–48 h earlier (Table 1). The relative influence of temperature during the previous 24 h and rainfall 48–24 h earlier were <5% so we do not interpret them further.

The final model explained 98% of the deviance using training data and 80% of the deviance using CV, and the AUC was 0.98 in training and 0.81 in CV. The partial response plots for the six variables (Figure 6) indicate that A. spurrelli are more likely to call after days with more rainfall (over ~50 mm) and higher humidity (over ~97%) (Figure 6a and b). The rug plot for rainfall during the previous 24 h shows that 92% of days in 2021 had insufficient rainfall to elicit calling (Figure 6a). Males were more likely to call 13–20 days after the full moon—that is, between new and waxing crescent lunar phases (Figure 3b). The probability of calling also varied strongly with DBO, increasing after 11 days without breeding activity (Figure 6D). Lastly, the probability of calling was higher between ~19 July and 7 October (i.e., DOY ~200–280) (Figure 6e).

There were three interactions between predictors of calling, with the most important being the interaction between DFM and humidity during the previous 24 h (Table 2). When humidity was below ca. 96%, the probability of calling was low and there was a little effect of DFM (Figure 4b). The effects of humidity and DBO on calling were dependent on rainfall during the previous 24 h; when rainfall was low (ca. <25 mm), the probability of calling was low, with little or no effect of humidity and DBO.

DISCUSSION

We determined the effects of several environmental factors on the probability of reproductive activity in gliding treefrogs, A. spurrelli, on Costa Rica’s Osa Peninsula. Our field surveys and machine-learning analyses indicate that A. spurrelli reproduce during multiple short explosive-breeding events distributed across several months at the beginning of the rainy season. This work improves our understanding of the reproductive phenology of explosive-breeding Neotropical frogs, and it highlights the value of combining natural history observations with advanced statistical methods to broaden our understanding of animal behavior and reproductive ecology.

The primary determinants of reproductive activity in A. spurrelli are rainfall during 1) the previous 24 h, particularly the afternoon and evening and 2) rainfall 48–24 h earlier. The positive relationship between rainfall and reproductive activity (i.e., calling and breeding) is consistent with associations found in other tropical explosive-breeding species where rainfall elicited reproductive activity, with 1–2-day lags (Aichinger 1987; Roberts 1994; Bevier 1997; Gottsberger and Gruber 2004; Greenberg and Tanner 2004; Saenz et al. 2006; Schalk and Saenz 2016; Ulloa et al. 2019). A lagged response to rainfall may reflect the time it takes frogs to migrate to breeding sites. However, our observation of individuals arriving to the breeding site as early as 20:00 h, just hours after heavy afternoon and early evening rain suggests that A. spurrelli can migrate to breeding sites relatively rapidly. Their parachuting ability may help them to arrive within hours after heavy rainfall (Scott and Starrett 1974). Rapid migration after afternoon and evening rainfall (Figures 1b and 5) might also explain why A. spurrelli begin to aggregate at night and continue breeding into the early morning, rather than the following night. Breeding 1 to 2 days after heavy rainfall might help to ensure the accumulation of enough water to create a persistent habitat for hatchlings. Nonetheless, we observed A. spurrelli breed over dry land early in the rainy season before the pond was full, in both 2019 (González et al. 2021) and 2023 (BAG personal observation).

It is possible that A. spurrelli breeding only occurs after the accumulation of enough rain to surpass some response threshold. Our analyses revealed that less than 10% of days across our three survey periods had enough rainfall for calling and/or breeding to occur (Figures 2 and 6). Breeding in response to heavy rainfall 48–24 h earlier may reflect a 2-day accumulation of enough rain to elicit breeding and/or extended durations of some explosive-breeding events. In total, we observed six 2-day breeding events (Figure 1a) and there are previous reports of 2–4-day explosive-breeding aggregations of A. spurrelli at our site and a nearby location on the Osa Peninsula (Scott and Starrett 1974; Thompson et al. 2016). Two-day breeding events also explain the high probability of breeding the day after breeding occurred and considerable decrease in probability two days later (Figure 2g). Breeding was also slightly more likely to occur when there was little-to-no rainfall 72–48 h earlier (Figure 2f), suggesting that breeding may be particularly triggered by heavy rainfall after dry spells. Whether explosive-breeding events last more than 1 day likely depends on several factors, including the amount and consistency of rainfall, or lack thereof, over the course of several days.

The strong negative relationship between DOY and breeding provides clear quantitative support that A. spurrelli primarily breed early in the rainy season (Figure 2c), but DOY is only important when rainfall during the previous 24 h is particularly high (Figure 4a and Table 2). Dropping DOY from our model had little effect on its performance or the distribution of relative importance among the other predictor variables (Table 3). The altered model explained 57% of the deviance using training data (vs. 58%), 47.7% of the deviance using 10-fold CV (vs. 48%) and had an AUC of 0.96 (vs. 0.97). In contrast, dropping the strongest predictors of breeding—accumulated rainfall during the previous 24 and 48–24 h—but including DOY, had a large impact on model performance (55% vs. 48% deviance explained using CV) and the distribution of relative importance among the other predictor variables (Table 3). These reduced models further highlight the critical importance of rainfall for A. spurrelli breeding and the secondary (interaction) role of DOY.

We presumably captured the end of the breeding season during our surveys, particularly in 2021 when our acoustic monitoring extended into December. However, breeding may have occurred before we started our surveys. For instance, in 2022 ENSO (La Niña) conditions brought about a very early start to the rainy season and A. spurrelli bred on 1 May (BAG personal observation). Their apparent high dependence on heavy rainfall and annual rhythm for reproduction suggest A. spurrelli may be particularly susceptible to climate change. Indeed, the increase in rainfall variability and dry spells in the tropics (Hulme and Viner 1998; Christensen et al. 2007; Trenberth 2011) is already negatively affecting egg survival in A. spurrelli (González et al. 2021) and may also affect the timing of future breeding events.

Temporal reproductive synchronization and predator avoidance are the dominant aspects of amphibian biology hypothesized to be affected by the moon (Grant et al. 2012). However, despite the recent finding that there are more cases of amphibians being affected by lunar cycles than unaffected, it is unclear whether these effects are due to lunar periodicity or moonlight intensities (reviewed in Grant et al. 2012). This distinction is beyond the scope of our study as we did not measure light levels. We did find that 54% of explosive breeding occurred between waxing crescent and waxing gibbous lunar phases (20–24 DFM) and only ~8% occurred during the full moon (Figure 3a). These findings are congruent with previous reports of breeding aggregations of A. spurrelli at our site in 2015, which also occurred between 20 and 24 DFM (Thompson et al. 2016), and a breeding event elsewhere on the Osa Peninsula on 10–11 August 1970, during the first quarter lunar phase (21–23 DFM; Scott and Starrett 1974). There is no indication that A. spurrelli explosive-breeding events are more likely during the full moon, in contrast to at least six species—including several explosive breeders—where migration to breeding sites and/or reproductive events are concentrated around the full moon (reviewed in Grant et al. 2012). Interestingly, A. spurrelli calling was more common during the darker new and waxing crescent lunar phases (13–20 DFM) and DFM had higher importance in predicting calling than breeding (12% vs. 8%, Figures 2d, 4 and 6c). This may indicate predator avoidance, since calling on darker nights may help male A. spurrelli avoid visual predators (Taylor et al. 2007; Pena et al. 2008); however, very few studies have empirically tested this hypothesis (see Grant et al. 2012).

Photoperiod has been identified as an important predictor of anuran activity, in some cases more important than rainfall or temperature (Both et al. 2008; Canavero and Arim 2009; Medeiros et al. 2016; Schalk and Saenz 2016). It is a common and unambiguous cue that provides information about time of year and can help predict environmental conditions in the future or at distant locations (Bradshaw and Holzapfel 2007). We cannot entirely reject that photoperiod affects reproductive timing in A. spurrelli. However, photoperiod is generally more important in temperate zones (above ca. 15°) where changes in day length are more pronounced and provide a more reliable seasonal cue (Bradshaw and Holzapfel 2007). For instance, Ulloa et al. (2019) monitored explosive breeding of several anuran species in French Guiana—at ca. 4° latitude—and found no association with photoperiod. Conversely, prior studies conducted at higher latitudes—in the Gran Chaco ecoregion of Bolivia at ca. 18°, southern Uruguay at ca. 34°, and southern Brazil at ca. 29°—found photoperiod to be an important factor in calling and breeding activity (Both et al. 2008; Canavero and Arim 2009; Medeiros et al. 2016; Schalk and Saenz 2016). The results of our reduced BRT model are consistent with this interpretation; excluding DOY from our model did not reduce its performance or predictive ability (Table 3), suggesting that photoperiod does not strongly influence A. spurrelli breeding activity at our study site in Costa Rica at ca. 8° latitude.

In addition to external environmental cues, reproductive timing in most frogs is influenced by biotic factors, such as endogenous reproductive rhythms. Most species show a predictable annual cycle of gametogenesis that closely matches their breeding season, and the ability to breed multiple times per season depends on physiological constraints as well as environmental conditions (Rastogi et al. 2011). In A. spurrelli, females likely impose stronger endogenous constraints on breeding phenology, since most males fail to mate during highly male-skewed aggregations where male–male competition is intense (Scott and Starrett 1974; Güell et al. 2019), and some call on non-breeding days (this study). Females presumably require more time to produce their next clutch of eggs than males do to replenish their potentially depleted sperm stores (Owens and Thompson 1994). In captivity, female A. spurrelli can breed 2–3 times per year at roughly 4-month intervals (Bland, personal communication; also see Bland 2013). These interclutch intervals may reflect the time needed for females to accumulate enough energy for another clutch, and they contrast with the population-level relationship between DBO and reproductive activity we found, with calling and breeding probability reduced for about 11 and 17 days after explosive breeding, respectively. It thus seems unlikely that individual females breed in consecutive events. Rather, interbreeding intervals may reflect the mean period of time between optimal environmental conditions for breeding—days with ≥100 mm of rain in 2019 (N = 11) occurred at 11 ± 9-day intervals (median = 10 days; mode = 11 days; range = 1–33 days).

Excluding the biotic factor DBO from our model did not change its performance or the distribution of relative importance of other variables (Table 3). The reduced model explained 57.6% of the deviance using training data (vs. 58%), 47.4% of the deviance using 10-fold CV (vs. 48%) and had an AUC of 0.97 (vs. 0.97). The fact that our models performed similarly with or without our only biotic predictor variable of DBO (Table 3) suggests that reproductive timing in our study population is mostly influenced by abiotic environmental factors. It also suggests A. spurrelli reproduction is somewhat asynchronous since intervals between breeding aggregations are shorter than likely interclutch intervals of females. The breeding phenology we found may reflect differences in when individual females are ready to reproduce and/or variation in female responsiveness to rain. If most females are ready and responsive at the beginning of the rainy season, it would explain why the first breeding events each year are the largest (BAG personal observation). Subsequent, smaller aggregations may include only females that were previously unready or unresponsive. Thus, even optimal environmental conditions shortly after an explosive breeding event might not elicit reproductive activity if too few females are ready to breed. Whether and how A. spurrelli interclutch intervals interact with rainfall patterns is unknown, but comparisons between gametogenic cycles, environmental conditions, and breeding cycles could elucidate endogenous constraints on their reproductive phenology.

This study provides the first description of the reproductive phenology of A. spurrelli and is the most comprehensive empirical analysis to date of the factors that determine explosive breeding in any Neotropical anuran across multiple breeding seasons. These data provide important baselines both for future comparisons and to inform conservation efforts, considering the current climate crisis. Across taxa, climate change is affecting amphibians the strongest, with many species showing a phenologically earlier onset of breeding activity (Blaustein et al. 2001; Parmesan 2007). Clearly, large-scale shifts in weather patterns (e.g., drought) can have negative effects on amphibian reproduction (Jansen et al. 2009); however, our results also highlight the importance of fine-scale changes in environmental conditions (e.g., daily and hourly rainfall) on reproductive activity. Longer-term data sets like ours are crucial to understanding how we can manage vulnerable populations, particularly during infrequent yet critical life events such as breeding. Finally, our work highlights the utility of long-term natural history observations and automated data collection strategies in conjunction with machine-learning techniques to provide new insights about the ecology of elusive tropical species.

FUNDING

This work was supported by the National Science Foundation (DGE-1247312 to BAG and IOS-1354072 to KMW), Sigma Xi (G2018031596022314 to BAG), and Boston University.

We thank our Egg Science Research Group, Sarah Davies, and Michael Caldwell for helpful comments on the manuscript, and Peter Buston for discussions about the research, implementation of BRTs and feedback that improved the manuscript at multiple stages. We thank Paulina Rodriguez and Priscilla Monserrat Solano for assistance in obtaining weather station data; Manuel Sánchez Mendoza, Elena K. Gomez, and Katherine González for field assistance; and Osa Conservation for laboratory space and logistical support. This research was conducted under Boston University IACUC protocol PROTO201800595 and permits from the Costa Rican Ministerio de Ambiente y Energía (MINAE) and Sistema Nacional de Áreas de Conservación (SINAC) (ACOSA-INV-048-18, ACOSA-INV-033-19, ACOSA-INV-058-19, ACOSA-DASP-OI-R-019-2021).

AUTHOR CONTRIBUTIONS

Brandon Güell (Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Validation, Visualization, Writing – original draft, Writing – review & editing), and Karen Warkentin (Conceptualization, Funding acquisition, Methodology, Project administration, Resources, Supervision, Validation, Writing – review & editing)

CONFLICT OF INTEREST

The authors have no conflict of interest to declare.

DATA AVAILABILITY

Analyses reported in this article can be reproduced using the data provided by Güell and Warkentin (2023).

REFERENCES