Citizen science reveals meteorological determinants of frog calling at a continental scale

Here we investigate the strength of the relationships between meteorological factors and calling behaviour of 100 Australian frog species using continent‐wide citizen science data. First, we use this dataset to quantify the meteorological factors that best predict frog calling. Second, we investigate the strength of interactions among predictor variables. Third, we assess whether frog species cluster into distinct groups based on shared drivers of calling.


| INTRODUC TI ON
The predictability of phenological events can be vital to our understanding of species distributions and population trends. For example, basic information, such as the timing of breeding seasons, is often used to guide biodiversity monitoring, with many species only detected reliably during the breeding season (e.g., Roth et al., 2014;Wilson & Bart, 1985). Accurate information on drivers of phenology is therefore critical to help inform conservation measures and priorities, such as well-timed, cost-effective population monitoring (Canavero et al., 2019;Visser et al., 2010).
Understanding the timing of phenological events with environmental conditions is necessary to understand the impact of climate changes to the species' phenology (e.g., Kellermann & van Riper, 2015). However, the cues driving circannual rhythms are difficult to assess as separate factors, as they often covary with each other over time and space. Seasonal trends in temperature are thought to be the main determinants of phenology for many species (Ficetola & Maiorano, 2016;Mazaris et al., 2013;Simonneau et al., 2016;Visser et al., 2006), though day of year (season-specific day length) may be an overlooked, noise-free factor that best explains temporal patterns. Day of year is intrinsically linked with other factors known to influence a species' phenology (e.g., temperature or rainfall) (Adole et al., 2019;Canavero & Arim, 2009). Understanding how meteorological conditions influence breeding is a critical issue, considering how essential breeding is to the survival of a species, the rapid changes in environmental conditions and stability throughout most of the world (Hällfors et al., 2020;Wadgymar et al., 2018), and -crucially -that historically coupled weather features, such as photoperiod and rainfall, may become uncoupled under future climate scenarios (Peñuelas et al., 2004;Thackeray et al., 2016).
One of the main challenges in understanding drivers of breeding phenology has been the lack of empirical data across the range of a species (Duursma et al., 2017;Elmore et al., 2016;Hurlbert & Liang, 2012). Robust population monitoring is often dependent on detecting breeding behaviour, and understanding optimal survey windows is a critical first step. This is particularly true for frogs, as many species are cryptic when not breeding. Their breeding phenology (Gomez-Mestre et al., 2012;Oseen & Wassersug, 2002;Ospina et al., 2013;Yoo & Jang, 2012) and consequently detection probability (MacKenzie, 2006;MacKenzie et al., 2003;Mackenzie & Royle, 2005) are strongly tied to meteorological cues. These limitations affect our ability to synthesize the results of multiple studies and in turn our ability to monitor meaningful trends in phenology and across populations. With so many logistical challenges to longterm and broad-scale monitoring, there are few monitoring programs of large magnitude.
Specific associations between frog calling and a range of abiotic variables including rainfall, temperature, humidity, vapour pressure, and moonlight have been investigated (see Table 1 for a global summary of previous research). These studies show that meteorological influences, either alone or in combination, frequently stimulate calling behaviour. However, the strength of the associations between calling and these factors varies among species and across studies (Heard et al., 2015;Oseen & Wassersug, 2002;Saenz et al., 2006). Broadly, rainfall has been found to be the strongest predictor of calling in the tropics, where photoperiod appears largely uninformative (Bradshaw & Holzapfel, 2007). In temperate zones, a combination of rainfall and temperature has been found to be most influential (Duellman & Trueb, 1994). In contrast, photoperiod was the most significant factor associated with frog calling in a long-term study in temperate Uruguay (Canavero & Arim, 2009).
The relationship of frog calling to meteorological factors varies among species and is not well documented across taxa, despite frog species being often categorized in the literature based on cues assumed to predict their breeding. The term "explosive breeder", for example, is used for species reliant on ephemeral ponds and impacted by rainfall, but potentially less influenced by abiotic factors once calling has commenced. "Seasonal breeders" are thought to be reliant on permanent ponds and influenced by temperature, while "generalists/prolonged breeders" call throughout the year or the rainy season (Heard et al., 2015;Lemckert & Grigg, 2010;Oseen & Wassersug, 2002;Saenz et al., 2006).
Increasingly, data at the broad spatial and temporal scales that could address these limitations are being collected by citizen science (Bird et al., 2014;Hochachka et al., 2012), though remarkably few citizen science programs (5%) have focused on frogs (Lloyd et al., 2020).
Frogs are one of, if not the most, imperilled vertebrate groups due to compounding threats including disease, habitat loss, competition from invasive species, and climate change (Gillespie et al., 2020).
Critically, evidence shows that frogs are already responding to the direct and indirect impacts of climate change (Blaustein et al., 2010;Cohen et al., 2018;Parmesan, 2006). Here, we use continental-scale citizen science data that document the calling of 100 Australian frogs (42% of species known in Australia) to investigate the strength of the relationships between meteorological factors and calling behaviour, an imperfect but common proxy for breeding (Crouch & Paton, 2002;Dorcas et al., 2009;Pellet & Schmidt, 2005). While many interacting external and internal cues inform phenology (i.e., hormones stimulated by perceived photoperiod), as a large, macroscale analysis, we consider only external meteorological cues. First, we use this dataset to quantify the meteorological factors that best predict frog calling. Second, we investigate the strength of interactions among predictor variables. Third, we assess whether frog species cluster into distinct groups based on shared drivers of calling (i.e., explosive, seasonal, or generalist breeders).

