Climate change is shifting the phenological timing, duration, and temporal overlap of interacting species in natural communities, reshaping temporal interaction networks worldwide. Despite much recent progress in documenting these phenological shifts, little is known about how the phenologies of species interactions are tracked across different life-history stages. Here we analyze four key phenological traits and the pairwise interaction potential of nine amphibian species for the adult (calling/breeding) and subsequent larval (tadpole) stage at eight different sites over six years. We found few strong correlations among phenological traits within species, but the strength of these correlations varied across species. As a consequence, phenological trait combinations of both stages varied substantially across species without clear signs of multidimensional clustering, indicating a distinct and diverse range of species-specific phenological strategies. Despite this considerable variation in the phenologies across species, the temporal overlap between species was largely preserved through the two life history stages. Further, we also detected significant correlations among the duration and temporal overlap of interactions with other species across stages in five species, demonstrating that temporal patterns of species interactions are mirrored across life history stages. For these species, these results indicate a strong tracking of phenologies and species interactions across life history stages even in species with complex life cycles where stages occupy completely different environments. This suggests that phenological shifts in one stage can impact the temporal dynamics and structure of interaction networks across developmental stages.

Study system

Studying the phenology of amphibian communities has several advantages and permits a highly sensitive and reliable approach to linking phenological distributions to competitive interactions. First, timing and duration of amphibian phenophases (life history stages) respond to tractable climatic cues like rainfall, temperature, and air pressure, but which cues are most important varies across species (Blankenhorn 1972, Pechmann et al. 1989, Kopp and Eterovick 2006, Saenz et al. 2006). Between-year fluctuations in weather thus precipitate some of the strongest species-specific phenological shifts of any taxa (Blaustein et al. 2001, Forchhammer et al. 1998, Oseen and Wassersug 2002, Parmesan 2007, Saenz et al. 2006). Secondly, amphibian species exhibit diverse life-history strategies, the phenologies of fall and winter-breeding amphibians are distinct from those of summer or spring-breeding species (Hocking et al. 2008, Parmesan 2007, Saenz et al. 2006). Finally, spatial competition among amphibians of different species is known to occur locally in multiple life history stages. Adults calling during the breeding phenophase produce audio interference that can interfere with heterospecifics’ search for potential mates (Schwartz 1987). Larvae and adults also compete over common resources such as food and space, yet the two stages experience different selective pressures between their distinct habitats and differ in which season they are most active (Alford and Wilbur 1999, Dayton and Fitzgerald 2001).

Data collection

We collected data on tadpole abundance and adult calling activity of 12 amphibian species over six years in eight ponds in Southeast Texas. We obtained estimates of tadpole abundance from funnel traps that were placed underwater in each pond and checked weekly. Two traps were deployed at each pond for the entire study, regardless of size of pond. Traps were placed in the littoral zone of each pond, where tadpoles tend to concentrate (Porej and Hetherington 2005). The traps were placed in shallow water so that the tops of the traps were exposed out of the water to allow any accidental bycatch species to reach the water surface to breathe. The entrance to the traps was always completely under the water’s surface. Abundances of all species were recorded and individuals returned to ponds, so recapturing individuals later was possible. We tracked adult calling activity and estimate abundance from audio recorders installed at each pond that recorded amphibian calling for one minute at six timestamps (9:00 pm, 10:00 pm, 11:00 pm, 12:00am, 1:00am, and 2:00am) each day in accordance with the times when species were actively calling (Bridges et al 2000). Manual processing of audio data involved using a sonogram reference to identify the number of calling species in a one-minute period. Up to five individuals were distinguishable from each minute-long recording, so the maximum number of individuals was bounded at 30, in line with the procedure outlined in Saenz et al. 2006. We summed observations of adults across all six recording times in a day to acquire measurements of total abundance for each day from May 2000 to December 2015. The weekly observations of funnel trap count data for tadpoles spanned May 2001 through December 2006, therefore this range was analyzed for both data types. Using 6 years of data has been shown to accurately capture typical phenological dynamics (Brown et al. 2016), but long-term trends are out of the question given that perceived declines may represent between-year fluctuations or site-selection bias (Blaustein et al. 1994, Bridges et al. 2000, Fournier et al. 2019). We obtained records of 12 species of amphibians (Hyla versicolor, Hyla cineria, Bufo valliceps, Bufo woodhouseii, Rana catesbeiana, Rana clamitans, Rana sphenocephala, Gastrophryne carolinensis, Pseudacris crucifer, Pseudacris triseriata, Acris crepitans and Rana palustris). However, some rare species (H. cineria, G. carolinensis, P. triseriata, A. crepitans, and R. palustris), with fewer than 15 days of observation across all ponds and years, were not included in our analysis of among-species phenological dynamics. In this community, the average larval duration preceding metamorphosis typically ranges from 20–35 days (B. valliceps, G. carolinensis, H. versicolor), 55-70 days (A. crepitans, H. cinerea, P. trisereata, R. sphenocephala), and 90 days (P. crucifer, R. clamitans), with R. catesbeina having a duration of 365 days (see Fig. 14 in Saenz 2004). However, tadpole larval period can be strongly influenced by biotic and abiotic factors and thus can be much shorter or longer depending on specific environmental conditions. Climatic metadata, specifically records of temperature, rainfall, and photoperiod length for each pond, could not be included because we did not track among-site abiotic data. Therefore, we could not investigate correlations between phenological shifts in interactions and climatic changes.

