Shifts in amphibian population dynamics in response to a change in the predator community

: Predation can affect prey behavior, demography, abundance, and distribution, particularly in lentic freshwater ecosystems. Fish are predators known to reduce the abundance of their prey and to restrict the distribution of species. Using time series which spanned 43 and 22 yr, respectively, we analyzed the effect of a change in the fish predator community on the dynamics of two pond‐breeding amphibian populations (Rana temporaria and Rana dalmatina). Specifically, we used a state‐space time series model which allows for density dependence and observation error, to ask whether the change in predation risk affects population growth rate and the return point around which the populations fluctuate. The results showed that the type of observation error assumed did not affect the biological parameters. We found evidence for density dependence in both populations. The effect of the change in fish predation on population growth rate and the return point was strong in the population where fish invaded a previously fish‐free pond. The effect was weaker in the population where the change was from cyprinid fish to pike. The results showed that fish predation can have strong effects on amphibian population dynamics. The observed population dynamical pattern is phenomenologically similar to alternative stable states. a change in the Abstract. Predation can affect prey behavior, demography, abundance, and distribution, particularly in lentic freshwater ecosystems. Fish are predators known to reduce the abundance of their prey and to restrict the distribution of species. Using time series which spanned 43 and 22 yr, respectively, we analyzed the effect of a change in the ﬁ sh predator community on the dynamics of two pond-breeding amphibian populations ( Rana temporaria and Rana dalmatina ). Speci ﬁ cally, we used a state-space time series model which allows for density dependence and observation error, to ask whether the change in predation risk affects population growth rate and the return point around which the populations ﬂ uctuate. The results showed that the type of observation error assumed did not affect the biological parameters. We found evidence for density dependence in both populations. The effect of the change in ﬁ sh predation on population growth rate and the return point was strong in the population where ﬁ sh invaded a previously ﬁ sh-free pond. The effect was weaker in the population where the change was from cyprinid ﬁ sh to pike. The results showed that ﬁ sh predation can have strong effects on amphibian population dynamics. The observed population dynamical pattern is phenomenologically similar to alternative stable states.