| Overview
We used frog occurrence data across 3 years to quantify the relationship between frog calling -a proxy for breeding activity -and meteorological variables. Our objective was not to predict the "timing" Note: For each study listed, the temporal and spatial scale of sampling is noted, the number of frog species studied, and the predictor variables explored in relation to calling behaviour. A checkmark indicates the study considered the variable indicated by that column. Lagged rainfall indicates the sum of rainfall over the preceding days (ranging from 1 to 7 days). If a latent temporal variable was included, details are provided. of peak breeding activity, but rather to assess the influence of meteorological variables on frog calling behaviour. Because we performed a cross-species analysis, our objective was focused on looking at the strength of relationship for given predictors, as opposed to the direction of the relationship (i.e., positive or negative). To accomplish this goal, we performed the following three overall steps ( Figure 1): (1) aggregated frog occurrence data from a popular citizen science dataset in Australia and spatially filtered these data to quantify whether a frog was calling in a particular grid cell on a given day; (2) integrated these data with meteorological variables (see Table 1); and (3) used boosted regression tree models to quantify which meteorological variables best predicted the likelihood of a species calling. We treat each of these steps in detail in the following sections.

| Frog occurrence data
We used frog occurrence data from FrogID, a national citizen science project led by the Australian Museum (Rowley et al., 2019;. Since its inception in November 2017, FrogID has collected over 600,000 validated observation records from 211 species -84% of frog species known in Australia. These data cover all of Australia (see Figure S1), with a bias towards the more highly populated east coast, particularly the state of New South Wales (10.4% land area of Australia yet 47% of all FrogID records). In addition to some spatial bias, there are temporal biases in submissions, with highest number of frog records in spring/summer, but this peak corresponds with peak breeding times for the majority of frog species (see Liu et al., 2021). A large part of its success is because auditory (call) surveys are one of the most common survey methods used to detect breeding frogs (Crouch & Paton, 2002;Da Silva, 2010;Lepage et al., 1997;Pellet & Schmidt, 2005). Participants submit 20-60-second audio recordings of calling frogs using a smartphone app, and the app adds associated metadata (time, date, latitude, longitude, and an estimate of precision of geographic location) to each submission. After a recording is submitted, a team of experts at the Australian Museum independently identifies any frog species heard calling in the recordings. Recordings with identifiable frog calls F I G U R E 1 (a) Methods illustrated using a single species as an example. Frog call data from FrogID was fitted to 10km 2 grid cells across the known range of each species and (ii) paired with meteorological variables from the same day and grid cell. This allowed us to (iv) pair call data with predictor variables within each species range, as well as infer pseudo-absences from other species detected in each spatiotemporal subsample. (v) We used boosted regression trees to assess the relationship between calling and our predictor variables. (vi) The scaled predictor importance from each analysis output allowed us to (b) compare both the predictor variables' importance to each other on an individual species level, as well as the strength of multiple species relationships to the same set of environmental drivers We define a "submission" as a submitted recording and an "observation" as a single record of a frog species originating from a submission for a particular site/date/time combination.
We used FrogID data from 11 November 2017 through 30 November 2020 (~36 months) and included observations of all species with a minimum of 100 validated observations, to reasonably represent each species calling behaviour and the environmental conditions it experiences. Observations were aggregated to 10 km 2 grid cells across Australia (Chase et al., 2019;Field et al., 2009).
First, for every species, we extracted all records within that species' geographic range (Rowley et al., 2019). And within each grid cell in a species' range, for each day, a species was either recorded as present or "absent." In this way we account for both the spatial and temporal biases commonly present in citizen science datasets (Bird et al., 2014). We inferred "absence," or pseudo-absence, if a species was not recorded in a grid cell on a particular day, but other species were detected (Rowley et al., 2019) (see Appendix S1 for raw data jitter plots displaying in-range presence and absence patterns by variable). We also calculated the total number of submissions for each grid cell to be used as a proxy for sampling effort (see modelling below). After applying these inclusion criteria, we used data from, on average, 22% of the grid cells within a species' known range (SD 0.15) (see Table S1 for the total number of grid cells, observations, and cv AUC for each species model).

| Meteorological variables
We collated data for the following environmental and meteorological variables: maximum temperature, mean of 10-day maximum temperature, minimum temperature, mean of 10-day minimum temperature, humidity, rainfall, cumulative rainfall from the previous 3 days, cumulative rainfall from the previous 10 days, and moon phase. We chose these variables based on a literature review of previously investigated variables, and their significance (see Table 1).
We downloaded Bureau of Meteorology (BOM) data for every day from 11 November 2017 to 30 November 2020 (BOM, 2020) to align with each cell of our aggregated frog occurrence data. Weather variables were aggregated and averaged from BOM 5 km 2 to 10 km 2 grid cells to create a database that paired call observations and weather data to uniform grid cells across the known range of a species for each grid cell and day of the 3-year dataset (see Figure 1 for a visual aid). Moon phase was calculated using the R package Suncalc (Thieurmel & Elmarhraoui, 2019). Aggregate rainfall over the previous 3 and 10 days was calculated as the cumulative total rainfall from BOM daily rainfall data. Although wind has also been investigated as a factor in call probability (Oseen & Wassersug, 2002;Penman et al., 2006), it was excluded because the available data are not accurate enough to associate with daily call activity (Jakob, 2010) as a grid cell scale. We assigned a numerical day of year to each day as an explanatory factor representing day length and harmonic regression (Chatfield & Xing, 2019;Weir et al., 2005) and distinguishing spring and autumn days from one another, unlike photoperiod (see Figure S2 for a plot displaying the correlation between day of year and photoperiod).

| Data analysis
Our objective was to quantify whether meteorological variables (described above) predict frog calling at a continental scale across Australia. We used boosted regression trees to assess the relationship between our binomial response variable of calling (i.e., presence/absence) and the meteorological traits used as predictor variables. Boosted regression trees are an additive regression model used for explanation and prediction (Elith et al., 2008). They are a combination of decision tree algorithms and boosting methods. Like Random Forest models, they fit many decision trees to improve the accuracy of the model. It is a modelling method growing in popularity for ecological and citizen science data Fink et al., 2020;Hochachka et al., 2012). The model structure suits our data because boosted regression trees allow modelling of nonlinear relationships that vary based on the nature of the relationship among different groups (i.e., observed species) and variables (i.e., meteorological variables), and test whether interactions have been detected and report the relative strength of these among predictor variables (Elith et al., 2008).
Predictor variables included in our full annual cycle models (i.e., data included for the entire year, all 3 years) for each species were (1) daily maximum temperature, (2) daily minimum temperature, (3) daily humidity, (4) daily rainfall, (5) cumulative rainfall from the previous 3 days, (6) cumulative rainfall from the previous 10 days, (7)  and longitude with the predictor variables of interest are included as a part of the modelling process. Moreover, we were not making spatial predictions, only accounting for potential differences (i.e., interactions among variables) in space. We also included the number of observations across all species per grid cell and day in the model as covariates to account for the biases that may influence species presence and detection rate (Brodie et al., 2020;Johnston et al., 2021). Boosted regression trees were fit using the R package dismo (Hijmans et al., 2021) gbm.step function, which uses crossvalidation to estimate the optimal number of tree for each model.
We used a tree complexity of 5, a learning rate of 0.005, and a bag fraction of 0.5 based on exploratory analysis of our data and suggestions by Elith et al. (2008). Tree complexity determines the degree to which predictors may interact with each other in relation to the response variable. The learning rate determines the contribution of each tree to the model. The bag fraction is the portion of the data drawn at random and without replacement from the full training set with each iteration. We fit the model with a "Bernoulli" distribution, as our response variable is a binary presence/absence (see methods illustrated for an example species in Figure 1).
For each species with over 100 frog observations within a target species' range (N = 100 species, observations = 152,534; Table S1), we extracted the relative influence of each predictor variable. We scaled these values from 0 (the least important variable) to 1 (the most important variable), excluding our covariates (latitude, longitude, and number of observations per grid and day) included in the models. The scaled independence of this variable allowed us to compare both the predictor variables' importance to each other on an individual species level, as well as the strength of multiple species relationships to the same set of meteorological drivers, while remaining population-and scale-independent across species. We also extracted predictor interactions from each model and scaled them from 0 to 1 to identify and compare relevant interactions among predictor variables. To test the robustness of our modelling process outlined above to spatial and temporal biases, for 10 of the most widely distributed species we re-ran our analysis (Figure 1). stratified by Australia's largely contiguous climate zones (i.e., temperate, subtropical, tropical, and desert) (BOM, 2021), and then stratified by year (i.e., 2018, 2019, and 2020).
To explore environmental predictors within a breeding season only, we repeated the above analysis for all species where day of year was the most important predictor (predictor importance = 1) (N = 67 species, observations = 80,774; Table S2). For each species, we identified the breeding season by creating a histogram of the observations grouped by day and clipping the annual data to consecutive days hosting 90% of the observations. For these models, day was not included as a variable in order to test the influence of other variables, aside from day.
A K-means cluster analysis (Aristeidou et al., 2017;Jain, 2010) was used to group frogs by their call patterns. K-means is an unsupervised learning algorithm used to identify patterns in data, and form groups based on those patterns. The iterative algorithm tests the Euclidian distance of each species to every group centroid. After a new species is classified, a new centroid is calculated as the mean of all species clustered in each group. The classification converges and the iterations stop when fewer species change their cluster assignment than in the previous iterations. We used the 10 predictor importance variables from the boosted regression tree analysis as inputs to a K-means cluster algorithm using the R packages cluster (Maechler et al., 2019) and factoextra (Kassambara & Mundt, 2020).
The output compiled categories of frogs based on the degree each variable explained calling behaviour in each species using the scaled predictor importance from the full annual cycle boosted regression tree models. We plotted the within cluster sum of squares for 1 to 20 groups and chose an optimal cluster size of seven by recognizing the asymptote change and performing a cluster validation using silhouette width.

