Autumn ungulate hunting in the Greater Yellowstone Ecosystem carries the risk of hunter–grizzly bear (Ursus arctos) conflict and creates a substantial challenge for managers. For Grand Teton National Park, Wyoming, USA, a key information need is whether increased availability of elk (Cervus canadensis) carcasses during a late autumn (Nov–Dec) harvest within the national park attracts grizzly bears and increases the potential for conflict with hunters. Using a robust design analysis with 6 primary sampling periods during 2014–2015, we tested the hypothesis that the elk harvest resulted in temporary movements of grizzly bears into the hunt areas, thus increasing bear numbers. We detected 31 unique individuals (6 F, 25 M) through genetic sampling and retained 26 encounter histories for analysis. Markovian movement models had more support than a null model of no temporary movement. Contrary to our research hypothesis, temporary movements into the study area occurred between the July–August (no hunt; N̄2014–2015 = 5) and September–October (no hunt; N̄2014–2015 = 24) primary periods each year, rather than during the transition from September–October (no hunt) to November–December (hunt; N̄2014–2015 = 15). A post hoc analysis indicated that September–October population estimates were biased high by detections of transient bears. Grizzly bear presence during the elk hunt was limited to approximately 15 resident bears that specialized in accessing elk carcasses. The late timing of the elk hunt likely moderated the effect of carcasses as a food attractant because it coincides with the onset of hibernation. From a population response perspective, the current timing of the elk harvest likely represents a scenario of low relative risk of hunter–bear conflicts. The risk of hunter–grizzly bear encounters remains, but may be more a function of factors that operate at the level of individual bears and hunters, such as hunter movements and bear responses to olfactory cues.
Grizzly bear (Ursus arctos) population recovery in the Greater Yellowstone Ecosystem has coincided with increases in both human populations on the periphery of the ecosystem and human visitation to public lands (Gude et al. 2006, Lawson 2014; U.S. Forest Service national visitor use monitoring statistics, https://apps.fs.usda.gov/nvum/results; National Park Service visitor use statistics, https://irma.nps.gov/Stats/). These combined trends are contributing factors to the increasing number of human–bear encounters and conflicts (Gunther et al. 2012). Whereas many direct encounters are chance events during non-hunting recreational activities, autumn elk (Cervus canadensis) hunting increases the probability of human–bear conflict because grizzly bears exploit resources provided by offal and other remains left behind by successful hunters and wounding losses (Haroldson et al. 2004). A successful elk harvest typically results in 30–40 kg (14% of body weight, excluding rumen) of edible offal left at kill locations (Wilmers et al. 2003a). Additional biomass is available from wounding losses, which may be 10–30% of harvested elk (Unsworth et al. 1993). Whereas grizzly bears generally avoid humans, even to the extent of foregoing foraging opportunities (Coleman 2012), the increase in elk carcasses and gut piles (hereafter, carcasses) during autumn ungulate-hunting seasons can create temporally and spatially predictable resources on the landscape at a time when caloric demand and intake is greatest (autumn hyperphagia). In effect, these areas may serve as a seasonal attractant to bears, particularly considering the ability of grizzly bears to locate and use carcasses (Cristescu et al. 2015, Ebinger et al. 2016). Haroldson et al. (2004) found that grizzly bears were 2.4–2.7 and 2.3–4.4 times more likely to be outside Yellowstone National Park's northern and southern boundaries, respectively, following the opening of the September elk season, thus increasing the risk of hunter–bear conflicts and grizzly bear mortality. Gunther et al. (2004) found that grizzly bears killed in self-defense by ungulate hunters (n = 27) represented 36.5% of all human-caused mortality during the period 1992–2000. This proportion declined but remained high at 27.5% in the following 15 years (2001–2015; Interagency Grizzly Bear Study Team, unpublished data).
With a management objective of 11,000 elk, the Jackson elk herd is large and a significant portion travels though Grand Teton National Park during annual autumn migrations to wintering areas on the National Elk Refuge and 3 state feeding grounds (U.S. Department of the Interior 2007). Supplemental winter feeding, which began at the refuge during the early 1900s, ensures strong site fidelity, high survival, and little winter loss of elk; thus, harvest management is used to help meet state objectives for the elk herd (U.S. Department of the Interior 2007). The national park is immediately adjacent to the National Elk Refuge (Fig. 1) and is used by elk for both summer and transitional range during the autumn migration. Under its 1950 establishing legislation, Grand Teton National Park is authorized to conduct a joint Elk Reduction Program with the State of Wyoming when necessary for conservation of the Jackson elk herd (U.S. Department of the Interior 2007). Using private citizens deputized as rangers, elk hunts have occurred in Grand Teton National Park every year since 1950 except two (1959–1960). Over the period 2008–2017, an average of 242 (range = 132–361) elk were harvested annually, representing approximately 20% (range = 14–25%) of total harvest for the Jackson herd (S. Dewey, Grand Teton National Park, unpublished data). The elk hunt is usually held during late October–early December and <20% of the national park is open to hunting (Fig. 1). Most elk are harvested as they migrate to the National Elk Refuge south of the national park, but there is a marked peak in the timing of harvest, with approximately 90% of elk harvested after the first week of November (S. Dewey, unpublished data). This pulse of spatially clustered elk carrion is fundamentally different from the scattered and scarce number of elk carcasses provided by gray wolves (Canis lupus) and other predators throughout the year (Wilmers et al. 2003b). The elk hunting season in Grand Teton National Park is typically open later (Nov–Dec) than ungulate hunting season on adjacent non-park lands, so there is potential for elk carcasses to attract grizzly bears to the park and increase hunter–grizzly bear conflicts. Several high-profile hunter–grizzly bear conflicts associated with the Elk Reduction Program have occurred, including the mauling of an elk hunter in 2011 and the death of a grizzly bear in 2012, which received substantial attention in local, regional, and national media.
Fig. 1.
Study area to determine population responses of grizzly bears (Ursus arctos) to elk (Cervus canadensis) hunting as part of the Elk Reduction Program in Grand Teton National Park, Wyoming, USA, 2014–2015. Wyoming Hunt Areas that are open during the Elk Reduction Program are shown, with the majority of harvest in Hunt Area 75. Elk autumn migration routes start in the northern portions of the map and converge along the north–south axis of the DNA sampling grid before entering the National Elk Refuge (U.S. Department of the Interior 2007).
Autumn elk hunting, in conjunction with increasing grizzly bear numbers, creates a unique and substantial challenge for national park managers. Several provisions for mitigating hunter–grizzly bear conflicts in Grand Teton National Park are in place, including requiring hunters to carry bear spray, prohibiting artificial elk calls, and providing hunters with a bear safety education packet. Additional measures may increase the safety of hunters, but little scientific information exists upon which to base such decisions. Therefore, park managers are seeking new, science-based information to help reduce conflict potential. A key information need is whether the autumn elk harvest attracts grizzly bears into the areas open for hunting. We tested the hypothesis that increased availability of elk carcasses during the elk harvest resulted in temporary immigration of grizzly bears, potentially increasing the probability of conflict between grizzly bears and elk hunters.
Study area
Our study area was located in the eastern portion of Grand Teton National Park and adjacent lands of the Bridger–Teton National Forest, including all of Wyoming Hunt Area 75 (120 km2) and portions of Hunt Areas 79 (53 km2) and 81 (34 km2; Fig. 1). The study area was within a high-elevation valley that included the upper Snake River drainage, bounded by the Teton Range to the west, the Gros Ventre and Absaroka mountains to the east, the Yellowstone Plateau to the north, and the town of Jackson, Wyoming, to the south. Elevations range from 1,890 m in the valley to 4,197 m atop surrounding peaks. Climate was characterized by long, cold, and snowy winters and short, cool summers. Approximately 70% of precipitation typically falls as snow.
Our genetic sampling was restricted to a 500-km2 grid (Fig. 1), which was bounded by the Snake River floodplain on the west, consisting of riparian forest (e.g., cottonwood [Populus spp.], willow [Salix spp.], and aspen [P. tremuloides]). The remainder of the sampling area consisted of terraces rising above the floodplain covered by sagebrush (Artemisia spp.) and grasses, and was occasionally interrupted by glacial moraines and forested buttes consisting of lodgepole pine (Pinus contorta), Douglas-fir (Pseudotsuga menziesii), and aspen (Jean et al. 2005). The highest elevation in the sampling area was Blacktail Butte (2,343 m). The sampling area was relatively small compared with estimates of annual home ranges of female (130 km2) and male (475 km2) grizzly bears in the ecosystem (Bjornlie et al. 2014), but focused on the area where paths of migratory elk moving to the National Elk Refuge converged (Fig. 1).
Methods
General approach
Genetic sampling can provide accurate and precise estimates of population size and have been widely applied in studies of bears and other large vertebrates (e.g., Woods et al. 1999, Mulders et al. 2007, Kendall et al. 2009). We conducted genetic sampling in 2014 and 2015 using DNA extracted from the roots of hair samples. Grizzly bear detections based on genotypes served as marks for our capture–recapture analysis. Given the large home ranges of grizzly bears in relation to the study area, temporary movements in and out of the study area are an important consideration and a focus of our study. These temporary movements violate assumptions of closed mark–recapture models as well as those of open models, which assume that emigration or immigration is permanent (Jolly 1965, Seber 1965, Otis et al. 1978). However, enhancements to the robust design framework of Pollock (1982) provide a mechanism to estimate temporary emigration and immigration, along with other parameters of interest that can help inform management regarding grizzly bear activity relative to the Elk Reduction Program (Kendall 1999, Kendall and Bjorkland 2001, Lind-berg and Rexstad 2002). Our research hypothesis was that grizzly bears temporarily immigrated into areas open to elk hunting during the park's Elk Reduction Program. With a robust design approach, we combined open and closed population models using 6 primary sampling periods across 2 years (2014–2015), within each of which we sampled 4 secondary periods. In each year, the 3 primary sampling periods spanned the approximate dates of 10 July–15 August for the first period, 5 September–15 October for the second period, and 1 November–10 December for the third period, during which the elk hunt was held. Although the elk hunt typically starts in late October, our third primary sampling period was timed to coincide with the peak in elk harvest previously mentioned.
Field sampling design
Active sampling. Previous studies showed that a spatial sampling unit of approximately 25 km2 satisfies the assumption that individual grizzly bears within the grid cell are exposed to olfactory cues of lures (Kendall and McKelvey 2008). We set up barbed-wire hair corrals distributed throughout the study area using a 5- × 5-km grid of 20 cells, with 1 corral trap/grid cell (Fig. 1). Each corral trap was paired with a motion-sensitive infrared camera to identify if and when individual bears left multiple hair samples and to document sampling of family groups, as subsequently discussed. The distribution of the corral traps within each of the 20 grid cells covering the study area followed sampling protocols outlined by Kendall and McKelvey (2008). We checked corral traps every 7–10 days and moved them between the primary sampling periods to reduce the probability of a trap response. We used cattle blood as a lure and refreshed lures at the onset of each secondary sampling period. For human and bear safety, all hair corrals were located >100 m from trails and roads and >500 m from human developments (Kendall and McKelvey 2008).
Passive sampling. Passive sampling of grizzly bear hair involves bear rubs on trees or human structures such as sign posts, cabins, power poles, fence posts, and gates (Kendall et al. 2009). We located and mapped rub sites using observations from National Park Service biologists and surveys of trails, roads, power lines, and other rub features (Kendall et al. 2009). To facilitate hair collection, we attached 4, 40-cm lengths of barbed wire to each rub object with 2–3 fencing staples each, using a “Z” pattern to cover the rubbing area. We tagged each rub object with a unique identification number nailed to its base, out of sight from the trail or road. We removed staples and nails once surveys were completed. We did not use a lure at these sites because of the habitual nature of bear visitation to rub sites. We also opportunistically sampled hair at elk carcasses. We collected passive samples in association with the field surveys for the hair sampling corrals, and thus linked samples to the same study area and primary and secondary sampling periods.
Sample collection and microsatellite analysis. We treated hair collected from individual barbs as separate samples. We sealed each hair sample in a paper coin envelope and stored it in a climate-controlled room upon field collection. We assigned unique identifiers to each hair sample and printed a bar code that was adhered to each sample envelope. Camera records indicated that hair on adjacent barbs were typically multiple samples from a single individual. Therefore, we subsampled by randomly selecting 1 sample when ≥2 neighboring barbs contained hair samples. Wildlife Genetics International (Nelson, British Columbia, Canada) conducted the genotyping and used 9 microsatellite markers (G10B, G1D, G10H, G10J, G10P, G10L, MU59, MU51, MU23; Paetkau et al. 1995, Taberlet at al. 1997) plus a sex marker (ZFX/ZFY; Durnin et al. 2007) to identify unique individuals. For new individuals that were not already in the genetic database of grizzly bears in the Greater Yellowstone Ecosystem, analyses were extended to 21 markers to allow additional genetic analyses (e.g., parentage). DNA was purified from the hair samples using QIAGEN DNeasy Blood and Tissue Kits, following the tissue protocol (Paetkau 2003). Genotyping was completed following a 3-phase process. After the first pass of the 10 markers, samples that failed to produce high-confidence genotypes (based on a combination of objective [peak height] and subjective [appearance] criteria) for ≥5 markers were removed. These samples were marked by removing the leading digit from the 3-digit allele score; 2-digit scores were treated as equivalent to missing data. Additionally, samples were removed that amplified at ≥3 alleles at ≥3 loci (mixed samples) and those that represented American black bear (U. americanus) genotypes, based on odd-numbered alleles for G10J. Data points that were weak or initially difficult to read were reanalyzed using 5 µL of DNA/reaction, instead of the 3 µL used on first pass. At the end of this process, samples with <10-locus genotypes were removed from further analysis. In the final phase, error-checking protocols involved evaluation of mismatching (MM) markers in 1-MM, 2-MM, and 3-MM pairs that fit the pattern expected of allelic dropout (Paetkau 2003) and effectively prevents the recognition of false individuals. Extensive blind testing by Kendall et al. (2009) indicated these protocols reduced genotyping error rates to inconsequential levels.
Population estimation techniques are based on the assumption that detections are independent among individuals. This assumption was reasonable for our study, except for family groups consisting of a female and her dependent young (<2 yr old). We had parentage information on approximately half of the individuals we detected (D. Paetkau, Wildlife Genetics International, unpublished data), and combined with frequent observations of family groups and telemetry data of bears captured in the study area, we were able to identify all genotypes that represented adult females with their dependent offspring during all or a portion of the study. We removed the encounter histories of known, dependent young.
Statistical analyses
Overview. We defined the study population (Ni) as those bears available for detection within the sampling area during a primary sampling period (i). The superpopulation (N0) represents not only the study population but the larger number of bears that ever enter the sampled population between the initial and last secondary periods and thus could be detected (Williams et al. 2011). The conditional probability of detecting a grizzly bear given that it is present (apparent detection probability p*i) and the probability that it is alive but outside the study area and unavailable for detection (γi*), together influence the probability that a bear from the superpopulation is detected during a particular sampling period. Kendall (1999) showed this relationship is the product (1 – γi*) × p*i , and referred to this as the effective detection probability.
We focused our inference on temporary movements, as measured by the availability for detection (observability) of bears in the study area. The robust design analysis allows estimation of 2 relevant parameters (Kendall and Nichols 1995, Kendall et al. 1997):
γi′ = the probability that a bear is outside the study area in primary period i, given that it was not present in the study area during primary period i – 1;
γi″ = the probability that a bear moves outside the study area in primary period i, given that it was present during primary period i – 1.
Parameterization of γi′ and γi″ can be used to estimate different types of movements. For our study, the term 1 – γi′ is of particular interest because it is the probability that an individual not on the study area during time i – 1 enters the sample between time i – 1 and time i (i.e., temporary immigration; Fig. 2). We focused our analyses on a Markovian movement model. Markovian movement refers to a scenario where the probability of an animal being in the sampling area and available for detection during primary period i is dependent on whether it was in or out of the sampling area at primary period i – 1. In other words, in the context of our study, bears would ‘remember’ they left the study area and return during the elk hunt because of the putative availability of a high-quality food resource (Kendall et al. 1997).
Fig. 2.
Schematic of the robust design for a single year and analysis results to test the hypothesis that grizzly bears (Ursus arctos) move into hunt areas of Grand Teton National Park, Wyoming, USA, during the autumn elk (Cervus canadensis) hunt, 2014–2015. (A) Circles show the study population nested within the superpopulation and mechanisms of temporary emigration under the Markovian movement model. Dashed lines indicate state (observable or unobservable) transitions and solid lines represent states remaining the same. Red lines indicate transitions for individuals that were within the study area (observable) at time i – 1 and black lines indicate transitions for individuals that were outside the study area (unobservable) at i – 1. Our primary parameter of interest was 1 – γ′ (black dashed line) for transition from the no-hunt to the hunt period. Gamma parameters apply to both years, because there was no evidence of a year effect. (B) Robust design framework with primary periods (PP) and secondary periods (SP), population closure assumptions (open vs. closed), and treatments (hunt vs. no hunt). Dates are approximate. Estimates of the movement parameters are shown next to the arrows for the top model in panel A (see Table 2).
Goodness-of-fit. No formal goodness-of-fit test exists for robust design analysis. Therefore, we considered an alternative by collapsing encounter histories across the primary sampling periods and assessing goodness-of-fit of a time-variant Cormack–Jolly–Seber model (Cormack 1964, Jolly 1965, Seber 1965). We used Program RELEASE (Burnham et al. 1987) to test the assumption that all marked animals in the population have 1) the same probability of being detected at any of the primary sampling periods (TEST 2 in Program RELEASE), and 2) have the same probability of surviving, regardless of the sampling period during which they were marked (TEST 3).
Robust design analysis. A key tenet of robust design is that sampling is integrated over 2 different time frames using primary sampling periods that are composed of multiple secondary periods (Pollock 1982). The time interval between secondary periods is short enough so that population closure can be assumed (i.e., no births, deaths, emigration, or immigration) and estimates of capture probability and population size can be obtained. By collapsing data from the secondary periods, estimates of apparent survival and temporary emigration are obtained across the primary periods with open models. Estimates from robust design analyses are generally more accurate and precise compared with those from separate closed or open models because they integrate the estimation of survival and abundance while accounting for temporary emigration (Kendall et al. 1995). In typical robust design applications, primary sampling periods are separated by a sufficient amount of time, such that the sampled population is expected to change through gains (birth and immigration) and losses (death and emigration; Clavel et al. 2008). In our study, we leveraged the high survival rate of grizzly bears and deviate from the typical robust design by selecting a relatively short interval between primary periods, thus focusing our inference on gains and losses to (temporary) immigration and emigration.
We used the genotypes of sampled bears to construct individual encounter histories for Program MARK using the Huggins closed-capture model with encounter parameters p (capture [detection] probability) and c (recapture probability) within the robust design framework (Huggins 1989, White and Burnham 1999). We used the closed population models to estimate capture (pij) and recapture (cij) probabilities and Ni, where i and j index primary and secondary sampling periods, respectively. Our original design for the 2-year study resulted in 24 encounter occasions (6 primary periods times 4 secondary periods each). However, a key assumption is that the population is closed to additions and deletions across secondary sampling periods, which may have been violated for the 4 secondary sampling periods associated with the last primary period each year by the fact of bears entering their dens (Haroldson et al. 2002). Kendall (1999) presented several scenarios where estimation of capture probabilities would still be unbiased, including one scenario in which the population is available for detection during the initial secondary sampling period, but starts leaving the study area, equivalent to entering a den, thereafter. This scenario may apply because our sampling for that primary period started in the first week of November, when almost all bears were still active. Based on telemetry data from 2014 to 2015 and using the methods of Haroldson et al. (2002), we estimated den entry dates for 17 radiocollared grizzly bears that frequented Grand Teton National Park and John D. Rockefeller, Jr. Memorial Parkway and found that 88% of those bears entered dens after the first week of November (Interagency Grizzly Bear Study Team, unpublished data). Consequently, we used the initial secondary period and collapsed the last 3 remaining secondary periods for the third primary period each year, resulting in 20 rather than 24 encounter occasions (i.e., 3 primary periods consisting of 4, 4, and 2 secondary periods, respectively, each year).
We followed a 2-stage approach to model development. We first identified which model structure provided the best fit for the encounter parameters p and c, followed by a model set to test our hypotheses regarding temporary emigration and Ni. For the first stage, we explored biologically plausible frameworks for the encounter parameters. We considered sex as an individual covariate to account for capture heterogeneity. We used the Markovian model (time-variant γi′ and γi″) as the basis for the stage-1 models because this provided the most parameterized formulation for the temporary movement parameters (MacKenzie et al. 2012, Pederson et al. 2012). Given that we excluded known encounter histories of dependent young, we assumed that genotypes from sampled hair primarily represented independent-age bears (≥2 yr old). The interval between primary sampling periods lasted only several weeks and mean annual survival rates are high (S = 0.95) once grizzly bears reach the age of 2 (Haroldson et al. 2006, Interagency Grizzly Bear Study Team 2012, van Manen et al. 2016); therefore, we fixed survival between primary periods to 1 (i.e., S = 1.0). The only exception to the short intervals was the period between the last primary sampling period in 2014 and the first in 2015, which included the denning period and early spring. Documented grizzly bear mortalities during the denning period are rare, so we assumed S = 1.0 for the denning months. However, survival is lower during April–June, so we used monthly survival estimates (S = 0.993) for those 3 months based on 2006–2015 data to set S between primary sampling periods 3 and 4 to 0.979. To determine the best model structure for the encounter parameters, we considered different combinations of p and c: 1) time-variant p and c across primary sampling periods, constant but different by study year, or constant across both years; 2) equal or different p and c probabilities to assess trap responses; and 3) p and c as a function of the individual covariate sex. We used Akaike's Information Criterion, corrected for small sample sizes (AICc), to examine support among competing models (Hurvich and Tsai 1989, Burnham and Anderson 2002).
For the second stage, we developed a simple model set to test our hypotheses regarding temporary emigration and Ni. Using the most supported model structure for the encounter parameters and individual covariates, we used different parameterizations for γ″ and γ′ and while fixing S as discussed previously. To test our research hypothesis, we modeled γi′ and γi″ as a function of a categorical covariate that designated the transition between primary periods as 1) no hunt to no hunt (i.e., first to second primary period each year), 2) no hunt to hunt (second to third primary period each year), and 3) a single period of hunt to no hunt transition between the last primary period in 2014 and the first period in 2015 (Fig. 2). If the availability of elk carcasses in the areas open to hunting during the Elk Reduction Program attracts grizzly bears from outside areas where ungulate hunting seasons have closed, we predicted greater support for the Markovian movement model compared with a null model of no movement. For the null model, we fixed parameter estimates for γi′ and γi″ to 1 and 0, respectively. In other words, animals that are available or unavailable for sampling remain in their respective states over all primary sampling periods. If the Markovian model had most support, we predicted that immigration (1 – γi′ ) would be greatest for the transition from the no hunt to hunt period. As an alternative model, we also examined the presence of a year effect on the gamma parameters associated with the hunt period.
Results
Microsatellite analysis resulted in successful grizzly bear genotypes for 209 of 274 samples for which DNA was purified; microsatellite analysis failed for 41 (15.0%) samples, whereas 20 (7.3%) were American black bear samples, and 4 (1.5%) were mixed samples. We identified 31 individuals (6 F, 25 M), with 111 detections of 22 individuals in 2014, and 98 detections of 21 individuals in 2015. Relative to a reference database of >1,000 grizzly bear genotypes sampled in the Greater Yellowstone Ecosystem from the mid-1980s through 2016, our sampling detected 11 new individuals. The total number of detections per individual averaged 7.7 and ranged from 1 to 32. The number of unique individuals detected at sampling sites ranged from 1 to 4. Seventy-seven of the 209 (37%) samples were collected from passive sampling, resulting in 21 of 87 detections (24%) among the 20 secondary sampling periods and 4 of 31 individuals (13%) that were not detected from hair corral samples. The proportion of passive samples was greater for males (62%) than females (47%).
Using auxiliary data (parentage, live captures, telemetry), we identified 4 female genotypes as adults and 2 females and 4 males as dependent (i.e., <2 yr old; cubs or yearlings) offspring. We excluded the entire encounter histories of 4 of these 6 dependent bears and the 2014 encounters for 2 bears that were dependent yearlings in 2014, but reached the age of independence (≥2 yr) in 2015. However, 1 of those 2 bears was not detected in 2015, so we conducted our analysis based on 26 encounter histories. Additionally, there was one instance in which we detected a dependent-age bear without a concurrent detection of its mother. We updated the detection history of the mother accordingly, based on visual observations confirming the family had not separated. The 26 individuals and family groups were detected 1–4 times during the 6 primary sampling periods and 1–11 times within the 20 secondary sampling periods (total detections = 74). Collapsing the last 3 secondary periods for the third primary period each year resulted in 72 detections across all secondary periods. During the 3 annual primary sampling periods, we detected 4, 15, and 4 unique individuals in 2014 and 3, 17, and 8 individuals in 2015, respectively (Table 1).
Table 1.
Summary of DNA-based detections of individual grizzly bears (Ursus arctos) for a robust design analysis of population size and temporary movements associated with the Elk Reduction Program in Grand Teton National Park, Wyoming, USA, 2014–2015.
Goodness-of-fit tests of the Cormack–Jolly–Seber model indicated adequate fit of the data (TEST 2: χ2 = 0.48, 3 df, P = 0.925; TEST 3: χ2 = 3.34, 4 df, P = 0.503). A model with sex as an individual covariate and in which p and c were equal and constant (p., sex = c., sex; AICc = 319.72) had the most support, but only slightly more than a model in which p and c were different and a function of year (pyr, cyr; ΔAICc = 0.31; Table S1 (ursu-30-01-04_s01.pdf)). All other models had less support (ΔAICc > 3.55). Sex is a known source of detection heterogeneity among bears, so we modeled the encounter parameters as p., sex = c., sex (model 1; Table S1 (ursu-30-01-04_s01.pdf)). Given that males were more likely than females to be sampled passively, we assumed that the sex covariate accounted for differences in detection probabilities due to sampling technique. Capture probabilities were 0.42 for females and 0.21 for males.
The model we used to test our research hypothesis had more support than the null model (ΔAICc = 10.06), suggesting transitions between unobservable and observable states between the primary periods (Table 2). Contrary to our hypothesis, however, estimates for γ′ and γ″ showed that bears from outside the study area were more likely to move into the study area from primary period 1 (no hunt) to period 2 (no hunt) each year (1 – γ2′ = 0.80), rather than during the period 2 (no hunt) to period 3 (hunt) transition (1 – γ3′ = 0.00; Fig. 2). Additionally, bears in the study area were likely to remain in the observable state regardless of the transition (1 – γ2″ = 0.69, 1 – γ3″ = 0.73). Derived estimates of population size based on the top model indicated population abundance was greatest during the second and fifth primary periods, but precision was poor (Fig. 3A). There was little evidence of a year effect (ΔAICc = 5.27; Table 2).
Table 2.
Model-selection results for robust design analysis of grizzly bear (Ursus arctos) population size and temporary movements associated with the Elk Reduction Program in Grand Teton National Park, Wyoming, USA, 2014–2015. We fixed survival (S) between primary periods 3 and 4 (denning period and spring) to 0.979 and for all other transitions between primary periods to 1. We modeled capture (c) and recapture (p) probabilities as equal and constant, with an individual covariate for sex to account for capture heterogeneity.
Fig. 3.
Estimated population abundance of grizzly bears (Ursus arctos) for 6 primary sampling periods in areas open to autumn elk (Cervus canadensis) hunting within Grand Teton National Park, Wyoming, USA, 2014–2015. Estimates of population abundance based on (A) robust design analysis, and (B) Clavel et al. (2008) approach, which accounted for transients in the population (dashed line), with Cormack–Jolly–Seber estimates (solid line) for comparison (no estimates can be derived for the first and last primary periods). Primary periods 1–3 were in 2014 and primary periods 4–6 in 2015. Error bars show the 95% confidence interval. For the Clavel et al. (2008) approach, we could not estimate confidence intervals for periods 3, 4, and 5 because several estimates for resident rates were 0 or 1 with standard errors of 0 or near 0.
Post hoc analysis
Examination of the encounter histories revealed potentially transient bears that may cause violation of the closure assumption within primary sampling periods. Therefore, we performed a post hoc analysis to investigate the potential effect of transiency on our results. Building on the transient model of Pradel et al. (1997), Hines et al. (2003) developed techniques to estimate the survival of residents that is not influenced by the existence of an unknown number of transients. Clavel et al. (2008) extended the approach to estimate population abundance while accounting for the transient rate, or τ.
We first classified each bear as a resident or unknown status, according to Hines et al. (2003). We defined residents as bears that were detected ≥2 times within the primary period in which they were first detected (i.e., at least approx. 1 week between detections); we classified remaining bears in the unknown group, which included residents and transients. We pooled detections over secondary periods to a single indicator for each primary period and estimated survival rates in Program MARK with the Cormack–Jolly–Seber routine, using residents and unknown bears as groups (Clavel et al. 2008). We estimated apparent survival of residents with the Cormack–Jolly–Seber model and apparent survival of the unknown group using Pradel et al.'s (1997) model, as described in Clavel et al. (2008). We estimated the resident rate among the unknown group (1 – τ) for each primary period based on the corresponding ratios of apparent survival for the unknown and resident groups.
There were 8 (1 F, 7 M) predefined residents among the 26 bears detected in our sample, with 1 initial detection in primary period 1, 4 in period 2, 1 in period 4, and 2 in period 5. We classified the remaining 18 bears in the unknown group. All bears in the unknown group were estimated to be residents for those detected in primary periods 1, 3, and 4. For bears in the unknown group detected in primary periods 2 and 5 (i.e., the period prior to the elk hunt), 64% and 100% were estimated to be transients, respectively. Applying the methods of Clavel et al. (2008), we estimated that resident abundance for periods 2 through 5 (abundance cannot be estimated for the first and last period because one of the estimated parameters is missing), averaged 15 and varied between 11 and 18 bears (Fig. 3B). Comparison with Cormack–Jolly–Seber estimates showed that abundance estimates were particularly affected by transients in periods 2 and 5 (Fig. 3B).
Discussion
Our findings suggest that temporary movement into the study area did occur, but primarily in the time period prior to the elk hunting season, rather than during the Elk Reduction Program. Results of our post hoc analysis suggest these temporary movements prior to the hunt were driven by detections of transient bears, which represented a large proportion of the abundance estimate. Once we accounted for transients and estimated the number of resident bears, there was little evidence of major shifts in abundance. These results indicate that bears using the study area during the hunt elk hunt are almost exclusively resident bears. Therefore, we conclude that the data did not support our research hypothesis.
Several factors likely contributed to this finding. First, the effect of elk carcasses as a food attractant may be moderated by the late timing of harvest associated with the Elk Reduction Program. Haroldson et al. (2002) found that 90% of female bears typically have entered their dens by the end of November, whereas 90% of males have entered their den by the second week of December. By the time elk carcasses have accumulated in significant numbers, only a small number of bears may remain active in the areas open to hunting. Given the lack of other food resources, these remaining bears specialize on elk carcasses, a notion that is supported by telemetry data and observations (Interagency Grizzly Bear Study Team, unpublished data; D. Gustine and K. Wilmot, National Park Service, personal observation). We examined putative family relationships among the resident bears for additional insights into the dynamics of bears present in the study area during the elk hunt. Parentage analyses and live-capture data provided information on known or putative parents or offspring for 6 of the 8 bears we initially defined as residents using the Hines et al. (2003) technique. Five of those bears were closely related (parent–offspring or siblings), involving 3 generations. Based on genetic, remote camera, and telemetry data, we documented almost all of these resident bears scavenging on ungulate carcasses during the elk reduction hunt. Although this foraging strategy may primarily be a function of the lack of other available resources during this time of year, learning may play an important role as well. In American black bears, for example, Mazur and Seher (2008) concluded that social learning was the dominant mode of transmission for foraging on anthropogenic foods. Similarly, Morehouse et al. (2016) found that offspring of female grizzly bears involved in conflicts associated with agricultural activities in southwestern Alberta, Canada, were more likely to be involved in conflicts themselves through social learning. In our study, several females with dependent young remained active into early December, when substantial snow accumulation had already occurred and food resources were limited to ungulate carcasses.
Overall estimated abundance was greatest during the September–October period, prior to the elk hunt, and our post hoc analysis indicates this was mostly due to transient bears. This time period coincides with the peak period of hyperphagia, during which grizzly bears increase their food consumption to store fat and protein in preparation for hibernation (Nelson et al. 1983, Schwartz et al. 2014). Grizzly bears may cover different portions of their large home ranges in search of, or en route to, high-calorie foods during this time (e.g., whitebark pine [Pinus albicaulis] seeds), possibly explaining the single detections that occurred during this period, as they moved through the sampling grid.
Data from telemetry studies during the 1980s and 1990s showed grizzly bears moving outside Yellowstone National Park once elk hunting began near the park's northern and southern boundaries, particularly during years of poor whitebark pine cone production (Ruth et al. 2003, Haroldson et al. 2004, Schwartz et al. 2010). Even when cone production was good, grizzly bears were >2 times more likely to be outside the park during the hunt, suggesting that scavenging carcasses provided a nutritional benefit (Haroldson et al. 2004). These previous findings do not necessarily contradict those of our study because a key difference between the studies is the timing of the elk hunting season. Elk hunts in the earlier telemetry studies involved early rifle seasons, starting on 15 September in Montana and 10 September in Wyoming, and extending into mid to late November. In our study, the elk hunt started 4–6 weeks later and extended into early December, with the peak of harvest occurring from mid- to late November. As discussed previously, the later starting dates likely contributed to the lack of a similar response of bears moving into the hunt areas during the Elk Reduction Program, which is a key consideration for managers.
Detection probabilities (M: p = 0.21; F: p = 0.42) were comparable to, or greater than, those reported for other DNA sampling studies on grizzly bears (e.g., Boulanger et al. 2004, Sawaya et al. 2012, Rovang et al. 2015) and above the threshold of p ≥ 0.20 suggested by several researchers to obtain reliable estimates of population size and trend in bear studies (Otis et al. 1978, Boulanger et al. 2002, Ebert et al. 2010, Laufenberg et al. 2013). Additionally, live captures conducted during a portion of the DNA sampling periods in both years indicated that only 1 out of 6 captured individuals was not detected from the hair samples.
A period effect was central to our research hypothesis; therefore, we considered whether seasonal patterns in successful extractions could have biased how many individuals were detected in the different primary periods. For example, extraction of DNA can be affected by weather conditions, particularly moisture due to rain or snow (Stetz et al. 2015), and the quality of hair samples may be affected by patterns of molting and hair growth. The proportions of failed, insufficient, and mixed samples combined did not vary much across 5 primary periods (0.20–0.24), but it was small (0.11) for period 1 in 2015. Given that the population estimate for period 1 was low in both years, this cannot account for the distinct shifts in the number of detected individuals we observed over the 3 primary sampling periods. Therefore, there was no evidence that seasonal differences in proportion of successfully extracted samples affected our findings.
We provide several important caveats to our study. As we noted in the methods, the robust design methods lack an adequate goodness-of-fit test and our use of an alternative approach based on primary period data only partially addressed this. We addressed a possible violation of the closure assumption in the last primary period each year (due to onset of denning behavior) by using an initial and second encounter history, the latter by collapsing the 3 final secondary sampling periods (Kendall 1999). Nonetheless, a few bears may have already started denning during the initial secondary sampling period of the third and sixth primary periods, which could have biased the estimates of N negatively. However, given that we estimated 88% of monitored bears were still active in the first week of November, such a bias would likely be limited and not change our conclusion regarding a lack of a population-level response to the Elk Reduction Program. An important assumption of the robust design method is that the probability of survival is the same for all animals in the sampled population, regardless of availability for detection. Estimates of survival from known-fate studies for the entire Greater Yellowstone Ecosystem for the period 2002–2011 showed no statistical differences among sex and age groups for independent-age bears (i.e., ≥2 yr old, S = 0.948), suggesting the data were robust to this assumption (Interagency Grizzly Bear Study Team 2012). Additionally, we have no evidence that survival would be different for bears that were or were not available for detection. Finally, the detection of transients in primary periods 2 and 5 is a violation of the closure assumption in the robust design analysis. We addressed this in our post hoc analysis and were able to estimate the abundance of resident versus transient bears separately, which further supported our conclusion that the data did not support our research hypothesis.
Inference from our study was limited by small sample sizes, model uncertainty, lack of precision among parameter estimates, and relatively short duration of 2 years. We probably could not fully account for heterogeneity in capture probabilities because of these factors. Despite these limitations, our inference was enhanced by similar patterns across both study years and the fact that we established 2 primary sampling periods prior to hunt. The first primary sampling period in each year allowed us to detect, particularly after adding the post hoc analysis, that there was temporary movement into the study area due to transients between early summer and late summer. Without these early sampling periods, we would not have known that the second primary sampling period in each year was unique in that regard. Given that the primary signal in the data was of transient grizzly bears entering the study area prior to the elk hunting season, but not remaining for the hunt period, the probability that we falsely rejected our research hypothesis is low.
Management implications
The primary motivation for this study was to better understand how grizzly bears respond to the elk harvest in Grand Teton National Park, so that managers can make science-based decisions to reduce the risk of encounters between elk hunters and grizzly bears. The data from our study suggest that, given current circumstances and environmental conditions, the late timing of the elk hunt helps to limit the use of elk carcasses to a small number of resident bears. We speculate that the familiarity of resident bears with people may increase their tolerance to accommodate the greater human presence in the study area during the hunt. Although continuation of the Elk Reduction Program with the current timing likely represents a scenario with a low relative risk, elk hunters should be aware that encounter risks remain real, as they are anywhere within occupied grizzly bear range. Thus, maintaining the status quo regarding the timing of the elk hunt would not diminish the importance of current strategies that are in place to reduce the risk of hunter–bear encounters, such as the requirement to carry bear spray and closure of areas near the Snake River bottoms. The timing and location of the Elk Reduction Program are unique, so we caution that our study findings may not apply elsewhere in the ecosystem. Finally, our analysis focused on population-level responses, but encounter risk is also a function of other factors that operate more at the level of individual grizzly bears and elk hunters, primarily the spatial and temporal distribution of elk carcasses, grizzly bear detection and use of elk carcasses, and the relative risk of hunter–bear encounters. Those factors are being evaluated as part of a broader study investigating grizzly bear responses to elk hunting in Grand Teton National Park.
Acknowledgments
This study was funded through the Natural Resource Protection Program, a joint program of the U.S. Geological Survey and National Park Service. Additional funding was provided by the Grand Teton National Park Foundation. We thank S. Cain for his support of this project and biological technicians T. Ritter, D. Weller, and V. Villalobos for assistance with the collection of hair samples. We thank B. Karabensh and S. Dewey for compiling auxiliary data. This study benefited substantially from previous data collections by the partner agencies of the Interagency Grizzly Bear Study Team: U.S. Geological Survey, National Park Service; U.S. Fish and Wildlife Service; U.S. Forest Service; Wyoming Game and Fish Department; Montana Fish, Wildlife and Parks; Idaho Department of Fish and Game; and the Eastern Shoshone and Northern Arapaho Tribal Fish and Game Department. We appreciate the constructive comments from the Associate Editor and reviewers, and thank J. D. Clark for review of an earlier draft of this manuscript as part of the U.S. Geological Survey's Fundamental Science Practices. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.