The eights ponds are divided equally between two neighboring national forests in South East Texas, with four located in Stephen F. Austin Experimental Forest (SFA), a disjunct unit of the Angelina National Forest under the control of the Southern Research Station, and four in the Davy Crockett National Forest (DC) (see Fig. 1 of Perez et al. 2021 for a map). Ponds in DC were built in 1992 with surface area ranging from 900–2000 m² and with a maximum depth of 2.5 m, while SFA ponds were built in 2000 and range in surface area from 500–600 m² with a maximum depth of 1 m (Saenz et al. 2006). Ponds in both forests support diverse communities, including many aquatic invertebrates, but only the larger DC forest ponds can support fish such as mosquitofish, green sunfish, and largemouth bass. There was not a significant difference in weather conditions among both forest types (temperature: p=0.155, mean ± SEM 19.08 ± 8 C° at DC, 18.09 ± 8 C° at SFA rainfall: p=0.212, mean ± SEM, 103.8 ± 26 cm/year at DC, 104.2 ± 29 cm/year at SFA), nor did these metrics shift significantly between years (temperature: p=0.777, R2=0.034, rainfall: p=0.401, R2=0.013). Within years, total rainfall levels differed by at most 91 cm between areas (in 2002, SFA saw 141 cm and DC had 50 cm). Thus, we did not include “forest” as a predictor in our final model (without any loss in model performances). The ponds represented 8 independent experimental units because an average distance of 17km between ponds made the dispersal of juveniles or adults between ponds highly unlikely within a season.

The shape of each species’ phenological distribution is highly unique, creating large variation in the potential for between-species interactions (Fig. 1a). Out of 9 species with sufficient observations to apply smoothers, 6 displayed one distinct larval activity peak in the spring or summer (Fig. 1a). We also detected the bimodal life histories of R. sphenocephala, R. catesbeina, and R. clamitans, which showed two distinct activity peaks separated from each other in time, with a primary peak in the spring and summer, for the first two, and a secondary peak in the fall for all three. These peaks were preceded by similar adult calling peaks in most species, but R. catesbeina and R. clamitans demonstrated a singular, broad, and cryptic peak over several spring and summer months, reflecting a highly flexible life history strategy that can yield multiple cohorts and the possibility of recovering tadpoles throughout the entire year. Larvae of R. clamitans and R. sphenocephala overwinter in ponds due to relatively warm year-round temperatures (Saenz 2004).  Adults of all species do not enter a period of diapause or dormancy during the year, given the mild winters of the study region.  Differential preferences in habitat choice among species caused some not to co-occur in the same forests, H. versicolor and R. sphenocephala tadpoles were common at SFA but rarely observed in DC ponds, while the opposite was true for R. clamitans (Fig. 1b). Pond 7 was removed from the analysis of phenological traits due to extremely low tadpole abundance over the years.

Response variables