| RE SULTS
We used a total of 152,534 frog occurrence records for 100 species, with a mean ± s.d. of 1510 ± 3664 records per species. The predictive power of the models was relatively high (mean cv AUC = 0.888 ± 0.01) (see Table S1 for species-specific scores). And unsurprisingly, the covariates included in the models were often of the highest relative influence (e.g., latitude, longitude, number of observations per grid and day).
Day of year was the strongest predictor of frog calling overall: in 67 out of 100 species, day of year had a predictor importance of 1, and the overall mean predictor importance among all species was 0.878 ( Figure 2a). Mean of the maximum temperature over a 10-day period was the second most important variable across all species with a mean influence of 0.532, followed by mean of the minimum temperature over a 10-day period temperature with a mean influence of 0.473. The relationship of calling to the four temperature variables was highly heterogeneous across species.
For some species (e.g., Crinia subsignifera, Cyclorana verrucosa, and Mixophyes fasciolatus), maximum and minimum temperature over a 10-day period were the most important variables, while mean minimum temperature on the day of calling had a mean influence of 0.299, and mean maximum temperature on day of calling had a mean influence of 0.298. Humidity had an overall importance of 0.322 and had a predictor importance of 1 in only one species (Litoria serrata). All rainfall variables showed a low relationship to calling across the majority of species, generally trending down from rainfall accumulated in the last 10 days (mean = 0.334), rainfall accumulated in the last 3 days (mean = 0.123), and rainfall the day of calling (mean = 0.008). However, some species (e.g., Austrochaperina pluvialis, Pseudophryne raveni, and Uperoleia tyleri) showed strong responses to previous rainfall. Moon phase was more significant than two of the rainfall variables (mean importance = 0.253), but not of high predictor importance for any species (see Figure S3 for species specific results).
We found the strongest model interactions among the most significant variables. Day and mean of maximum temperature over a 10-day period (mean = 1), followed by interactions between day and mean of minimum temperature over a 10-day period (mean = 0.61), and interactions between mean of maximum and minimum temperature over a 10-day period (mean = 0.61) (Figure 3).
When we repeated our analysis for those species with a strong dependence on day of year ("seasonal breeders"), restricting data to within the breeding season, we found notable differences between overall predictor importance compared to the annual species model (Figure 2b). For the 67 seasonal breeders assessed, the most important variables in predicting calling within their breeding season were minimum temperature over a 10-day period (mean importance = 0.772) and maximum temperature over a 10-day period (mean importance = 0.716). The next most significant variable was rainfall accumulated over the past 10 days with a mean predictor importance of 0.651. Fourth, humidity had an overall variable importance of 0.557, followed by mean maximum and mean minimum temperature on the day of calling (mean importance = 0.505 and 0.487, respectively). Again, moon phase was of moderate significance (mean = 0.429), and showed a stronger relationship to calling than rainfall over the previous 3 days and rainfall on the day of calling (mean importance = 0.238 and 0.018, respectively) (see Figure 4 for species specific results).
We found no distinct grouping of frog species based on their calling with respect to the predictor variables. Frogs did not cluster into groups of species based on the strength of their relationship to predictor variables (i.e., combined effects of recent rainfall and day of year or combined effects of temperature and daily rainfall).
Rather, the frogs fell along a spectrum of shared predictor importance patterns (Figure 4 [for an interactive figure with species names see Appendix S1]). The variance explained by the clustering of frogs based on predictor variables indicated the relationship between calling and the other explanatory variables produced groups with as much in common within clusters as between clusters.

