Elk population dynamics when carrying capacities vary within and among herds

Study areas

Time series of population survey data were used from nonmigratory elk populations in three different locations along the West Coast of the USA (Fig. 4). Five of the populations were in the Prairie Creek drainage (Davison), the Lower Redwood Creek drainage (Levee Soc), the Stone Lagoon area, the Gold Bluffs region, and the Bald Hills region of Redwood National and State Parks (RNSP), Humboldt County, California (41.2132° N, 124.0046° W). These populations occupy an area of about 380 km². The climate in this region was mild, with cool summers and rainy winters. Annual precipitation was usually between 120 and 180 cm and most of the precipitation fell between October and April. Snow was rare since average winter temperatures rarely dropped below freezing and ranged from 3 to 5 °C. Average summer temperatures ranged from 10 to 27 °C, depending on the distance inland. Elk in RNSP were not legally hunted, and displayed strong social bonding between females, juveniles, and sub-adult males⁷.

An elk population in the Point Reyes National Seashore inhabited part of the Point Reyes Peninsula in Marin County, California (38.0723° N, 122.8817° W). The elk were restricted to an area of 10.52 km² on the northern tip of the peninsula by a 3-m-tall fence. The climate of this study area was Mediterranean, with an average annual precipitation of 87 cm²⁷. Most of the precipitation fell from autumn to early spring. Temperatures averaged about 7 °C in winter and 13 °C in summer^27,35.

Another elk population was in the Arid Lands Ecology (ALE) Reserve and occupied a 300 km² area within the U.S. Department of Energy’s Hanford Site, Washington (46.68778° N, 119.6292° W). The climate in this area was semi-arid with dry, hot summers and wet, moderately-cold winters. Average summer temperatures were around 20 °C and average winter temperatures were around 5 °C with an average annual precipitation of 16 cm, half of which fell in the winter as rain³⁶.

Population surveys

In RNSP, females, juveniles, and subadult males were often in the same group and tended to use open meadow habitat more frequently than adult males^37,38. These behavioral patterns likely explain why females, juveniles, and subadult males were sighted more frequently than adult males⁷. Moreover, in size-dimorphic ungulates such as elk, recruitment was strongly correlated with female abundance and weakly correlated with male abundance^7,13,39. In RNSP, the abundance of groups of females, juveniles, and subadult males drove the dynamics of the group and of adult males⁷. Therefore, for the RNSP populations, we used herd counts where a herd was comprised of females, juveniles, and subadult males. We also used herd counts for the Point Reyes and ALE Reserve populations to remain consistent.

Systematic herd surveys of elk were conducted during January from 1997 to 2019 in RNSP. Surveys in the Davison meadows, the Levee Soc area, the Stone Lagoon area, the Gold Bluffs region, and the Bald Hills region were conducted by driving specified routes 4 to 10 times on different days throughout the month of January. The time series for these five herds ranged from 19 to 23 years of data. The elk were counted and classified by age and sex as adult males, subadult males, females, and juveniles. Females could not be visually differentiated into adult and subadult age categories³⁷. The highest count of females, juveniles, and subadult males from the surveys conducted each year was used as an index of abundance of each herd since the detection probabilities were high both on an absolute basis (> 0.8) and relative to variation in detection probabilities (CV_sighting = 0.05)^7,40. For the Bald Hills herd, which is the only herd in RNSP where harvests occurred, we added hunter harvests to the highest count of each year to account for this source of mortality. These harvests occurred only when elk from the Bald Hills herd left RNSP.

Elk population surveys were conducted at the Point Reyes National Seashore from 1982 to 2008. Weekly surveys were conducted after the mating season. Surveys were conducted on foot or horseback of female elk that were ear-tagged or had a collar containing radio telemetry^32,35. Individuals counted were classified as females, juveniles, subadult males, and adult males. Data were not available for the years 1984 to 1989 and 1993, so the time series included 20 years of data. We used the highest count of females, juveniles, and subadult males in each year in our analyses. This herd was also not hunted.

Elk population surveys in the ALE Reserve were conducted in winters after hunting and before parturition. From 1982 to 2000, biologists used aerial telemetry studies, in which they located all collared elk during each survey and classified them by sex and age. We used the total counts of females, juveniles and subadult males. For years in which multiple surveys were conducted, we used the highest count in each year as an index of abundance for that year^25,41. We omitted population survey data collected in 1982 from our analysis because individuals were not classified by sex and age in this year. Consequently, the time series included 18 years of data. For all years of data used, we added hunter harvests to the highest count of each year to account for this source of mortality. The count in 2000 was much lower than in the previous year, likely due in part to a large wildfire which occurred in the summer of 2000, which probably had an immediate effect of reducing available elk forage in the reserve and caused elk to spend more time outside of the ALE Reserve^42,43. In addition, the highest recorded number of elk (about 291) were harvested that year⁴³.

