DOI: 10.5061/dryad.280gb5mrn
title: Responses to natural gas development differ by season for two migratory ungulates

description: 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.
description: 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. <strong>Analysis</strong> - 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 &gt; 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 &gt;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 &gt;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 &gt; 0.01%
 surface disturbance. Further, we only included individual winter sequences
 in the RSF analysis if they had &gt;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).
description: 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).