| DISCUSS ION
Using more than 150,000 citizen science observations from 100 species, we demonstrate the importance of day of year and temperature thresholds as predictors of calling behaviour in Australian frogs. Day of year was by far the most important variable predicting calling behaviour at the scale examined, with a maximum predictor importance (PI = 1) in 67% of species examined. Conducted at a continental scale over multiple years, our analysis revealed strongly day-driven seasonal trends in calling across Australian frogs, but also unique species-specific responses to meteorological variables.
One of the reasons for the overriding importance of day of year, a proxy for photoperiod, may be that animals have evolved to respond to photoperiod as a harbinger of other important conditions (i.e., seasonal shifts in temperature and rainfall) (Bradshaw & Holzapfel, 2007). Indeed, for ectotherms, photoperiod tends to be more important to phenology than temperature because it is more consistent (Gotthard, 2001). Our results suggest that this may also be the case for Australian frogs. Although temperature and photoperiod have both been considered significant to calling (Duellman & Trueb, 1994), the importance of photoperiod may be higher when investigating calling on broad temporal scales such as in this study (36 months). In an 18-month study period in Uruguay, photoperiod was similarly recorded as the main predictor of frog calling (Canavero & Arim, 2009). Within a season, some environmental variables may F I G U R E 2 Scaled importance of predictor variables influence on calling (a) for all 100 species, and (b) for the 67 species with a strong dependence on day of year ("seasonal breeders"), restricting data to within each species' breeding season Annual all-species model Seasonal model (day species only) be important (see Oseen & Wassersug, 2002;Saenz et al., 2006;Weir et al., 2005), but across multiple seasons, photoperiod may be the predominant driver of calling behaviour (Both et al., 2008;Canavero & Arim, 2009;Schalk & Saenz, 2016). And the relative strength of photoperiod may depend on the interaction of photoperiod with other variables, as suggested by our analysis where we found the importance of day of year was strongly correlated with temperature ( Figure 3).
While we did not explicitly model the "timing" of frog calling, but rather the presence/absence of calling on a given day, our results suggest that the timing of frog calling (i.e., a proxy for the breeding season) is strongly seasonal in most frog species.
Therefore, these results can be used to quantify approximate breeding seasons for Australian frogs, which are often necessary for biodiversity monitoring (e.g., Roth et al., 2014;Wilson & Bart, 1985). Additionally, while we present the strength of Surprisingly, we found very little influence of rainfall on the probability of calling. The low significance of all three rain variables underscores the findings in some regional Australian studies, where the probability of calling was uncoupled with rainfall (Lemckert & Grigg, 2010), and more related to lagged rainfall (over the previous week) than rainfall on the day of calling (Heard et al., 2015).
Likewise, we found some evidence that cumulative rainfall within a 10-day period was more important than cumulative 3-day rainfall, which, in turn, was more important than mean rainfall on a given day F I G U R E 3 Heatmap of the scaled strength of interactions among predictor variables for the annual all-species model. Interactions were strongest among the variables with the strongest relationship to frog calling ( Figure 2). A negative association between rainfall and calling has even been demonstrated in some temperate environments (Heard et al., 2015;Oseen & Wassersug, 2002). Possible reasons for minimal correlation with rainfall in frog species include increased risk of eggs washing away for lotic breeding frogs (Heard et al., 2015), noise interference from the precipitation (Dorcas & Foltz, 1991;Henzi et al., 1995;Saenz et al., 2006), and competition avoidance with other frog species (Duellman & Pyles, 1983;Heard et al., 2015). Similarly, a global meta-analysis of amphibian phenology also found phenological shifts in calling across taxa were associated more strongly with temperature than precipitation (Ficetola & Maiorano, 2016).
Overall, there is evidence both for and against the positive correlation of rainfall and lagged rainfall on calling patterns -varying by species, time scale, and ecoregion (Canavero & Arim, 2009;Heard et al., 2015;Lemckert & Grigg, 2010;Oseen & Wassersug, 2002;Shalk and Saenz 2006). Indeed, although we found overall weak support for rainfall, this varied among species, and rainfall over the past 10 days was the most important predictor for 15% of species during the breeding season (i.e., Limnodynastes dumerilii, Litoria tornieri, and

Neobatrachus sudellae).
The coarse spatial scale of our study is also likely to have played a role in the strength of the effect of photoperiod over meteorological conditions in our results. BOM data are already interpolated to 5 km 2 grid cells (BOM, 2020) before we smoothed them further in our analysis. Weather is highly variable over space and time and short-duration extreme precipitation (lasting an hour or less and covering only a few km 2 ) cannot be well documented by rain gauge networks (Lengfeld et al., 2020). Weather variables can be difficult to reliably interpolate from weather station observations, especially where weather stations can be thousands of km apart, such as throughout much of arid and semi-arid Australia (Peña- Arancibia et al., 2013). Conversely, at a finer scale, pairing on-site data loggers with observations may result in stronger associations between meteorological variables and calling (Dorcas et al., 2009; F I G U R E 4 PCA plot showing K-means clusters of all species grouped by the similarity of their response to predictor variables. Clusters were not recovered as distinct groups, illustrating the continuum of variable importance across frog species Weir et al., 2005). Together, these results suggest the importance of considering temporal and spatial scale in predicting phenological patterns among species. Improvements in meteorological interpolation may allow more accurate, large-scale analyses in the future. The strength of the macroscale, aggregate analysis presented here does not lie in revealing the influences of extremely localized factors that surely do impact frog survival (Scheffers et al., 2014) and successful reproduction (Blaustein et al., 1999;Harkey & Semlitsch, 1988;Kiesecker & Blaustein, 1998;Watkins & Vraspir, 2006). Our results trend towards common species and broad scale patterns.
Sampling biases are common in citizen science data. In this study, urban and temperate areas were disproportionately well sampled compared with dry and remote regions. Participants most often record in the evening and near areas of high population density Callaghan et al., 2020;Liu et al., 2021). In both model outputs, burrowing and terrestrial frogs (such as Austrochaperinia pluvialis, Austrochaperina robusta, and Heleioporus eyrei) comprised the bulk of the outliers (high and low predictor importance) across all variables. As a result of small ranges or infrequent detection (i.e., large periods of aestivation and often remote locations), these frog species also had fewer grid cells and observations included in analysis and, in some cases, were disproportionately impacted by our spatiotemporal subsampling (Tables S1-S2). For example, while citizen science participants sample frogs all year round in many parts of New South Wales, this is less true in more remote parts of our Australia, potentially influencing our results. When we stratified a subset of species to climate zone, predictor importance varied across climate zones. For species with records in both zones, rainfall variables were often more important in subtropical than temperate zones, while humidity was more important in temperate than subtropical zones ( Figure S5), adding to some of the variation. Additionally, when we stratified the same subset of species by observation year, predictor importance varied somewhat from year to year ( Figure S6), suggesting the importance of using robust multi-year datasets to investigate phenology. While our work was largely focused on temporal differences, future work should test how breeding cues differ in space. Our exploratory analysis suggests that macro-evolutionary constraints (i.e., different evolutionary responses to different climate zones) may influence breeding cues in frogs.
Understanding how meteorological conditions influence the onset of phenological events, such as breeding, is particularly important considering the rapid changes in environmental conditions and stability throughout most of the world, and how important breeding is to species survival, population dynamics, and resilience. While we focus on species-specific responses, community-level data (e.g., species richness) could also better inform risk assessment models. For example, by understanding how the underlying species diversity in space leads to co-occurrence and competition, we can uncover frog communities facing shared risks. Indeed, future work should look to test how the effect of meteorological determinants investigated here influence spatiotemporal co-occurrence. To paraphrase Gwinner and Helm (2003): circannual rhythms are intimately involved in the seasonal organization of breeding behaviour, providing the substrate onto which seasonal environmental factors act. The large volume of data across broad spatial and temporal scales necessary for elucidating phenological patterns is rapidly becoming available through citizen science (Bird et al., 2014;Hochachka et al., 2012;Sullivan et al., 2014). Indeed, citizen science data are increasingly being used to document phenological changes, including spatiotemporal changes in butterfly (Soroye et al., 2018) and bird migration (Hurlbert & Liang, 2012), flowering time (Gonsamo et al., 2013), and patterns in bird call phenology (Dickinson et al., 2010;Sullivan et al., 2014) -and now, frog calling behaviour. To the best of our knowledge, our research represents significantly more data, from a greater number of species, and over a greater timespan than any frog call phenology study to date. While the availability of freshwater breeding sites is vital for frogs, at the scale examined, calling may not be as tightly linked to rainfall for all frog species as often assumed. The correlation between calling and temperature recovered, particularly within breeding seasons, may result in a shift in breeding with climate change (Gibbs & Breisch, 2001), and this has the potential to further affect many frog species breeding success and survival. Our results illustrate the importance of day of the year as a strong, but not isolated, predictor for frog calling behaviour at a macroecological scale.

ACK N OWLED G EM ENTS
We would like to thank the Citizen Science Grants of the

CO N FLI C T O F I NTE R E S T
The authors declare no conflicts of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available here: https://doi.org/10.5281/zenodo.7042152.