Study area

We conducted our research at the Beaver River Wildlife Management Area (WMA) which is located in Beaver County, OK (36°50′21.62″N, 100°42′15.93″W). The total area of the WMA is approximately 11 315 ha, however for our research, the WMA was divided into two separate units (Beaver River Unit 6824 ha; McFarland Unit 4491 ha), in which the Beaver River Unit was exposed to hunting pressure while the McFarland Unit was not. Unit boundaries were separated by barbed wire fencing and had signs indicating whether or not bobwhite hunting was permitted based on our study design. Both units primarily consist of upland areas dominated by tivilo fine sand soils and a floodplain dominated by lesho silty clay loam. Areas bordering the WMA consisted of private property with land use practices that included grazed rangelands, Conservation Reserve Program (CRP) enrollment and row crops. The only nonprivate property that bordered the WMA was a city park on the eastern boundary, which was used for off-road vehicle recreation and hiking. We were unable to determine if adjacent private properties experienced hunting pressure. However, hunting on these properties would not have occurred without consent from landowners and no hunting was permitted within the boundary of the city park.

Anthropogenic features on the WMA consisted of roads, power lines, buildings and oil/gas structures. Power lines were not included in our analyses as very few of these features existed within our study area and bobwhite rarely encountered them (density of 0.08/km²). Overall density of roads was 2.12 km/km². Additionally, there were six buildings (0.05 structures/km²) and 94 oil/gas structures (0.83 structures/km²) on our study area. A study area map illustrating the spatial orientation of roads, buildings, oil/gas structures, and hunting units at Beaver River WMA is given in the supplementary material (Supplementary material Appendix 1 Fig. A1).

During our study, quail hunting season began on 10, 9 and 8 November in 2012, 2013 and 2014 respectively and ended on February 15 of the following year for all three years. Hunting activity was recorded for all two week periods during our study except weeks 8–9 of each year (Christmas and New Year holiday). However we suspect there was hunting activity occurring during these periods though technicians were not available to verify.

During the course of the study, annual precipitation was 34.44, 50.29 and 39.42 cm in 2012, 2013 and 2014 respectively, while the long term (1895–2014) average annual precipitation for this region is 49.63 cm. Climate data were obtained from the Beaver Mesonet station (Brock et al. 1995, McPherson et al. 2007). At no time were our two study units out of drought conditions (The National Drought Mitigation Center, Lincoln, Nebraska, USA).

Vegetation map development

An unsupervised maximum combined vegetation classification method was used to develop our vegetation map from 2 m resolution satellite imagery using ArcMap 10.1. Aerial imagery was collected in July 2013 when cloud cover was minimized. This method resulted in 65 different classes which were reclassified into 10 ecologically meaningful cover types based on field observations and 214 ground-truthed points. The primary cover types that comprised both units were: mixed shrub (consisting of sand plum Prunus angustifolia, fragrant sumac Rhus aromatic, and sand sagebrush Artemisia filifolia), sand sagebrush, mixed grass (little bluestem Schizachyrium scopariu, switchgrass Panicum virgatum, bromes [Bromus spp.; non-native]), short-grass/yucca Yucca glauca, sparse vegetation/exposed soil, bare ground, salt cedar (Tamarix spp.; non-native), open water, developed housing, and food plots (primarily winter wheat Triticum aestivum). A more detailed description of plants found within these cover types is described in Tanner et al. (2015).

Radio-telemetry

Bobwhite were captured between 2012–2015 using walk-in funnel traps (Stoddard 1931). We attached a necklace-style radio transmitter weighing 6 g if a bird met a minimum body mass requirement of 130 g. We located radio-marked individuals a minimum of three times per week. Locations of individuals were determined using the homing method (White and Garrot 1990). We homed in on individuals to within 15 m and recorded the distance and azimuth to the actual bird location as well as recording the Universal Transverse Mercator (UTM) coordinates of the observer with a GPS. Individuals and coveys were located at different times across days to capture any variability of diurnal patterns throughout the non-breeding season. When a mortality signal occurred (12-h signal), we located the collar and confirmed the mortality. Mortalities were categorized as avian, mammalian, unknown, or other based on evidence at the mortality site. A detailed explanation of radio-tracking methods is described in Tanner et al. (2015). All trapping and handling methods were approved by the Oklahoma State University's Institutional Animal Care and Use Committee Permit (no. AG-11-22).

