Estimating the survival of unobservable life stages for a declining frog with a complex life history

. Demographic models enhance understanding of drivers of population growth and inform conservation efforts to prevent population declines and extinction. For species with complex life histories, however, parameterizing demographic models is challenging because some life stages can be dif ﬁ cult to study directly. Integrated population models (IPMs) empower researchers to estimate vital rates for organisms that have cryptic orwidely dispersing early life stages by integrating multiple demographic data sources.For a stream-inhabit-ing frog ( Rana boylii ) that is declining through much of its range in Oregon and California, USA, we collected egg-mass counts and capture – mark – recapture data on adults from two populations in California to ﬁ t IPMs that estimate adult abundance and the survival rate of both marked and unobserved life stages. Estimates of adult abundance based on long-term monitoring of egg-mass counts showed that study populations ﬂ uctuated greatly inter-annually but were stable at longer timescales (i.e., decades). Adult female survival during 5 – 6 yr of capture – mark – recapture study periods was nearly equal in each population. Survival rate of R. boylii eggs to the subadult stage is low on average (0.002) but highly variable among years depending on post-oviposition stream ﬂ ow. Population viability analysis showed that survival of adult and subadult life stages has the greatest proportional effect on population growth; the survival of egg and tadpole life stages, however, is more mal-leable by management interventions. For example, simulations showed head-starting of tadpoles, salvaging stranded egg masses, and limiting aseasonal pulsed ﬂ ows could dramatically reduce the threat of extirpation. This study demonstrates the value of integrating multiple demographic data sources to construct models of population dynamics in species with complex life histories.