INTRODUCTION
Predation is a biotic interaction that is well known to affect prey populations. Predators canhavebothdirectandindirecteffectsonprey populations. They may reduce the performance of prey species and alter competition within prey communities (Morin 1981, Peacor and Werner 2001, Anholt et al. 2005. Ultimately, the many effects of predation can reduce abundance, restrict the distribution of (potential) prey species, or lead to niche contraction (Paine 1966, Jeffries and Lawton 1984, Sih et al. 1985, particularly in aquatic ecosystems (McPeek v www.esajournals.org 1 May 2021 v Volume 12(5) v Article e03528 1990, McPeek 1994, Resetarits 2005). Predation by fish is one of the main forces that structures freshwater aquatic communities (Wellborn et al. 1996) and that restricts and drives species distributions across the freshwater habitat gradient (Bradford et al. 1993, Eaton et al. 2005, Werner et al. 2007); predation by fish can even affect adjacent terrestrial communities (Epanchin et al. 2010, Rudman et al. 2016. The presence of fish often excludes prey species from otherwise suitable habitats or prey species actively avoid water bodies inhabited by fish (Brönmark and Edenhamn 1994, Resetarits 2005, Winandy et al. 2017. The introduction of fish into previously fish-free habitats can have devastating effects on prey communities (Knapp et al. 2001, Nyström et al. 2001, Kats and Ferrer 2003. The distribution and abundance of amphibians are well known to be strongly affected by fish predation (Kats et al. 1988, Hecnar and M'Closkey 1997, Hartel et al. 2007, Werner et al. 2007. A likely reason is that predation by fish has comparatively stronger effects on amphibians than predation by other predators such as invertebrates, with fish being able to completely eliminate larval amphibian populations (Semlitsch 1993, Leu et al. 2009). It is therefore no surprise that the introduction of fish into habitats that were naturally fish-free is one of the drivers of the global decline of amphibians (Houlahan et al. 2000, Kats and Ferrer 2003, Knapp et al. 2016. The introduction of fish has often led to the local extinction of amphibian populations and the disruption of metapopulation processes (Bradford et al. 1993, Denoel et al. 2005, Orizaola and Brana 2006, but many other effects, including sublethal ones, on individual performance and populations have been documented (Kats and Ferrer 2003). For example, several studies have reported reduced abundance of amphibian populations in the presence of predatory fish, with the strength of the effect depending on the species (Denoel and Lehmann 2006, Orizaola and Brana 2006, Roth et al. 2016. Some amphibian populations, however, may show resilience to fish introductions and may quickly recovery after the removal of nonnative fish (Vredenburg 2004, Pope 2008, Knapp et al. 2016, Tiberti et al. 2019).
Here, we use time series analysis to study the effect of a change in the fish predator community on the population dynamics of pond-breeding frogs. There are only few studies which have documented the temporal response of amphibian populations to a change in predation (Gamradt and Kats 1996, Eaton et al. 2005, Knapp et al. 2007, Tiberti 2018. Specifically, we address the question whether a change in the community of predatory fish leads to a change in the growth rate of amphibian populations and to a change in the carrying capacity (i.e., the return point of a population regulated by density dependence; see Materials and Methods section for details).
We analyzed a 42-yr (1975-2017) time series of egg mass of a common frog (Rana temporaria) population and 21-yr (1997-2018) time series of egg mass counts data of the agile frog (Rana dalmatina), where counts of egg masses serve as a proxy for the abundance of female breeding population size. Meyer et al. (1998) and Bȃncilȃ et al. (2016) analyzed data from the years 1975 to 1997 for the common frog and 1997 to 2011 for the agile frog, respectively, and provided evidence for direct and delayed density dependence in the two populations. During the time that the data sets were collected, fish predation changed in the community of each frog species. In the pond used by the agile frog population, the fish predator community changed from dominated by cyprinid fish to a community dominated by pike (Esox lucius). The change occurred in 2009 (Bȃncilȃ et al. 2016). The common frog population experienced fish predation after 1989 when nonnative goldfish (Carassius auratus) were introduced into the pond (Meyer et al. 1998). While pike prey on both larvae and adult amphibians, goldfish and other cyprinids prey primarily on tadpoles (Glandt 1985, Semlitsch 1993, Bauer and Laufer 2007.
We expected to find similar patterns of population regulation in both populations. We expected that the change in the fish predator community would lead to reduced population growth rates and lower return points in the later years; the latter would lead to reduced abundance. We expected a stronger effect in the common frog population than in the agile frog population because the introduction of fish is expected to have a stronger impact than a change in the fish community.

Study species and study areas
Both the agile frog and the common frog are widely distributed over most of Europe. Females typically lay a single egg mass in early spring (February-March). Tadpoles metamorphose after a larval period of 2-3 months and juveniles reach sexual maturity after about two years (Günther et al. 1996).
The agile frog population is located in central Romania, north of the town Sighişoara. The pond has a surface area of 2.2 ha. The shorelines of the pond are covered by reeds (Typha sp. and Phragmites sp.). For additional information on this population and study area, see Hartel (2008a, b) and Bȃncilȃ et al. (2016). Four fish species occur in the pond: C. auratus, Pseudorasbora parva, E. lucius, and Leucaspius delineatus. The presence of C. auratus, P. parva, and L. delineatus was recorded during the entire period of the study whilst E. lucius was first detected in 2009 and occurred until the end of the study (Hartel et al. 2007, Bȃncilȃ et al. 2016. The common frog population breeds in a pond in the Bermoos, a fen north of Bern, Switzerland. The pond has a surface area of 0.96 ha. Frogs breed in shallow parts of the pond in Carex stands. For additional information on the study site and population, see Meyer et al. (1998). There were no fish in the pond at the start of the study but C. auratus were first observed in 1989. The presence of this fish species was confirmed regularly between 1989 and 2003. We stopped looking for fish in the later years of the study, but we have no reason to assume that fish did not persist. Populations of C. auratus and other fish species generally persist for a long time (B. Schmidt, personal observation), go extinct only under exceptional conditions (Eaton et al. 2005), and are not easy to eradicate (Rytwinski et al. 2019). The pond dried on 2018, and this event may have eliminated the fish population (K. Grossenbacher, personal observation).

Population monitoring
We collected data from 1997 to 2018 (agile frog) and from 1975 to 2017 (common frog). Every year we counted the number of egg masses. The number of egg masses serves as a proxy of the female breeding population (Griffths and Raper 1994, Meyer et al. 1998, Crouch and Paton 2000, Grant et al. 2005, Hartel 2008a. The egg masses were counted along the shoreline because the frogs deposit the egg masses in shallow water. The egg masses were counted each year from mid-February until early-mid-April. The pond was visited and egg masses counted up to six times a year, depending on the frog breeding phenology (Crouch andPaton 2000, Grossenbacher et al. 2002). The maximum count (per year and population) was used in the analysis.

Description of the modeling approach
A simple model for density-dependent population dynamics is Taper 1994, Meyer et al. 1998) N t is population size at time t, the constant a is the population growth rate in the absence of density dependence (i.e., b 1 = 0) and b 1 describes density dependence of the population growth rate. If this model describes population dynamics well, then it can be shown (Dennis and Taper 1994) that the population will fluctuate around a return point (ñ)ñ If population size is greater thanñ, then, on average, the population size will become smaller. Conversely, if the population size is smaller thañ n, the population will increase in size. The Eq. 1 can be extended to include the effect of delayed density-dependence (b 2 ) and an effect of the change in the predator (i.e., fish) community (like any other environmental covariate; Dennis and Otten 2000, Pellet et al. 2006, Pasinelli et al. 2011: The variable "fish" is a binary variable taking value 0 before the change in the fish community and 1 afterward. If the variable "fish" affects the population growth rate (i.e., the regression coefficient b 3 ≠ 0), then it will also affect the return point (Pellet et al. 2006).
Thus, the model shows how a change in the fish community may change the growth rate of an amphibian population and the return point around which population size fluctuates (ñ t ). Population growth rate can take two values. It is a in the first period and a + b 3 in the period after the change in the fish community. Depending on the strength of the effect of fish predation a change in the predator community may not only reduce population growth rate but may lead to the exclusion of an amphibian species from a habitat patch where the predator occurs. In doing so, fish predation may restrict the distribution of amphibians (i.e., whenñ t = 0 in the presence of predators).
Because there is generally observation error in amphibian time series (Schmidt 2003), we explored how different models for the observation error affect inference. We had no a priori expectations as to how the model for observation error would affect estimates of the parameters which are of biological interest.

Time series analysis
We used a state-space model (de Valpine and Hastings 2002, Dennis et al. 2006, Lebreton and Gimenez 2013, as extended by Dennis and Otten (2000) and Pasinelli et al. (2011) and previously used for amphibian populations (Tobler et al. 2012, Buckley et al. 2014, Bȃncilȃ et al. 2016, to analyse the time series data. The state process equation describing the unknown true population size, N t, at time t is the same model as developed above (Eq. 3), but arranged in slightly a different way and including process variance: where a is the intercept, b 1 and b 2 are the estimates of the strength of direct and delayed density dependence, b 3 describes the effect of change in the fish community and σ 2 proc is the process variance. The model assumes an immediate change in population dynamics after the change in the fish predator community. In reality, the response may be slower but we feel that this simplification is justified because it makes the model more conservative. The second equation links the state process with the observations (X t ), that is, counts of egg masses. We adopted three different descriptions of the observation process: 1. Normal observation error: where σ 2 obs is the observation variance. 2. Binomial observation error: where p is the detection probability. There is no information in the data that can be used to estimate detection probability. Therefore, we used an informative prior (p~Beta(12,2)). The prior was chosen such that its mean (0.86) is close to observed detection probabilities for frog egg mass counts (Grant et al. 2005).
3. No observation error. To fit this model to the data, we fixed p at 1 in Eq. 7.
Based on the previous published analyses (Meyer et al. 1998, Bȃncilȃ et al. 2016), a preliminary analysis using a partial rate correlation (Berryman and Turchin 2001), and visual inspection of the data, we selected different models for the two frog species. For the agile frog, we selected a model with direct and time-lagged density dependence (i.e., Eq. 5) while for the common frog, we dropped the time-lagged effect from Eq. 5. For both species, we allowed process variance (σ 2 proc ) to differ between the two periods (before/ after the change in the fish community).
We fitted the models to the data using R and JAGS as described in Kéry and Schaub (2012). We used normal priors with a mean of zero and a precision of 0.0001 for a, b 1 , b 2 , and b 3 . We used uniform priors in the interval 0-500 and 0-10 for the standard deviations of σ 2 obs and σ 2 proc , respectively. Note that population size in the first and second year is not defined in Eq. 5. For N 1 and N 2 we specified a normal prior with mean equal to the log-transformed count and a variance of 100. We ran three parallel MCMC chains with 200,000 iterations and a burn-in of 100,000. Chains were thinned by a factor 2. Convergence was assessed using the Gelman-Rubin diagnostic R-hat. Convergence was satisfactory when R-hat values were smaller than 1.1 (Roy 2020). The v www.esajournals.org time series data and the JAGS code are available in Appendix S1 and Data S1.

RESULTS
Figs. 1 and 2 show the counts, the estimates of population size, and the estimated return points for the two time periods, that is, before and after the change in the fish community. The parameter estimates, the 95% credible intervals (CI), and the proportion of the posterior distribution of the parameters, which has the same sign as the mean (f; Wade 2000) are shown in Table 1.
In the common frog population, there was strong evidence for direct negative density dependence (as judged by the f values;  Fig. 1 shows that the population fluctuates more widely). The type of observation error did not have a strong effect on the parameter estimates (Table 1).
In the agile frog population, there was strong evidence for both direct negative density dependence and for delayed negative density dependence (i.e., an effect of N t−2 on N t , Table 1). While the estimate of delayed density dependence was negative, the 95% CI slightly overlapped zero but the f were greater than 0.96. Even though the 95% CI of the estimate of the fish effect (b 3 ) slightly overlapped zero, the f values were greater than 0.95 (Table 1). Fish had a negative effect on abundance (Fig. 2). The return point (ñ) changed from 400 (median; 95% CI 307, 505) to 328 (median; 95% CI 206, 465; values are based on the model with normal observation errors). Process variance was much lower in the first years of the study before the fish community was dominated by pike ( Table 1). The type of observation error did not have a strong effect on the parameter estimates (Table 1). Effects were generally weaker in the agile frog than in the common frog. The reason is most likely that (1) the time series is shorter, (2) the model is more

Rana temporaria
Year Number of egg masses

DISCUSSION
Predation is well known to affect the distribution, abundance, and behavior of prey (Kerfoot and Sih 1987). Here, we documented changes in the population dynamics of two species of pondbreeding frogs in response to a change in the fish predator community. The amphibian populations were regulated, that is, there was densitydependent population growth. The change in the fish predator community negatively affected the population (i.e., the estimate of b 3 in Eq. 5 was negative) and led to a lower return point, and to a lower mean abundance. In addition, the magnitude of population fluctuations changed. The population of the common frog fluctuated less after the reduction in population size caused by the change in the predator community whereas the agile frog population had higher process variance after the change in the fish community (but this was caused by the exceptionally high abundance in 2017). Overall, the pattern is similar to alternative stable states which are commonly observed in aquatic ecosystems (Scheffer 2009).
The biological results we report were not affected by how we modeled observation error. Observation error is ubiquitous in animal abundance data (Schmidt 2003) and should be accounted for in the analysis because it can affect inference (Dennis and Taper 1994, Link and Nichols 1994, Linden and Knape 2009, Lebreton and Gimenez 2013. We modeled three types of observation error that are commonly used for time series of population counts (Pellet et al. 2006, Hostetler and Chandler 2015, Amburgey et al. 2018. The results showed that the model used to account for observation error had little effect on estimates of the parameters of biological interest. We believe that this was the case because there is generally little observation error in counts of egg masses of pond-breeding frogs, which breed early in the season and are explosive breeders. Several studies have found that detectability of egg masses was very high and did not vary much among years (Griffths and Raper 1994, Grant et al. 2005, Scherer 2008, . The horizontal solid lines show the return points (ñ), and the dashed lines represent 95% CI for the two time periods, that is, before and after the change in the fish community. The arrow shows when the fish community changed. Fellers et al. 2017). Thus, observation error in this type of count data may be negligible. Other types of count data may have much higher levels of observation error (Schmidt 2003), making data analysis more challenging.
While there is ample evidence for density dependence in amphibian populations from experimental studies (Alford 1999), there is much less evidence for density dependence at the population level derived from individual mark-recapture data or time series data (Berven 1990, Salvidio 2009, Cayuela et al. 2019). Here, we found strong support for density dependence in the populations of both species. This is important, because it shows that amphibian populations can be regulated by intrinsic processes despite the erratic fluctuations in population size, which are often observed (Meyer et al. 1998, Pellet et al. 2006, Salvidio 2009). The regulation suggests that populations may be buffered, at least to some extent, against external perturbations (Lebreton and Gimenez 2013). Here, the buffering affected the responses of the populations to changes in the fish community. Predation by fish reduced the size of the population but it continued to fluctuate in a density-dependent manner around a lower return point. The lowered return point was likely set by the densely vegetated shallow edges of the ponds where tadpoles are likely to suffer reduced predation rates by fish (as opposed to open water). Shallow slopes with dense vegetation may mitigate the effects of fish introductions into ponds and may avert local extirpation of amphibian populations (Hecnar and M'Closkey 1997, Hartel et al. 2007, Shulse et al. 2012, Hartman et al. 2014, Kloskowski et al. 2020, the most commonly observed effect of fish introductions (Knapp et al. 2001, Kats and Ferrer 2003, Denoel et al. 2005. It is important to note, however, that the time series models cannot distinguish between density dependence at the larval or adult stages and that delayed density We also provide the probability that the parameters differ from zero (f). a is intrinsic rate of increase, b 1 and b 2 are the strength of direct and delayed density dependence, respectively, b 3 is strength of the effect of fish, σ 2 obs is the observation variance, σ 2 proc is the process variance, and p is detection probability.
[1] and [2] refer to the two time periods (before and ofter the change in the fish community).
dependence should not be accepted as evidence for density dependence at the larval stage (Meyer et al. 1998). Our conclusions regarding the underlying mechanisms of the reduced return point are therefore tentative. Identifying the stages at which density dependence occurs in natural amphibian populations remains an important challenge for the future. Despite evidence for regulation, we observed years (2009 for common frog and 2017 for agile frog) with high abundance in the later years of the time series when return points were low. We suspect that the peak years occurred after years with exceptionally high recruitment but we cannot pinpoint the mechanism which may have led to high recruitment despite fish predation (apart from stochastic events; Eaton et al. 2005). Despite the occurrence of single years with high abundance, the lower return point was also associated with reduced variance in population size in the common frog. This contrasts with the results of Green (2003) who reported higher variance in smaller amphibian populations. Reduced variability in population size may reduce extinction risk after the change in the fish predator community (Inchausti and Halley 2003).
The two populations showed the same general response to the change in the fish predator community but the magnitude of the response was different. The effect was very strong in the Swiss common frog population where nonnative fish were introduced into a formerly fish-free pond while the effect for the Romanian agile frog population was weaker. The posterior distribution strongly suggested an effect of the change in the fish predator community but the 95% CI of the parameter estimate overlapped zero. It is worth mentioning that the effect of fish predation would be strong without the year 2017, when the number of egg masses counted was exceptionally high. However, the agile frog population did not experience a sudden change from no fish predation to fish predation but rather a shift from predation by cyprinid fish to predation by (mostly) pike. In contrast to the other fish species, pike are large enough to prey on both tadpoles and adults. The population dynamic change seen in this population may be equivalent to the effects of the transition from no fish to non-predatory to predatory fish on amphibian distributions at the landscape level (Hecnar and M'Closkey 1997).
From this observation, we tentatively conclude that the factors that cause patterns in the distribution of species may also cause changes in the abundance of species.
Most studies have previously reported the extinction of amphibian populations after the introduction of fish. Here, we showed that a change in predation pressure by fish can affect amphibian population dynamics in a way that is phenomenologically similar to alternative stable states. The change in the return points can also be viewed as a niche contraction (Scheele et al. 2017) where heterogeneity in prey species, predator, and habitat can lead to variable outcomes, as seen here by the differences between the two populations. Our results lead us to conclude that responses to perturbations can be equally population specific. This corroborates the result of Amburgey et al. (2018) who reported population-specific responses to environmental responses. Nevertheless, since the overall response-al o w e rr e t u r np o i n t -was the same in both populations and only the magnitude of the response differed, we may begin to explore the context-dependency of the response. Ultimately, understanding context-dependency may lead to a more general and predictive ecological science (Nichols et al. 2020).

ACKNOWLEDGMENTS
We thank Claudio Bozzuto, Rebecca McCaffery, and anonymous reviewers for comments on the manuscript andÁrpád Szapanyos who helped monitoring the Rana dalmatina population since 2014. We did not receive funding for this study. BRS and RIB conceived the ideas; BRS and MS designed the analysis; TH and KG collected the amphibian data; BRS analyzed the data and led the writing of the manuscript. All authors contributed critically to the drafts and gave final approval for publication. We declare no conflicts of interest. No permits were necessary for this study, because no animals were captured.