Harvest limits are an important component of wildlife management but evaluating whether realized harvest rates are close to established regulations is often not done in practice. Also, managers typically only partially control harvest rates via changes to metrics such as season length and hunting party size, but the relationships between these proxies and harvest rate are often unknown. We studied a population of northern bobwhite Colinus virginianus over a three-year period to estimate harvest, survival and relationships between harvest rate and proxies of hunting pressure. We captured and marked bobwhite with only bands (n = 479) or bands and radio-transmitters (n = 592). We monitored radio-marked individuals weekly until radio failure, lost signal or death. A comparison of annual survival rates indicated birds marked with only bands survived at similar rates (14.5%, 95% CI: 7.6, 21.4%) compared to those marked with both bands and radio-transmitters (10.6%, 95% CI: 7.9, 13.3%). A cumulative incidence function for cause-specific mortality produced an annual harvest rate estimate of 12.0% (95% CI: 8.6, 15.4%). However, this estimate was likely negatively biased in the presence of undocumented crippling loss. We also derived a 14.4% annual harvest rate using average observed levels of hunting pressure and a joint live–dead model that was not dependent on classification of mortality sources. We predicted annual harvest rates over a range of population sizes and hunting pressures to inform managers of possible thresholds to target. Hunting pressure was an important predictor of harvest rate. For instance, we estimated that adding approximately two more hunters per week would push seasonal harvest rates past the targeted harvest threshold, indicating managers should carefully regulate hunting pressure to avoid exceeding set harvest limits. Estimating relationships between weekly survival and hunting pressure can provide information to evaluate different management scenarios and help develop hunting regulations.
Harvest guidelines are used to ensure game populations are exploited sustainably. Gamebirds are commonly managed by setting limits on the number of individuals that can be harvested per hunter during a defined period of time (Ellison 1991). The goal is to limit harvest to a percentage of the population that will not incur additive effects on natural mortality (Pollock et al. 1989a, Rolland et al. 2010). Directly estimating the threshold above which harvest mortality becomes additive to natural mortality requires experimental studies that manipulate harvest rates (Nichols et al. 1984) or long-term studies of marked populations subjected to varying levels of annual harvest (Burnham and Anderson 1984, Sedinger et al. 2010). The results of such studies are commonly used to set harvest regulations, yet realized harvest rates are rarely evaluated. Furthermore, understanding how hunting pressure contributes to variation in survival within the hunting season is important because harvest rates are not easily manipulated, i.e. harvest is only partially controllable through management actions (e.g. number of hunts per season, number of hunters and bag limits).
Hunted populations that are comprised of marked individuals offer opportunities to directly estimate harvest rates. Information on individual fates can be collected using electronic devices (e.g. very high frequency [VHF] transmitters), band-recovery data reported by hunters, capture–recapture data from live reencounters, or some combination (Williams et al. 2002). These data allow for the application of time-to-event models (White and Garrott 1990), band-recovery models (Brownie et al. 1978) and joint band-recovery and capture–recapture models (Burnham 1993). These models provide different information regarding the influence of hunting on wildlife populations; time-to-event models are practical for direct estimates of harvest rates (Robinson et al. 2014), while band-recovery and capture–recapture models can provide ways to examine hunting dynamics and their influence on survival through the incorporation of time varying covariates.
Careful evaluation of model assumptions and parameter interpretation can help quantify bias in survival estimates. For example, there is an explicit assumption that marks (e.g. radio transmitters) do not affect the vital rate(s) of interest, otherwise estimates produced will directly reflect bias introduced by the marker used (Winterstein et al. 2001). This assumption has been previously investigated for several VHF transmitter studies of gallinaceous species, with most studies of grouse concluding that radio-collar transmitters do not influence survival (Thirgood et al. 1995, Hagen et al. 2006). However, some have suggested that radio-marks might produce a severe negative bias when applied to smaller gallinaceous species, such as quail (Osborne et al. 1997, Cox et al. 2004). Two prior studies found no support for a radio-transmitter effect on quail (Palmer and Wallendorf 2007, Terhune et al. 2007). Nonetheless, Terhune et al. (2007) recommended the need for additional investigation given local conditions may interact with transmitter marks to influence survival.
Here, we studied a species that is regularly subjected to exploitation by hunting, the northern bobwhite Colinus virginianus (hereafter, bobwhite) – a gallinaceous gamebird distributed throughout the eastern half of the United States. Bobwhite populations have been in decline at the continental scale since at least the 1960s (Hernandez et al. 2013), despite management efforts throughout much of their range (Williams et al. 2004). There are general guidelines to help managers set harvest policies, such as focusing hunting over short periods early in the winter and harvesting a fixed percentage of the fall population (Morgan et al. 2016). Annual harvest rates are largely dependent on geographic location and are variable. For example, Morgan et al. (2016) reported that state agencies targeted harvest rates for bobwhite between 20 and 40%, with higher rates generally occurring in the northern portion of the bobwhite range. In the southeastern United States, harvest targets are often set at 15% because this rate has been suggested to avoid additive mortality (Sisson et al. 2017). Consequently, a 15% harvest rate is currently considered a conservative policy to target by some for sustainable population management in the southeast.
We assessed the effects of hunting on a population of bobwhite at Di-Lane Wildlife Management Area in Georgia, USA. Our first objective was to determine if Di-Lane was meeting target harvest rates of 15%. For this first objective, we produced estimates from known-fate models and joint live–dead models using radio-marked birds. Our second objective was to estimate annual cause-specific mortality probabilities, including both harvest and non-harvest sources of morality, and whether cause-specific hazard rates changed through the year. Our third objective was to inform management at Di-Lane and similar sites by estimating effects of hunting pressure (i.e. the numbers of hunters, hunter hours and dogs that were present during hunts), sex and season on survival. Our fourth objective was to test for an effect of radio transmitters on survival to understand if our estimates based on radio-tagged birds were likely to be biased. Taken together, these research objectives allowed us to directly assess if hunting management guidelines were being met, identify the primary sources of mortality for the population, test an important model assumption, and predict the weekly and annual harvest rates likely to occur given changes in hunting pressure.
Material and methods
Study area
We conducted our research on the 3278 ha Di-Lane Wildlife Management Area located in Burke County in the Coastal Plain of east-central Georgia. Approximately 2023 ha at Di-Lane are managed for bobwhite. Forest stands are typically treated with prescribed fire during the dormant season and early growing season to stimulate herbaceous growth and maintain canopy openings. Forest stands at Di-Lane are primarily composed of mixed hardwood forests (oak species – Quercus spp., sweetgum – Liquidambar spp.) and loblolly pine Pinus taeda (Moseley et al. 2003). The site is generally flat with little topographic variation. Wildlife on Di-Lane are managed by the State of Georgia's Wildlife Resources Division. Management includes predator control through targeted removal of meso-mammals (primarily raccoons Procyon lotor and opossums Didelphis virginiana). Additionally, supplemental feeding occurs at various locations within the managed area.
Bobwhite are hunted on several designated dates during December and the first week of February every Wednesday and Saturday (10-hunt days total over the winter period falling within 6–7 weeks). The daily quotas during our study were 6 bobwhite per person and 12 per party (1–3 hunters per party), with a maximum of 24 hunters allowed per hunt day. There were no reserve areas at Di-Lane, other than a small safety zone located around the property headquarters. Hunters were required to check-in prior to hunting and check-out after they finished hunting at a single check station, regardless of whether their hunt was successful. Exit surveys conducted by Georgia's Wildlife Resources Division occurred following hunts, which provided information to calculate metrics of hunting pressure, including the number of hunters, dogs used and hours hunted.
Capture and monitoring
We monitored bobwhite on Di-Lane from October 2016 through April 2019. Fall trapping occurred primarily during early October – early November for a total of 26 daily fall trapping occasions each year. Spring trapping primarily occurred from late February to early March, for a total of 14 daily spring trapping occasions each year. We captured bobwhite using ‘walk-in’ funnel traps baited with sorghum (Stoddard 1931). We set between 127 (Spring 2019) and 396 (Fall 2016) traps within dense cover across the study site. Traps were spaced an average of 146m (Fall 2016) and 397m (Spring 2018) apart. This configuration met our targeted trap spacing of no more than 400m to optimize spatial coverage across the property and meet operational constraints for checking this number of traps.
We weighed captured individuals to the nearest gram and banded them with a unique numbered leg band (National Band & Tag Company, Newport, Kentucky). We assigned the sex and age of captured individuals based on plumage characteristics and the shape of the outer primary feathers (Petrides and Nestler 1943). VHF radio transmitters equipped with 8-h motion delayed mortality signals (Holohil Systems, Carp, Ontario, Canada and American Wildlife Enterprises, Monticello, FL, USA) were fit to a subset of captured bobwhite, which we refer to as birds belonging to the radio-marked group. Birds that did not receive radio-marks are referred to as belonging to the band-only group. Birds remained in their respective group throughout the study (i.e. birds did not have their mark type changed if subsequently recaptured). Transmitters weighed ∼6 g and individuals that were ≥130 g (<5% body weight) and deemed healthy (based on physical appearance and behavior) received a radio-mark. We tested for differences in mass between the radio-marked and band-only groups using a two-sample t-test. The estimated battery life on each transmitter was 10–12 months. We monitored bobwhite approximately daily during weekdays and some weekends to determine their locations and live–dead status via homing telemetry (White and Garrott 1990). We maintained ≥30 live birds with radio-marks during each week of the study.
We attempted to identify the cause of death for all individuals with a mortality signal that we were able to locate. Hunting mortalities were recovered and reported directly by hunters. We assigned birds that were not recovered by hunters but subsequently found dead and intact (i.e. not scavenged) within three days post-hunt as hunter mortalities, similar to methods used by Haines et al. (2009). Predation events were classified as either unknown, avian, mammalian or snake. We assumed mortalities classified as unknown were primarily from predation since human related causes outside the hunting season rarely occurred. We based these classifications on the condition of the radio-collar (e.g. presence of bite marks, bent antenna, nearby tracks), condition of the bird (e.g. a pile of feathers versus scattered remains), or the location of the transmitter (e.g. burrow, ground, tree). Unknown classifications occurred when no clear evidence of the mortality cause was present (e.g. only a few remaining feathers or wing found near collar). Snake classification was straightforward as the transmitter was typically consumed along with the carcass.
Known-fate analysis
We fit time-to-event models to data from radio-marked birds to estimate harvest rate (objective 1), cause-specific hazard rates and whether cause-specific hazard rates varied over the 52-week annual cycle (objective 2).
Before fitting models to data, we standardized the week number so that week 1 began on 1 October and week 52 ended on 30 September. A subset of bobwhite was monitored in >1 year, in which case the bird received a separate row entry in the data file for each year it was alive. The preceding year was then censored at week 52, and the first entry for the following year was set to week 1 (following Sandercock et al. 2011). We left-truncated data because individuals entered the at-risk group during different weeks (Pollock et al. 1989b). We right-censored data in cases of transmitter failure, dead transmitter batteries or birds that were still alive at the end of the study. If a bird went missing for more than a week before subsequently being found dead, we assigned the date of the mortality event as the midpoint between the date last observed alive and the date found dead (Mayfield 1975), with the exception of known mortality dates for reported harvested birds. Band number was unique to each bobwhite, and it was used as a random effect to account for lack of independence of individuals occurring in multiple years (i.e. having multiple data rows). Differences between survival of males and females were assessed using Cox proportional hazards models and the Wald Z-test. All analyses were conducted using the survival package in R (ver. 3.4.4; < www.r.project.org>).
We calculated cause-specific mortality using the competing-risks nonparametric cumulative incidence function estimator (NPCIFE) described by Heisey and Patterson (2006). We implemented the NPCIFE using the wild1 R package (Sargeant 2011), which allowed for left-truncation and right censoring. We partitioned mortality sources in two ways. First, we summarized mortality events into two competing risks, and estimated cumulative hunting mortality (i.e. harvest rate; coded as ‘1') and cumulative natural mortality (all sources of mortality other than hunting pooled and coded as ‘2'). Second, we summarized mortality events into five competing risks, and estimated cumulative mortality for avian predation (‘1'), mammalian predation (‘2'), snake predation (‘3'), hunter harvest (‘4') and unknown causes (‘5'). In both cases right censored individuals were coded as ‘0'. We tested for differences in cause-specific mortality rates between males and females by fitting Cox proportional hazards models with a sex effect to each mortality source separately (i.e. mortality source of interest coded as ‘1’ and all other records censored as ‘0'). Again, the Wald Z-test was used to assess statistical significance for each model.
To quantify how hazard rates changed throughout the annual cycle, we modeled the smoothed instantaneous hazard functions separately for hunting mortality, cause-specific natural mortality and all sources of natural mortality combined. We used R and the gss package to estimate smoothing splines of the hazard rates across the annual cycle (‘gss' package, Gu 2014), as demonstrated by Sandercock et al. (2011).
Joint live–dead analysis
We estimated survival on the weekly scale (radio-marked birds) and between fall and spring capture periods (radio-marked and band-only birds) by fitting the joint live–dead model developed by Burnham (1993) to provide additional estimates of harvest rate (objective 1). This model also provided a way to estimate weekly survival as a function of sex, season and hunting pressure (objective 3). In this way, we were able to quantify how changes in metrics that managers can control (i.e. the number of hunters, hunter hours and dogs that were present during hunts) would affect bobwhite survival rates. Additionally, the joint live–dead model provided a way to test for evidence of a transmitter effect on survival between capture occasions (objective 4).
Model description
After initial capture, bobwhite could be reencountered alive through radio-tracking (radio-marked birds only) and recapture, or dead through radio-tracking or hunter-harvested recoveries. We compiled capture, monitoring and recovery data for bobwhite into encounter histories using a live–dead (LD) coding to represent the status of a bird during each encounter period (White and Burnham 1999). The joint model includes probability parameters for survival (S), live reencounter (p), dead encounter (r) and fidelity (F). The definitions for probabilities S and p are standard from general capture–recapture models and include the probability of survival between interval j and j+ 1 (S), and the probability of being reencountered alive at occasion j given a marked individual is alive in the study area and available for reencounter (p). The dead encounter probability r is defined as the probability that an animal that died between occasion j and j+ 1 is both found and reported. Reporting rates at Di-Lane were assumed to be close to 100% because Georgia Department of Natural Resource's personnel administered surveys during hunting days and hunters were required to report harvest and effort before leaving the area. Hence, between hunting periods for the band-only group, r is an estimate of the proportion of birds that died and were retrieved as a result of hunting, and 1 – r represents the proportion of mortality due to natural causes or unretrieved hunter kills. The parameter F is defined as the probability a bird at risk of capture at occasion j is again at risk of capture at time j+ 1. We used program MARK (White and Burnham 1999) and the R package RMark (Laake 2013) to fit the joint live–dead model to encounter histories.
Weekly survival and derived harvest rate
We constructed encounter histories for radio-marked birds on the weekly scale to estimate survival as a function of sex effect and time-varying covariates. We used Akaikie's information criterion for small sample sizes (AICc) to select the most parsimonious model structures (Burnham and Anderson 2002) based on the following two-step procedure. First, we reduced the number of possible models to compare by selecting only a limited set of structures for parameters r, p and F (i.e. the ‘nuisance’ parameters). We expected a priori that r and p would vary temporally due to variable tracking effort throughout the year, so we allowed weekly variation by year-month (31-level factor, average weekly parameter estimate constrained to be same within each unique year-month period). We fixed F to 1 because there was no evidence our radio-marked birds permanently emigrated. Model structures for S included a null effect (i.e. intercept only), hunting occasion (‘hunt’, 2-level factor covariate; 0 = non-hunting week and 1 = hunting week), season and sex (2-level factor, 0 = female, 1 = male). Covariates for S were considered individually or in additive combinations with other covariates.
For the second stage of model selection, we carried the best base model structure of S forward to evaluate hunting pressure covariates on weekly survival. If the hunt covariate was in the top ranked model, it was replaced with time-varying hunting covariates. These included the number of hunters (‘nhunters', continuous covariate), the cumulative number of hunter hours (‘nhours', continuous covariate), and the number of dogs used during hunts (‘ndogs', continuous covariate). We divided each covariate by the fall population size, which was estimated from a spatial capture–recapture model applied to the Di-Lane population during each year of our study (Howell et al. unpubl.). Scaling the covariates in this way allowed us to interpret changes in survival as a function of changes in the covariates per bird in the fall population. We also retained the best base structure (i.e. {S(season + hunt)}) which produced four candidate models. Again, we used AICc to rank candidate models and choose the most parsimonious for inference on how hunting pressure affects weekly survival rates. We used the ‘RMark' package (‘covariate.predictions’ function) to model average seasonal survival predictions for each group across the entire candidate model set for the purposes of presenting estimates while accounting for model selection uncertainty. We used a function in RMark (‘deltamethod.special') to obtain standard errors of the product of weekly survival estimates to calculate 95% confidence intervals for annual survival (i.e. the delta method).
We used coefficients estimated from the joint live–deadsurvival models {S(season+nhunters)} and {S(season+nhours)} to predict annual harvest rates (ha). We predicted ha under different management scenarios by varying levels of average weekly number of hunters (nhunters=0–40) and abundance prior to the start of the hunting season, i.e. fall population size during week i (N[i=1]=800–1600). We produced 400 predicted annual harvest rates for different combinations of number of hunters and population size. The annual harvest rate depended on the cumulative effects of harvest during each week i of the season (hw[i]). Weekly harvest rate was estimated as the difference between the predicted weekly survival rate with hunters equal to 0 and hunters >0:
The number of birds harvested (H[i]) during each week of the hunting season was then calculated as N[i]hw[i], and the population at i+1 was calculated as N[i]–H[i].
Thus, if n represents a vector of weekly indices (i) during which hunting occurs, the derived annual harvest rate ha is: 1 – N[i = last(n) + 1]/N[i = 1]. A benefit of this derived rate is that it should not be susceptible to negative bias in the presence of incorrect classification of cripples attributable to other causes, because the joint live–dead model does not distinguish between sources of mortality. However, it may overestimate the annual harvest rate when losses between hunting weeks (i.e. non-hunting weeks are interspersed with hunting weeks during the hunting season) are not subtracted from population abundance in the preceding time step. This bias can be accounted for by multiplying the weekly harvest rate hw[i] by the predicted proportion of surviving birds between each hunting week, i.e. logit–1(β0 + βseason)I, where I is the number of weeks between hunting periods. We did not develop a variance estimator for the derived harvest rate. Thus, only point estimates are presented.
Survival between capture periods and transmitter effect
We tested for the presence of a transmitter effect by comparing survival rates of radio-marked and band-only birds between capture periods using the joint live–dead model. We lacked data for the band-only group to make this comparison at the weekly scale because recoveries and live reencounters were not possible during all weeks. Unlike weekly survival, intervals between the mid-dates of fall and spring capture occasions were roughly five and seven months in length, respectively. Thus, we constructed new encounter histories with six encounter occasions (i.e. season each year) that also included a band-only group. Similar to our analysis of weekly survival, we reduced the number of possible models by limiting the number of structures for parameters r, p and F, and evaluating different structures of S. We expected a priori that parameters r and p should differ between mark groups. Therefore, structures for r and p included a mark type effect (‘mtype', 2-level factor covariate; 0 = band-only and 1 = radio-mark). We also included a variable time effect (‘time', 6-level factor allowing unique differences between each capture occasion interval) for r, and a capture season effect (‘capseason', 2-level factor allowing differences between fall and spring capture periods) for p. For the band-only group, we fixed r to 0 during periods between spring and fall because recoveries could not occur then. We fixed F to 1 for both mark groups because no radio-marked birds permanently left the study area and no bands were recovered outside the managed hunting area. We considered several possible model structures for S, including separate models for mtype, sex, capseason and time, and additive combinations of each. The same method applied to weekly survival estimates was used to model average survival predictions between capture occasions for each group in the candidate model set. The delta method and RMark were again used to derive annual survival rates for birds in the band-only and radio-marked groups.
Results
We analyzed data from 1071 bobwhite captured at Di-Lane from October 2016 through April 2019 (Table 1). Our sample included slightly more males than females (560 versus 511, respectively) and more juveniles than adults (782 versus 289, respectively). Of those captured, 479 received bands, and 592 received both bands and a VHF radio collar. The average mass at first capture of birds in our sample did not differ between those in the band-only group (x = 167.2 g, SD = 17.3 g) and radio-marked group (x = 166.7 g, SD = 14.1 g; two-sample t-test: df = 1008, Tcrit = 2.0, p = 0.6).
Table 1.
Summary of captured northern bobwhite Colinus virginianus at the Di-Lane Wildlife Management Area in Burke County, Georgia. The numbers of newly captured birds each year are presented by category. Categories are marker type (band-only or radio-marked), age and sex.
Hunting pressure was highly variable among weeks during the study. Hunters spent an average of 2.2 h hunting per hunt (SD = 0.4). The average number of hunters was 14 per week (SD = 5.1), and the average number of hunting parties was 5.6 per week (SD = 2.0). Nearly all hunting parties used dogs (there were only three recorded instances of a party hunting without them), with an average of 18.2 dogs used per week (SD = 8.3). Over the study period hunters harvested 443 birds, 119 of which were marked. Sex and age classifications of marked harvested birds indicated juvenile males were most frequently harvested (39%), followed by juvenile females (27%), unknown sex juveniles (20%), adult males (8%) and adult females (5%).
Known-fate model: cause-specific mortality
Over the study period we right-censored 120 birds due to lost radio contact (potentially due to radio failure, dead batteries, undocumented permanent emigration or radio destruction during predation events). Score tests from Cox proportional hazards models indicated the hazard rate of males was slightly higher than females (hazard ratio [HR] sex = 1.221; 95% CI: 1.00, 1.50; Z = 4.21; p = 0.06). The assumption of proportional hazards was met (χ2 = 0.22; p = 0.64).
Sources of mortality estimated from the competing risk models indicated unknown causes were the biggest contributor (0.256; 95% CI: 0.205, 0.307), followed by avian predation (0.245; 95% CI: 0.190, 0.299), mammalian predation (0.174; 95% CI: 0.114, 0.234), hunter harvests (i.e. harvest rate, 0.120; 95% CI: 0.086, 0.154) and snake predation (0.031; 95% CI: 0.000, 0.063) (Fig. 1). Wald Z-tests produced from Cox proportional hazards models applied to cause-specific mortality sources indicated no differences in the hazard rates between sexes for any of the mortality sources (all p>0.1). Therefore, neither males nor females were more or less susceptible to any of the classified causes of mortality.
Smoothed instantaneous hazard rates indicated the risk of mortality from hunting peaked in December and February. The steepest increase in instantaneous hazard occurred from mammals in late August through September. However, there was a gradual increase for all other mortality sources from October through September (Fig. 2).
Figure 1.
Cumulative incidence functions of different sources of mortality for northern bobwhite Colinus virginianus studied at the Di-Lane Wildlife Management Area in Burke County, Georgia. Solid lines represent the predicted step function and dashed lines represent the associated 95% confidence intervals. Natural mortality excluded hunter harvest but may over estimate its effects in the presence of a compensatory mechanism.
Joint live–dead model: weekly survival, derived harvest and transmitter effect
The base structure for weekly survival of radio-marked birds receiving the most model support included effects for season and hunting occasion {S(season + hunt)} (ω = 0.64, Supplementary material Appendix 1 Table A1). Model-averaged predictions indicated average weekly survival probability was lowest during winter hunting periods followed by summer and winter periods, respectively (Fig. 3a). Model-averaged predictions indicated no influence of sex on survival (Fig. 3a). The derived mean annual survival probability of radio-marked birds taken from weekly survival estimates was 0.108 (95% CI: 0.084, 0.132).
In the second stage of model selection the number of cumulative hunter hours was the best supported time-varying hunting metric (Table 2), receiving the majority of AICc weight (ω = 0.61; Fig. 4a). The second and third ranked models included the number of hunters (Fig. 4b) and the number of dogs used (Fig. 4c), respectively. The base structure received no support and had a much higher AICc score (ΔAICc = 19.1), indicating that time-varying hunting covariates substantially improved explanatory power. Encounter probability estimates for parameters r and p from the top model varied considerably by year-month (Supplementary material Appendix 1 Table A2).
Figure 2.
Smoothed instantaneous hazard rates for different causes of mortality of radio-marked northern bobwhite Colinus virginianus studied at the Di-Lane Wildlife Management Area in Burke County, Georgia. Natural mortality included the combined effects of avian, mammalian, snake and unknown sources of mortality.
Our derived annual harvest rates were slightly higher than the mean estimate from the analysis of radio-marked birds only, but within the confidence bounds. A mean annual abundance of 1199 birds and a mean annual accumulation of 327.6 hunter hours (46.8 hunter hours per week times seven weeks) produced a predicted annual harvest of 14.5%. Similarly, using the same abundance and a mean annual accumulation of 149.1 hunters (21.3 hunters per week times seven weeks) produced a predicted annual harvest of 14.4% (Fig. 5).
The best structure for seasonal survival (between fall and spring and spring and fall capture periods) included a capture interval effect (ω = 0.36; Table 3). The second ranked model included the addition of a marker effect (relative to the top-ranked model), which increased its ΔAICc by 0.83 (ω = 0.24). Similarly, the third ranked model included the addition of a sex effect (relative to the top-ranked model), which increased its ΔAICc by 1.73 (ω = 0.15). The accumulated weights for models that included a mark type effect was 0.321 compared to 0.745 for models including capseason effect. Model-averaged survival predictions indicated that while mark-type and sex were included in top models, only capture interval appeared to influence survival (Fig. 3b). A comparison of the derived annual survival probabilities between mark type groups indicated they were similar (band-only: 0.145, 95% CI = 0.076, 0.21; radio-marked: 0.106, 95% CI = 0.079, 0.133). A comparison of the derived annual survival probabilities between males (10.8%, 95% CI: 7.5, 14.1%) and females (11.9%, 95% CI: 8.6, 15.2%) also indicated similar survival. Encounter probability estimates for parameters r and p from the top-ranked model are presented in Supplementary material Appendix 1 Table A3.
Figure 3.
Model-averaged survival probability estimates for a northern bobwhite Colinus virginianus population at the Di-Lane Wildlife Management Area in Burke County, Georgia. Joint live–dead models were fit to radio-marked birds to obtain average weekly survival estimates during three seasonal periods (a). Survival comparisons were also made between capture periods for birds marked with very high frequency (VHF) transmitters and birds that received only bands (b). Capture periods occurred twice each year in the fall and spring, and survival was estimated between these periods (i.e. ‘spring' in the figure corresponds to survival from spring capture to fall capture, and ‘fall’ correspond to fall capture to spring capture).
Discussion
Marked populations exposed to hunting provide important opportunities to critically evaluate harvest guidelines. We studied a population of bobwhite in east-central Georgia using band recovery data and radio telemetry monitoring. Realized harvest rates were close to harvest limits set in the management area's hunting guidelines. We also found that weekly survival was strongly correlated with several measures of hunting pressure, indicating that managers should carefully monitor hunting pressure to avoid overharvest.
Harvest rate and hunting pressure
The mean annual harvest rate based on radio-marked birds monitored at the Di-Lane Wildlife Management Area was 12%, produced from the known fate competing-risks nonparametric cumulative incidence function estimator. This estimate was similar to a harvest rate of 9% reported by Terhune et al. (2007) in a southwest Georgia population, but substantially lower than estimates produced in other studies, which ranged from 17% (Pollock et al. 1989a) to 48% (Cox et al. 2004). We consider our estimate from radio-marked birds to represent a minimum realized harvest rate, because in the case of crippled birds that were not retrieved by hunters, its accuracy depended on those crippled birds being correctly classified as hunter harvested. Any crippled birds that were subsequently scavenged by wildlife prior to retrieval by our field technicians may have been incorrectly classified as predator killed. A worst case scenario for bias due to undocumented crippling loss might be if all mortalities attributed to unknown causes were in fact the result of hunter-caused injuries. The cumulative incidence function for unknown causes indicated that 5.8% of mortalities occurring during hunting weeks were classified to an unknown cause. We also included one additional week after hunting ended to this period to account for potential delays in mortality attributed to hunting. Under this scenario, 17.8% (12% + 5.8%) of annual mortalities would be attributed to hunting. Birds included in our time-to-event models that were lost to radio-telemetry monitoring (i.e. signals disappeared) were censored, but future studies might consider modeling endpoint hazards using hierarchical models to infer the proportion of unknown mortality attributed to cripple loss, while also accounting for the possibility that censored individuals were lost due to this cause (Stenglein et al. 2018).
Table 2.
Model selection results for comparison of hunting pressure covariates on northern bobwhite Colinus virginianus weekly survival at the Di-Lane Wildlife Management Area in Burke County, Georgia. All covariate models build off the best-supported base structure model {S (season + hunt)p(year – month)r(year – month)F(1)} by replacing ‘hunt’ with the hunting pressure covariate. Covariates were the cumulative number of hours hunted each week (nhours), number of hunters (nhunters) and number of dogs used (ndogs). Each covariate was scaled by fall density and is represented on the logit scale. Effects of season (1 = winter; 0 = summer) and hunt (1 = hunt week; 0 = no hunt week) were represented as dummy variables. The number of parameters (K), ΔAIC (difference between AICc of 17799.3 for top-ranked model and current model), model weight (ω) and deviance are reported for each model.
The annual harvest rate based on the joint-live dead analysis of radio-marked birds (14.4%) was higher than our analysis of radio-marked birds using the competing-risks nonparametric cumulative incidence function estimator, although within its 95% confidence interval (9–16%). Survival during the hunting period was strongly responsive to several hunting pressure metrics (Fig. 4a–c). For our model set, the cumulative number of hunter hours was the most well supported hunting pressure effect, but it is probably easier for managers to control the number of hunters allowed per hunt than hours hunted. If we assume an abundance of 1199 bobwhite (the mean abundance estimated at Di-Lane across the three hunting seasons, Howell et al. unpubl.) and observe an average of 21.3 hunters per week, the weekly harvest rate is expected to range from 2.3 to 2.9% of the population. However, if the number of weekly hunters averaged 23, the predicted annual harvest rate after seven weeks would equal 16.1%, exceeding harvest guidelines. Naturally, as population abundance increases, the predicted harvest rates will be lower for the same number of hunters (Fig. 5). Similarly, reducing the number of hunting weeks (or hunts that occur within a week) will reduce the predicted annual harvest rate. For example, following from the above with 21.3 hunters and a population of 1199 bobwhite, hunting only six weeks lowers the predicted annual harvest rate to 12.0%. The joint-live dead modeling approach offers the benefit of being unbiased with respect to crippling loss since the model does not depend on cause-specific classifications of recorded mortalities. Indeed, higher mean predicted annual harvest rates obtained from this model were consistent with our expectations. In short, if fall abundance is estimated, managers could use the model coefficient estimates and check station/hunter exit survey data to predict how the number of hunters in a given week influences overall harvest rate and make decisions about harvest regulations on a weekly time scale. This may provide opportunities for managers to use weekly hunting pressure information to actively manage the population within a hunting season. We caution that our models assume a linear relationship between effort and kill rate but non-linear or threshold efforts may exist at low or high densities.
Management at Di-Lane targets 15% of the fall population size as the maximum allowable harvest. We did not assess the appropriateness of this targeted harvest rate or whether additional harvests above this threshold will lead to additive mortality, and we acknowledge its usefulness will depend on its accuracy. Regardless, from the standpoint of our modeling approach, the actual defined harvest limit is arbitrary; the realized harvest rates we calculated provide a useful comparison to any targeted rates that are determined based on sound management principles. Therefore, these methods are likely to be a useful management tool whenever paired with targeted harvest rates that do not lead to additive mortality.
Survival and cause-specific mortality
Our annual survival estimate derived from weekly survival for radio-marked bobwhite (10.6%) was within the range of other published bobwhite studies (Brennan et al. 2020, Sandercock et al. 2008). However, our seasonal and sex-specific survival estimates generally differed from those reported in the literature. Extrapolating our average weekly winter (including hunting periods) and summer survival estimates – pooled across demographic classes – over each period would produce a winter survival rate of 40% (95% CI: 36, 55%) and summer survival rate of 28% (95% CI: 23, 33%). These findings contradict most bobwhite studies reporting estimates of seasonal survival, which indicate most mortality occurs during the winter period. For example, Sandercock et al. (2008) compiled an extensive summary of studies reporting seasonal survival for many bobwhite populations. Based on those summaries (presented in Table 4 and 5 of Sandercock et al. 2008), the average winter and summer survival rates were 26% and 39%, respectively. It is worth noting that previous average seasonal survival estimates reported from studies in Georgia, however, were similar (37% winter versus 42% summer) (Sandercock et al. 2008). Our derived annual sex-specific survival estimates were similar for males (10.8%) and females (11.9%). These results are in contrast to differing annual survival estimates reported for males (19%) and females (14.3%) in a previous study (Pollock et al. 1989a). However, Burger et al. (1995) found no evidence for sex-specific differences in annual survival. Therefore, reports in the literature appear to be equivocal.
Figure 4.
The relationship between the estimated weekly survival probability of northern bobwhite Colinus virginianus and different measures of hunting pressure at the Di-Lane Wildlife Management Area in Burke County, Georgia. The black line represents the mean predicted survival and light gray shading represents the 95% confidence interval of the mean predictions. Covariates included number of cumulative hunter hours (a), number of hunters (b) and number of dogs used by hunters (c). Covariates were divided by the fall abundance estimate for each hunting season to scale per bobwhite in the population. All probabilities are on the weekly scale.
Figure 5.
Relationships between fall population size and number of hunters permitted to hunt northern bobwhite Colinus virginianus at the Di-Lane Wildlife Management Area in Burke County, Georgia, over three hunting seasons (2016–2019). Contour lines represent derived annual harvest rates predicted from the covariate model {S(season + nhunters)}. Red vertical and horizontal lines represent the average observed number of hunters per season, and fall population abundance, respectively.
Predation was responsible for most recorded mortality events during our study. Based on our natural mortality classifications, unknown predation accounted for the majority of predator kills, followed by avian, mammalian and snake predation, respectively. These results are consistent with previous studies conducted on bobwhite using radio telemetry. Burger et al. (1995) attributed 29% of mortality to avian predation, and 26% to mammalian predation. Madison et al. (2002) attributed 55% of overwinter mortality to avian and mammalian predators. Our smoothed hazard rates indicated that the instantaneous hazards for all cause-specific natural sources increased from October through the end of September (Fig. 2). Thus, rising predation throughout the breeding season and early fall appears to be the driver behind lower summer survival. Past work has implicated avian predation as a source of high predation during the fall period when raptor migration commences (Holt et al. 2012, Turner et al. 2014), but our work indicates mammalian predation rates have the steepest rise during the summer to fall period (Fig. 2), similar to results reported in Burger et al. (1995).
Table 3.
Model selection results for comparison of survival between capture events of northern bobwhite Colinus virginianus studied at the Di-Lane Wildlife Management Area in Burke County, Georgia. Effects of capseason (0 = fall to spring, 1 = spring to fall), mtype (0 = band-only, 1 = radio-marked) and sex (0 = female, 1 = male) were coded with dummy variables. Additional model structures included time (unique survival between each capture interval) and an intercept-only model (null). All models included the same structures for live reencounter {p(mtype + capseason)}, dead encounter {r(mtype + time)} and fidelity {F(1)}. The number of parameters (K), ΔAIC (difference between AICc of 3892.9 for top-ranked model and current model), model weight (ω) and deviance are reported for each model.
Transmitter effect
The accuracy of our estimates of weekly survival probability were dependent on the assumption that the transmitters used did not influence this rate. Prior studies on bobwhite (Palmer and Wallendorf 2007, Terhune et al. 2007) did not find evidence to refute the validity of this assumption. Our results indicated birds belonging to the band-only group had similar annual survival (14.5%, 95% CI: 7.6, 21.4%) compared to radio-marked birds (10.6%, 95% CI: 7.9, 13.3%). Thus, we did not find strong evidence for a transmitter effect, consistent with previous bobwhite studies examining marker effects. We did not attempt to evaluate other possible influences of transmitters on vital rates other than survival, such as reproductive rates (e.g. nest and brood survival) or behavior; this was beyond the scope of our work. As discussed by Terhune et al. (2007), our findings should be viewed in the context of our population, which was provided supplemental food. If supplemental food alters the ability of bobwhite to withstand a transmitter that it otherwise could not tolerate (e.g. improved body condition, Corteville et al. 2000), then our results should not be extended to other populations given such a mitigating effect. To our knowledge, transmitter effects have not been evaluated for bobwhite populations where supplemental food is absent.
Conclusions
The development of hunting guidelines is challenging, even for populations that are highly regulated and monitored. Herein we demonstrate how evaluating realized harvest rates can be beneficial in the absence of long-term marking programs or large-scale experiments, given an informed and targeted harvest limit exists. We found that realized bobwhite harvest rates at Di-Lane Wildlife Management Area in east-central Georgia were within targeted guidelines and that hunting pressure strongly influenced survival and annual harvest. We provide examples of how our models can be used to predict a threshold of hunting pressure after which harvest is predicted to exceed targeted rates. This approach may be useful in other exploited populations where animals are marked.
Acknowledgements
We thank Todd Arnold for providing suggestiosn that greatly improved the quality of this manuscript. We also thank the field technicians and graduate students: Jamie L. Goethlich, Melissa L. Maleckar, Ravin R. Grove, Collin J. Richter, Joseph D. King, Lacey L. Clarke, Jake G. Spears, Jessica L. Mohlman, Rachel Gardner and Nathan Wilhite.
Funding – This study was funded by Georgia Department of Natural Resources and Warnell School of Forestry and Natural Resources for the financial support over the course of this project (grant no. W-GA-F16AF00189).
Author contributions – GTW analyzed the data, JMY and PEH quality checked the program code used in the analysis, and JAM and IBP conceived of the study design and oversaw data collection. GTW, JMY, PEH and JAM contributed equally to writing.
References
Appendices
Supplementary material (available online as Appendix wlb-00667 at < www.wildlifebiology.org/appendix/wlb-00667>). Appendix 1.
