Drivers of variation in the population dynamics of bighorn sheep

,


INTRODUCTION
The successful management and conservation of wild populations, particularly those at risk, requires information on sources of variation in vital rates and the contribution of those vital rates to variation in demography and population trajectories (Nichols and Williams 2006). Where management resources are limited, such understanding allows the optimal allocation of resources to restoration efforts (Johnson et al. 2010, Mills 2012. Though the importance of understanding the relative influence of vital rates for shaping population trajectories has long been recognized, such assessments are still comparatively rare due the challenges of estimating the vital rates and population structure required (Heppell et al. 1996, Johnson et al. 2010, Koons et al. 2017. Consequently, information is frequently borrowed from species and systems for which there is sufficient information. For ungulates, the often-referenced expectation for population dynamics is that population growth rates are most sensitive to variation in adult survival; however, variation in adult survival is small enough to contribute little to actual variation in population growth rates ). In contrast, substantial variation in offspring recruitment is the primary driver of variation in population growth rates, even if growth rates are theoretically less sensitive to changes in this vital rate than changes in adult survival. This paradigm was developed for, and intended to be applied to, specific kinds of ungulate populations, that is, from temperate zones, non-harvested and with stable or increasing population growth rates (Coulson et al. 2005, Nilsen et al. 2009). Evidence suggests that the dynamics of populations that are small, in decline, subject to stochastic and substantial variation in vital rates, or in non-temperate ecosystems defy these expectations (Owen-Smith and Mason 2005, Nilsen et al. 2009, Johnson et al. 2010, Lee et al. 2016, and the evolving understanding of ungulate population dynamics has broadened the perspective to demonstrate how elasticities and variances of vital rates integrate to shape trajectories (Hilde et al. 2020).
The population dynamics of bighorn sheep (Ovis canadensis) may or may not operate according to this paradigm. Populations of bighorn sheep suffered significant declines near the turn of the 20th century in response to a series of pressures including disease, competition with livestock, and over-harvest (Buechner 1960, Berger 1990, and restoration efforts have demonstrated mixed success (Singer et al. 2000, Picton and Lonner 2008, Hedrick 2014. Major areas of research required to provide the information for bighorn sheep restoration have been identified for decades (Buechner 1960), and a substantial body of work has developed characterizing the important vital rates for bighorn sheep. Although adult survival is generally high with limited among-year variation, stochastic events such as predation by specialist predators and respiratory disease epizootics, as well as hunter harvest, can induce variation in this key vital rate and generate very different population dynamics and trajectories compared to populations not subject to such pressure (Festa-Bianchet 1989, 2006, Ross et al. 1997, Cassirer and Sinclair 2007, Brewer et al. 2014, Manlove et al. 2016, Parr et al. 2018). Moreover, a recent comprehensive estimation of sources of variation in ewe survival suggests that, apart from the influence of stochastic events and harvest, ewe survival varies among years in response to environmental drivers (Proffitt et al. 2021). Lamb survival demonstrates substantial among-year variation in response to predation, disease, and environmental variation, which is consistent with the empirical and theoretical expectation for high variation in this rate for ungulates in general ) and bighorn sheep in particular (Douglas and Leslie 1986, Hass 1989, Portier et al. 1998, Smith et al. 2014. However, unlike the ephemeral decline of adult survival in response to respiratory disease epizootics, lamb survival can be depressed for years following a disease event Sinclair 2007, Cassirer et al. 2013). Pregnancy rates are typically high (Festa-Bianchet 1988, Singer et al. 2000. However, pregnancy is the result of a complex set of metabolic processes that integrate environmental variation and state processes (e.g., previous year's reproductive success) such that pregnancy rates can also vary among years (Parker et al. 2009).
Given the potential for wide variation in vital rates, it is unsurprising that there is no consensus on which rate is most important to the growth of bighorn sheep populations. Previous work has suggested that population trajectories are driven by recruitment (Bender andWeisenberger 2005, Manlove et al. 2016), adult survival (Rubin et al. 2002, Singer andSchoenecker 2004), or a combination of the two (Parr et al. 2018). The most comprehensive evaluation to date of the relative contributions from vital rates to population growth for bighorn sheep strongly suggests that these drivers can differ between populations over comparatively small spatial scale (Johnson et al. 2010), a conclusion supported by recent work assessing variation in population trajectories for populations of bighorn sheep in the v www.esajournals.org Rocky Mountain west (Donovan et al. 2020). Rather than being contradictory, these disparate conclusions point to the diversity of challenges facing bighorn sheep populations, where the aggregated pressure from predation, disease, small population size, and environmental drivers varies spatially and temporally (Buechner 1960, Brewer et al. 2014. A better understanding of the mechanisms through which population trajectories of bighorn sheep transition from increasing or stable to declining (or vice versa) is important for informed management and restoration efforts (Coulson et al. 2005, Johnson et al. 2010, Koons et al. 2016. Prior work on bighorn sheep and other ungulate species has implicated declines in adult survival as the primary driver of declining population growth rates, a pattern seen in declining populations in both temperate (Wittmer et al. 2005, Nilsen et al. 2009, Johnson et al. 2010) and tropical systems (Owen-Smith andMason 2005, Lee et al. 2016). In contrast, improved offspring recruitment is heavily associated with the recovery of populations of ungulates across a similar range of ecosystems such that understanding limiting factors of offspring recruitment is crucial (Beissinger and Peery 2007, Mitchell et al. 2009, Manlik et al. 2016. For bighorn sheep, there are comparatively few studies of how lamb recruitment varies in response to ecological drivers, and the extant work suggests a potentially complex interplay of effects from environmental, predation, and disease-related factors (Hass 1989, Portier et al. 1998, Manlove et al. 2016, Butler et al. 2018. Here, we use multiple long-term time series of survey data on 17 populations of bighorn sheep in Montana and Wyoming in combination with a recent comprehensive estimation of bighorn sheep pregnancy and survival rates (Proffitt et al. 2021) to estimate the parameters of a population model for bighorn sheep to address three objectives: (1) to characterize the variation in key vital rates (lamb survival and adult survival) and population growth rates, (2) to estimate the contributions of changes in these vital rates to observed changes in population trajectories, and (3) to identify potentially important sources of variation in lamb survival including environmental sources of variations (e.g., precipitation, primary production, and winter severity), predation, and disease.

Study area and populations
The study was conducted in Montana and Wyoming and included data for 17 bighorn sheep populations (Fig. 1). In all but two cases, populations were defined based on historic delineations of management units used by the management agencies in the two states thought to (generally) represent demographically closed populations. In two cases, we combined data from separate management units in Wyoming to create single demographically closed populations (Whiskey Mountain [Mtn]-West and Whiskey Mtn-East) that each were the combination of two underlying management units. Sixteen of the 17 populations were located in mountainous areas, and one population (Middle Missouri) was located in a prairie breaks landscape. Although the composition of the communities of pathogens that are associated with bighorn sheep respiratory disease are not known in these populations, the most comprehensive assessment to date suggests that all of these populations host Pasteurellaceae bacterial pathogens, and all but Paradise, Petty Creek, Targhee, and Middle Missouri host Mycoplasma ovipneumoniae (Butler et al. 2018). All 17 populations had a management history that included removals from the population due to translocations of animals to other herds or harvest. Every population had ram removals, with the number of rams taken in any given year ranging from 0 to 69 (median = 7 removed). However, only eight of the 17 populations had ewe removals; in contrast to ram removals, the number of ewes removed in any given population varied from year to year ranging from 0 to 79 (median = 0 removed). The management histories indicate that nine of the populations experienced all-age die-off events at least once during the study period.

Survey data
Survey data for the 17 populations were compiled based on agency records that contained information obtained using multiple methods, including fixed-wing and helicopter-and ground-based survey efforts to count animals and classify animals into age and sex classes. The timing of surveys differed among and within bighorn sheep populations (median = March, v www.esajournals.org range = December-May). Ideally, a survey was conducted annually, and during each survey, a sample of the total animal inventory was classified as lambs, adult males, and adult females. A count of all observed animals was recorded if the survey was thought to be representative of the population in the consistent area of core seasonal range. However, the nature of the data and the practical challenges associated with annual surveys resulted in a discontinuous time series for nearly every population such that some years were missing a count (but had classification data), some years were missing classification data (but had a count), some years were missing all data, and some years had both count and classification data (complete) (Fig. 2, see Appendix S1 for complete population observation histories). We left-truncated the time series for each population at the first year for which there was a representative count. v www.esajournals.org 4 July 2021 v Volume 12(7) v Article e03679 Fig. 2. Structure of the bighorn sheep survey data (a) and variation in covariate values used as potential sources of variation in lamb survival (b). The survey data are largely discontinuous time series, where data in any year can be missing: both total count data and classification data (missing all data), classification data, or total count data. Covariates used for lamb survival index growing conditions (primary production [NDVI] and precipitation [PREC]) during the spring (early) and summer (late) and winter severity (PREC winter ) (colors represent populations). NDVI, normalized difference vegetation index.

Modeling approach
Our modeling approach was to characterize the drivers of variation in the population dynamics of bighorn sheep using the time series of management count and classification data. Prior work has demonstrated that statespace approaches have desirable characteristics for modeling population dynamics (de Valpine and Hastings 2002, Buckland et al. 2004, Newman et al. 2006, Schaub and Abadi 2011. A state-space model in this context consists of two processes: a biological model that connects changes in the population structure and size through time via key vital rates such as survival and fecundity and an observation model that handles the stochastic and imperfect observation process (Buckland et al. 2004). One of the key features of such state-space models for population dynamics is that they separate the observation process from the biological process and, conditional on the model, improve precision on inference for population dynamics (de Valpine and Hastings 2002). Such hierarchical models naturally fit within the Bayesian paradigm of statistical inference, and the flexible nature of model specification in the popular and freely available software packages accommodate such models well (e.g., JAGS, STAN, OpenBUGS) and also have two key additional benefits. First, by specifying a biological process model that connects vital rates to size of age/sex classes through time, it is straightforward to derive a number of biologically relevant quantities such as population growth rates and sex/age class ratios that are the integrated result of variation in multiple vital rates. Second, it is easy to introduce informed priors for vital rates that are not directly estimated in the model.
Biological process model.-We developed a stage-based model for the population dynamics of bighorn sheep (Caswell 2001). We defined a biological year from June 1 in year t to May 31 in year t + 1 to account for the pre-birth pulse surveys. At the time of the surveys in the late spring, each population was comprised of lambs (<1 yr old), yearlings (>1 and <2 yr old), and adults (≥2 yr old). Recent work in this system demonstrated that pregnancy rates for yearling ewes can be substantially lower than that for older ewes, suggesting that it is important to treat yearlings as a separate class in a biological model to avoid bias in estimates of offspring survival (Proffitt et al. 2021). Therefore, our biological process model used six age-sex classes (lambs, yearlings, and adults by sex).
Bighorn sheep females typically have a single offspring, and fecundity in a pre-birth pulse model is therefore the product of pregnancy rates and lamb survival to the end of their first year. For each population, we separately modeled the number of lambs (P l ) produced from yearling ewes (ye) and adult ewes (ae) in biological year t as a function of the size of these two reproductive classes in the previous year and age class-specific pregnancy rates (τ ye and τ ae ), and a binomial process to incorporate demographic stochasticity (or, variation in fates at the population level due to chance, May 1973, Lande 1993, Kendall and Fox 2002: such that the total number of lambs produced at the start of biological year t is simply The survival of lambs from the beginning of biological year t to the end of biological year t (immediately prior to transitioning to the yearling class in year t + 1) was modeled using a similar binomial process based on lamb survival (S l ): The number of yearling ewes (N ye t ) and yearling rams (N yr t ) at the end of biological year t was modeled using a binomial process based on the number of lambs in year t − 1 (N l tÀ1 ), the assumption of an equal sex ratio of lambs (0.5), and sex-specific adult survival (ewes: S e , rams: S r ): Finally, the number of adult ewes (N ae t ) and rams (N ar t ) at the end of year t were modeled using a binomial process based on the number of yearling ewes and rams in year t − 1 (N ye tÀ1 , N yr tÀ1 ), the number of adult ewes and rams in year t − 1 (N e t , N r t ), and sex-specific adult survival: Our biological model relied on three simplifications of bighorn sheep vital rates: (1) The survival of male and female lambs are equal, (2) there was no difference in the sex-specific survival of yearlings and adults, and (3) there was no age-related variation in either survival or pregnancy rates (e.g., reproductive or actuarial senescence). Our biological model also assumed density independence of vital rates, a reasonable assumption for the generally smaller, recovering populations used here for whom densitydependent resource limitation was unlikely to be limiting.
Observation process model.-To connect the expected number of animals in each sex and age class generated by the biological process model to the observed data, we used a model that could accommodate survey data where only a fraction of the total number of animals counted were classified into the management age/sex categories (lamb, adult female, adult male) (Paterson et al. 2019). Frequently, populations were surveyed for counts, but only a fraction of the total number of animals counted were classified into the management age/sex categories (lamb, adult female, adult male). Therefore, we modeled the total count for each population in each year t (Count t ) using a Poisson process where the expected value was the total population size (the sum of all age and sex classes, where ζ t was an observation-level random effect used to account for overdispersion in the observation process. We then connected the total number of animals classified in each survey (Classified t ) to the age and sex classes using a multinomial distribution. During these management surveys for population trends, animals were classified only as lambs, ewes, and rams, for example, yearling ewes and adult ewes are not separately observable. However, inference on an unobserved age class is possible in our model due to the structure of the biological process model. We incorporated the yearling class into the classification process in year t by relating the number of classified lambs (Count l ), ewes (Count e ), and rams (Count r ) to their proportion in the population: This observation model connected the survey data to the underlying vital rates, a central assumption of which is that these populations were demographically closed. Where this assumption was violated in the presence of immigration/emigration, it would be expected to result in biased estimates of vital rates, that is, survival biased high in years with immigration. Where animals were removed by harvest or for translocations, it would bias survival rates low if removals were additive to other sources of mortality. Additionally, the presence of incomplete data (years with completely missing data, or years missing a total count or a classification) would be expected to inflate the variance of estimates of vital rates from our model and/or generate bias in sex-specific survival rates.
Accommodating inconsistent survey timing.-A common complication in our survey data is that surveys were not always conducted at the end of the biological year in late spring. Thus, if we ignored differences in survey timing and naively used the biological and process models above, we would have introduced bias into our estimates of vital rates by conflating the timing of surveys with the survival processes, that is, early surveys that missed late winter mortality would tend to cause year-specific survival rates to be inflated. We therefore adapted the above model to the timing of the surveys with a model for how survival changes within a year and used the timing of the observed data to help estimate the population size at the end of year t. Specifically, we introduced another latent state for the size of each age and sex class at the time of the surveys and connected it to the model for the size of each age and sex class at the end of year t. For example, the number of lambs in the population in year t at the time of the survey (N l t * ) was modeled as v www.esajournals.org and the number of lambs at the end of the biological year (N l t , the ideal timing) was based on this latent state where f(S l , d * ) was a function that connected the relationship between the timing of the observation (d * ) and the expected survival of lambs to that point, and g(S l , d * ) was a function that connected the difference in timing between the surveys and the end of the year to the expected number of lambs at the end of the year. We chose to use simple functions that factored the survival rate into monthly survival and then calculated survival to the time of observation as monthly survival raised to a power equivalent to the number of months since the beginning of year t. For example, if the survey was conducted in March, it was in month 10 (d * = 10) of biologi- cal year t then f S l , 10 À Á ¼ S l À Á 1 12 10 , and g S l , 10 . We used these same functions to model the latent states of all age/sex classes. Such survival functions are a considerable simplification of a reality in which the shape of the survival curve likely differs between age and sex classes and among years. However, we considered this approach to be a reasonable approximation in the face of sparse information for this species and a considerable improvement over naively using survey data without accounting for the timing.

Modeling goals
Our state-space model for the time series of counts is a hierarchical model with a large set of latent states required to accommodate the structure of the data. However, these models for populations are driven by a small number of vital rates: pregnancy rates of yearling and adult ewes (τ ye and τ ae ), the survival rate of lambs (S l ), and the survival rates of ewes (S e ) and rams (S r ). Our goals were to estimate these vital rates from the count data, to characterize their associations with population trajectories (a derived parameter), and to understand sources of variation in lamb survival. To meet these three goals, we used two versions of the above model: (1) a version wherein S l was estimated for each populationyear using a random effects structure (REmodel) and (2) a version that used logistic regression to assess potential sources of variation in lamb survival in each population using covariates as proxies for environmental variation (COVmodel). All other model structures were the same between models.
Common model structures. -We used informed priors for adult survival and yearling and adult ewe pregnancy rates, constructed from a recent comprehensive assessment of pregnancy and survival rates in bighorn sheep in this same area and for many of these populations (Proffitt et al. 2021). We took the empirical distribution of the estimated probabilities for yearling pregnancy and adult pregnancy and fit a beta distribution to each to construct the priors for the populationspecific probabilities of pregnancy in the statespace model (τ ye p ∼ Betað29:737, 9:731Þ and τ ae p ∼ Beta ð384:49, 33:25Þ). Similarly, for ewe and ram survival, we used the results from Proffitt et al. (2021) for ewe survival to construct an informed prior to use for sex-specific, time-varying survival rates for each population (S e p,t ∼ Betað10:452, 3:292Þ and S r p,t ∼ Betað10:452, 3:292Þ). Although this previous work was specific to ewes, we considered this distribution (with an approximate mean of 0.76 and standard deviation of 0.11) to be sufficiently vague so as to serve as a general prior for adult survival. For the initial sizes of each age and sex class in each population, we used a vague uniform prior (Uniform(0, upper p )), where upper p was an upper limit for each population defined as the maximum total count in the time series, that is, the number of individuals in each age and sex class in the initial population had to be less than the maximum aggregate count of all classes). Observation-level random effects were drawn from a normal distribution (ζ t ∼ Normal 0, σ 2 ζ ) with a vague prior for the standard deviation (σ ζ~U niform(0, 10)).
Distributions of vital rates and their association with demographic performance.-The first of our three goals was to estimate the key vital rates underlying our process model, and our second was to characterize the associations between vital rates and the observed trajectories of bighorn sheep populations. However, unlike survival and pregnancy rates, we did not have an informed prior to place on lamb survival. To provide a moderate amount of shrinkage and improve estimation for years with missing data while not overly influencing estimates from years with v www.esajournals.org complete data, we modeled all population-year lamb survival rates, S l p,t , using a random effects structure on the logit scale (REmodel), where logit S l with a Logistic(0, 1) prior for μ and a vague prior for the standard deviation σ ξ~U niform(0, 10). It is important to note that lamb survival and sexspecific adult survival were assumed to be independent (i.e., no explicit modeling of the potential covariation between rates), after initial model runs suggested estimation problems for this more complex model using count data and generally vague priors for survival and pregnancy (Koons et al. 2017).
Rather than a prospective analysis such as lifestage simulation analysis (Wisdom et al. 2000) that estimates how hypothetical changes in vital rates may impact population growth rates, we wanted to decompose the observed variation in population growth rates into contributions from the vital rates. Life-table response experiments (LTREs) provide a powerful method of investigating these drivers of observed variation in population growth rates, and recent work has derived a series of transient LTREs that do not assume a stationary environment (Cooch et al. 2001, Koons et al. 2016, 2017. The population growth rate in year t (λ t ) can be expressed as a function of underlying vital rates and the population structure in year t − 1, and sensitivities of population growth rates to demographic parameters derived using partial derivatives (Appendix S2). Our study focused on: (1) estimating the contribution of each of the vital rates to the total observed variation in λ t (contribution varðλ t Þ θ i for vital rate θ i from Koons et al. 2016) and (2) understanding the drivers of changes in λ t between successive time steps (contribution Δλ t θ i for vital rate θ i from Koons et al. 2016) (details in Appendix S2).
Sudden, precipitous declines in vital rates that are associated with epizootics, shifts in predation risk, or severe environmental conditions are a persistent source of variation in many populations of bighorn sheep (Ross et al. 1997, Portier et al. 1998, Festa-Bianchet et al. 2006, Cassirer and Sinclair 2007, Manlove et al. 2016). For such circumstances, our use of informed priors for adult survival rates and a single, common distribution of lamb survival rates could have been too restrictive. The priors for adult survival might have been too informative to capture these dramatically lower rates, and rarity of these events resulted in the estimation of the distribution of lamb survival rates being dominated by "normal" years and masking the signal of these events. To evaluate the sensitivity of our results to our prior specification, we constructed an alternative version of the model wherein both adult and lamb survival were modeled using mixture distributions that allowed for more variation in vital rates (Appendix S2).
Correlates of lamb survival.-Our third main goal was to assess a series of covariates as potential sources of variation in lamb survival. Our a priori expectation was that lamb survival should be positively related to favorable conditions during the growing season (e.g., indexed by the normalized difference vegetation index [NDVI] or precipitation) and negatively related to winter conditions (e.g., indexed by indices of winter severity). We used a derived metric of NDVI as an index of primary production during the growing season (Pettorelli et al. 2011). Using the AVHRR (Advanced Very High Resolution Radiometer) time series of NDVI values (Vermote and NOAA CDR Program 2019), we performed a series of smoothing steps on each pixel in our study area. First, we calculated the monthly median of the daily NDVI values produced by the AVHRR time series and then applied a smoothing spline on the monthly time series of values. Second, we estimated the start of the growing season as the point in the smoothed signal at which the value of NDVI first attained 50% of the seasonal maximum, and the end of the season as the first point in time past the seasonal maximum at which the smoothed signal that the NDVI value was below 50% of its seasonal maximum (Bradley et al. 2007). Third, to ameliorate the impact of canopy cover on the NDVI signal, we calculated the integrated NDVI from the start of the season to the end of the season as the sum of the differences between each NDVI value and the value of NDVI at the start of the growing season (Ruimy et al. 1994, Jönsson andEklundh 2004). Prior work has suggested that population responses to NDVI might depend on the timing (Paterson et al. 2019), that is, high NDVI during the spring could indicate substantial primary production and favorable growing conditions conducive to neonates that could be coupled to a late season characterized by drought and poor growing conditions such that a single NDVI metric for an entire growing season could be misleading. Therefore, to investigate differential impacts of the seasonal timing of NDVI (e.g., early-season vs late-season growing conditions), we defined the integrated NDVI over two periods: from the start of the growing season to June 30 (NDVI early ), and from July 1 to the end of the growing season (NDVI late ). Additionally, we derived a second metric of environmental conditions during the growing season. Using the PRISM data set (PRISM Climate Group, Oregon State University, http://prism.ore gonstate.edu), we derived precipitation metrics similar to those for NDVI for each pixel in our study area as the sum of precipitation from May 1 to June 30 (PREC early ) and the sum of precipitation from July 1 to September 30 (PREC late ). Moreover, we used the sum of precipitation from December to March as an index of winter severity (PREC winter ). To connect the values of each of these five metrics to the populations in our study, we relied on estimated seasonal ranges for these populations using either GPS collar data from a regional research project on these populations (n = 16, details in Proffitt et al. 2021), or seasonal ranges derived from professional opinion (n = 1, Dubois Badlands). Environmental conditions in each biological year were approximated by using the median value of all the pixels contained in each seasonal range for that year.
We evaluated potential sources of environmental variation in lamb survival using logistic regression to connect covariates in population p and year t (e.g., NDVI early p,t ) to survival using the logit link (COVmodel): where the population-specific mean, μ, was given a Logistic(0, 1) prior and regression coefficients were given independent Normal(0, 4) priors that are vague on the logit scale (Northrup and Gerber 2018).
There were two additional potential sources of variation in lamb survival that were not amenable to assessment using our full data set. First, predation on bighorn sheep lambs can have a dramatic impact on survival rates (Berger 1991, Rominger et al. 2004, Festa-Bianchet et al. 2006, Brewer et al. 2014) and be critical to consider when evaluating potential sources of variation in population dynamics. Although the larger region has seen the recovery of populations of multiple predator species during our study period, mountain lions (Puma concolor) are implicated as the primary predator affecting bighorn sheep population dynamics (Rominger 2018). However, there is no direct measure of predation pressure from mountain lions or an index of it (e.g., mountain lion abundances) for a region the size of our study area and over the timespan of the data collection. Therefore, to assess the potential impact of mountain lions on lamb survival, we developed a spatial covariate by extrapolating a previously published winter resource selection function (RSF) for Montana (Robinson et al. 2015) over the extent of our study area and then used the median value of the RSF within each seasonal range as an index of risk from mountain lions (details in Appendix S2). Due to the underlying assumption of stationarity, we limited the data set for this analysis to the last 10 yr of observations for each population. We then evaluated the strength of evidence for a relationship between lamb survival and the median value of the extrapolated mountain lion RSF for each population by modifying Eq. 15 to include the population-year value of the derived covariate (LION). Second, disease processes are thought to play a pivotal role in the population dynamics of bighorn sheep. To explore the potential role of disease as a limiting factor in lamb survival, we assessed the evidence for a relationship between pathogen communities in these populations and lamb survival. Using the results of the most comprehensive and rigorous study to date on the composition of communities of pathogens associated with respiratory disease in bighorn sheep (Butler et al. 2018), we compared the distributions of lamb survival for those populations that hosted both M. ovipneumoniae and Pasteurellaceae bacteria (Bibersteinia trehalosi, Pasteurella multocida, and Mannheimia haemolytica) (n = 12) to those populations that only hosted v www.esajournals.org Pasteurellacea (n = 4). Pathogen data for the Spanish Peaks population were not published as part of this study, and data were censored from this population for this portion of our analysis. Similar to the mountain lion covariate, our comparison of lamb survival values between the two groups assumed stationarity, that is, the presence of the pathogen(s) was considered a fixed characteristic that did not vary through time. Therefore, we only used the last 10 yr of estimates of lamb survival for each population from the first version of our model (Eq. 14, REmodel) and used the approximate posterior distribution of estimated survival values to estimate the median lamb survival for all years as well as annual values for both groups.
The covariates used to assess sources of variation in lamb survival showed substantial variation both across and within populations (Fig. 2, Appendix S3). Prior to being used in analyses, all covariates were screened for collinearity (with a threshold of 0.5) and were centered and scaled using the mean and one standard deviation. To facilitate interpretation for estimated covariate effects, for each covariate for which the 90% highest posterior density interval of the estimated regression coefficient did not overlap zero, we made predictions for the populationspecific relationship between the covariate and lamb survival by holding all other covariates to their mean. For each covariate, we then selected a set of example populations to demonstrate how the predicted relationship translated into variation in lamb survival from the 5th percentile of the covariate values (S l (covariate 0.05 )) to the 95th percentile of the covariate values (S l (covariate 0.95 )).

Model estimation and evaluation
Models were fit in the R programming environment (R Development Core Team 2019) using the runjags package (Denwood 2016) as an interface to the JAGS program (Plummer 2017) for Markov chain Monte Carlo (MCMC) estimation of our model. For both versions of our model (REmodel and COVmodel), we ran two parallel chains for 150,000 iterations and discarded the first 50,000 as the burn-in and adaptation phase. Due to the large number of derived parameters saved from these simulations, we kept every 20th iteration (i.e., a thinning interval of 20), which resulted in a combined total of 10,000 samples used for inference. Convergence of these chains was graphically assessed using traceplots for the underlying vital rates. For the large number of derived parameters (population-, sex-, age-, and time-specific sizes of populations), a graphical assessment was not practical and convergence was evaluated using the Gelman-Rubin statistic and convergence was assumed for values <1.05 (Gelman and Rubin 1992). We summarized the approximate posterior distributions for all of the model parameters using the median and 90% highest posterior density interval.
There are no standard test statistics for evaluating the adequacy of the fit of state-space models to time series data such as these. Therefore, we evaluated the goodness of fit of both versions of our model to each population using two posterior predictive checks (Gelman et al. 1996) that compared the fit of the models to the data using two biologically relevant characteristics of the data: the observed population growth rates and the observed ratio of lambs to ewes (both yearling and adult) (details in Appendix S2). For each iteration of the MCMC chain for each model, we derived a replicated data set to compare to the observed data. We calculated a Bayesian p value as the proportion of iterations where the value of the statistic for the observed data exceeded the value of the statistic for replicated data. Bayesian p values close to 0.5 indicate the model is faithfully replicating the variation of the observed data, whereas values close to 0 or 1 indicate a severe lack of fit (Gelman et al. 1996).

RESULTS
The survey data were comprised of (typically) discontinuous time series of counts and age and sex classifications for each population from 1983 to 2018 ( Fig. 2; for population-specific summaries, see Appendix S1). In total, the data set included 470 population-years of observations (including years missing observations within the time series). Across populations, there was considerable variation in the median total count from Targhee (68 animals) to Wapiti Ridge (716 animals) coupled to similar variation in range of counts from Targhee (maximum-minimum = 59) to Whiskey Mtn-East (maximum-minimum = 640). The ratio of lambs to ewes across populations demonstrated similar variation with the median ratio ranging from a minimum in Lost Creek (0.14) to a maximum in Middle Missouri (0.47) and the range of ratios from Francs Peak (maximum ratio-minimum ratio = 0.24) to Whiskey Mtn-West (maximum-minimum = 0.81).
Convergence was achieved for parameter estimates in both models; bR values for all latent states (i.e., sizes of sex/age classes) were <1.05, and bothbR values and visual inspection of traceplots for underlying vital rates of both models indicated convergence and adequate mixing. All of the posterior predictive checks (two posterior predictive checks for the REmodel and COVmodel) suggested adequate fit for the models, given that the one-sided Bayesian p values were far from 0 or 1 for both population growth (REmodel: p = 0.67; COVmodel: p = 0.52) and age ratios (REmodel: p = 0.26; COVmodel: p = 0.81).

Variation in vital rates
Using the REmodel for inference, we found substantial variation in the vital rates (Fig. 3). Finally, we estimated the geometric means for each population and found that nine populations had estimated geometric means with the upper limit of a 90% highest posterior density interval <1 (indicating long-term population decline; Castle Reef, Clarks Fork, Francs Peak, Lost Creek, Targhee, Trout Peak, Wapiti Ridge, Whiskey Mtn-East, and Younts Peak), two populations had estimated geometric means with a lower 90% highest posterior density interval limit >1 (indicating long-term population growth; v www.esajournals.org Middle Missouri and South Madison), and the remainder had 90% highest posterior density intervals that included 1 (Fig. 3).
We found evidence for constraining relationships between lamb survival, ewe survival, and the population growth rate for the female component of the population ( bλ f ) (Fig. 4). Using the medians of the approximate posterior distribution from the REmodel for each annual parameter as a point estimate, we found that lamb survival rates corresponding to estimates of bλ f ≥ 1 had a similar distribution (median = 0.46; range =-0.21-0.78) to those corresponding to bλ f < 1 (median = 0.28; range = 0.05-0.61). There was similar ambiguity in the relationship between ewe survival and λ f . The distribution of ewe survival rates corresponding to estimates of bλ f ≥ 1 (median = 0.86; range = 0.77-0.90) was essentially the same as the distribution of values corresponding tobλ f < 1 (median = 0.83; range = 0.60-0.90). We found that the association between these vital rates depended on their values, such that the impact of lamb or ewe survival on λ f depended on the value of the other rate (Fig. 4). For example, at values of ewe survival ≥0.85, we found λ f ≥ 1 over a wide distribution of lamb survival values (median = 0.43; range = 0.21-0.68). Once ewe survival dropped to ≤0.80, however, we found that λ f ≥ 1 over only a narrow and high distribution of lamb survival values (median = 0.58; range = 0.54-0.68). If ewe survival dropped below 0.75, no value of lamb survival was capable of yielding λ f ≥ 1. Conversely, when lamb survival rates were ≥0.40, we found that λ f ≥ 1 over a relatively wide distribution of ewe survival rates (median = 0.85; range = 0.77-0.90). When lamb survival rates were ≤0.20, no value of ewe survival resulted in λ f ≥ 1. Finally, we noted a relationship between the growth rate of the female population (λ f ) and the number of females removed from each population due to either harvest or translocation, which was incorporated into the annual survival rates in the model. We calculated the ratio of the number females removed in year t + 1 to the estimated size of the population from the model in year t and compared it to λ f in year t. We found that when the ratio of female removals to female population size exceeded approximately 0.18, all point estimates of λ f were <1.

Transient life-table response experiments
We decomposed the observed variation and changes in the population growth rates of the female portion of each population into contributions from changes in constituent demographic parameters using the REmodel for inference and found substantial among-population differences in drivers of population trajectories. Partitioning the total variation in population growth rates (see Fig. 3 for the distribution of point estimates for annual growth rates) into the contributions from lamb survival, ewe survival, and changes in the age structure revealed that, in general, variation in ewe survival accounts for a large amount of the variation in bλ f (Fig. 5). For six of the 17 populations (Clarks Fork, Jackson, Lost Creek, Paradise, Petty Creek, and Targhee) the 90% highest posterior density interval for the estimated proportion of variation inbλ f attributed to variation in ewe survival did not overlap the similar estimate for the proportion of variation inbλ f due to variation in lamb survival, strongly suggesting that variation in ewe survival dominated variation in bλ f for these populations. For two of the populations (Middle Missouri and Trout Peak), overlap between the credible intervals prevented strong inference but suggested a similar relationship. For seven of the 17 populations (Clarks Fork, Dubois Badlands, Francs Peak, Spanish Peaks, Wapiti Ridge, Whiskey Mtn-East, and Younts Peak), the largely coincident credible intervals suggest roughly equal contributions from lamb and ewe survival. In only one case Fig. 4. Relationship between key vital rates and the growth rates of the female component of the population (λ f ). Panel (a) represents the univariate relationships between lamb and ewe survival and estimated λ f . Panel (b) represents the multivariate relationship between these two key underlying vital rates and λ f . For panel (a), the dots are the point estimates of rates. For panel (b), the colored and contoured surface represents a simple linear interpolation of the underlying estimated values (shown by the rug on the xand y-axes).
(South Madison) did we find suggestive evidence that lamb survival accounted for more variation in bλ f than did ewe survival. Variation in the age class structure never accounted for more than 3% of the variation inbλ f . We also decomposed the actual changes in population growth rates between years into the contributions from demographic parameters. In all cases, the changes in population growth rates were driven by variation in lamb survival and ewe survival with only a trivial contribution from changes in the age structure (unsurprising given the simple age structure used, that is, only yearlings and adults), and we focused on variation in those two vital rates as drivers of population growth rates. Our results suggest a set of potential mechanisms associated with changes to the trajectory of a population. When a population's trajectory changed from increasing or stable ( bλ f t ≥ 1) to decreasing ( bλ f tþ1 < 1) (n = 63, upper right quadrant of Fig. 6), declines in lamb survival were involved in 54 (86%) of the cases and declines in ewe survival were involved in 31 (41%) of the cases. Notably, in 32 of the 54 cases associated with declines in lamb survival (59%), the population trajectory was reversed despite increases in ewe survival. In nine of the 31 cases associated with declines in ewe survival (29%), the trajectory was reversed despite increases in lamb survival. When a population's trajectory changed from decreasing ( bλ f t < 1) to increasing or stable ( bλ f tþ1 ≥ 1) (n = 60, lower left quadrant of Fig. 6), increases in lamb survival were involved in 57 (95%) of the cases, and increases in ewe survival were involved in 40 (67%) of the cases. A small number (three, or 5%) of cases involved Fig. 6. Across-population results of the life-table response experiment analysis of the relationship between demographic parameters and the changes in the estimated population growth rate of the female component of the population λ f . The columns represent λ f in year t (<1 or >1), the rows represent λ f in year t + 1, the xand y-axes represent the contribution of changes in ewe survival and lamb survival to changes in λ f (Δλ f ), and the color represents Δλ f . For example, the upper right quadrant is the case where the population was increasing in year t and declining in year t + 1, and the dots denote how ewe survival and lamb survival contributed to that change in λ f . increases in ewe survival and decreases in lamb survival, whereas 20 (33%) involved increases in lamb survival and decreases in ewe survival. For successive years in which population growth rates were below 1 (λ f t < 1, bλ f tþ1 < 1) (n = 159 cases, upper left quadrant of Fig. 6), there was considerable variation in the contribution of ewe survival and lamb survival to λ f t , including positive changes that were insufficient to reverse the negative trajectory of the population. Similarly, for successive years in which population growth rates were at or above 1 (λ f t ≥ 1, bλ f tþ1 ≥ 1) (n = 154 cases, upper left quadrant of Fig. 6), there was no clear pattern as the positive trajectory of the population absorbed wide variation in cases of negative contributions from lamb and ewe survival. Together, these patterns suggest that these populations are subject to a complex interplay of changes in vital rates, and changes in single vital rates are typically poor indicators of population trajectories.
Finally, we note that our inference on the distributions of vital rates and sources of variation in population growth rates was essentially identical to results from the mixture model (Appendix S4). Given the mixture model was an attempt to capture dramatic changes in vital rates due to the large-scale mortality events in the data set, the lack of any clear difference in inference between the two modeling approaches suggests that our population model, which was estimated using management data alone, may not be adequate to capture these die-off events.

Correlates of lamb survival
Our second major goal was to understand correlates of lamb survival, and our first analysis using the COVmodel revealed a diverse set of relationships between environmental conditions and lamb survival. We found strong support for our a priori hypothesis that lamb survival would be negatively associated with winter severity as indexed by cumulative winter precipitation. In eight out of 17 populations, the association between PREC winter and lamb survival was negative and the associated 90% highest posterior density intervals for bβ PREC winter did not overlap zero: The strongest effect was estimated for Lost Creek ( bβ PREC winter ¼ À1:63½À2:19, À 1:09, and the weakest effect was estimated for Wapiti Ridge ( bβ PREC winter ¼ À0:14½À0:23, À 0:04). These estimated effects translated into substantial changes in predicted lamb survival over the range of observed PREC winter values (Fig. 7). For Lost Creek, predicted lamb survival decreased over the range of PREC winter values from S l PREC winter 0:05 ð Þ ¼ 0:58 ½0:41, 0:79 to S l PREC winter 0:95 ð Þ ¼ 0:01 ½0:002, 0:02, and lamb survival for Wapiti Ridge decreased from S l PREC winter 0:05 ð Þ ¼ 0:26 ½0:24, 0:29 to S l PREC winter 0:95 ð Þ ¼ 0:17 ½0:13, 0:22. We found no evidence that the strength of the association between winter precipitation and lamb survival depended on the range of winter conditions experienced (Fig. 8).
In contrast to the strong and consistent support that we found for the relationship between lamb survival and winter severity, we found conflicting evidence regarding the relationship between lamb survival and growing conditions as indexed by summer precipitation. Our results suggested both positive and negative relationships between lamb survival and PREC early and PREC late covariates. Where the 90% highest posterior density interval for the regression coefficient for PREC late did not overlap zero (five out of the 18 populations), three populations had positive relationships with PREC late (Dubois Badlands, Francs Peak, and Lost Creek) and two populations had negative relationship (South Madison and Whiskey Mtn-East). The strongest positive effect was for Lost Creek ( bβ PREC late ¼ 1:25 ½0:74, 1:80, which translated into a predicted increase in lamb survival from S l PREC late 0:05 ð Þ¼ 0:03½0:01, 0:06 to S l PREC late 0:95 ð Þ¼0:81 ½0:61, 0:98 (Fig. 7). The strongest negative effect was for the South Madison population ( bβ PREC late = −0.60 [−1.01, −0.17], which translated into a predicted decline in lamb survival from S l PREC late 0:05 ð Þ¼ 0:55½0:40, 0:70 to S l PREC late 0:95 ð Þ¼0:23 ½0:12, 0:32 (Fig. 7). Where the 90% highest posterior density interval for the regression coefficient for PREC early did not overlap zero (four out of the 18 populations), two populations had positive relationships with PREC early (Castle Reef and Lost Creek) and two populations had negative relationship (Dubois Badlands and Wapiti Ridge). The strongest positive effect was for Lost Creek ( bβ PREC early = 1.26 [0.81, 1.70]), which translated into a predicted increase in lamb survival from S l PREC early 0:05 = 0.04 [0.01, 0.07] to S l PREC late 0:95 ð Þ¼0:66½0:47, 0:84. The strongest negative effect was for Dubois Badlands (bβ PREC early ¼ À0:56½À1:06, À 0:07, which translated into a predicted decline in lamb survival from S l PREC early 0:05 ¼ 0:46 ½0:30, 0:64 to S l PREC early 0:95 ¼ 0:13½0:04, 0:23. Although we found no strong evidence for a relationship between the strength and/or direction of the relationship between early-or late-summer precipitation and lamb survival (Fig. 8), the patterns for the coefficient estimates whose 90% highest posterior density intervals do not include zero provide some evidence that populations that experience higher late-season precipitation have a negative coefficient, whereas populations that experience higher earlyseason precipitation have a positive coefficient.
We found equivocal evidence for a positive relationship between lamb survival and growing conditions as indexed by NDVI. Similar to results for early and late-summer precipitation, we estimated both positive and negative relationships between lamb survival and NDVI early and NDVI late covariates. Where the 90% highest posterior density interval for the regression coefficients for NDVI late did not overlap zero (seven out of the 18 populations), two populations had positive relationships with NDVI late (Dubois Badlands and Francs Peak) and five populations had negative relationship (Jackson, Lost Creek, Middle Missouri, Wapiti Ridge and Whiskey Mtn-East) (Fig. 9). The strongest positive effect was for Dubois Badlands ( bβ NDVI late ¼ 0:68 ½0:02, 1:38), which translated into a predicted increase in lamb survival from S l NDVI late 0:05 ð Þ ¼ 0:10 ½0:01, À 0:25 to S l NDVI late 0:95 ð Þ ¼ 0:53 ½0:30, 0:79. The strongest negative effect was for the Middle Missouri population ( bβ NDVI late ¼ Fig. 7. Population-specific predicted relationships between precipitation and lamb survival. All predictions were constructed by holding other covariates to their mean value. The black line represents the median of the estimated response to change in that covariate (the gray ribbon represents the 90% highest posterior density interval). À0:62 ½À1:23, À 0:06), which translated into a predicted decline in lamb survival from S l NDVI late 0:05 ð Þ ¼ 0:84 ½0:68, 1:00 to S l NDVI late 0:95 ð Þ ¼0:43 ½0:28, 0:61. Where the 90% highest posterior density interval for the regression coefficients for NDVI early did not overlap zero (seven out of the 18 populations), three populations had a positive relationship with Fig. 8. Relationships between estimated regression coefficients for normalized difference vegetation index (NDVI) and precipitation (PREC) covariates and the median value of each covariate for each population. For each of our five covariates, we graphed the median value through time for each population (x-axis) against the estimated regression coefficient (y-axis). Estimated regression coefficients whose 90% credible interval did not overlap zero are in black; the remainder are in gray. The dot represents the median and the line the 90% highest posterior density interval. NDVI early (Castle Reef, Lost Creek and Spanish Peaks) and four populations had a negative relationship (Jackson, Paradise, Wapiti Ridge, and Younts Peak). The strongest positive effect was for Castle Reef ( bβ NDVI early ¼ 0:73½0:37, 1:11), which translated into a predicted increase in lamb survival from S l NDVI early 0:05 ¼ 0:04 ½0:01, 0:08 to S l NDVI early 0:95 ¼ 0:27 ½0:18, 0:35. The strongest negative effect was for Wapiti Ridge ( bβ NDVI early ¼ À0:32 ½À0:46, À 0:18), which translated into a predicted decline in lamb survival from S l NDVI early 0:05 ¼ 0:32½0:28, 0:37 to S l NDVI early 0:95 ¼ 0:15½0:11, 0:18. Similar to the results for precipitation, we found no strong evidence for a relationship between the magnitude and/or direction of these estimated effects and the range of conditions experienced (Fig. 8), but suggestive evidence for a relationship in which populations with a higher late-season NDVI had a negative coefficient, and populations with a higher early-season NDVI had a positive relationship. We note that the predicted relationships between growing season conditions, winter severity, and lamb survival have important implications for population growth rates. As Fig. 9. Population-specific predicted relationships between normalized difference vegetation index (NDVI) and lamb survival. All predictions were constructed by holding other covariates to their mean value. The black line represents the median of the estimated response to change in that covariate (the gray ribbon represents the 90% highest posterior density interval). demonstrated in the previous section, lamb survival plays a significant role in reversing population trajectories. The magnitudes of these predicted changes in lamb survival associated with environmental variation are large enough to suggest environmental drivers can play a significant role in altering population growth rates, for example, the predicted decline in lamb survival associated with winter precipitation for the Clarks Fork population (S l PREC late 0:05 ð Þ¼0:39 ½0:30, 0:46 to S l PREC late 0:95 ð Þ ¼0:15½0:10, 0:22) was sufficient to push lamb survival below the threshold level where no value of ewe survival could have resulted in λ f ≥ 1.
We used the last 10 yr of observations and the COVmodel to evaluate evidence for whether mountain lion densities indexed by the quality of mountain lion habitat might serve as a potential limiting factor for lamb survival. Our estimate was small and quite imprecise such that we found essentially no evidence for such a relationship ( bβ LION ¼ 0:09 ½À0:62, 0:79).
Our assessment of the relationship between pathogen communities and lamb survival for the last 10 yr of observations provided strong evidence that the median value of lamb survival for those populations that hosted both M. ovipneumoniae and Pasteurellaceae (0.24 [0.22, 0.26]) was lower than the median value of lamb survival for populations that hosted Pasteurellaceae alone (0.41 [0.37, 0.46]) (Fig. 10a), with an estimated difference between the medians of 0.18 [0.13, 0.23]. However, although the distribution of lamb survival values for populations that hosted both pathogens was lower than the distribution for populations that hosted Pasteurellaceae alone, we note a considerable amount of overlap (Fig. 10b). Approximately 48% of the point estimates for annual lamb survival values for populations that hosted both pathogens were equal to or exceeded the lowest estimated lamb survival value for populations that hosted Pasteurellaceae alone; moreover, approximately 18% were equal to or greater than 0.40, a threshold of lamb survival we previously associated with values of λ f ≥ 1 despite a wide range in ewe survival.

DISCUSSION
Results from our novel state-space model in a Bayesian framework highlight the benefits of using population models and information on vital rates to estimate vital rates and evaluate how variation in vital rates shapes the trajectories of populations. Although we note that populations of ungulates are successfully managed with a conventional use of age ratios and count data (Bender 2006, Harris et al. 2008, our model allows for a more comprehensive assessment of the mechanisms by which population trajectories may vary. We demonstrated a variety of mechanisms by which variation in vital rates affects changes in population growth rates, including combinations of lamb survival and adult survival that result in similar population growth rates, and suggest that these mechanisms may differ among populations. Our assessment of sources of variation in lamb survival rates indicated a potentially complex interplay of ecological and disease processes affecting lamb survival, including generally strong effects of winter climate, contrasting effects of summer growing conditions, and no effect of mountain lion predation or relationship with resident pathogen community. Together, these results offer a rigorous and comprehensive evaluation of bighorn sheep vital rates. Vital rates demonstrated substantial variation across years and populations such that a variety of combinations of vital rates could yield similar population trajectories. The variation in lamb survival that we observed is consistent with results from decades of empirical work for other ungulates that documents high variability in offspring survival rates , Gaillard and Yoccoz 2003, Raithel et al. 2007). The evidence for very low lamb survival rates in some years and populations highlights the concern that this rate could play a limiting role in conservation and restoration efforts for bighorn sheep in particular (Johnson et al. 2010). In contrast to the expectation of high adult survival with minor variation for many ungulate populations, our estimated adult ram and ewe survival rates were also highly variable. This is an unsurprising result for populations of ungulates like bighorn sheep that are subject to stochastic disease and severe weather events, harvest, or predation that all result in lower and highly variable adult survival rates (Owen-Smith and Mason 2005, Nilsen et al. 2009, Eacker et al. 2017. The integrated effects of this variation in vital rates on the population growth, particularly in small populations, is poorly understood. To our knowledge, our work is the first to empirically define the existence of a demographic "safe space" (de Silva and Leimgruber 2019) for bighorn sheep, that is, combinations of vital rates that result in positive growth rates. These results indicate that management can address declines in populations by mitigating either (or both) lamb survival and adult survival to yield positive growth rates (e.g., predator control, treatment for respiratory disease, and supplemental provisioning) and provide minimum thresholds of vital rate values required for positive population trajectories, below which removals from translocations or harvest may need to be reduced and/or augmentation may be necessary for population persistence (Johnson . Although our results suggest that a constraining relationship may exist between removals and the growth rate for the female population, this relationship was obscured by two sources of uncertainty that prevented strict inference: (1) uncertainty in estimates of growth rates and (2) uncertainty in the relationship between the larger population from which animals are removed and the portion of this larger population that is surveyed (uncertainty that arose from the distribution of animals on the landscape and the quality of the survey process).
The expectation for the dynamics of unharvested ungulate populations lacking substantial exposure to predation and/or disease is that low variation in adult survival (despite high elasticity) renders the high variation in offspring recruitment the dominant source of variation in population dynamics (despite lower elasticity) , Raithel et al. 2007. A key component of that expectation is the low variation in adult survival; where this rate is affected by harvest, disease, or predation, it is unsurprising that variation in adult survival can dominate population dynamics given its higher elasticity. Our approach found a diversity of relationships between vital rates and population growth rates among populations similar to results from other work on bighorn sheep that used asymptotic approaches (Johnson et al. 2010). Although variation in adult survival explained the most overall variation in population growth rates, improved lamb survival was found to be the dominant mechanisms by which populations actually reversed declines (Manlove et al. 2016). The relative importance of these mechanisms was not consistent between populations, and we add to the small number of studies suggesting the drivers of population dynamics for ungulates may vary spatially (Albon et al. 2000, Coulson et al. 2005, Garrott et al. 2008a, Nilsen et al. 2009, Johnson et al. 2010. These results also help further broaden the perspective on the drivers of population dynamics for ungulates. We suggest that any apparent conflict between these results and the expectation for ungulate population dynamics merely reflects the fact that populations live in different portions of the vital rate parameter space: Where adult survival is high and nonvariable, growth rates are driven by variation in offspring recruitment, but where adult survival is variable due to intrinsic/extrinsic drivers, it can overwhelm the role of offspring recruitment, particularly for small and declining populations (Nilsen et al. 2009). These results caution against a one-size-fits-all approach to conservation and restoration efforts for bighorn sheep and highlight how managing for a particular objective requires a careful assessment of the range of vital rates in the population so as to determine the most efficient strategy.
The diversity of mechanisms that underlie changes in population growth rates shown here is at apparent odds with two prevailing ideas about population dynamics of bighorn sheep. All-age die-offs due to disease events in bighorn sheep can potentially devastate populations, both during the initial epizootic and for years afterward as enzootic disease reduces lamb survival (Cassirer and Sinclair 2007). Although our results do provide some evidence for years in which survival rates of both lambs and adults precipitously declined (consistent with all-age die-offs), we also demonstrated that population trajectories can change due to a diversity of vital rate combinations, many of which are not consistent with all-age die-offs. We found evidence for both positive reversals in trajectory due to improved lamb survival even as adult survival declines, and negative reversals in trajectory due to decreased lamb survival even as adult survival improves. Although all-age die-offs undoubtedly play a dominant role in population dynamics when they occur, population growth rates are also shaped by a complex interplay of variation in other vital rates that are (likely) decoupled from the disease process as well as from one another, for example, declines in survival induced by harvest. Second, our work suggests that the presence of pathogens associated with respiratory disease in bighorn sheep is insufficient on its own to account for the realized level of variation in population dynamics. Consistent with prior work on many of these same populations (Butler et al. 2018), we found that the presence of M. ovipneumoniae was not a sufficient factor in limiting population growth. Despite our conclusion that populations hosting M. ovipneumoniae and Pasteurellacea had an overall lower median lamb survival than those that hosted Pasteurellacea alone, the ranges of lamb survival estimates for both types of populations were v www.esajournals.org very similar and included values that resulted in positive population growth rates and reversals of declining population trajectories, regardless of ewe survival rate values. If variation in offspring survival was not well-explained by the composition of pathogen communities, then ecological and environmental factors have to be the key players either independent of resident pathogen communities or interacting with pathogen communities to impact disease expression. Overall, our work strongly suggests that multiple ecological factors are important drivers of population dynamics with similar or stronger effects on bighorn sheep population trajectories as compared to disease epizootics (short of all-age die-offs), and thus, bighorn sheep management and restoration should focus on advancing ecological understanding of the species, rather than assume disease dynamics are predominately responsible for population trajectories.
Our work supports this supposition by demonstrating the significant impact that variation in environmental conditions can have on lamb survival. In contrast to results for other sympatric ungulates where correlates of offspring recruitment and survival have been welldocumented to include intrinsic factors such as body mass and extrinsic factors such as predation and environmental factors (Singer et al. 1997, Côté and Festa-Bianchet 2001, Garrott et al. 2003, 2008b, White et al. 2010, Creel et al. 2011, Lukacs et al. 2018, fewer assessments of similar drivers have been done for bighorn sheep lambs (Douglas and Leslie 1986, Hass 1989, Portier et al. 1998, Douglas 2001, Brewer et al. 2014. Results of those assessments suggest that bighorn sheep lamb survival can be affected by a similarly diverse set of processes and that lambs may be quite susceptible to predation from mountain lions (Rominger 2018), and our results are broadly consistent with those findings, with the notable exception that we found no relationship between lamb survival and our proxy for the risk of predation from mountain lions. We attribute this result to the difficulties inherent in developing a suitable proxy for predation in the absence of relationship between lamb survival and growing season conditions (as measured by spring and summer precipitation and primary production) can vary among populations and suggest that the nature of the relationship might differ depending on the range of conditions to which a population is exposed, a possibility that will require further work to elucidate. We speculate that these different responses may arise from the manner in which populations interact with their environment, for example, where lateseason precipitation is typically low, a high annual value may indicate greener conditions during the later season, whereas where lateseason precipitation is typically high, a high annual value may indicate late-season snow events that predisposes lambs to nutritional or other physiological stresses. There are two possible explanations for the ambiguity of these results. First, it is possible that there is a relationship between the range of environmental conditions experienced by a population and the response of lamb survival that is obscured by either our small sample size of populations (n = 17), or the indirect method we used to estimate lamb survival. Second, the broad scale of our covariate derivation (median values for entire ranges within a year) only captures a portion of the variation actually experienced by populations. These populations are made up of individuals that experienced different environmental conditions due to different movement strategies (Lowrey et al. 2019), and our derived covariates may be a crude average of more complex processes. Future work using individual collared animals and an appropriate sampling design that accounts for such a population structure is required.
There are two further limitations of our approach that impact the interpretation of our results and the broader applicability. First, uncertainty in the observation process can cripple inference, and although our approach yields improved inference on vital rates and population trajectories, we suggest that the quality of the inference is related to the quality of the survey data. One practical possibility for improving the observation process would be to keep some component of the population collared to continually reassess the distribution of the animals on the landscape to inform survey coverage and estimate survey-specific detection probabilities. Second, our method is not a substitute for simultaneously monitoring vital rates and population trends, and we likely suffer from low precision on estimates of vital rates due to v www.esajournals.org estimating vital rates based on count data alone. We speculate that relying solely on count data reduced our power to detect precipitous drops in vital rates using flexible mixture models and to estimate the relationships between lamb survival and environmental covariates. We acknowledge that some of the relationships identified here may be spurious and actually represent the conflated effects of unmodeled processes (e.g., stochastic die-off events, resulting multi-year lower lamb survival values and environmental processes). Future work that incorporates data at the individual level is needed to more clearly elucidate the relationships between disease, predation, environmental variation, and lamb survival.
In conclusion, we emphasize our broadscale finding that the majority of the bighorn sheep populations in this study had marginal population performance. Given the life-history characteristics of this species, their overall low density and assuming that resources are not limiting, a naïve expectation for population growth from broadly similar species suggest that they are capable of rates >1.2 per year (McCorquodale et al. 1988, Eberhardt et al. 1996. Although marginal population performance is the intended result of harvest and translocations for some bighorn sheep populations in this study (by reducing population growth [Jorgenson et al. 1993] and rendering populations less prone to epizootics [Monello et al. 2001]), we found that that 15 of 17 populations had geometric mean growth rates ≤1, and no point estimate for an annual value of λ >1.08, indicating a general disjunction between the species' potential and reality. Unfortunately, our work demonstrated multiple mechanisms by which variation in vital rates affects population growth rates and that these mechanisms may differ among populations, such that conservation and restoration efforts need to be populationspecific. We suggest that the first steps to improving our ability to understand the drivers of bighorn sheep population dynamics are population-specific enhanced monitoring programs. By maintaining a long-term, collared sample from each population, the observation process can be improved by providing more information on the locations of animals, and the biological model can be improved with direct estimates of vital rates and estimation of surveyspecific detection probabilities. The net effect would be improved precision on vital rate estimates, the ability to detect stochastic die-off events and to assess the impact of predation, an improved ability to assess the impacts of removals on population trajectories, and, as a result, a deeper understanding of populationspecific dynamics. Moreover, the detailed information would better be able to estimate the influence of density-dependent processes in population dynamics if/when these populations recover to a point where resources may become limiting. In the face of significant remaining uncertainty regarding the drivers of populations for bighorn sheep, we recommend such an enhanced monitoring program as a necessary tool for helping to develop management plans that will conserve and restore this iconic species.

ACKNOWLEDGMENTS
We are indebted to the regional biologists and game wardens who assisted in gathering and curating the historical data. We thank the many agency and university staff, capture crews, and volunteers that participated in all aspects of this bighorn sheep study and specifically to Kevin Monteith, Sarah Dewey, and Alyson Courtemanch for sharing survival data from collared animals from the Whiskey Mountain and Targhee populations. Primary funding for this work was provided by the Federal Aid in Wildlife Restoration Grant W-159-R to Montana Fish Wildlife and Parks and the annual auction sale of a Montana bighorn sheep hunting license, as well as the Wyoming Governor's Big Game License Coalition, the Wyoming Wild Sheep Foundation, and the Wyoming Game and Fish Department.