To estimate the accuracy of our telemetry locations, we conducted a telemetry error trial in which each researcher (n = 12) conducted 10 repeated telemetry trials. These trials consisted of individuals estimating the location of a radio transmitter with a known location via homing techniques across 10 random locations within the WMA. The estimated telemetry error during the course of our study was 8.97 m (95% CI = 6.48 to 11.46).

For all analyses in our study, if an individual or covey changed units (hunted or non-hunted) during our monitoring period, we censored that individual or covey from our study. With respect to our survival analysis, two, ten, and nine individuals changed units during the 2012–2013, 2013–2014 and 2014–2015 non-breeding seasons, respectively. For our space use analysis, only one covey had an estimated home range that overlapped between units.

Andersen—Gill models

We used Andersen—Gill (AG) models to estimate hazard rates for bobwhite across both units (Andersen and Gill 1982) using the survival package in program R (ver. 3.1.1, <  www.r-project.org >). We used AG models due to monitoring gaps and left-truncated data entry for individuals (based on staggered entry design, Pollock et al. 1989). This model is similar to a Cox proportional hazard model (CPHM), however it allows for time-varying covariates when estimating hazard rates (Fleming and Harrington 1991, Therneau and Grambsch 2000, Murray 2006, Fieberg and DelGiudice 2009). To estimate bobwhite hazard rates, we left-censored individuals if they entered the population after our initial time interval (1 October) and right-censored individuals if their fate was unknown (Johnson et al. 2004). We also censored individuals that did not survive the first seven days after they were released to control for any effects of capture myopathy (Guthery and Lusk 2004).