Ricker growth models

We fit linearized Ricker growth models simultaneously to the seven time series to estimate population growth parameters as well as temporal variation in r and β. We estimated K as the x-intercept of the Ricker growth model (i.e., when r = 0). Notably, preliminary analyses showed that not accounting for observer error did not bias our results (see Supplementary Information).

We used a Bayesian Markov Chain Monte Carlo (MCMC) algorithm with 3 chains, 150,000 iterations, a burn-in period of 75,000, an adaptation period of 75,000, and no thinning. We used Bayesian inference and MCMC because these methods offer advantages when fitting hierarchical models to model variation in ecological data^44,45. We conducted these analyses in the RJAGS program (JAGS Version 4.0.0; https://sourceforge.net/projects/mcmc-jags/files/JAGS/4.x/Windows/) in RStudio (R Version 3.5.0; https://cran.r-project.org/bin/windows/base/old/3.5.0/). We used uninformative priors for the y-intercept (i.e., r_max) and the slope (i.e., β) in order to allow solely the data to influence posterior estimates of these parameters. Informative priors were not necessary as long as parameter estimates from each chain converged. Convergence among chains was determined when the Gelman-Rubin diagnostic ((hat{R})) was less than 1.01, and through visual checks of trace and density plots⁴⁶.

The estimate of r_max borrowed information among herds because this parameter should be similar among populations within a species²². Therefore, we modeled r_max for each herd (j) as a random effect following a normal distribution with (mu_{{r_{max} }} sim Normalleft( {0, 0.001} right)) and (sigma_{{r_{max} }} sim Uniformleft( {0, 100} right)). To model temporal variation in r for each herd, we included a zero-centered random effect which was also modeled following the normal distribution (gamma_{t,j} sim Normalleft( {0,sigma_{{gamma_{j} }} } right)), where (sigma_{{gamma_{j} }} sim Uniformleft( {0, 100} right)). The estimate of β did not borrow information among herds because this parameter can vary widely among herds¹⁸. The prior for β for each herd (j) followed the normal distribution (beta_{j} sim Normalleft( {0, 0.001} right)). To model temporal variation in β for each herd, we modified how we modeled β by using a normal distribution ({beta_{{delta }_{t,j}}} sim Normalleft( {mu_{{beta_{{delta }_{j} }}}} , {sigma_{{beta_{{delta }_{j}} }} } right)), where ({mu_{{beta_{{delta }_{j}} }}} sim Normalleft( {0, 0.001} right)) and ({sigma_{{beta_{{delta }_{j}} }}} sim Uniformleft( {0, 100} right)). Thus, there were four possible Ricker growth models for each herd; (1) no temporal variation in r and β,

(2) temporal variation in r,

(3) temporal variation in β,

and

(4) temporal variation in both r_max and β

The residual variance was modeled as (varepsilon sim Uniformleft( {0,100} right)). We fit the model with no temporal variation (Eq. (3)) in either parameter to all seven time series simultaneously. All parameters except for r_max were estimated independently for each herd. For each time series of population survey data, we determined whether models with more parameters provided a better fit. We did so by fitting each possible growth model (Eqs. (4)–(6)) to each time series one at a time, while modeling all other time series with no temporal variation in r_max or β (Eq. (3)). The model with the lowest mean deviance from RJAGS by more than 2 was selected for that herd⁴⁷.

Environmental and demographic stochasticity

We estimated fluctuation in abundance which can be attributed to demographic and environmental stochasticity for herds with different K for each herd. The stochasticity model was outlined by Ferguson and Ponciano⁹;

where (Varleft( {N_{t – 1} } right)) was total population stochasticity, (Var_{dem} left( {N_{t – 1,} } right)) was population abundance fluctuation due to demographic stochasticity, (Var_{r} left( {N_{t – 1} } right)) was population abundance fluctuation due to changes in r (i.e., density-independent environmental stochasticity), and (Var_{beta } left( {N_{t – 1} } right)) was population abundance fluctuation due to changes in β. The model assumes density-dependent survival following the Ricker model. Demographic stochasticity was calculated as follows;

where (sigma_{dem}^{2}) was assumed to be equal to α⁹. Environmental stochasticity that is expressed as changes in r, otherwise known as density-independent or additive stochasticity, was calculated as follows;

and environmental stochasticity that is expressed as changes in β, otherwise known as density-dependent or multiplicative stochasticity, was calculated as follows;

Population growth parameters from the selected Ricker growth model for each herd were used in these equations to estimate each of these sources of stochasticity for each herd across abundances ranging from five to above K. The relative total population stochasticity was expressed as the total population stochasticity at K for each herd divided by that herd’s K.

Source: Ecology - nature.com