While migrating, animals make directionally persistent movements and may only respond to human-induced rapid environmental change (HIREC), such as climate and land-use change, once a threshold of HIREC is surpassed. In contrast, animals on other seasonal ranges (e.g., winter range) make more localized and tortuous movements while foraging and may have the flexibility to adjust the location of their range and the intensity of use within it to minimize interactions with HIREC. Because of these seasonal differences in movement, animals on seasonal ranges should avoid areas that contain any level of HIREC, however, during migration, animals should use areas that contain low levels of HIREC, avoiding it only once a threshold of HIREC has been surpassed. We tested this hypothesis using a decade of GPS collar data collected from migratory mule deer (Odocoileus hemionus; n = 56 migration, 143 winter) and pronghorn (Antilocapra americana; n = 70 migration, 89 winter) that winter on and migrate through a natural gas field in western Wyoming. Using surface disturbance caused by well pads and roads as an index of HIREC, we evaluated behavioral responses across three spatial scales during winter and migration seasons. During migration, both species tolerated low levels of disturbance. Once a disturbance threshold was surpassed, however, they avoided HIREC. For mule deer, thresholds were consistently ~3%, whereas thresholds for pronghorn ranged from 1-9.25% surface disturbance. In contrast to migration, both species generally avoided all levels of HIREC while on winter range. Our study suggests that animal responses to HIREC are mediated by season-specific movement patterns. Our results provide further evidence of ungulates avoiding human disturbance on winter range and reveal disturbance thresholds that trigger mule deer and pronghorn responses during migration – information that managers can use to maintain the ecological function of migration routes and winter ranges.

We used helicopter net-gunning to capture female mule deer (n=183) and pronghorn (n=90) on winter range. Captured animals were fit with GPS store-on-board collars (Telonics, Mesa, Arizona) that collected locations every 2-3 hours.We attempted to sample animals in proportion to their abundance within the study area, as determined by a pre-capture aerial survey(Sawyer et al., 2017). All captures were completed in compliance with The University of Wyoming Institutional Animal Care and Use Committee and recommendations of the American Society of Mammalogists(Sikes et al., 2011). Mule deer were captured across a 17-year period, from 2001 to 2017, and pronghorn were captured across an 8-year period, from 2009 to 2016. Mule deer were collared for 1-7 years and pronghorn for 1-2 years. Most mule deer (n=127) migrated along the western flanks of the study area and did not interact with development infrastructure. To evaluate how energy development affects mule deer and pronghorn movement behavior, our analyses wererestricted to 56 mule deer and 89 pronghorn that migrated through or within ~1 km of the study area boundary, resulting in 15 years of data from mule deer and 8 years for pronghorn.

Analysis - To assess the response of mule deer and pronghorn to HIREC during migrationand while on winter range, we evaluated movement and habitat selection patterns as a function of surface disturbance associated with energy infrastructure. We built on earlier work on non-linear responses to development by migrating mule deer (Sawyer et al., 2020) to consider species-specific behavioral responses relative to season. Movement and habitat selection analyses were conducted at 3 spatial scales. At the finest scale, we assessed how animals made fine-scale movements in response to disturbance. Second, at a broad scale, we assessed how much disturbance occurred within each animals’ migration and winter range contours. This broad scale analysis differed from the other two in that it was a use analysis instead of a selection analysis.Finally, at the broadest scale, we assessed where animals placed their migration routes relative to all disturbance within the gas field and how animals select habitat within their winter ranges relative to disturbance. Habitat use by animals is scale-dependent and thus conducting analyses across spatial scales provides strong inference on animal behavior across space and time (DeCesare et al., 2012).

We extracted migration sequences by identifying start and end dates for spring and autumn migration using plots of net squared displacement (Bunnefeld et al., 2011). We defined winter sequences by a fixed period from 15 December to 15 March. We obtained 277,350 deer locations in winter and 33,600 deer locations during migration. For pronghorn, we obtained 109,237 locations in winter and 35,803 locations during migration. We quantified surface disturbance as the proportionof native habitat that was converted to both well pads and roads within species and season-specific buffered GPS points (Sawyer et al., 2020).The amount of well pads and roads (surface disturbance) being added to our study area changed from year to year and leveled off around 2012. Therefore, we matched the amount of surface disturbance in each year with its respective years of GPS data.