We worked with time series count data of twelve species covering ~6 years for each of the eight ponds. Temporal resolution differed between adult calling data and tadpole data, since calling observations were made daily in contrast to the weekly checks of funnel traps. Thus, we converted daily adult calling data to a weekly time scale to match the tadpole sampling scale. We found that weekly resolution still showed the same general patterns and provided sufficient temporal resolution for model-fitting without significantly affecting model diagnostics. We analyzed four metrics of phenological traits for each species and stage in every pond and year where they were observed. These “single-species metrics” were first day and median day of observation, temporal range, and integrated phenophase area. We used Julian day of the year (1—365) rather than the calendar date to calculate these metrics. This allowed us to directly compare years, ponds, and species on the same temporal scale. First day was defined as the day when an individual of a given species and stage was first observed at a given pond in a year. We chose to use the absolute first date rather than another metric (e.g. third day of observation, 5% quantile) since these other metrics seriously truncated our tadpole distributions, many of which were quite narrow and were calculated on a sensitive, weekly grain. Median day represents the day when 50% of the total abundance of each stage of a species was reached in a given year and pond. These two traits captured the timing of phenologies, while the temporal range metric targeted their duration. The temporal range metric (duration) represented the total number of days in a year when at least one individual of a species was observed. For adults, daily observations allowed us to calculate a sensitive days-duration metric. However, tadpoles were sampled on a weekly basis, thus duration measures were in increments of 7 days (weekly measures of observation were multiplied by seven). No species in our system can complete development in under 12 days, however, in some extreme cases duration may be underestimated or overestimated. If tadpoles were only observed during one sampling period but not after, we assumed that they were present for 1.5 weeks (10.5 days) since we do not know whether they reached metamorphosis right after the first trap check or right before the second. This assumption did not bias our duration results since all measurements of larval duration were standardized to the same scale. We calculated phenophase area by integrating under a Lowess-smoothed curve (f = 1/50, iter = 3, delta = 4) of the count time series data in a given year and pond. Lowess smoothing was done following Carter et al. 2018 using the Lowess and integrate.xy functions from R packages sfsmisc in R version 4.0.3. This metric represents the duration of a phenophase weighted by abundance of individuals in a given life history stage.

Measuring pairwise temporal overlap of species phenologies provides a holistic, community-level assay of phenological interaction potential across years and ponds that is more accurate and powerful than single trait metrics to measure the duration of species’ co-occurrence (Post 2019, Carter et al. 2018). We calculated the area of temporal overlap between Lowess-smoothed phenological curves of all possible species-pairs for each pond and year and derived a metric, the proportional overlap, for both species in each pair. Proportional overlap is calculated as the overlap area divided by total the area under the phenology curve of a given focal species (A_overlap/A_focal) and represents how much of a focal species’ overall distribution is accounted for by the area of overlap between both species and fits a zero-one inflated beta distribution. An overlap of 1 means that the entirety of a focal species’ distribution overlaps with the other species, often implying that the focal species distribution is small and contained within the other species’ distribution. 164 species pairs showed non-zero temporal overlap with 20 having complete overlap. In 18 of those 164 cases, two species co-occurred on the same day, but those days were not consecutive. The integration function interpreted non-consecutive days of co-occurrence as zero temporal overlap. However, since these species still co-occurred, we did not eliminate these apparent observations of zero overlap from analysis. Additionally, a metric of “median difference” was calculated to account for differences in the timing of peak abundance between species pairs. For this, the median days of observation for the species in a pair were subtracted and then standardized by the total day range of the focal species.

The phenologies of some species spanned the winter season, December through January, bridging two years. Thus, a year with a cutoff of day 365 did not capture the full phenological distributions of both the tadpoles and adults of these species. We rescaled the data so that the time periods analyzed contained the full distributions of both adults and tadpoles for one cycle. These “relative years” began with the first date of adults calling and ended either 380 days from that point or on the last day of tadpole observation prior to the onset of a new adult calling period. We only considered instances where both tadpoles and adults of a species were found in a given year. As a result, some species had relative year metrics for the entire year range, while others had only two relative years, each relative year potentially disjunct in time from the other by several years. 

The same procedure was used to capture relative interaction envelopes that contained the full distributions of both species in a pair. The left bound of this interaction envelope was set as the first day of adult calling for the earlier-calling species in the pair. The right bound was set as the last day of tadpole observation for the later-hatching species, ending this relative interaction envelope before the initiation of another adult breeding peak for the earlier-calling species in the pair.

Phenological metrics such as duration of phenophase, interaction potential, and median day of tadpole observation were calculated relative to the day of first adult calls. However, as these metrics required the presence of both tadpoles and adults in a given year and pond, sample size was reduced relative to data considering only tadpole interactions: from n=129 to n=76 in the case of single species metrics and from n=164 to n=88 for species pairs.

Single-species metrics

We first examined the factors that drive the differences in phenological distributions among species to understand patterns in phenology between species and stages. To handle unequal variance across species, we fit linear models with generalized least squares (GLS) (nlme package in R) to investigate the influence of species identity, site, and year on single-species phenological metrics for both tadpoles and adults. We fit separate models for first date, median date, and day range, log-transforming the first date response variable and weighting variance by species identity. We fit a second set of models to combined adult and tadpole data to test for an interaction between species and stage identity on the same single-species metrics. Models and corresponding error structures were selected based on diagnostic residual plots, as well as Cooks-distance measurements, pseudo-R-squared calculations, and AIC comparisons. We used the Anova function in car package to determine significance of predictors. Significant predictors would indicate that timing and duration of each stage depends on some combination of differences between sites, years and the inherent life histories of species. The likelihood-ratio test was applied and estimated marginal means (emmeans function from package emmeans) were calculated for all models to compare typical timing and duration between species, ponds, and years.

