Reproductive success of mule deer in a natural gas development area

Natural gas development is increasing across North America and causing concern over the potential impacts on wildlife populations and their habitat, particularly for ungulate species. Understanding how this development impacts reproductive success metrics that are influential for ungulate population dynamics is important to guide management of ungulates. However, the influences of natural gas development on reproductive success metrics of mule deer Odocoileus hemionus have not been studied. We used statistical models to examine the influence of natural gas development and temporal factors on reproductive success metrics of mule deer in the Piceance Basin, northwest Colorado during 2012–2014. We focused on study areas with relatively high or low levels of natural gas development. Pregnancy and in utero fetal rates were high and statistically indistinguishable between study areas. Fetal survival rates increased over time and survival was lower in the high versus low development study areas in 2012 possibly influenced by drought coupled with habitat loss and fragmentation associated with development. Our novel results suggest managers should be concerned with the influences of development on fetal survival, particularly during extreme environmental conditions (e.g. drought) and our results can be used to guide development planning and/or mitigation. Developers and wildlife managers should continue to collaborate on development planning, such as implementing habitat treatments to improve forage availability and quality, minimizing disturbance to hiding and foraging habitat particularly during parturition, and implementing directional drilling to minimize pad disturbance density to increase fetal survival in developed areas.

Natural gas development is increasing worldwide, causing concern over the potential impacts on wildlife and their habitat (Northrup and Wittemyer 2013). Impacts on mule deer Odocoileus hemionus population dynamics and their habitat are of particular interest in North America due to the recreational, social, and economic importance of deer as a game species (Sawyer et al. 2009, Northrup et al. 2015. While a number of studies have assessed the impacts of natural gas development on mule deer behavior (Sawyer et al. 2006, Northrup et al. 2015, the influencess on reproductive success have not been studied. Specifically, accurate estimates of pregnancy rates (i.e. proportion of adult females carrying 1 fetus), in utero fetal rates (i.e. the number of fetuses per pregnant female), and fetal survival rates (i.e. survival of fetuses to birth) are needed to quantify fawn recruitment and population dynamics (Bonenfant et al. 2005, DeCesare et al. 2012) Thus, we quantified reproductive success parameters to assess if and how natural gas development may influence mule deer populations.
The mechanism by which natural gas development may influence ungulate reproductive success is through direct and indirect habitat loss. Direct habitat loss results from construction of well pads, roads, compressor stations, and pipelines. Conversely, indirect habitat loss may result from activity and noise associated with increased human presence and development, which could lead to a zone of avoidance around development that is greater than the footprint itself. Past studies suggest that deer generally decrease time spent near roads (Webb et al. 2011, Lendrum et al. 2012) and well pads (Sawyer et al. 2006, Northrup et al. 2015, suggesting indirect habitat loss. In addition, development disturbances may cause stress, alter behavior and habitat use, and decrease forage and habitat availability (Sawyer et al. 2006, Lendrum et al. 2012, Northrup et al. 2015. These processes all have the potential to influence body condition of maternal females by reducing available foraging habitat or by limiting their time spent foraging or causing increased energy expenditure (Frid and Dill 2002). Thus, reproductive success could be negatively influenced by development.
We examined the influence of natural gas development and temporal (e.g. year) factors on reproductive success metrics of mule deer in the Piceance Basin, northwest Colorado, during 2012-2014. We estimated reproductive success metrics in areas with relatively high (0.04-0.90 well pads km -1 ) or low (0.00-0.10 well pads km -1 ) levels of natural gas development. Our objectives were to test hypotheses that reproductive success metrics would be lower in the high development areas than the low development areas and vary by year with increased precipitation influencing vegetation availability and quality. Our study provides the first insights into reproductive success of mule deer in a natural gas development area, which is helpful to comprehend mule deer population dynamics to assist in future development planning.

Study area
During 2012-2014, we examined reproductive success metrics of migratory mule deer in the Piceance Basin, northwest Colorado. Our winter range study area included four study units in the Piceance Basin and are part of a larger research project (Anderson 2016). Deer occupied two winter range study units with relatively high levels of natural gas development (0.6-0.9 well pads km -1 ) and two winter range study units with relatively low levels of development (0.0-0.1 well pads km -1 ). We note that wells in our study area were directionally drilled from multiple well pads, which reduces the development density compared to coalbed methane and single-well drilling development. Winter range habitat was geographically diverse and comprised of two-needle pinyon pine Pinus edulis, Utah juniper Juniperus osteosperma, and mountain shrublands.
Summer range study units included parts of Garfield, Moffat, Rio Blanco, and Routt counties in northwestern Colorado. Deer from winter range units with relatively high levels of natural gas development generally migrated to the Roan Plateau summer range (Lendrum et al. 2013) where deer potentially encountered natural gas development (0.04-0.06 well pads km -1 ). Hereafter we refer to deer inhabiting the winter and summer range study units with relatively high levels of development as being in the high development study areas. Deer from winter range units with relatively low levels of natural gas development generally migrated towards the Flat Tops Mountain Range summer range (Lendrum et al. 2013) where deer encountered minimal natural gas development (0.00-0.01 well pads km -1 ). Hereafter we refer to deer inhabiting the winter and summer range study units with relatively low levels of development as being in the low development study areas. Summer range habitat was dominated by Gambel oak Quercus gambelii, quaking aspen Populus tremuloides, two-needle pinyon pine, Utah juniper, and mountain shrublands. Shrublands included alderleaf mountain mahogany Cercocarpus montanus, antelope bitterbrush Purshia tridentate, big sagebrush Artemisia tridentate, mountain snowberry Symphoricarpos oreophilus, rubber rabbitbrush Ericameria nauseosa, and Utah serviceberry Amelanchier utahensis. Drainages bisected the study units and most of the primary drainage bottoms have been converted to irrigated, grass hay fields. Shrubs, forbs and grasses common to the area are listed in Bartmann (1983) and Bartmann et al. (1992).

Adult female capture and handling
During December 2011-2013, adult (2.5 years old) female mule deer were captured in each of the four winter range study units using helicopter net gunning (Webb et al. 2008, Jacques et al. 2009). Deer were blindfolded, hobbled, chemically sedated with 0.5 mg kg -1 of Midazolam and 0.25 mg kg -1 of Azaperone given intramuscularly and ferried to a central handling location. We fit each captured deer with a global positioning system (GPS) radio collar (Model G2110D, Advanced Telemetry Systems).
During early March 2012-2014, radio-collared adult females were recaptured on winter ranges using helicopter net gunning. We performed transabdominal ultrasonography (SonoVet 2000, Universal Medical Systems) to determine pregnancy status and number of in utero fetuses (Stephenson et al. 1995, Bishop et al. 2007). If an adult female was pregnant, we inserted a vaginal implant transmitter (VIT; Model M3930, Advanced Telemetry Systems) following VIT insertion procedures described in detail by Bishop et al. (2011) and Peterson (2016). Helicopter net gunning and vaginal implant transmitters have been used without influencing reproductive success and are effective methods to capture females and neonates (Carstensen et al. 2003, Bishop et al. 2009, Monteith et al. 2014).

Adult female monitoring and neonate capture
During parturition (late May-mid-July), we checked for radio collar and VIT signals daily from a fixed-wing aircraft. When we detected a fast pulse rate (80 beats min -1 ) signifying parturition, ground crews used radio telemetry to simultaneously locate the radio-collared female and expelled VIT. Ground crews then searched for neonates around (400 m) the female and VIT for up to 1 h. If neonate(s) documented in utero were not captured during the initial attempt, crews located the female on the next two days and searched near the female to locate neonates.
Ground crews attempted to determine the fate of each female's fetus(es) documented in March as live or stillborn neonates, including scavenged remains. Unless evidence suggested a neonate was born alive at a birth site (e.g. milk in the abomasum), crews classified a dead neonate as stillborn. We submitted stillborn neonates to the Colorado Parks and Wildlife's Health Laboratory (Fort Collins, CO) for necropsy to confirm that a neonate had died before birth.
During 2012 and 2013, ground crews captured neonates in the high and low development study areas. In 2014, crews captured neonates predominantly in the high development study areas and sporadically in the low development study areas because VIT photo sensors malfunctioned. We focused our effort in the high development areas during 2014 because it was logistically more difficult to monitor deer from the ground and capture neonates in the low development areas due to geographic separation among birth sites. Each captured neonate was blindfolded and sexed. All individuals who handled neonates wore latex gloves to minimize transfer of human scent. All capture, handling, radio collaring, and VIT insertion procedures were approved by the Institutional Animal Care and Use Committee at Colorado Parks andWildlife (protocol no. 17-2008 andno. 01-2012).

Statistical methods
We modeled pregnancy and fetal rates of adult females as a function of winter range study area and year using PROC LOGISTIC and PROC MIXED (e.g. generalized linear models) in SAS (SAS Inst.), respectively. We modeled fetal survival from March to birth as a function of study area and year using PROC NLMIXED in SAS and a joint-likelihood described in Bishop et al. (2008). We were not able to determine fate of all fetuses detected in utero because neonates were challenging to detect and some VITs malfunctioned, thus we used the joint-likelihood with six nuisance parameters (relative to our interests in this paper) to estimate fetal survival probability (S 1 ). The six nuisance parameters are neonatal survival probability from birth to 5 days old (S 2 ), the probability of detecting a neonatal fawn 1 day old given that field crews conducted a search 1 day after birth (p 1 ), the probability of detecting a neonatal fawn 1 day old given that crews conducted a search 1 day after birth (p 2 ), the probability of detecting a stillborn fetus when a VIT was not expelled at a birth site (r), the probability of locating a radio-collared adult female and searching for her neonate(s) 1 day after birth (a), and the probability a VIT was expelled at a birth site (b). We modeled S 2 as constant or as a function of study area to account for survival differences between areas. We modeled p 1 , p 2 , a and b as constant or as a function of study area and year to account for temporal differences in detection probabilities. We constrained r to be constant because crews did not locate stillborns without the aid of a VIT during some years and in some study areas, thus we could not separately estimate r. We assumed fetal survival data were not overdispersed based on the recommendation of Bishop et al. (2008). Lastly, we fit the same model set for reproductive success metrics as Bishop et al. (2009), except age class, and that we hypothesized would influence reproductive success.
We used Akaike's information criterion adjusted for small sample size (AIC c ), ΔAIC c , and AIC c weights (w i ) for model selection (Burnham and Anderson 2002). We used model averaging to obtain model-averaged parameter estimates when more than one model was within 2 AIC c units of the top-ranked model (Burnham and Anderson 2002).

Results
We documented pregnancy status of 346 adult females, of which 204 produced 383 fetuses (31, 167, and 6 females with 1, 2 or 3 fetus(es), respectively). Seventeen females were not pregnant and we were unable to determine accurate fetal counts for 127 females for various reasons (e.g. denied access to private property, VIT malfunctioned) and we excluded these females from the fetal analyses. Ultimately, we documented sex of 195 fetuses (99 males and 96 females).

Reproductive success metrics
A model indicating constant pregnancy rates across years ranked highest (w i = 0.776; Supplementary material Appendix 1 Table A1). We found minimal support for a study area effect on pregnancy rates (Supplementary material Appendix 1 Table A1). Pregnancy rate for all adult females during the study was 0.948 (SE = 0.012). Pregnancy rate for females in the high and low development areas was 0.953 (SE = 0.016) and 0.942 (SE = 0.018), respectively.
A model indicating constant fetal rates across years ranked highest (w i = 0.766; Supplementary material Appendix 2 Table A2). We found minimal support for a study area effect on fetal rates (Supplementary material Appendix 2  Table A2). Fetal rate for all adult females during the study was 1.877 (SE = 0.029) and most females produced twins (0.819; Supplementary material Appendix 3 Table A3). In utero fetal rate for females in the high and low development areas was 1.849 (SE = 0.037) and 1.908 (SE = 0.044), respectively. The most parsimonious model for fetal survival from March until birth included an interaction between study areas and year (w i = 0.714; Supplementary material Appendix 4 Table A4). The same model for fetal survival, but without the study area variable had little support (ΔAIC c = 17.598, w i  0.0014). Fetal survival was higher in the low development areas than the high development areas in 2012, whereas we found no difference in 2013 and 2014 (Fig. 1) In the high and low development areas, respectively, females produced eight (11%) and zero stillborn fetuses in 2012, eight (12%) and three (4%) stillborns in 2013 and zero and zero stillborns in 2014.

Discussion
We found pregnancy and in utero fetal rates were high, showed little variation across years, and were similar in the high versus low development areas. Bishop et al. (2009) found no difference in pregnancy and in utero fetal rates when examining the effects of supplemental nutrition treatments versus a control group in a different area of Colorado and Anderson (2016) found no difference in body condition of adult females between the high and low development areas in our study area. Pregnancy and fetal rates were high in each area and in the upper range of previous estimates (0.860-1.000 and 1.650-2.010, respectively) across Colorado (Andelt et al. 2004, Bishop et al. 2009). Of note, deer abundance is trending upward in the Piceance Basin (Anderson 2016) after a decline during the 1990s (White and Bartmann 1998) and could be partly explained by high fetal rates coupled with fawn recruitment, which has largely been driven by relatively high overwinter fawn survival (Anderson 2016). Ultimately, high pregnancy and fetal rates seem to be the norm for deer despite a wide range of spatial and temporal differences across populations (Bishop et al. 2009, Hurley et al. 2011, Monteith et al. 2014 including our study area with natural gas development. Fetal survival from March until birth was higher in the low development areas than the high development areas, suggesting an influence of development, although survival varied annually. However, fetal survival rates exceeded previous estimates (0.747-0.983) measured on the Uncompahgre Plateau, southwest Colorado (Bishop et al. 2009). Annual variation in fetal survival could be the result of an interaction between environmental conditions and development. Annual variation in precipitation may alter the onset of spring green-up (Pettorelli et al. 2005), which can influence maternal condition (Parker et al. 2009) and possibly reduce fetal survival. Increased precipitation in arid environments is linked to forage availability (Derner et al. 2008) and quality and growth of forbs (Marshal et al. 2005), thus drought conditions may reduce forage availability and/or quality below levels needed for growth of fetuses (Parker et al. 2009). Access to high quality forage is necessary to meet the energetic demands of the last trimester when most fetal growth occurs (Armstrong 1950). Precipitation during the third trimester (1 April -15 June) was lower in 2012 (4 cm) than 2013 (12 cm) and 2014 (12 cm), suggesting reduced forage availability and growth of forbs, which may have contributed to lower fetal survival particularly in the high development areas during 2012. Further, dry weather likely reduced forage availability which may have been exacerbated by habitat loss and fragmentation associated with development, possibly contributing to lower fetal survival in the high development areas during 2012. Of note, eight stillborn neonates were produced in the high development areas versus zero in the low development areas during 2012. Stillborn fetuses were mostly small and lightweight suggesting reduced forage availability and quality contributed to increased stillborns (Verme 1969) and consequently decreased fetal survival in the high development areas during 2012. However, active natural gas development activity was minimal during our study because most wells were in production (a phase that is characterized by less human activity and construction than the more active drilling phase), thus the influence of development could be stronger with increased development intensity and associated disturbances. Overall, development coupled with extreme environmental conditions (e.g. drought) may have contributed to lower fetal survival we observed during 2012.
The probability of detecting a neonate 1 day old was low, but was highest in 2012 and similar in 2013 and 2014 because neonates were challenging to detect and some VITs malfunctioned particularly in 2014. The probability of detecting a neonate 1 day and 5 days old was also low, but increased from 2012 to 2014. Many VITs failed in 2014 providing minimal assistance in detecting neonates at birth sites, thus contributing to higher detection of older neonates as mothers and presumably neonates move farther from VITs and birth sites as they age (Vore andSchmidt 2001, Long et al. 2009).
Estimating reproductive success metrics from marked adult females is helpful to understand fawn recruitment and population dynamics of ungulates (Bonenfant et al. 2005). Our results suggest managers should not be concerned with the influences of natural gas development on pregnancy and fetal rates under existing development conditions during this study. However, we suggest that future research should be conducted in areas with increased development intensity to further evaluate the influence of natural gas development on pregnancy and fetal rates. Contrarily, managers should be concerned with the potential influences of development on fetal survival as our results suggest fetal survival was lower during 2012 from increased stillbirths in the high development areas when drought conditions also were present. If development indeed caused an increase in stillbirths during extreme environmental conditions, fetal survival could be increased if forage availability and quality is improved. Thus, developers and wildlife managers should continue to collaborate during development planning to avoid important habitats during critical time periods (e.g. parturition) and consider habitat treatments (e.g. hydro-ax, roller chopping, and seeding) and/or reclamation plans to improve forage availability and quality (Johnston andChapman 2014, Stephens et al. 2016) to enhance fetal survival and possibly fawn recruitment. Future research should focus on attempting to more precisely identify the mechanism underlying the documented differences in fetal survival.