At the finest scale, we assessed whether surface disturbance influenced fine scale movement through the study area during migration and while on winter range using a step selection function (SSF; Fortin et al.,2005). SSFs compare used steps (the endpoint of a linear segment between two consecutive animal locations) to randomly generated endpoints from steps taken from the same starting location. For each used step, we generated 25 available target steps drawn from the step-length and turning angle distributions from all animals in the population calculated separately for each species and season. We confirmed that we properly sampled availability in SSFs by resampling available points and re-fitting the models, and verifying that beta coefficients did not change. For both migration and winter range, surface disturbance was calculated at each end point of each step by creating a buffer around those points. The appropriate buffer size (analogous to the scale at which animals respond to surface disturbance) was determined, for each species and season separately, by comparing relative empirical support for SSFs with a linear effect of surface disturbance calculated at 11 different buffer sizes with the following radii: 25, 50, 100, 150, 200, 250, 300, 350, 400, 450, and 500 m. To focus only on movement when individuals are in or near the development, we only included a used step and its paired 25 available steps if at least 1 of the 26 steps intersected with > 0.01% surface disturbance. During migration, our sample size thus totaled 56 mule deer with 145 migration sequences, and 68 pronghorn with 136 migration sequences. For winter range, our sample size totaled 135 mule deer with 268 winter range sequences, and 84 pronghorn with 113 winter range sequences. All SSFs were fit using conditional logistic regression. Temporal and spatial autocorrelation were accounted for in SSFs by using generalized estimation equations (GEE) to calculate standard errors and 95% confidence intervals (Craiu et al., 2008). In the GEE analysis, we assigned a unique cluster for strata within a given individual’s seasonal migration sequence or an individual’s annual winter range (Craiu et al., 2008). Empirical support for SSFs was determined by comparing quasilikelihood under the Independence model Criterion (QIC)(Pan, 2001).We validated our SSFs using a simple modification of the cross validation methods outlined in Fortin et al. (2009). See Appendix S2 for details.

Piecewise regression has been shown to be an effective tool in identifying ecological thresholds (Toms & Lesperance, 2003). Thus, once the appropriate buffer size was identified, we tested for non-linear relationships between relative use and surface disturbance by conducting piecewise regression. In the piecewise regression analysis, one or two breakpoints (splines) were specified, and we used a grid search of 0.25% surface disturbance units to test at what surface disturbance value the breakpoints were most supported by the data. If a non-linear relationship was detected, we then plotted the response and determined visually whether thresholds occurred. We defined thresholds as stable or increases in relative probability of use of HIREC initially, followed by a decrease in relative probability of use as HIREC values increased.

At the broad scale, we estimated how much surface disturbance animals experienced along their migration routes and within their winter range contours, using the 99% contour estimated from Brownian bridge movement models (BBMM; Horne et al. 2007; Sawyer et al. 2009). BBMM estimate the movement path of an animal and provide an occurrence distribution. We used the 99% contour of the BBMM because it captures the maximum area of where an animal could have been between GPS fixes. This scale was essentially an analysis of habitat use or where animals placed their migration routes and winter ranges and how much disturbance occurred within them. For migration, we restricted the analysis to individuals who migrated through the study area and interacted with energy development infrastructure (n=56 mule deer with 143 migration sequences, and n=70 pronghorn with 148 migration sequences). All migration routes extended past the study area (up to 100 km beyond the study area). Thus, to isolate when animals migrate through the study area, we defined the study area as a minimum convex polygon (MCP) around all well pads installed from 2000-2017. We then cropped each animal’s migration route to the study area polygon and calculated the percent of surface disturbance that occurred within the cropped migration route. For winter range, to effectively capture how animals used their winter ranges in relation to surface disturbance, we only included individuals whose winter range contour contained >0.01% surface disturbance (n=130 mule deer with 242 winter range sequences and n= 74 pronghorn with 103 winter range sequences). We then simply calculated the percent of surface disturbance that occurred within each winter contour. To test for thresholds at this scale, we plotted histograms of surface disturbance across migration routes and winter range contours. These plots provided a distribution of surface disturbance that each animal experienced along their migration route within the study area or throughout their winter range each year. We then used the segmented package in Program R to assess empirical support for non-linear patterns in this distribution.