We conducted non-metric dimensional scaling (NMDS) to determine whether species showed distinct phenology strategies. We used the median date, first date, and days duration variables averaged across all ponds and standardized to the (0,1) interval by subtracting the mean of each trait from each individual observation and dividing by the range of the trait. We then applied metaMDS from vegan set to 3 dimensions (3 dimensions showed the lowest stress measurements) and plotted the coordinates for each species. This cohesive approach combined the information on timing and duration we had on both adults and tadpoles, therefore we expected any clusters to represent overarching differences in phenological strategies for breeding and hatching. We included vector loadings for each trait used in ordination to determine which traits were important to creating phenological differences among species.

Consistencies in phenological patterns across stages

Correlation coefficients between the phenological traits of each stage were calculated using repeated measures correlation coefficients in the rmcorr package, controlling for the species’ identity and using the built-in bootstrapping method to calculate confidence intervals (set to 100 resamples). We used GLMs to test whether adult phenology traits predict tadpole phenology traits, with a given response tadpole metric as the dependent variable and the same metric for adults as the predictor, along with species and pond identity and the breeding period being considered (all predictors were fixed effects). We only included species with at least five non-zero observations in both stages in the same pond and year to acquire power for an interaction in these models: H. versicolor, R. sphenocephala, B. woodhousii, R. clamitans, and P. crucifer. Relative median date, duration, and area for tadpoles were log-transformed to improve model performance, but we did not fit a model to the first day of observation, as the relative first day for adults was always 1 and thus no variability existed in that predictor. Interaction effects between adult metrics and species identity on tadpole metrics were initially included in GLM models but did not significantly determine any metric or help improve model performance and thus were dropped for final analysis. To test whether variation in the timing and duration of adults affects the same phenological traits in tadpoles, we used the Anova (type 2) function in car package to assess the correlation of phenology metrics between stages. Linear relationships for main effects and species-specific effects between adult and tadpole metrics were calculated and compared in post hoc contrasts using the emtrends function from R package emmeans. A significant correlation between timing or duration of the two stages would indicate that shifts in adult timing and duration are tracked in tadpole phenologies.

Correlations in temporal overlap between species and stages

To determine what factors explain the interaction potential, or temporal overlap, between species within each stage we fit a beta regression to proportional overlap between species pairs. This was done separately for tadpoles and adult stages using the betareg function in the R package betareg. Since proportional overlap data contained both zeros and ones (see section on Response variables), we applied a transformation recommended for zero-one-inflated beta distributions in betareg (Douma and Weedon 2019): (y(n-1)+.5)/n. We did not include year in this model because only 6 years of data were available and we therefore did not expect fluctuations among years to help explain species interactions. The final model included tadpole (or adult) proportional overlap, the identity of the species pair, and the pond, and the difference between median date of observation for the two species. Including median difference accounted for differences in timing between species. We used the joint.tests function to assess the significance of predictors, and calculated marginal means. A model was also fit to combined adult and tadpole data to test for differences in interaction potential based on stage, with the same checking process applied. A significant interaction would indicate marked differentiation in the interactions between species in each stage.

Temporal patterns of interspecific interactions across stages

We modeled how tadpole and adult interaction potential, or temporal overlap, is correlated. We fit a beta regression similar to the one described in the previous section (betareg in R) to non-zero observations (with the same transformation from Douma and Weedon 2019 described in the previous section) of the 10 species pairs with enough power (at least 5 observations) for a potential pair-specific relationship between overlap in age classes. We included adult proportional overlap, pair, and site identity as fixed effects. An interaction term between adult overlap and pair identity reduced model performance and did not significantly determine variation in tadpole overlap. Therefore, simple species-pair-specific slopes were calculated between adult and tadpole overlap to approximate the relationship between overlap in stages for each species pair. We used the L.R. Test and Anova functions to estimate the significance of model terms and cplot (effects in R) to plot model predictions. Predicted values were estimated on the scale of the response and compared to raw data. A significant relationship between overlap in the two age classes would indicate that tadpoles track shifts in the interspecific interactions between adults.

All files are in Microsoft Excel. Data analysis was performed in R Studio.