In semiarid ecosystems ungulate densities can compound the effects of drought on forage availability, which can feed back to regulate reproduction and offspring recruitment. Climatic changes in the frequency and severity of drought could thus exacerbate these effects. Here, we use long-term data sets to examine the dual influences of population density, cohort, and climatic variation on recruitment in the Henry Mountains bison (Bison bison bison) population. We found that 1-year lagged annual precipitation had a positive effect on recruitment (β = 0.032, SE = 0.009) and population density had a negative effect (β = −0.0023, SE = 0.0004), but did not detect cohort effects. Furthermore, the effect of population density seemed to be more severe in dry relative to wet conditions, indicating that competition for forage could be severe in drought conditions. These results could help managers gain insight into the effects of climate change on bison population dynamics, and help guide the management of bison abundance on rangelands shared with cattle.
Understanding the effects of climate change on vertebrate populations is currently of critical importance in conservation biology, wildlife management, and numerous other fields. These issues have thus been central in studies of large ungulates, which often play a key role in ecosystem dynamics (e.g., Berteaux and Stenseth 2006; Knapp et al. 1999). Density-dependent mechanisms often interact with density-independent climatic fluctuations to regulate the population dynamics of large herbivores. For example, deep snowpack or cold and wet conditions at northern latitudes can exacerbate the effects of competition for available food and negatively affect vital rates (Coulson et al. 2001; Creel and Creel 2009; Forchhammer et al. 2001; Garrott et al. 2003; Sauer and Severinghaus 1983; Solberg et al. 2001). Positive effects of summer precipitation on ungulate survival and reproduction can nevertheless counteract the deleterious effects of winter snowpack (Johnson et al. 2010). In arid and semiarid environments, precipitation often limits plant growth and productivity (Noy-Meier 1973; Rosenzweig 1968), thereby driving the strength of herbivore competition and demographic performance (Bender and Weisenberger 2005; Douglas 2001; McKinney et al. 2008; Owen-Smith and Ogutu 2003).
Reintroduced populations of plains bison (Bison bison bison; hereafter “bison”) in arid and semiarid environments offer unique opportunities to investigate the dual influences of climate and density dependence on demography. Between 1941 and 1942, 23 bison were transported from Yellowstone National Park (YNP) to the San Rafael Desert of south-central Utah. Forage in the red sand desert was scarce and the bison soon found their own way to the Henry Mountains (HM) where they have remained ever since (Popov and Low 1950). The herd grew rapidly and an annual hunt, initiated in 1960, provided a new opportunity for Utah sportsmen. The population experienced an all-time high in 2007 and abundance now appears to be approaching a stochastic equilibrium (which could be influenced by the forces of density dependence, harvest, and translocation tempered by environmental variability; Fig. 1).
Fig. 1.
Annual Utah Division of Wildlife Resources counts of total abundance of plains bison (Bison bison bison) in early August on the Henry Mountains and surrounding plateaus in southern Utah (left axis, closed circles with dashed line), and annual precipitation (cm) at the nearby Hanksville airport (right axis, open circles with dotted line).
Understanding the drivers of HM bison demography is important because 3 prominent management issues revolve around their abundance: hunting opportunities; removal for translocation within a metapopulation management program; and competition for available forage with domestic cattle. Because it is a free-ranging, genetically pure herd (Halbert and Derr 2007), the HM bison provide a valuable source for establishing other conservation herds while also offering a unique sport-hunting experience. Conversely, the local ranching community, which holds leases on public land grazing allotments in the HM, perceives bison as competitors with domestic cattle for available forage. Accordingly, resolving disagreement between ranchers, sportsmen, and natural resource agencies has become a focal issue in the management of bison and other herbivores in the southwestern United States (Garrott and White 1993). On the basis of input from multiple stakeholder groups, the current objective is to maintain 315 adults in the population after the annual hunting season (in 2012 the objective will move up to 325), but a better understanding of demographic mechanisms is needed to consistently meet this management objective (Utah Division of Wildlife Resources [UDWR]).
Hunting, density-dependent mechanisms, and climate may all play a role in regulating bison population dynamics (Fuller et al. 2007a, 2007b). Similar to cattle, climate can have a strong effect on bison forage quality and weight gain (Craine et al. 2009). Competition for available forage could thus exacerbate conflict among stakeholders of the range if the southwestern United States were to become drier (e.g., Seager et al. 2007). Yet little is known about the interactive effects of climate and population densities on bison demography (but see Fuller et al. 2007b for a study in YNP). Research on environmental drivers of bison population dynamics is nevertheless needed to develop adaptive management protocols in changing environments (Walters 1986).
Here, we use long-term data sets to examine the influence of population density and climatic variables on the demography of HM bison while controlling for possible cohort effects. We focus our attention on offspring recruitment because we expect it to be more responsive than adult survival to density-dependent biotic interactions (Bonenfant et al. 2009; Eberhardt 2002) and abiotic climatic variability (Gaillard et al. 1998). Because precipitation can modulate per capita forage availability (Holmgren et al. 2006), we hypothesize that population density has its most severe effects on recruitment following dry periods (i.e., an interactive effect between local precipitation and population density). Moreover, the dual influences of population density and precipitation on recruitment should be most severe when the population is large as opposed to the phase of rapid population growth from low abundance after introduction (Picton 1984). Our aim is to provide empirical results for guiding the multifaceted management of wild bison in arid to semiarid ecosystems amidst changing environmental conditions.
Materials and Methods
Study area and data collection
The HM in south-central Utah provide arid, semiarid, and alpine habitats for bison during their seasonal migrations from low to high altitudes. Apart from bison, cattle are the most common grazers in the region. A small herd of elk (Cervus canadensis) uses the HM, but given the low abundance (20–30; UDWR) they likely play a minor role in the herbivore community. Mule deer (Odocoileus hemionus) are common in the HM but their preference for forbs would suggest negligible levels of competition with the grazers (van Vuren and Bray 1983). Management efforts are actually aimed at keeping elk rare on the HM to favor the other herbivores (UDWR). Mountain lions (Puma concolor) and coyotes (Canis latrans) utilize the study area, but likely focus their foraging efforts on deer rather than risking the injuries that adult bison could impose on such predators (Lott 1991). A detailed description of the study area can be found in Nelson (1965) and van Vuren and Bray (1986).
Since 1949, the UDWR conducted surveys to count the number of bison at the end of July or in early August each year (Nt). Historically these surveys involved 12 or more personnel on foot, horseback, or in a vehicle counting bison in specified zones across the study area, but since 1990 an aerial helicopter survey was used (R. Hodson, pers. comm.). Since the 1941 and 1942 introductions, the bison herd has moved and used different parts of the HM in response to forage availability, hunting, and the 1962 brucellosis inoculation (Nelson 1965; van Vuren and Bray 1986; W. Paskett, pers. obs.). It was impossible to formally quantify population density per se, and, as in many studies, we assumed that the population count (Nt) was a good index.
At approximately the same time of year, the UDWR conducted ground surveys near roadsides to estimate calf-to-cow (hereafter, calf ∶ cow) and bull ∶ cow ratios. Calves, bulls, and cows were distinguished on the basis of their color, body size, horn shape, and conformation (e.g., lateral projection in calves; Fuller 1959; van Vuren and Bray 1986). As calving generally took place in May (Nelson 1965), most calves were 2–3 months old at the time of annual surveys, although it was not uncommon to see an orange calf <1 month old. Given the timing of the ground survey, the calf ∶ cow ratio encompassed rates of pregnancy, abortion, and early-life survival, and thereby provided a holistic measure of recruitment (i.e., the production and survival of bison offspring to August 1 of a given year t; van Vuren and Bray 1986). The population composition surveys did not begin after the initial introduction, but only took place between 1960 and 1964, and then from 1978 to 2010.
The total number of bison harvested and removed each autumn and early winter (Ht) was known because reporting of harvest (primarily adults) and hunting success to the UDWR was mandatory, and the number of immature animals removed and translocated was recorded directly. Surveyed abundance in subsequent summers (Nt+1) naturally included the subtractive process of Ht that occurred between surveys.
Precipitation data were recorded from 1946 to 2010 at the nearby Hanksville weather station (38°22′N, 110°43′W; 35 km N of the study area: X¯ = 14.08 cm of precipitation annually; data managed by the Western Regional Climate Center). This was the closest weather station with a continuous time series of climate data that overlapped the duration of the calf ∶ cow ratio surveys (Fig. 1). Although Hanksville was lower in elevation (1,300 m) and drier than the HM (which range up to 3,512 m), the temporal variability in precipitation likely reflected that which occurred on the HM. The upper slopes of the mountains received ∼3 times more precipitation than at the base, but from 1965 to 1973, precipitation on the upper slopes was highly correlated with the Hanksville data (R2 = 0.91; van Vuren and Bray 1986).
Explanatory covariates
We retrospectively examined the influence of population density, birth cohort, and climate variables on HM bison calf ∶ cow ratios. In this effort, we considered 2 different parameterizations of population density. We began by examining the influence of population density in August at t − 1 (Nt−1) on the observed calf ∶ cow ratio in August of year t (a bison survey year went from August through July). We also considered the effect of postharvest abundance in January (Npost,t−1 ≈ Nt−1 − Ht−1) because inter- and intraspecific competition could be at their highest during stressful winter conditions. Given that bison have high life expectancy (Peterson et al. 1991) and delayed age at maturity (2+, with primiparity being most frequent at age 3; Wolfe et al. 1999), we also considered lagged effects of each measure of population density at t − 2 and t − 3 on calf ∶ cow ratios in year t. However, we did not consider additive effects of abundance across seasons or years (e.g., Nt–1 + Nt−2) because of strong multicolinearity.
In species with delayed maturity, calf ∶ cow ratios could also exhibit strong cohort effects. For example, a bumper crop of calves born in year t could have led to a surplus of immature females that were counted as cows in years t + 1 and t + 2, which would have reduced the corresponding calf ∶ cow ratios. Similarly, production of a proportionally large calf cohort in year t could have increased intraspecific competition. Undernourishment during development could have caused such cohorts to be in poor physical condition upon primiparity, leading to reduced recruitment and a low calf ∶ cow ratio in year t + 3 (Descamps et al. 2008; Hamel et al. 2009; Pettorelli et al. 2002). To account for these possible cohort effects on the HM calf ∶ cow ratios, we considered autoregressive terms up to order 3 (i.e., the effects of calf ∶ cow ratios at t − 1, t − 2, and t − 3 on the calf ∶ cow ratio at time t).
Because of its controlling influence on primary productivity, precipitation is expected to be the dominant climatic driver of herbivore performance in arid and semiarid environments (Holmgren et al. 2006). We considered 7 precipitation variables that could have affected HM bison recruitment: total annual precipitation (cm) over the 12 months leading up to the calf ∶ cow survey in August of year t (annualt), 2 × 6-month seasons (August–January, including conception and early gestation; February–July, including late gestation, parturition and lactation), and 4 × 3-month seasons. We did not include nested precipitation variables in the same model (e.g., annualt + August–Januaryt) because of strong multicolinearity. However, correlations between precipitation variables in different years were low (all correlation coefficients <0.2), so we allowed immediate and lagged effects (up to 2-year lag) of a given precipitation variable to enter the same model (e.g., annualt + annualt−1). In preliminary analyses we also considered United States Geologic Survey streamflow data from rivers to the north and east of the HM (perhaps capturing the runoff of melting snow); however, Hanksville precipitation variables were much better predictors of calf ∶ cow ratios, likely because they were better indicators of local precipitation than stream runoff that included precipitation from distant sources.
Statistical analysis
Because bison rarely produce twins (Meagher 1986), calf ∶ cow ratios represented probabilistic reproductive outcomes that ranged between 0 and 1. We therefore examined the influence of the covariates described above on calf ∶ cow ratios using generalized linear mixed models with a binomial distribution and logit link function (McCulloch 2003). We computed the annual sample size of cows by correcting the aerial survey of total abundance with the ground surveys of calf ∶ cow and bull ∶ cow ratios (n = Nt/[1+ calf ∶ cow + bull ∶ cow]). Before analysis, we standardized all explanatory variables to X¯ = 0, σ = 1.
Rather than examining the hundreds of possible covariate combinations, we used a tiered approach to identifying the best explanatory model given the available data (Franklin et al. 2000). We began by comparing models with immediate (Nt−1) or lagged effects of population density (Nt−x), postharvest density analogues (Npost, each model containing a single explanatory covariate), and a null model using Akaike's information criterion adjusted for sample size (AICc—Akaike 1973; Burnham and Anderson 2002). Taking a pluralistic approach to scientific inference, we also assessed statistical significance of the individual parameters in models ranked best by AICc (Scheiner 2004). Using the same methodology, we then assessed possible cohort effects by comparing null, 1st-, 2nd-, and 3rd-order autoregressive models with one another (described above). Next, we compared models with immediate and/or lagged effects of the various precipitation covariates described above.
Last, we compared all possible combinations, including additive and interactive effects, of the population density, cohort, and climate covariates that performed best in the preceding analyses. In addition to these fixed effects, we included a random year-effect in each candidate model to account for temporal stochasticity caused by variables other than the measured covariates (using the “glmmML” procedure in R 2.13.0). Using AICc to compare the final candidate models, we were able to assess the influence of population density, cohort, climate, interactive effects, and temporal stochasticity on recruitment of the HM bison herd.
Results
Using 38 years of available UDWR survey data (1960–1964 and 1978–2010), we found that HM bison calf ∶ cow ratios were generally lower and less variable (X¯ = 0.37, σ = 0.10, range: 0.17 to 0.56) over the course of our study than observations collected by van Vuren and Bray (1986) during the late 1970s and early 1980s when precipitation was above average (X¯ = 0.52, σ = 0.19, range: 0.10 to 0.65). We also found that calf ∶ cow ratios were density dependent over the long term, and best explained by abundance in the previous survey (
= −0.0022, P < 0.01) as opposed to abundance after the hunting season or lagged abundance (all had higher AICc values). Without accounting for other variables, the calf ∶ cow ratios exhibited positive 2nd-order cohort effects (
= 0.65, 
= 0.93, both P < 0.05) as opposed to the predicted negative effects. Furthermore, additive effects on calf ∶ cow ratios of precipitation in the 2 years leading up to the survey (
= 0.015, 
= 0.038, both P < 0.01) were better supported than effects of precipitation in specific seasons (models with seasonal precipitation had higher AICc values).
After using the top-ranked covariates from above to develop a final list of candidate models, we found that population density and precipitation simultaneously affected recruitment of HM bison (Table 1; Fig. 2). Abundance in the previous survey had a strong negative effect on calf ∶ cow ratios (
= −0.0023, SE = 0.0004, P < 0.01), and lagged annual precipitation between t − 2 and t − 1 had a positive effect (Fig. 2; 
= 0.032, SE = 0.009, P < 0.01). In addition, a significant amount of random variation existed around the predicted relationship (σ = 0.21, SE = 0.04), suggesting that unmeasured variables also affected fluctuations in recruitment. Inclusion of an interaction between Nt−1 and annual precipitation between t − 2 and t − 1 indicated that the effect of population density was more severe in dry relative to wet conditions (Fig. 3). However, the interaction term increased AICc (Table 1) and was not significant (P = 0.15). All other models in the final analysis, including those with effects of cohort and nonlagged annual precipitation (annualt), were less competitive (ΔAICc > 2, all P > 0.05; Table 1).
Fig. 2.
Predicted recruitment (calf ∶ cow ratio) of plains bison (Bison bison bison) in the Henry Mountains of southern Utah (1960–1964 and 1978–2010) across the observed range of lagged annual precipitation between t − 2 and t − 1 (in cm) and total abundance at t − 1 (top-ranked model in Table 1).
Fig. 3.
Predicted recruitment (calf ∶ cow ratio) of plains bison (Bison bison bison) in the Henry Mountains of southern Utah (1960–1964 and 1978–2010) across the observed range of interacting lagged annual precipitation between t − 2 and t − 1 (in cm) and total abundance at t − 1 (2nd-ranked model in Table 1).
Table 1.
The top 10 generalized linear mixed models (GLMMs) for plains bison (Bison bison bison) calf ∶ cow ratios in the Henry Mountains of southern Utah (1960–1964 and 1978–2010) as a function of total abundance in the previous survey (Nt−1), a 2nd-order cohort effect (AR2, which contains 2 parameters), total annual precipitation between t − 1 and t (annualt), or between t − 2 and t − 1 (annualt−1). Evidence for a model is inversely related to its change in Akaike's information critierion adjusted for sample size (ΔAICc), with ΔAICc = 0 having the best fit among considered models. All models contained an intercept term and a random year effect; K represents the total number of estimated parameters in a model.
Discussion
Our findings for bison in the HM add to a growing body of evidence indicating that large herbivore populations are not driven by either density-independent or density-dependent mechanisms, but rather reflect both factors (Sæther 1997; Sibly and Hone 2002). Consistent with our hypothesis, we found that lagged annual precipitation (between August of year t − 2 through July of t − 1) and population density at t − 1 collectively influence the calf ∶ cow ratio in year t. Similar to ungulates in Montana (Picton 1984) and mule deer on the Kaibab Plateau (Elton 1927; Leopold 1943; McCulloch and Smith 1991), we moreover found that the effect of drought on recruitment was much stronger during phases of high relative to low population density (Fig. 3) Yet, statistical support for this interactive effect was not significant (P = 0.15), and observation over a longer time period with high population densities may be needed to improve precision in this estimate.
On the tallgrass prairie, the timing of precipitation affected the nutritional quality of forage and subsequent rates of growth in bison calves. For example, late-summer precipitation increased forage quality and calf growth rates. Conversely, midsummer precipitation increased stem biomass, reduced forage nutrients, and consequently reduced calf growth (Craine et al. 2009). A historical short-term study on the HM found that precipitation in winter and spring had a positive effect on calf ∶ cow ratios, presumably through its effects on forage quality (van Vuren and Bray 1986). Although we found that precipitation was an important driver of long-term variation in bison calf ∶ cow ratios, we did not find the seasonal timing of precipitation to be especially important. Nor did we detect cohort effects after accounting for population density and precipitation. Ongoing vegetation research will allow us to more closely examine whether or not seasonal changes in precipitation affect nutritive quality, and the abundance and diversity of forage throughout the year. Biologists could then assess how precipitation affects the development of a calf cohort, which could later affect their age at maturity and ensuing survival and reproductive success (e.g., in mammals see Descamps et al. 2008; Hamel et al. 2009; Pettorelli et al. 2002; e.g., in birds see Nevoux et al. 2010).
Calf ∶ cow ratios of bison in YNP were also affected by population density and climate, primarily snowpack, whereas the survival of prime-aged adults was robust to variation in both parameters (Fuller et al. 2007b). Even on our own study area, survival estimated from a life table during the late 1970s and early 1980s was stable across an array of precipitation regimes, whereas calf ∶ cow ratios varied greatly (van Vuren and Bray 1986). This result should nevertheless be reassessed because life tables often yield biased estimates of survival for mobile and elusive animals (Williams et al. 2002). For now, our analysis of long-term data readily indicates that HM calf ∶ cow ratios fluctuate greatly with both population density and precipitation (Figs. 2 and 3).
Although population growth in long-lived herbivores is more sensitive to changes in adult survival than to equivalent changes in calf recruitment (as measured with elasticities), fluctuations in recruitment can nevertheless have important effects on population dynamics (Gaillard et al. 1998). In many organisms, the deleterious effects temporal variation can have on fitness may have selected against the lability of highly influential demographic parameters (e.g., adult survival; see Boyce et al. 2006; Koons et al. 2009; Pfister 1998). Consistent with this view, bison and other long-lived herbivores have evolved life-history strategies where adult survival is much more resistant to (i.e., buffered against) environmental variability than reproductive parameters (Fuller et al. 2007b; Gaillard et al. 1998, Gaillard and Yoccoz 2003). The evolutionary mechanisms that drive selection for “demographic buffering” may also help explain the Eberhardt pattern (2002), where survival probabilities of prime-aged adult herbivores are less likely to be affected by changes in population density than fecundity and offspring survival (Bonenfant et al. 2009).
Managers should remain cognizant of the large influence that culling and changes in survival can have on bison population dynamics (Buhnerkempe et al. 2011; Fuller et al. 2007b), but the strong response of bison recruitment to changes in population density and precipitation might contribute more to population growth on an annual basis than subtle changes in natural mortality of adults (Gaillard and Yoccoz 2003). Changes in precipitation regimes could thus have strong effects on the management of bison populations in arid and semiarid environments (Holmgren et al. 2006), especially in areas like the HM where cattle and bison share forage resources (van Vuren and Bray 1983). Recent climate projections for the southwestern United States predict a drying trend in the future, but suggest that such a trend has yet to be realized (Cayan et al. 2010; Seager et al. 2007; Seager and Vecchi 2010). Over the past 65 years, we did not detect a trend in precipitation at the Hanksville weather station close to our study area (P = 0.33 and R2 = 0.02 for a trend regression). Precipitation is nevertheless highly variable (Fig. 1, coefficient of variation = 0.35), and in addition to other sources of environmental variability (e.g., cattle stocking rates), has an important influence on HM bison recruitment (Table 1; Fig. 2).
Combined with continued monitoring of population density, precipitation, and calf ∶ cow ratios, the functional relationships presented here could help management agencies adjust their actions to changing environmental conditions. For example, our top-performing model can be used in combination with already-collected monitoring data to predict whether there will be a bumper crop of calves, which is likely to occur following a wet year when population density is low, or a bust, which is likely to occur following a dry year when population density is high (Fig. 3). To meet bison abundance objectives, management agencies could thus use dynamically updated predictions of calf recruitment (i.e., calf ∶ cow ratio) to plan the timing and scale of hunting and translocation management efforts (similar to disease management strategies in African buffalo, Syncerus caffer, populations—du Toit 2010). For example, if bison numbers were above the management objective (325 posthunt adults in 2012), then the removal of animals timed to coincide with poor recruitment following droughts would be the most effective way to meet objectives. Conversely, if bison numbers were below objective during and following a time of drought and low recruitment, minimal harvest and removal might be advised. Updated information on age at maturity and adult survival will nevertheless be needed to develop empirically based population models for guiding sound management of the HM herd (e.g., in YNP see Fuller et al. 2007b), which must also be balanced with genetic viability of the Utah metapopulation.
Acknowledgments
We thank the Bureau of Land Management, United States Fish and Wildlife Service, and the Utah Division of Wildlife Resources for funding and the latter for the bison data.