At the broadest scale, we assessed where animals placed their migration routes relative to surface disturbance and how animals select habitat within their winter ranges using a resource selection function (RSF) framework(RSF; Manly et al. 1993). RSFs are similar to SSFs, however, they differ in that they compare used GPS locations to randomly generated available locations within a broader availability domain (Fig 1; Manly et al. 1993). For migration, we defined availability as the area of intersection between a MCP around all well pads installed from 2000-2017 (i.e., the study area polygon) and a MCP around the distribution of all observed GPS locations during migration. For migration, we again restricted this analysis to individuals who migrated through the study area and interacted with at least 0.01% energy development infrastructure, and keeping only individuals in the RSF analysis if they had >5 GPS points in the study area per migration (n=49 mule deer with 87 migration sequences, and n=49 pronghorn with 80 migration sequences). For winter range, we defined availability as individual animal winter ranges, calculated using the 99% contour estimated from the occurrence distribution from the BBMM. To ensure our analysis was based only on animals that clearly interacted with surface disturbance, we only included individual winter ranges in the RSF if they overlapped with > 0.01% surface disturbance. Further, we only included individual winter sequences in the RSF analysis if they had >10 GPS points during winter (n=143 mule deer with 291 winter range sequences, and n=89 pronghorn with 133 winter range sequences).

For all RSFs, we chose the number of available points by starting with a 1:1 used to available ratio. We then increased the number of available points and reran the analysis until beta coefficients did not change. In the migration analyses, for each used GPS point within the availability domain within a given year and season, we randomly sampled 10 and 100 available locations for mule deer and pronghorn, respectively, within the study area polygon. In the winter range analysis, for every used GPS point within each individual’s winter range, we sampled one and three available locations for pronghorn and mule deer, respectively. We accounted for psuedoreplication in the winter data by reducing each animal’s used points to one randomly used point per day. We validated our RSFs using a simple modification of the cross validation methods outlined in Boyce et al. (2002). See Appendix S2 for details.

For both migration and winter range, surface disturbance was calculated at each used and available point by creating a buffer around those points. The appropriate buffer size was determined, for each species and season separately, by comparing relative empirical support for RSFs with a linear effect of surface disturbance calculated at 12 different buffer sizes with the following radii: 50, 100, 150, 200, 250, 300, 350, 400, 450, 500, 750, and 1000 m. All RSFs were fit using mixed effects logistic regression. We accounted for psuedoreplication within individuals and years, and variation among individuals and years by including a random intercept for each individual’s migration and winter range in each year (Gillies et al., 2006). Potential non-linear models were identified using the same piece-wise regression methods described for the fine-scale analysis. Relative empirical support for linear and non-linear RSFs was assessed by using Akaike’s Information Criterion (AIC; Anderson, Burnham, and White 1994).

For the migration data (mig in file name), the column id_mig_yr identifies data from one individual in one year and in either spring or fall (indicated by either sp or fa followed by a two digit year). Similarly, for winter range data (WR in file name), the column id_wint_year identifies data for one individual in one year for the winter season (animal id is first followed by four digit year). The column case includes an animal's used (1) and randomly generated available points (0). In the step selection function (SSF) analyses, there is a strata column that groups used steps with their randomly generated available steps. Columns named surfdist_res100 show the amount of surface disturbance calculated in used and available points. To calculate surface disturbance, we buffered used and available points by a given diameter, in this case 100m (50m radii in our manuscript). 

For the Brownian Bridge Movement Model datasets (BBMM in file name), the column area_bb_km2 gives the area of the animal's migration route or winter range, which was defined using the BBMM. The column area_surfdist_km2 gives the area (square kilometers) of the migration route or winter range that is covered by surface disturbance. The column density_surfdist gives the percent of the migration route or winter range that is covered by surface disturbance. The column acres_per_sqrmile gives how many acres per square mile are covered by surface disturbance in the migration or winter range.

All datasets include an animal id column (AID).