INTRODUCTION
Species extinction is occurring at an unprecedented rate (Barnosky et al. 2011, IPBES 2019. Actions to protect whole ecosystems can be effective to slow biodiversity loss, especially when a master variable, such as the flow regime of a river, has profound influences on survival of many taxa (Poff et al. 1997, 2010, Palmer and Ruhi 2019. Still, extinction can be idiosyncratic (Colchero et al. 2019, Howard et al. 2019. Demographic knowledge at a fine-grained level of detail is essential (Akçakaya and Sjögren-Gulve 2000), because species may face disparate stressors across their ranges, and life stages within a species can display incongruent risks from similar stressors (Otero et al. 2011, Sturrock et al. 2020. Structured population models, which condense demography into a few key age-or stage-specific vital rates, have proven useful for identifying the environmental factors, biological traits, and anthropogenic stressors driving growth or decline (Brook et al. 2000, Caswell 2001). Here, we develop an integrated population model (IPM) to quantify the efficacy of measures, including the management of water flow in rivers with dams, designed to minimize extinction risk of an imperiled frog. We also address the difficulties of vital rate estimation for hard-to-study life stages.
Because of complex life histories and cryptic life stages, vital rates can be challenging to quantify for a wide range of taxa that are emblematic of the current biodiversity crisis, from lepidoptera to amphibians (Dirzo et al. 2014, Scheele et al. 2019. The life history of a free-living larval stage followed by a dispersing post-metamorphic stage can complicate parameterization of models (Williams et al. 2002, Oberhauser et al. 2017, Earl 2019. In particular, survival of recently metamorphosed individuals can be difficult to estimate in the wild when recapture rates for marked individuals are low (Kupferberg et al. 2009). This has pushed researchers to use in situ enclosures that might not fully reflect natural conditions (Altwegg and Reyer 2003). When vital rates are highly uncertain and models are parameterized using independent data sources, demographic models, such as matrix models, can produce unrealistic estimates of population growth rates (Schaub and Abadi 2011). Despite the challenges posed by high uncertainty in estimates of vital rates for early life stages, many authors have used stage-based population models to study the demography of amphibians (Earl 2019). Integrated population models can partially resolve the aforementioned challenges by explicitly integrating different data types into a joint likelihood and can produce demographic estimates with greater precision and accuracy than independent models (Abadi et al. 2010). Integrated population models can enable the estimation of survival and potentially can uncover environmental covariates of survival for unobserved life stages (Abadi et al. 2017).
Integrated population models hold promise for improving conservation of species with breeding migrations and complex life cycles, such as the Foothill Yellow-Legged Frog (Rana boylii). This stream and river dweller, endemic to Oregon and California (USA), has a life cycle synchronized with the region's winter flood/summer drought flow regime and punctuated by spring migration of adults from resident habitats to lek sites where they mate and lay eggs. Some populations inhabit tributary streams too cool and shaded to support timely larval development; wider sun-lit channels where tadpoles graze on abundant periphyton are necessary for metamorphosis prior to the next flood season (Catenazzi and Kupferberg 2013). Thus, the quantity of clutches of eggs laid in mainstem channels reflects adult female abundance in adjacent tributaries. Several factors have driven the extirpation of R. boylii from almost two thirds of its historical range (Lind 2005). Flow regulation by dams creates stressors for early life stages (Kupferberg et al. 2012, Catenazzi andKupferberg 2013), while conditions for post-metamorphic life stages in un-dammed tributaries remain unimpaired. Other stressors include invasive species that thrive under artificial low flow or drought conditions (Kupferberg 1997), chytridiomycosis (Adams et al. 2017b), habitat destruction (Davidson et al. 2002), and synergies among all these factors (Adams et al. 2017a). Artificial fluctuation in stream flow released from dams is a particularly potent cause of early life stage mortality because it can scour or strand a substantial portion of an annual cohort of egg masses and larvae (Lind et al. 1996, Kupferberg et al. 2012. Cold water from the hypolimnion of upstream reservoirs released in the summer shifts periphyton communities toward taxa inedible by tadpoles (Furey et al. 2014), reinforcing the effects of low temperatures on reduced tadpole growth and higher mortality Kupferberg 2013, Wheeler et al. 2015). Despite such mechanistic knowledge of the causes of this species' decline, critical information gaps about survival of R. boylii remain, in part due to dispersal of young of the year away from the mainstem lek sites into complex fluvial networks. Although monitoring studies have tracked trends in adult R. boylii abundance using egg-mass counts (Kupferberg et al. 2012), and manipulative experiments have provided insight into mortality of pre-metamorphic life stages (Kupferberg et al. 2011, Catenazzi andKupferberg 2013), capturemark-recapture (CMR) data suitable for estimating post-metamorphic survival are limited, and survival of pre-reproductive life stages remains largely unknown. Quantifying these vital rates in systems without anthropogenic impairment is a necessary step to project the viability of populations when stressors are added and to assess the need to protect R. boylii under the federal Endangered Species Act (United States Fish and Wildlife Service 2015).
In this study, we use an IPM to model the demography of robust populations of R. boylii inhabiting free-flowing streams at two northern California sites where the aforementioned anthropogenic stressors are largely absent (Hurdygurdy Creek in the Smith River watershed and Fox Creek and its confluence with the South Fork Eel River). Starting with that reference model, we then conduct scenario analysis to project the viability of populations in rivers with dams or in un-dammed rivers where increasing flow variability associated with climate change and disease presents future threats. Our IPM integrates long-term monitoring of egg-mass counts (7 yr at one site and 27 yr at the other), CMR data from adult frogs, fecundity, and survival of eggs to hatching to estimate recruitment to the subadult stage and survival of adult and subadult frogs, and to quantify trends in adult female abundance for these two populations. Importantly, the integration of multiple data sources, and the linkage between population trends and demographic rates, produces precise and realistic estimates of survival, including survival from egg to subadult, which is not directly observed. Using these demographic vital rates to parameterize a stochastic population model allowed for an examination of how different conservation and water management scenarios could affect population viability. By simulating the effects of increasing survival to adulthood (as would occur in a zoo rearing and reintroduction effort to bolster small populations) or decreasing recruitment from the egg stage (as when ill-timed releases from dams occur), we quantify anthropogenic influences on the relative probability of persistence vs. extirpation for R. boylii populations.

Data sources
We compiled data on the survival of each life stage of R. boylii and the fecundity of adult females from published literature and unpublished data (Appendix S1: Table S1). The number of eggs in individual egg masses (fecundity) and embryo survival rates (proportion of eggs that hatched) were counted from the South Fork Eel River (hereinafter "SF Eel River") in 1992 (n = 3) and 2007 (n = 20; Kupferberg et al. 2009).
We collected CMR data from adult and subadult R. boylii at Fox Creek, a tributary of the SF Eel River, from 2014 to 2019 (Fig. 1). We sampled for one or two primary periods each summer, between June and early September. Each primary period consisted of 2-6 secondary periods of sampling. We sampled frogs using visual encounter surveys both during the day (between 09:00 and 18:00 hours) and at night (between 20:00 and 01:00 hours). We marked frogs with passive integrated transponder (PIT) tags and determined sex based on the presence or absence of nuptial pads.
We marked and counted egg masses along an approximately 5.2-km reach of the SF Eel River each spring/summer from 1993 to 2019 using the census protocol described by Kupferberg et al. (2012). We performed two to six egg-mass surveys per year (average = 3.8 surveys) between the beginning of April and mid-June; egg masses were marked to prevent double counting. For the Fox Creek IPM, we included egg masses counted at the lek sites immediately upstream and downstream of the confluence of Fox Creek (where CMR surveys took place) and the SF Eel River (reach length = 507 m). We documented eggv www.esajournals.org mass mortality from stranding and scouring on the SF Eel River from 1992 to 2013 and 2017 to estimate the relationship between egg-mass survival and post-oviposition stream flow (following Kupferberg et al. 2012

Integrated population model
The IPM links the survival rate of each life stage and the fecundity of adult females to create a demographic model for R. boylii (Fig. 2). The survival rate of adult frogs is informed by both the CMR data and the annual egg-mass counts (Fig. 2). Egg-mass counts are also used to estimate the number of adult females present in the population each year. The inclusion of egg-mass counts as an estimate of adult female abundance constrains all vital rate parameters to values that match the temporal trend in the abundance of adult frogs in the breeding population. Likewise, the estimate of adult survival from CMR data constrains estimates of trends in adult female abundance.
The IPM is a female-based population model with three stages: egg (age 0), subadult (age 2), and adult (ages 3+). The number of eggs produced in year t is the product of the number of adult females, N ad [t], and the average number of eggs per egg mass, f, divided by 2 (assuming a 1:1 sex ratio of eggs; Eq. 1). The survival rate from the egg stage to subadult stage (ϕ rec ) is the product of the survival rate of egg masses in the face of stranding and scouring (ϕ em ), the  (Twitty et al. 1967) and observed near residential roads up to 331 m from a natal stream (Cook et al. 2012). For scale, the size of the adult female, the bottom frog in (D), is 55 mm (snout-to-vent length), and the size of the juvenile in (F) is 26 mm.
v www.esajournals.org survival rate of embryos hatching to become tadpoles (ϕ h ), and the survival rate from the tadpole stage until the 2-yr-old subadult stage (ϕ t . s ; Eq. 2). Individuals take two years to grow from the egg stage to the subadult stage (Eq. 3) and reach the adult stage (sexual maturity) at 3 yr of age. Empirical data are available for each survival parameter except ϕ t . s , which is estimated by leveraging the relationships between the other parameters in the IPM (Fig. 2). The number of female recruits that survive to two years old is the product of the recruit survival rate (ϕ rec ) and the expected number of eggs produced two years earlier, N egg [t − 2]. The survival rate of subadult frogs from 2 yr old to 3 yr old is assumed to be equal to the annual survival rate of adult frogs (ϕ a ). The number of adult frogs at time t, N ad [t], is the sum of adults present at time t − 1 that survive to time t and the number of subadults present at time t − 1 that survive to time t (Eq. 4).
We analyzed the CMR data for adult and subadult R. boylii using a Bayesian implementation of a Cormack-Jolly-Seber (CJS) model (Cormack 1964, Jolly 1965, Seber 1965, Kéry and Schaub 2012). The CJS model conditions on first capture and estimates the recapture rate (p) of marked individuals and their survival rate (ϕ a ) between primary sampling periods. We used a binomial  (IPM). Solid rectangles represent data sources, and circles represent parameters in the IPM. Dashed rectangles represent the two sub-models that make up the IPM. Parameters in the IPM are as follows: ϕ a , adult survival; ϕ t.s , tadpole to subadult survival; ϕ h , hatching probability (embryo survival); ϕ em , egg-mass survival; N ad , number of adult female frogs; N sub , number of subadult frogs; N egg , number of eggs; and N tad , number of tadpoles. Survival from egg to subadult, ϕ rec , is the product of ϕ em , ϕ h , and ϕ t.s v www.esajournals.org indicator of captured (1) or not (0) for each year of sampling; if we sampled frogs during multiple primary periods within a year, we collapsed data into a single binomial response. Because we cannot separate true survival from emigration, ϕ a is an estimate of apparent survival of adult and subadult frogs. We tested for an effect of frog sex on p and ϕ a in each CMR dataset by fitting sexdependent models and calculating the difference between females and males for each parameter.
We used a state-space model to analyze trends in egg-mass counts (Kéry and Schaub 2012). Female R. boylii, like many ranid frogs, lay one clutch per breeding season, so the number of egg masses is effectively a count of the number of adult females present and breeding in a given year (Zweifel 1955). The egg-mass count results from an observation process that samples the true underlying number of breeding adult females. The observation process is modeled using a Poisson distribution, and both over-and under-counting are possible, although overcounting is unlikely given our method of marking egg masses. The IPM does not explicitly include immigration or emigration, and we assume that all adult females breed each year. These assumptions can be violated for amphibian populations (Earl 2019) but were necessary because we had no data to estimate these parameters. We calculated the long-term log stochastic growth rate (log λ s ) for each population using the estimated abundance of adult females each year (Caswell 2001). To assess how annual population growth rates were related to time-varying vital rates, we calculated the correlation between λ[t] and adult survival that year (ϕ a [t]) and egg-mass survival three years prior (ϕ em [t − 3]), following the methodology of Kéry and Schaub (2012).
We summarized the fit of each sub-model with a Bayesian P value, where P equals the proportion of iterations in which the fit is better for replicate data generated by the model than the observed data. In general, Bayesian P values near 0.5 indicate good fit of the model to data, and values far from 0.5 indicate the model fits the replicate data better (near 1) or much worse (near 0) than the observed data (see Appendix S1; Kéry and Schaub 2012).
We used uninformative Beta(1,1) priors for mean ϕ a and mean p. We allowed adult survival and recapture probability to vary among years using a logit-normal random effect (Table 1). Survival from tadpole to subadult stage (ϕ t.s ) was constant because we lacked data to estimate annual variation in this parameter. We used an uninformative Beta(1,1) prior for ϕ t.s . The prior for hatching survival (ϕ h , the proportion of eggs that hatch into tadpoles) in the Hurdygurdy Creek IPM was a Beta(3.8, 0.75) distribution based on empirical data from the SF Eel River study reach (Table 1).
We modeled egg-mass survival at the SF Eel River as a function of the natural log of the ratio of maximum to minimum post-oviposition stream flow (hereinafter "flow magnitude index"; Kupferberg et al. 2012) using a binomial regression with the observed number of egg masses that survived stranding and scouring as the response. For Hurdygurdy Creek, egg-mass survival was modeled using an informative Beta(6,1) prior based on observed egg-mass survival rates at the SF Eel River study reach (Table 1).
We ran IPMs in JAGS version 4.3.0 (Plummer 2003) accessed through R version 3.6.3 (R Core Team 2020) using the runjags package (Denwood 2016). We ran IPMs on five independent chains of 2,000,000 iterations each after a burn-in of 100,000 iterations, with chains thinned by a factor of 100 to produce a final posterior of 100,000 samples. We inspected trace plots and calculated the Gelman-Rubin diagnostic,R (Brooks and Gelman 1998), to assess mixing of chains and convergence. All parameters hadR < 1.01 and exhibited no signs of failure to converge. For parameter estimates, we report the median of the posterior distribution followed by the 95% highest posterior density interval (HPDI) in parentheses. R and JAGS codes to replicate the IPM are available on ScienceBase (Rose et al. 2020).

Population viability analysis
We used vital rate estimates from the IPM to parameterize a stochastic population model for a population viability analysis (PVA). The goal of the PVA was to compare the sensitivity of population viability to changes in the survival rate of early life stages (eggs and tadpoles) and survival of adult and subadult frogs, not to estimate actual probabilities of focal populations going extinct. We projected the future trend of populations under different scenarios over a time period of 20 yr, similar to the study period for the SF Eel v www.esajournals.org River population. Shorter time horizons are preferable for PVA projections because extinction probabilities over longer time periods (e.g., 50 or 100 yr) are sensitive to initial population growth rates and environmental conditions over longer time horizons are highly uncertain (Fieberg and Ellner 2000). We incorporated both temporal variance in vital rates (i.e., process variance, Morris and Doak 2002) and parametric uncertainty (i.e., sampling variance or variance in estimates of survival rates from the IPM) to fully account for uncertainty in projections of population viability under different management scenarios (McGowan et al. 2011). We ran all PVAs for 10,000 simulations. For each simulation, we drew survival and fecundity parameters from the posterior distribution of the Fox Creek IPM to reflect uncertainty in vital rates. We then used the value drawn for that simulation for every year of the projection interval. We drew the adult survival value for each year (ϕ a [t]) from a normal distribution with a mean equal to the estimate from the posterior draw and a standard deviation based on a random draw from the posterior for the standard deviation of the temporal random effect on adult survival (σ ϕ a ).
We evaluated scenarios representing unregulated and regulated (i.e., dammed) streams. For the unregulated and regulated stream scenarios, the flow magnitude index was drawn randomly each year from a gamma distribution based on empirical stream discharge (m 3 /s) data from the SF Eel River (1990-2019) and Alameda Creek below Calaveras Reservoir (1996-2019), respectively (Appendix S1: Fig. S1; Kupferberg et al. 2012). We based egg-mass survival each year on the estimated relationship between ϕ em and flow magnitude (Fig. 3); heavy flows after oviposition (e.g., from late spring rains) can reduce survival of egg masses by scouring them from their attachment sites (Lind et al. 1996, Kupferberg et al. 2012. Pulsed flows released from upstream dams during Notes: Parameters for temporal variation in survival (σ ϕ a ) and recapture probability (σ p ) are presented on the logit scale. Lower and upper limits from the posterior distribution represent the highest posterior density interval. SD is standard deviation.
v www.esajournals.org the summer reduce tadpole survival by reducing algal resources, slowing growth, and washing tadpoles downstream (Kupferberg et al. 2009, 2011, Catenazzi and Kupferberg 2013. For the regulated stream scenario, the number of summer pulsed flows from upstream dams was drawn each year from a Poisson (rate = 1) distribution. The survival rate of tadpoles during each pulsed flow was 0.61 (Kupferberg et al. 2009). No summer pulsed flows occur in the v www.esajournals.org unregulated scenario. Adult and subadult survival rates are equal for the regulated and unregulated stream scenarios.
For both the regulated and unregulated scenarios, we simulated the effect of altering (1) survival from egg to subadult stage (ϕ rec ) and (2) survival of subadult and adult life stages (ϕ a ) by a set percentage (from +50% to −50% in 10% increments) on the probability of quasi-extirpation. We chose these simulated reductions in survival to reflect empirical estimates of the effect of known stressors for each life stage. Recent evidence indicates that declines in some populations of R. boylii are related to drought-associated outbreaks of chytridiomycosis caused by Batrachochytrium dendrobatidis (hereinafter Bd; Adams et al. 2017a). During outbreak events, mortality rates of post-metamorphic frogs in the closely related Rana muscosa and Rana sierrae can be very high (Rachowicz et al. 2006, Briggs et al. 2010). In anuran populations in the western United States where Bd is endemic, the effects of chytridiomycosis on post-metamorphic survival ranged from a 6% reduction to as high as 55% (Russell et al. 2019). We designed the PVA scenarios with reduced ϕ a to represent populations where Bd is endemic and causes a persistently lower adult survival rate. We based the PVA scenarios with reduced ϕ rec on observed fluctuations in eggmass survival as a function of post-oviposition flow magnitude at the Fox Creek population.
We simulated the effect of management interventions targeted toward early, pre-metamorphic life stages and subadult and adult R. boylii by increasing the survival rate of these life stages at equal levels to the reductions described above. Interventions that could increase the survival of early life stages have been demonstrated from observation and experimental studies. For example, R. boylii tadpoles reared within in-stream enclosures exhibited survival rates that were likely much higher than tadpoles in natural conditions (Kupferberg et al. 2009), and limiting extreme flows after oviposition can increase eggmass survival (Kupferberg et al. 2012). Interventions that increase the survival of adult and subadult frogs are not as well established, but could include active removal of invasive predators and competitors and restoration of natural flow regimes, which keep non-native predators in check (Moyle 1973, Hayes et al. 2016).
For each of the above scenarios, we ran simulations with the initial abundance of adult females (N 0,ad ) set to 25, 50, 100, and 250 individuals/km based on the range of egg-mass densities (no./ km) observed at SF Eel River, Hurdygurdy Creek, and other rivers in central and northern California (Catenazzi and Kupferberg 2017). We set a quasi-extirpation threshold of 10 adult females per river km; if the abundance of adult females dropped below this threshold, we considered the simulated population functionally extinct. We calculated the probability of quasi -extirpation as the proportion of simulations that dropped below the quasi-extirpation threshold during the 20-yr projection interval. We ran the PVA using custom R code that is available on ScienceBase (Rose et al. 2020).

We made 571 captures of 340 frogs at Fox
Creek from 2014 to 2019, including 145 females, 182 males, and 13 individuals of unknown sex. Of these 340 frogs, 270 frogs were captured in only one year, 53 were captured in two years, 15 were captured in three years, and two were captured in four years. The CJS (P = 0.36) and eggmass count sub-models (P = 0.56) of the Fox Creek IPM showed no evidence of a lack of fit to the observed data. There was no evidence for differences in survival (ϕ f − ϕ m = 0.014, 95% HPDI = −0.14 to 0.17) and recapture rates (p f − p m = −0.001, 95% HPDI = −0.24 to 0.23) between female and male R. boylii at Fox Creek. Therefore, we focused our inference on an IPM where we pooled captures of both sexes to estimate adult survival (Table 1). The estimated mean annual adult survival rate at Fox Creek was 0.46 (95% HPDI = 0.37-0.56). When transformed back to the probability scale, the standard deviation of annual variation in adult survival was 0.08 (0.01-0.14) or approximately 17% of the average annual survival rate of adults at Fox Creek.
The median estimate of survival from the egg to subadult stage (ϕ rec ) for a year with average post-oviposition flow was 0.0016 (95% HPDI = 0.0007-0.0026; Appendix S1: Fig. S2). There was a negative relationship between egg-mass survival and the post-oviposition flow magnitude index (β flow = −0.98, −1.16 to −0.81); egg-mass v www.esajournals.org survival was lower in years with a greater difference between maximum and minimum postoviposition flows (Fig. 3). Given the relationship between egg-mass survival and the magnitude of post-oviposition flows, the proportion of eggs that survive to reach the subadult stage (ϕ rec ) can vary over time. For the lowest post-oviposition flow magnitude measured at the SF Eel River, survival from egg to the subadult stage would be 0.0018 (0.0011-0.0027), whereas for the highest measured post-oviposition flow magnitude, survival to subadult would only be 0.0006 (0.0004-0.0010). With an average of 1526 eggs per egg mass, assuming a 1:1 sex ratio and no immigration or emigration of recruits, only 1.3 (0.8-1.8) females per egg mass are expected to reach two years of age on average; recruitment of two-year-old frogs could vary from 0.5 (0.3-0.8) females per egg mass in a bad flow year to 1.4 (0.9-2.0) females in a good flow year.
The abundance of adult female R. boylii at Fox Creek fluctuated from 1993 to 2019 (Fig. 3  A) and closely mirrored fluctuations of eggmass counts in the larger SF Eel River study reach (Appendix S1: Fig. S3). The average population growth rate, log λ s , indicated a stable population at Fox Creek over this time period (median = 0.000, 95% HPDI = −0.012 to 0.011). Despite the apparent stability in the long term, the abundance of adult females fluctuated substantially from year to year. The annual population growth rate, λ[t], at Fox Creek was positively correlated with ϕ a [t] (r = 0.36; 95% HPDI = −0.05 to 0.68) and with egg-mass survival three years earlier, ϕ em [t − 3] (r = 0.32; 0.14-0.49). Although the median correlation between λ[t] and ϕ em [t − 3] was lower, this relationship was less uncertain than the relationship between λ[t] and ϕ a [t] (Fig. 4). The effect of low ϕ em on λ[t] is apparent for two observed declines in N ad . High post-oviposition flow magnitude resulted in low ϕ em in 2005, which was followed by a decline in N ad from 2008 to 2009 with estimated λ = 0.74 (0.5-1.02). Similarly, low ϕ em in 2009 was followed by a decline in N ad from 2012 to 2013, with estimated λ = 0.72 (0.51-0.97; Fig. 3). Stream discharge at the SF Eel River study reach during the post-oviposition period appeared more variable from 1990 to 2019 than during the mid-20th century, with years of both higher and lower post-oviposition flow magnitude recorded from 1990 to 2019 compared to data from 1946 to 1970 (Appendix S1: Fig. S4). for each parameter. r is the median correlation coefficient between λ t and the survival parameters, followed by the 95% HPDI in parentheses. P(r > 0) is the posterior probability that the correlation coefficient is greater than zero.
We made 1052 captures of 319 frogs at Hurdygurdy Creek from 2002 to 2008, including 121 females and 198 males. Of these 319 frogs, 258 were captured in only one year, 47 in two years, 12 in three years, and 2 in four years. The CJS (P = 0.19) and egg-mass count sub-models (P = 0.57) of the Hurdygurdy Creek IPM showed no evidence of a lack of fit to the observed data. Sex differences were apparent in the survival and recapture rates of adult frogs from Hurdygurdy Creek; females had higher survival and lower recapture rates than males ( Table 1). The estimated survival rate of adult females at Hurdygurdy Creek (0.48, 0.32-0.67) was similar to that of those at Fox Creek. The standard deviation of annual variation in female adult survival was 0.10 (<0.01-0.23) or approximately 21% of the average annual survival rate of adults at Hurdygurdy Creek. The estimated survival rate from egg to subadult at Hurdygurdy Creek was also similar to that at Fox Creek (0.0019, 0.0008-0.0035), although the median estimate for ϕ rec was slightly higher at Hurdygurdy Creek (Appendix S1: Fig. S2). The abundance of adult female R. boylii at Hurdygurdy Creek was stable from 2002 to 2008 ( Fig. 3B; median log λ s = 0.019, −0.034 to 0.072).
Of the scenarios evaluated in the PVA, the biggest influence on quasi-extirpation probability was whether the simulation represented a regulated or unregulated stream. The probability of quasi-extirpation over a 20-yr period was much higher for simulated populations in regulated streams, which experienced pulsed summer flows that reduced tadpole survival, compared to unregulated streams without pulsed flows (Fig. 5). Within a given flow scenario, reducing ϕ a produced a rapid increase in the probability of quasi-extirpation. In comparison, reducing ϕ rec had a weaker effect on quasi-extirpation, particularly for the unregulated scenario. Increasing the survival rate of adults and subadults resulted in much lower probabilities of quasi-extirpation compared to equivalent increases in the survival rate from egg to subadult.

DISCUSSION
In contrast to the pernicious loss of amphibian biodiversity globally (Scheele et al. 2019), the population growth rate for both focal populations of R. boylii was stable when assessing the study period as a whole, despite inter-annual fluctuations in abundance and recruitment, as well as multi-year periods of apparent decline. Measuring stability of populations requires a long time series relative to a species' generation time (Connell and Sousa 1983); indeed, had we sampled the Fox Creek population for only a fiveor 10-yr period we might have concluded it was declining (e.g., from 2011 to 2015) or growing (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011). The dynamic equilibrium of the populations is notable in light of the survival rate of adult female R. boylii (0.46-0.48), which is lower than estimates of adult survival for three declining Rana species in California: 0.73 in R. sierrae (Fellers et al. 2013), 0.80 in R. muscosa, and 0.58 in Rana draytonii (Russell et al. 2019), but similar to reintroduced Rana pretiosa in Oregon (0.40; Duarte et al. 2017). A potential explanation for the apparent stability might be that episodes of recruitment into the adult population in some years compensate for periods of low recruitment and low adult survival. This interpretation is supported by: (1) the estimated threefold variation in survival from egg to subadult, driven by fluctuations in the magnitude of stream flow following oviposition and (2) the positive correlation between annual population growth rates for and the estimated survival rate of egg masses laid three years earlier. An alternative explanation for the observed stability is that immigration of adult females from outside of our sampling area offset losses due to death and emigration from the sampling area. Because the model estimated apparent survival, true survival could be higher if females emigrated from our study area and were unavailable for capture. It is unlikely that immigration and emigration were the main drivers of the observed dynamics at Fox Creek, because we observed similar fluctuations in egg-mass counts for a much longer reach of the SF Eel River over the same time period. The congruence in vital rate estimates from the two models suggests similar demographic processes govern population growth for these two populations inhabiting unregulated streams in northwestern California. Population densities of R. boylii in northwestern California are robust relative to those in Oregon (Olson and Davis 2009), the Sierra Nevada, and Coast Ranges of central and southern California (California v www.esajournals.org Department of Fish and Wildlife 2019). Therefore, demographic patterns observed in this study might not be found in other parts of the species' range where hydrologic regimes differ (e.g., snowmelt driven in the Sierra or xeric conditions in the south).
The importance of recruitment to the population dynamics and persistence of R. boylii heightens the need for precise estimates of the survival rate of early life stages. We did not have data that directly informed the survival rate from the egg stage to the subadult stage. However, based upon the fecundity of R. boylii, adult survival rates (estimated from CMR data), and trends in the adult female population (from egg-mass counts), we were able to obtain estimates of the proportion of eggs that survive to two years of age. Because we had concurrent egg-mass counts as an index of adult female abundance and CMR data to estimate adult survival, the IPM could constrain vital rate estimates to realistic values that matched the observed relative stability in abundance averaged over the time series for each population. Leveraging multiple data sources that provide information about the same vital rate yields more precise and realistic estimates and is one of the great advantages of IPMs (Schaub and Abadi 2011), although estimates of additional parameters can be sensitive to model assumptions (Riecke et al. 2019).
The average survival of early life stages is not the ultimate driver of the population dynamics of R. boylii; however, mortality of eggs and tadpoles varies greatly from year to year due to hydrologic conditions (Kupferberg et al. 2012) and is likely to increase in variability as climate changes (Grantham et al. 2018). In some years, survival from egg to subadult is higher than our average estimate, which could result in an episode of high recruitment, whereas in other years, survival to subadult is lower, potentially resulting in little recruitment. Such episodic recruitment events are common among streambreeding fish, where environmental variability can be strongly linked to larval survival and productivity (Sakaris andIrwin 2010, Lobón-Cerviá 2014). Therefore, predicting the viability of any single population of R. boylii is fraught, because its fate will be inextricably linked to the joint effect of temperature and precipitation on hydrological conditions in the future. Even these robust populations may not be immune to the effects of climate change. Our retrospective analysis of post-oviposition flow fluctuation in the SF Eel River revealed greater variability from 1990 to 2019 than during the mid-20th century, and California streamflow predictions for the mid-21st century  forecast that trend will continue. Projections informed by global climate models indicate that May (typically when most R. boylii breeding occurs) will experience the largest decline in monthly flow relative to a baseline  in coastal northern California, and change in May flows will be the greatest in the xeric central and southern regions (Grantham et al. 2018), where R. boylii is largely extirpated (Adams et al. 2017b).
Population models for organisms as disparate as butterflies and lizards show that modeling the response of each life stage is important for evaluating the outcome of future environmental scenarios (Radchuk et al. 2013, Levy et al. 2015. Creating detailed population models can be challenging for species with cryptic life stages, however. Rather than forecasting the trajectories of our study populations, the best use of PVA is to evaluate the relative risk of extirpation under different hypothetical management scenarios (Lindenmayer and Possingham 1996, Akçakaya and Sjögren-Gulve 2000). Our approach of combining multiple data sources in an IPM to estimate a vital rate that was not directly observable and using vital rate estimates to evaluate relative extinction risk under different management scenarios could be applied to many species with difficult-to-study life stages. Two key considerations for applying our methods are (1) the life cycle model of the focal species must reflect the data available to inform each vital rate parameter (i.e., one cannot estimate multiple vital rates that are each lacking empirical data) and (2) simulated scenarios should reflect natural or humaninduced perturbations of the estimated vital rates to explore which interventions and vital rates influence population viability the most.
For R. boylii, PVA scenarios could include efforts to mimic ecologically relevant aspects of a river's natural flow regime, to augment survival of eggs and tadpoles via in situ rearing, or other types of conservation translocation. Each of these management methods is increasingly employed to address biodiversity loss and reverse population extirpation (Poff et al. 1997, Harding et al. 2016, Bubac et al. 2019. Population viability analysis results demonstrated that reduced survival rates of early life stages (eggs, tadpoles), designed to simulate the effect of ill-timed releases of water in regulated rivers, could have a large effect on the probability of extirpation. This result reinforces the finding from a regional study of population trends in R. boylii, which concluded that altered flow regimes lead to declines in abundance and extirpation downstream of large dams (Kupferberg et al. 2012), as well as v www.esajournals.org the global finding that reversing the degradation of riverine biota depends on imitating natural flow regimes (Palmer and Ruhi 2019).
In our simulation of disease effects, decreasing the survival of subadult and adult frogs by 10-20% led to a much-increased probability of extirpation in modeled populations, even those in unregulated streams. A recent die-off of post-metamorphic R. boylii at a stream in central California, likely caused by chytridiomycosis, raises concerns that conditions favorable to Bd outbreaks, such as drought that crowds animals together or the presence of Bd-tolerant bullfrogs, could cause population-level declines (Adams et al. 2017a). Bd epizootics in other Rana species can lead to high mortality in post-metamorphic life stages, resulting in rapid extirpation (Briggs et al. 2010). In populations of western Rana species with enzootic Bd, adult female survival was reduced by 6-12%, but for western toads (Anaxyrus boreas), survival was reduced by up to 55% (Russell et al. 2019). Based on our PVA, even small reductions in adult survival could be enough to drive population declines and increase the risk of extirpation. We detected Bd on frogs from Fox Creek in 2018 and 2019 (USGS, unpublished data), but the potential effect of the pathogen on the survival of post-metamorphic frogs in this population is not known.
Although altering the survival of adult and subadult frogs had a greater influence on population viability than an equivalent change to the survival from egg to subadult, pre-metamorphic life stages are likely more amenable to management intervention. Head-starting of early life stages is a popular strategy to rescue amphibian species on the verge of extinction (Griffiths and Pavajeau 2008). In addition to the aforementioned effects of flow on survival of egg masses, it has been shown that R. boylii tadpoles can survive to metamorphosis at a higher rate when reared in an in situ enclosure and fed high-quality food (Catenazzi and Kupferberg 2013). The difference in egg-mass survival between benign and highly volatile post-oviposition flow regimes can be >50%, illustrating the benefits of avoiding high post-oviposition flows in regulated streams. In contrast, it is unlikely that management actions could increase adult survival even 1.5-1.75 times greater than its average value of 0.46-0.48. Compared to the survival of early life stages, the natural variability of adult survival was low, even though the 2014-2019 CMR period at Fox Creek included both the most severe multi-year drought of the last 1200 yr (Griffin and Anchukaitis 2014) and extreme winter flooding. This stability in adult survival could mean population growth is driven more by the variation in early life stage survival (Gaillard et al. 2000). It is valuable to consider not only the proportional effect of management actions, but also the range of vital rates that are possible in natural populations (Mills et al. 1999).
There are gaps in the data we have available to estimate demographic rates for R. boylii, and future studies could improve on this work by using novel methods to collect empirical data that could fill these gaps. For example, in the PVA we assumed that the survival rate of postmetamorphic life stages was the same in regulated and unregulated streams. It would be valuable to collect concurrent egg-mass counts, adult CMR data, and fecundity data from a population of R. boylii inhabiting a regulated stream to see whether the demography of post-metamorphic life stages differs from unregulated streams. Also, we assumed that all females breed each year and therefore that the number of egg masses could represent the abundance of adult females in the population. Female ranid frogs might not always breed annually; for example, female R. pretiosa are known to skip breeding in some years (Turner 1960). If adult female R. boylii skip years of breeding, then our estimated vital rates could be biased, because only changes in survival rates were able to account for increased numbers of egg masses in certain years in our model (e.g., 2006, 2011, 2017). It would be valuable to link the CMR data for adult females to the egg-mass count data by tracking the identities of which females return to the oviposition sites to breed each year. Because we conducted our CMR sampling at Fox Creek in the non-breeding, summer habitat, our captures of adult females are not affected by which individuals breed or not in a given year. Installing a monitor array that tracks individuals going to and from the breeding site would provide useful data to estimate the proportion of adult females that breed each year and potentially to link that vital rate with v www.esajournals.org environmental factors that drive females' reproductive decisions.

SUMMARY AND CONCLUSIONS
Collecting concurrent egg-mass count and CMR data and analyzing these data in an IPM enabled us to obtain precise estimates of adult survival for R. boylii and to estimate the low survival rate of unobserved early life stages. Researchers studying species with complex life histories can employ IPMs to gain inference into the vital rates of life stages that are challenging to study directly. A relatively small increase in effort expended, to collect data that complement an existing study, can have great benefits for modeling population dynamics. Still, our understanding of the population dynamics of R. boylii would benefit from more data on the survival of unobserved life stages, such as recently metamorphosed frogs, along with potential drivers of temporal variation in adult survival. Given the increased probability of extirpation of R. boylii in regulated streams due to low survival of premetamorphic life stages, threats that could increase the mortality of adult frogs, such as drought-induced outbreaks of chytridiomycosis, could greatly threaten extant populations. Active conservation efforts, such as habitat enhancement, flow management, and head-starting of early life stages, may be needed for the continued persistence of populations.