Estimating sex-ratio, survival, and harvest susceptibility in greater sage-grouse: making the most of hunter harvests

We analyzed banding (3259) and recovery (six years) data from a hunted population of greater sage-grouse Centrocercus urophasianus in northwestern Moffat County, Colorado to examine vulnerability to hunter harvests and annual survival of adult and hatch-year (juvenile) birds. Additionally, we combined the recovery data with hunter-harvested wings and applied the Lincoln estimator to provide unbiased estimates of tertiary sex-ratio. Our results yielded the following findings: juveniles were harvested at twice the rate of adults, but harvest vulnerability was similar between adult males and females. Annual survival of juveniles was highly variable but similar between sexes. Sex ratios of adults and juveniles largely conformed to previously assumed proportions despite having adjusted those estimates with harvest rates. We suggest there is potential to effectively model populations of game birds using a combination of band recovery and hunter-harvested samples of species that can be readily captured, marked and have reasonable harvest rates.

The greater sage-grouse Centrocercus urophasianus (hereafter sage-grouse), a lekking game bird species of the sagebrush (Artemisia) ecosystem, was once widespread throughout the Intermountain West and harvested in every state in which it occurred (Reese and Connelly 2011). Long-term contraction of its distribution (56%) and more recent declining trends have led to unprecedented conservation efforts to maintain populations and their habitats (USFWS 2013). Research efforts to increase our understanding of the species' demography have been commensurate with increased conservation and management (Garton et al. 2011, Taylor et al. 2012. However, there are limited data on survival and harvest rates of hatch-year (HY) birds (Taylor et al. 2012), and unbiased estimates of tertiary sex ratio (Atamian and Sedinger 2010). Generally, tertiary sex ratio is reserved to describe that of the breeding cohort; we use it here to refer to the sex ratio of both juveniles transitioning to breeding age class, and adults. Recruitment of hatch-year birds (i.e. production) and tertiary sex ratio are foundational to the understanding of avian population dynamics. Sage-grouse are relatively long-lived for a gallinaceous bird, and population dynamics are often most sensitive to adult survival (Taylor et al. 2012). However, recruitment and, perhaps more specifically, survival of juveniles is frequently the principal driver of regulating population growth. Estimates of these demographic rates remain elusive for this species in some cases (Taylor et al. 2012). Population parameters are important for modeling population dynamics in response to natural (e.g. wildfire) or anthropogenic perturbations (e.g. energy development) and conservation (e.g. hunting).
Understanding age-and sex-specific susceptibility to harvest is vital to establishing hunting regulations and adaptively managing populations. Currently, there is a paucity of data on juvenile vulnerability to hunter harvest (Beck et al. 2006, Caudill et al. 2014, but the body of information on adult sage-grouse is becoming well founded (Connelly et al. 2000, Zablan et al. 2003, Broms et al. 2010). Adult harvest rates vary between 5 and 10% and appear to be similar between males and females (Zablan et al. 2003. Survival of juveniles to first breeding (i.e. recruitment) is poorly understood. However, the few studies that have examined juvenile survival from autumn (~Sep-Nov) to first breeding generally indicate that overwinter survival is similar to that of adults (Beck et al. 2006, Battazzo 2007, Caudill et al. 2014, Apa et al. 2017. Given the short-term (3-year) duration of these studies, we are lacking information on annual variation in recruitment (but see Blomberg et al. 2014).
Hunter-harvested wildlife affords an opportunity to collect data on vital rates that may otherwise be logistically or cost prohibitive (Broms et al. 2010, Dusek et al. 2014). We used six years of band-recovery and composition of age and sex information from sage-grouse wings collected from hunters in Moffat County, Colorado (1978-1983. We used a novel application of the Lincoln estimator (Lincoln 1930, Alisauskas et al. 2014 to generate an unbiased estimate of tertiary sex-ratio and production using band recoveries and wings collected from hunter-harvested birds. Our primary objective was to examine age-specific patterns of survival and harvest rates of male and female sage-grouse. Specifically, we sought to provide estimates 1) of juvenile survival to adulthood, as this parameter is not well described in the literature, 2) of harvest rates for juvenile and adult sage-grouse, and 3) of the tertiary sex-ratio.

Study area
We studied sage-grouse in Moffat County which is in the northwest corner of Colorado adjacent to Utah on the west and Wyoming on the north. Our efforts were concentrated at Blue Mountain (adjacent to Utah on the west), Cold Spring Mountain (adjacent to Utah on the west and Wyoming on the north), northcentral Moffat County (adjacent to Wyoming on the north), and eastern Moffat County including a small portion of Routt County on the east (both adjacent to Wyoming on the north). However, for brevity we refer to the study area as Moffat County. These areas varied in elevation with the highest at Blue Mountain and Cold Spring Mountain = ~ 2785 m, lowest = ~ 1525 to 2000 m in the northcentral portion of the county . All areas had extensive areas dominated by big sagebrush Artemisia spp. with understories of native forbs and grasses with some aspen Populus tremuloides, mountain shrubs (Amelanchier spp., Purshia spp., Rosa spp.), and small trees (Juniperus spp.). Some agricultural fields (barley, oats, wheat, native hay) occurred at lower elevations.

Capture
Sage-grouse were located and captured using spotlights and long-handled nets (Giesen et al. 1982) near lek sites during the breeding season and in known summer ranges during July and August using, walk-in drive traps (Giesen et al. 1982) and bumper-mounted (Latcher and Latcher 1964) and stationary cannon nets (Braun 1976). Each bird was classified to an age and sex category based on body size and diagnostic plumage characteristics (Braun and Schroeder 2015). Once age and sex were ascertained, each bird was marked with a uniquely numbered aluminum band provided by the Colorado Division of Wildlife (CDOW) and released where captured. The CDOW did not have an IACUC protocol prior to about 1990, however, guidelines of Gaunt et al. (1987) were followed for appropriate care while capturing and handling sage-grouse.

Data recovery
Wings were obtained from hunter-harvested sage-grouse at three hunter check stations and 16 wing barrels (Hoffman and Braun 1975) located at major access points to areas where sage-grouse were commonly hunted. They were labeled as to site of collection and stored frozen until thawed and classified to age and sex (Braun and Schroeder 2015). Bands were obtained at check stations, in wing barrels, and through return of hunter questionnaires as well as personal contacts such as walk-ins and reported through the mail or via the phone. All hunters were required to obtain a permit prior to hunting during 1978-1980 which allowed a 100% questionnaire survey about hunting activities including: area of hunting, dates of hunting, number of grouse harvested/lost, a question about recovery of banded birds, etc. (Braun 1981). Data sheets were prepared for each band reported and were later transferred to digital files for storage and manipulation. Hunting seasons generally varied between 9 (1978), 16 (1979, 1981-1983), and 25 days (1980), beginning on the second Saturday of September each year. Hunters were allowed three birds per day and six in possession .

Statistical approach
We used Brownie parameterization of age-structured recovery models to estimate annual recovery rates (f ) of male and female sage-grouse (Brownie et al. 1985). The Brownie parameterization provides estimates that can separate birds reported from hunter harvest (f ) from other sources of mortality (m = 1 -S -f ) (Cooch and White 2014). We included age-class at banding (age; adult or juvenile) and timedependence (t) to investigate the effects of age and annual variation on survival (S). We modeled recovery probability f as a function of time (t), age and sex. Because of the relatively short duration of our study, and the absence of capture and release of juvenile birds in the final year, we limited timedependence in our model set. Specifically, we constrained S to estimate temporal variation in f and vice versa. This approach enabled us to elucidate some of the temporal variation in survival and harvest, and potential changes in sexratio over time. However, this approach limits inference of S and f, especially under the Brownie parameterization because f contains some information about mortality and constraints on either parameter and necessarily makes assumptions about the pattern of variation in the unconstrained parameter (Cooch and White 2014). Mark-recapture analyses were conducted in Program MARK (Cooch and White 2014).
Model selection was used to identify the most parsimonious model from which to estimate age and sex specific parameters of survival and harvest. After examining the fit of the global model, models with fewer parameters were fit to the data. The number of potential models was large and we used a hierarchical procedure to guide model fitting (Lebreton et al. 1992). In both S and f, we started with additive models with main effects only (e.g. sex + age + t).
Model selection was based on Akaike's information criterion adjusted for small sample sizes (AIC c , Burnham and Anderson 2002). Models where ∆AIC c  2 from the best fit model (∆AIC c = 0) were considered equally parsimonious. The ratio of AIC c weights between two models was used to quantify the relative degree that a pair of models was supported by the data (Burnham and Anderson 2002). All parameter estimates (including means) are presented with 95%CIs otherwise noted.
We applied the capture-recapture Lincoln (1930) estimator for closed populations to provide unbiased annual (t) estimates of abundance (N) for each age (i) and sex (j) of a harvested population (Alisauskas et al. 2009). The Lincoln-Peterson capture-recapture estimator of abundance for closed populations is where n 2 is the number of animals captured in the second of two sampling occasions, and p  2 is the ratio of recaptures to animals captured on the first occasion. The estimator of population size, N, from hunter harvest was where H t was an estimate of harvest for year t, and h t , was the estimated probability that a bird alive at the time of banding in year t was harvested in year t. In the context of our work, the number of wings reported was the equivalent of H t and f was the equivalent of h t , assuming a band reporting rate of 1.0. Thus, from N tij we then estimated annual sex-ratio and age-ratios (i.e. production, HY: AHY females) in the harvest as the ratios of male and female population size or juveniles per adult. We note that the assumption that band reporting rates = 1.0 could be violated but not as long as band-reporting rates were the same for both males and females or for juveniles and adults. We used the delta method (Powell 2007) to estimate confidence intervals for sex and age-ratios estimated from the Lincoln calculations. For count data from wings, we used Wilson's confidence intervals as an estimate of variance for each ratio as it provides excellent coverage of the ratio estimators (Hagen and Loughin 2008).

Results
In total, 3259 sage-grouse were captured and banded from 1978-1983. Sixteen birds were removed from the study because of recovery not associated with harvest. Hunters harvested and reported 368 sage-grouse bands (3243 available for harvest) to the Colorado Division of Wildlife (Table 1). Hunters provided 10,231 wings from which age and sex of individual birds could be classified (Table 2).
Model selection based on AIC c indicated the best fit model for estimating temporal variation in survival (S Ad sex + t, S Juv c + t, f age) was one with age specific recovery rates, but similar rates between sexes in each age class ( Table 3). The model also indicated sex specific survival of the adult class, but juvenile survival was best estimated as one probability for both sexes. The same model but with sex specific survival of juveniles was not well supported (∆AIC c = 5.56). The best fit model for estimating temporal variation in recovery rate (S Ad sex, S Juv c, f age + t) was one with age-specific recovery rates, but similar rates between sexes in each age class (Table 3). However, a competitive model (∆AIC c = 1.82) was one in which sex-specific recovery rates were estimated for each age class, but this model was 2.47 times less likely than the best model. Because of our temporally limited data and relatively small effect sizes, we defaulted to make inference from the simpler model rather than use model averaging.
Annual variation in survival and recovery rates was greatest among the juvenile age class (Table 4). Generally, recovery rates were twice ( f  = 0.112, 95% CI: 0.097-0.131) that of adult sage-grouse ( f  = 0.063, 95% CI: 0.053-0.073), and survival of juveniles was 13 to 40% less than that of adult  Age and sex  1978  1979  1980  1981  1982  1983  Total   Adult male  49  764  597  11  21  36  1478  Adult female  36  168  150  31  45  14  444  Juvenile male  217  136  154  81  72  0  660  Juvenile female  206  177  167  70  57  0  677  Total  508  1245  1068  193  195  50  3259   Table 2. Number of greater sage-grouse wings reported from hunter harvests in Moffat County, CO, 1978CO, -1983 Age and sex 1978 1979 1980 1981 1982 1983 Table 3. Recovery models used to examine the effects of age survival of male and female greater sage-grouse in Moffat County, CO, 1978CO, -1983. Model fit is described with deviance (Dev), the number of parameters (K), the difference in Akaike's information criterion corrected for small sample size from the best fit model (ΔAICc), and AIC c weights (wi). Model structure estimated annual survival (S) for adult (S Ad ) and juveniles (S Juv ), and probability of recovery (f). The letter c denotes that a parameter was estimated as constant within a particular age or sex of birds, and t indicates an additive time effect. Additive models (+) included only main effects. males and females, respectively. Because we constrained our temporal trend to additive models only, we are limited in our inference about temporal variation. These additive models indicated that survival was greatest at the beginning of the study and steadily declined until the penultimate year (Table 4) Our population estimates averaged 29 611 over the study period (Table 5). From those estimates we derived tertiary sex-ratio for adult and juveniles. Sex ratio was generally skewed toward females in both adults (range = 1.352-2.607 females per male, Fig. 1A-B.) and juveniles (range = 1.016-1.320 females per male, Fig.  1C-D). However, statistically the tertiary sex ratio of juveniles was not measurably different from 1:1 based on Lincoln estimates (Fig. 1C). Adult sex ratios estimated from counts of wings alone were generally 11% greater than those based on Lincoln estimates (Fig. 1B). Juvenile sex ratio estimates were 6% greater in wing count data than those based on Lincoln estimates (Fig. 1C). Because of the differences in age-specific band recovery rates, estimates of production based solely on wings were on average 1.2 juveniles per adult female greater than those adjusted for the band recovery rate (Fig. 1E-F).

Discussion
Our objectives were to provide insights to harvest rates of juvenile and adult sage-grouse and unbiased estimates of production and tertiary sex-ratios. Our study provided several key findings, 1) juvenile harvest rates were twice those of adults, 2) harvest vulnerability was similar between adult males and females, 3) survival of juveniles was highly variable, and 4) sex ratios of adults and juveniles largely conformed to previously assumed proportions despite having based those estimates on harvest rates .
Our analyses yielded novel insights about harvest mortality, indicating that juveniles were nearly two-times more likely to be harvested than adults, and adult males and females shared a similar vulnerability to harvest. The latter finding was contrary to anecdotal evidence from hunters which suggested that adult females exhibited substantially greater susceptibility to harvest than males . Beck et al. (2006) reported a harvest rate of 0.02 from 1 juvenile of 56 susceptible to hunter harvest in southeastern Idaho. Our findings were similar to harvest rates in Nevada and elsewhere in Colorado, with no differences in harvest rate between the sexes of adults . However, sub-adult (after-hatch year) birds did exhibit differences in harvest susceptibility between sexes, and compared to Table 4. Survival probabilities (S) and recovery rates (f) estimated from greater sage-grouse band recovery data using Brownie parameterization, Moffat County, CO, 1978-1983 adults (Zablan et al. 2003). Averaged across age and sex class, and year, our harvest susceptibility ( f  = 0.083 SE = 0.005) was nearly equivalent to that found in Nevada ( f  = 0.09 SE = 0.008). Sedinger et al. (2010) determined that harvest rates 0.10 were compensatory to natural mortality. If that estimate is correct then harvest was compensatory for all cohorts in our study except for juveniles in 1979 and 1980 (Table 4). Caudill et al. (2014) reported greater harvest rates (0.09 and 0.27 in 2008 and 2009, respectively) from radiomarked juvenile sage-grouse at Parker Mountain, Utah. Although sample sizes were limited (50 birds available for harvest) these higher rates may be reflective of susceptibility to harvest in geographically isolated populations (Gibson et al. 2011). Alternatively, these higher rates could be a result of increased vulnerability of radio-marked birds as has been observed in some waterfowl studies (Parker 1991). Figure 1. Annual estimates (95% CIs) of greater sage-grouse tertiary sex ratio for adults estimated from Lincoln estimator (A) and from counts of wings (B), juveniles estimated from Lincoln estimator (C) and from counts of wings (D), and production (juveniles to female ratio) estimated from Lincoln estimator (E) and from counts of wings (F), Moffat County, CO, 1978CO, -1983 Variation in post-fledging survival has been linked to body condition and climatic variation (Blomberg et al. 2014). Sage-grouse juveniles appear to have improved survival in years with cooler and wetter growing season conditions (Blomberg et al. 2014). Given the limitations with available climatic data, we were unable to link survival to such covariates. However, we found survival of juveniles to be highly variable and lower than studies of radio-marked birds (Battazzo 2007, Beck et al. 2006, Blomberg et al. 2014. Direct comparisons of our juvenile survival estimates are limited because telemetry studies followed juveniles to the first breeding season (Mar-Apr) and our estimates follow them to the subsequent hunting season. However, we transformed our annual estimates to monthly survival estimates ( S  month = S 1/12 ) and then projected those out seven months to provide comparable estimates. We estimated juvenile survival to first breeding (Oct-Mar) and it varied from 0.131 to 0.708. These estimates were considerably lower than those from radiomarked birds in Idaho ( S  = 0.80, Beck et al. 2006) and Montana ( S  = 0.91, Battazzo 2007), but similar between 2005 and 2008 for the same region of Colorado as our study ( S  range of means 0.477 to 0.657, Apa et al. 2017). One limitation to these comparisons could be unaccounted for seasonal variation in survival with reductions during the breeding season. Thus, our study will tend to underestimate mean monthly survival compared to telemetry studies that monitored birds until the nesting season. While our data did not support sex specific survival rates of juveniles, recent work by Apa et al. (2017) demonstrated consistent measurable differences. Juvenile female survival was between 7 and 13 percentage points greater than for males, which was similar to our estimates (male S  = 0.258, female S  = 0.388), albeit from a model not well supported in the data set. These data points lend support to previous hypotheses regarding differential physiological demands between the sexes that may result in the observed lower survival of juvenile male sage-grouse (Swenson 1986).
Our estimates of adult annual survival were comparable to those reported previously , Taylor et al. 2012, Dahlgren et al. 2016 with the exception of 1981 when adult male and female survival was  0.14 and 0.30, respectively. Reduced annual survival has been linked to severe winter weather conditions (Anthony and Willis 2009). As measured by Anthony and Willis (2009), the winter of 1980 would have been considered severe in our study. Per the weather station at Craig, Colorado there were 31 days with temperatures  -18ºC and a snowpack of 917 cm. Although these were commensurate with lower survival, the winter of 1979 was more severe but adult survival was seemingly unaffected. Juvenile annual survival was only directly comparable to that of harvested birds in Nevada at S  = 0.40 ) and our estimates were considerably less ( S  0.098) in 3 of 5 years, and overall was 23% lower. Extrapolating estimates to first breeding from Idaho ( S  = 0.68) and Montana ( S  = 0.80) to annual survival suggests greater survival of juveniles in those populations (Beck et al. 2006, Battazzo 2007. However, such extrapolation assumes constant survival during the breeding season which typically is a period of high mortality and our estimates may be biased low. Nonetheless, if adjusted for breeding season mortality we suggest juvenile survival estimates from other systems would remain considerably greater than our lowest annual survival estimates that were  10%. Our estimates of tertiary sex-ratio of juveniles and adults are the first to be correctly adjusted for harvest susceptibility (Atamian and Sedinger 2010). Because our modeling indicated little variation in harvest rates within age classes it did not result in measurable differences in our estimates of sex-ratio from counts of wings alone. The juvenile sex ratio was essentially 1:1, and adults exhibited strong bias towards females in the population 2.04 (95%CI: 1.45-2.64, Braun et al. 2015). This pattern is reasonably supported by the differential in survival between the sexes , Taylor et al. 2012). However, our adjusted estimates of production were considerably less (1.2 juveniles: female) than those generated from counts alone (Fig. 1D-E). While the trends of each data set tracked one another, the adjustment for harvest susceptibility may provide for a more accurate estimate of production. Additional work, perhaps a simulation study could assist in more clearly depicting the bias associated with using only count data.
We are encouraged by our approach and application of mark-recovery data to the Lincoln estimator in combination with hunter harvest information as a potential tool for reconstructing population sizes and age ratios of harvested populations (Alisauskas et al. 2014). Although contemporary harvest of sage-grouse may be more limited than in our study, there is potential for future application of these methods to hunted populations. Additionally, there appears to be some promise for other gallinaceous birds with greater harvests such as sharp-tailed grouse Tympanuchus phasianellus, gray partridge Perdix perdix, chukar Alectoris alectoris or quail Callipepla sp. that can be readily captured, marked and have reasonable harvest rates.