DOI: 10.25338/b8g92p
title: Tracking phenological distributions and interaction potential across life stages

description: 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.
description: <strong>Study system</strong>
 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).
 <strong>Data collection</strong>
 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 (<em>Hyla versicolor, Hyla cineria, Bufo valliceps, Bufo
 woodhouseii, Rana catesbeiana, Rana clamitans, Rana sphenocephala,
 Gastrophryne carolinensis, Pseudacris crucifer, Pseudacris triseriata,
 Acris crepitans</em> and <em>Rana</em>
 <em>palustris</em>). However, some rare species (<em>H.
 cineria, G. carolinensis, P. triseriata, A. crepitans</em>, and
 <em>R. palustris</em>), 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
 (<em>B. valliceps, G. carolinensis, H. versicolor</em>), 55-70
 days (<em>A. crepitans</em>, <em>H. cinerea</em>,
 <em>P. trisereata, R. sphenocephala</em>), and 90 days
 (<em>P. crucifer, R. clamitans</em>), with <em>R.
 catesbeina </em>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<sup>2</sup>
 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<sup>2</sup> 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 <em>R. sphenocephala, R.
 catesbeina</em>, and <em>R. clamitans</em>, 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 <em>R. catesbeina</em> and
 <em>R. clamitans</em> 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
 <em>R. clamitans</em> and <em>R.
 sphenocephala</em> 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, <em>H.
 versicolor</em> and <em>R. sphenocephala</em> tadpoles
 were common at SFA but rarely observed in DC ponds, while the opposite was
 true for <em>R. clamitans</em> (Fig. 1b). Pond 7 was removed
 from the analysis of phenological traits due to extremely low tadpole
 abundance over the years. <strong>Response
 variables</strong> 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
 <em>sfsmisc</em> 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<sub>overlap</sub>/A<sub>focal</sub>) 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. <strong>Single-species
 metrics</strong> 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) (<em>nlme</em> 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 <em>Anova</em> function in
 <em>car</em> 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 (<em>emmeans</em>
 function from package <em>emmeans</em>) 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
 <em>vegan</em> 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.
 <strong>Consistencies in phenological patterns across
 stages</strong> Correlation coefficients between
 the phenological traits of each stage were calculated using repeated
 measures correlation coefficients in the <em>rmcorr</em>
 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: <em>H.
 versicolor, R. sphenocephala, B. woodhousii, R. clamitans</em>, and
 <em>P. crucifer</em>. 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 <em>Anova</em> (type 2) function in
 <em>car</em> 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 <em>emtrends</em>
 function from R package <em>emmeans</em>. 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. <strong>Correlations in temporal
 overlap between species and stages</strong> 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 <em>betareg</em> function
 in the R package <em>betareg</em>. 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 <em>betareg</em> (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.
 <strong>Temporal patterns of interspecific interactions
 across stages</strong> 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
 (<em>betareg</em> 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 <em>Anova</em>
 functions to estimate the significance of model terms and
 <em>cplot</em> (<em>effects</em> 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.
description: All files are in Microsoft Excel. Data analysis was performed in
 R Studio.