Our dataset consisted of 26 time intervals, which were the number of weeks during the non-breeding season (1 October — 31 March). To estimate the effects of anthropogenic features on bobwhite survival, we estimated the mean weekly Euclidean distance (m) to a feature for each individual. This consisted of distance to: oil/gas structures, buildings, and the four different road types (county road, primary WMA roads, restricted access WMA roads [truck and all-terrain vehicle {ATV} traffic], and restricted access WMA roads [ATV traffic only]). We included both functioning and non-functioning oil/gas structures in our analysis as we were primarily interested in responses to the structure (and associated roads) rather than actual activity. Only 6% of oil/gas structures within our study units were considered non-functioning. Furthermore, there were no new oil/gas structures developed during the course of our study, and direct human activity related to oil/gas development was limited to periodic maintenance checks throughout the year. All oil/gas structure locations used in our study consisted of pumpjacks, condensate tanks and meter stations and were within 15 m in height. Spatial oil/gas structure data were obtained through the Oklahoma Corporation Commission in 2013 and were verified through ground-truthing efforts. County road data were obtained from the Oklahoma Department of Transportation (< http://okmaps.org/ogi/search.aspx>), while buildings and all other road data were hand digitized and confirmed via ground-truthing. To determine if the presence of hunting affected survival in our population, we also included a categorical variable based on the unit in which an individual occurred. Other categorical variables included in our analysis were age (adult or juvenile) and year (2012–2013 (year 1), 2013–2014 (year 2) and 2014–2015 (year 3)). Sex (male or female) of individuals was not included as a covariate in our survival analysis as we expected no difference in harvest rates (Shupe et al. 1990) or survival (Cox et al. 2004, Seckinger et al. 2008, Tanner et al. 2012) between sexes for the nonbreeding season.

The primary assumption of the AG model, like the CPHM, is that hazards from covariates are proportional over time (Johnson et al. 2004). To test this assumption, we plotted Schoenfeld residuals and assessed significant deviances of residual plots from 0 (Therneau et al. 1990, Fox 2002). We stratified our third road category (restricted access roads [truck and ATV traffic]) values into three distance categories (< 500 m, 500–1499 m, and ≥ 1500 m) as this variable did not meet the proportional hazard assumption (Fox 2002). Finally, we included a global model in our survival analysis, which included the additive effects of all variables of interest.

We used Akaike's information criterion adjusted for small sample sizes (AIC_c) to rank models relating covariates to hazard rates for bobwhite over the non-breeding season. We considered models with a ΔAIC_c ≤ 2 to be plausible models and determined those to be the most parsimonious based on model weights (w_i) and ΔAIC_c values (Burnham and Anderson 2002). We built models that we found biologically meaningful or models that specifically addressed our research questions. We considered parameters with hazard ratios that had 95% confidence intervals overlapping one to be statistically uninformative to our survival analysis. This is because a hazard ratio of one indicates no difference in proportional hazards (Fox 2002).

Resource utilization functions

We used resource utilization functions (RUFs; Marzluff et al. 2004, Millspaugh et al. 2006) to estimate the relationships between covey space use and environmental variables. We estimated RUFs for coveys rather than individuals as space use by individuals within coveys has been shown to be nonindependent (Janke and Gates 2013, Brooke et al. 2015). We estimated 95% fixed-kernel densities (Worton 1989, Seaman et al. 1999) for coveys having ≥ 30 radio-telemetry locations during the hunting season (Seaman and Powell 1996, Seaman et al. 1999, Lohr et al. 2011) using the geospatial modelling environment (GME; Spatial Ecology LLC, USA). A likelihood cross-validation bandwidth estimator was used to obtain kernel density estimates (KDE), which has been shown to outperform other bandwidth estimators when sample sizes are small (Horne and Garton 2006).

Along with the anthropogenic variables included in our survival analysis, we also incorporated vegetation cover types and theoretical hunting pressure variables into our RUF analysis. As we did not have a direct measure of hunting pressure within our hunting unit over the course of the study, we incorporated hunter behavior data discussed in Richardson et al. (2008) to estimate areas of potentially high, medium, low, and no hunting pressure. The data presented by Richardson et al. (2008) incorporated vegetation cover, distance from roads (< 500 m, 500-< 1500 m, 1500-< 2500 m, and ≥ 2500 m), and % slope (≤ 3% and > 3%) data and used GPS data from hunters at Packsaddle WMA (Ellis County, Oklahoma, USA) to determine selection indices in specific areas by bobwhite hunters (Supplementary material Appendix 1 Table A1). They separated slope categories so that both categories contained ∼50% of the WMA (Richardson 2006). We used these data to model potential hunting pressure on our study site because road density (Packsaddle WMA: 1.86 km/km², Dunkin et al. 2009; Beaver River WMA 2.12 km/km²) and slope (< 3% slope: 50.76% of the area; > 3% slope: 49.24% of the area) were similar between WMAs. We incorporated these data into a model of theoretical hunting pressure on our study site through the use of the weighted overlay tool in ArcGIS 10.2. Taking into consideration the selection indices of hunters provided by Richardson et al. (2008), we used vegetation cover, distance from roads, and % slope in our model, with each variable having equal weight. We assigned values (1–4) to each category within these variables, where one represented the highest level of theoretical hunting pressure and four represented the lowest. Table 1 shows the values assigned to all categories within our variables and Supplementary material Appendix 1 Fig. A2 shows the spatially explicit theoretical hunting pressure model for the hunted unit during our study.

We extracted values for space use and all environmental variables to points centered on every cell within each covey's home range. Cells in our analysis consisted of a square 10 m spatial resolution. We used the Ruf.fit package in program R (ver. 3.1.1) to estimate coefficients of resource use for each variable and for every covey. Our response variable (utilization distribution) was right skewed, thus all values of space use were log-transformed (Hooten et al. 2013, Winder et al. 2014). Because vegetation cover type and theoretical hunting pressure were categorical variables, we removed a class in each variable to serve as a reference class in our analysis (Jachowski et al. 2014). Therefore, we used the sand sagebrush cover type and the highest level of theoretical hunting pressure as the reference class for the vegetation cover and hunting pressure variables, respectively. The sand sagebrush class was used as a reference because it is the most abundant vegetation type on our study site (Jachowski et al. 2014). To directly address the question of whether bobwhite were altering their space use in relation to higher hunting pressure, we used the highest theoretical hunting pressure class as a reference class to compare bobwhite space use in relation to other hunting pressure categories. Mean standardized β coefficients () and conservative estimates of variance were calculated for each environmental variable to estimate overall population responses to these variables across the hunted and non-hunted units (Marzluff et al. 2004). Standardized coefficients with 95% confidence intervals overlapping 0 were considered non-significant. Finally, we estimated the number of individual coveys that had significant positive, negative, or non-significant relationships to our environmental variables to indicate differences among coveys. We were unable to evaluate all levels of theoretical hunting pressure (high, medium, low, or no pressure) for 10 coveys. As such, low hunting pressure and no hunting pressure variables were not contained within the KDEs for all possible coveys in our analysis and population level estimates were limited to the coveys that encountered these variables (Jachowski et al. 2014).

Table 1.

Variable weights¹ and assigned values given to vegetation cover types, distance from road categories, and slope (%) categories in estimating potential hunting pressure for northern bobwhite across the hunted unit of Beaver River WMA, Beaver County, Oklahoma, USA, 2012–2015. Values were derived from data presented by Richardson et al. (2008) where 1 represents the highest potential hunting pressure and 4 represents the lowest, and were incorporated into a weighted overlay analysis in ArcGIS 10.2.

Movement analysis

To compare estimates of covey movement across hunted and non-hunted units, we calculated a conservative estimate of average daily movement across the non-breeding season for coveys with ≥ 10 locations (Brøseth and Pedersen 2010). Coveys with ≥ 10 locations, rather than those with ≥ 30 locations, were used in movement analysis because we were not estimating KDEs for this stage of our analysis. We conservatively estimated average daily movement to be the Euclidean distance between a covey's locations across consecutive days (Williams et al. 2000, Brøseth and Pedersen 2010, Lohr et al. 2011, Unger et al. 2012), which we considered an index of daily movement patterns. The time between daily consecutive locations varied from 24 h, thus this metric was considered an index of average daily movement rather than true average daily movement. The median time between covey locations for average daily movement estimates were 25.07 h (range: 24 to 32.22 h), 25.5 h (range: 24 to 34.45 h), and 24.27 h (range: 24 to 29.13 h) for the 2012–2013, 2013–2014 and 2014–2015 hunting seasons, respectively. Because the time between consecutive covey locations varied from 24 h, we included velocity as the dependent variable in our movement analysis to account for this variability. Thus, velocity was considered a proxy for the conservative estimate of average daily movement. Linear mixed effect models (Pinheiro and Bates 2000) were used to assess the influence of hunted/non-hunted units, years, week time, and all possible additive and interactive combinations between these variables on covey movement (Brøseth and Pedersen 2010). To meet the assumption of data normality, we used a Box—Cox transformation (Box and Cox 1964) approach to determine the most appropriate transformation for our movement data. Based on this approach, we used x0.106 to transform our data. Covey identity was included as a random effect to account for interdependence of movement data within each covey (Brøseth and Pedersen 2010). We used an AIC_c approach, and used model weights (wi) and a ΔAIC_c ≤ 2 to determine the most parsimonious model (Burnham and Anderson 2002). Finally, we used a restricted maximum likelihood (REML) approach to obtain parameter estimates for fixed effects in our models (Brøseth and Pedersen 2010) and considered any parameters with 95% confidence intervals overlapping 0 to be non-significant in explaining average daily movement between coveys.

Table 2.

Model selection of Andersen—Gill hazard models of survival for northern bobwhite during the non-breeding season at Beaver River WMA, Beaver County, Oklahoma, 2012–2015.

Bobwhite survival

After censoring individuals, at total of 475 bobwhite were included in our survival analysis. This resulted in 192 (72 adults and 120 juveniles), 164 (74 adults and 90 juveniles) and 119 (69 adults and 50 juveniles) individuals during the 2012–2013, 2013–2014 and 2014–2015 non-breeding seasons respectively. A total of 342 (149 adults and 193 juveniles) and 133 (66 adults and 67 juveniles) individuals were located on the hunted and non-hunted units, respectively. For adults, mortalities were categorized as 28% mammalian, 25% unknown, 23% avian, 5% harvested and 1% other. For juveniles, mortalities were categorized as 30% mammalian, 18% avian, 16% unknown, 2% harvested and 1% other. Collar failures (slipped collars or dead batteries) accounted for 18% and 33% of bobwhite losses in adults and juveniles, respectively.

Based on AICc values, the global model was the most parsimonious model when explaining bobwhite survival in relation to anthropogenic features and disturbance during the non-breeding season (Table 2). However, of the variables in this model, only two were considered significant. These included an effect for the second year season (year two season hazard rate [HR] = 0.55, 95% confidence interval [CI] = 0.34 to 0.87) and distance to primary WMA roads (HR = 1.0008, 95% CI = 1.0003 to 1.0013). Individuals alive during the year two non-breeding season were 45% less likely to experience mortality compared to birds during the year one non-breeding season, while they were only 10% less likely to experience mortality when compared to individuals alive during the year three non-breeding season (Fig. 1A; year three season HR = 0.65, 95% CI = 0.38 to 1.11). Finally, every 10 m decrease in distance from primary WMA roads was associated with a 0.08% increase in probability of mortality.

Figure 1.

Northern bobwhite Colinus virginianus non-breeding season survival as determined from Andersen—Gill hazard models. Survival curves are separated by year (A) and our overall best performing model (B) for bobwhite on Beaver River WMA, Beaver County, Oklahoma, USA, 2012–2015. Week numbers correspond to the non-breeding season beginning on 1 October of each year. Vertical lines indicate the beginning and end of the quail hunting season in Oklahoma. Survival curves: panel (A): n = 192, 164 and 119 bobwhite for the 2012–2013, 2013–2014 and 2014–2015 seasons, respectively; panel (B): n = 475 bobwhite.

There were no differences in survival for individuals on hunted versus non-hunted units in our top model (non-hunted HR = 0.81, 95% CI = 0.49 to 1.35) nor in our model selection results (unit model ΔAIC_c = 6.47). Furthermore, based on our hazard rate curves, there is no indication that once the hunting season started, hazard rates for birds increased significantly (Fig. 1B). However, survival did consistently decrease across weeks for both hunted and non-hunted units during the non-breeding season with ∼20% of individuals surviving through the season (Fig. 1B).

Covey resource selection

Across all three years, a total of 28 coveys representing 62 radio-marked birds were used in estimating RUFs where locations occurred only during the quail hunting season. In general, there was little difference in space use of hunted versus non-hunted coveys in relation to our variables of interest. Furthermore, a majority of the variables included in our analysis had a non-significant relationship to covey space use based on () estimates with 95% confidence intervals that overlapped 0 (Table 3). Of all the variables analyzed, only the exposed soil/sparse vegetation and bare ground cover classes had significant differences between hunted and non-hunted coveys (Table 3) when compared to use of sand sagebrush. Coveys on the hunted unit significantly avoided the exposed soil/sparse vegetation cover type (() = -0.01, 95% CI = -0.02 to -0.002) and bare ground (() = -0.01, 95% CI = -0.02 to -0.002) when compared to non-hunted coveys. However, if pooled across all coveys (both hunted and non-hunted coveys), these relationships were not considered significant (exposed soil/spare vegetation pooled () = -0.005, 95% CI = -0.012 to 0.002; bare ground pooled () = -0.003, 95% CI = -0.012 to 0.005).

Covey movement

There were no differences in covey average daily movement between hunted and non-hunted units across all three seasons during our study (hunted unit β = 0.002, 95% CI = -0.01 to 0.03; non-hunted unit β = -0.01, 95% CI = -0.03 to 0.01). Time-related variables (week and year) were the best explanatory variables included in our analysis (Table 4). However, the parameter estimate for the week variable (β = 0.001, 95% CI = -< 0.001 to 0.002) was not significantly different from 0 and thus was not considered a strong explanatory variable for covey average daily movement. When compared to year one, only year three was significantly different based on 95% confidence intervals (β = -0.03, 95% CI = -0.05 to -0.01), indicating that average daily movement for coveys during year three was lower than years one and two.

Implications

Relatively low levels of harvest pressure appear to have no negative impact on bobwhite populations, as illustrated by our results. We were also not able to provide evidence that anthropogenic features increased hunting pressure across our study site. However, primary WMA roads appeared to increase mortality risk due to some unknown cause. We emphasize that low densities of anthropogenic features such as roads and oil/gas structures are compatible with bobwhite management within the context of landscapes already providing large areas of usable space. However, negative compounding impacts related to interactions between anthropogenic features and hunting pressure may exist in other regions with higher densities of developed features or higher levels of hunting pressure within the bobwhite distribution. Future studies should seek to empirically quantify hunting pressure to determine if there is a threshold in which these potential compounding impacts begin to negatively influence a population. As Williams et al. (2004) discussed, regional efforts should be made to assess whether anthropogenic development may be increasing hunting pressure so that harvest management is scaled appropriately based on local landscape configuration.