Data

Harvest statistics in Denmark have been compiled from hunter submissions since 1941 (Strandgaard and Asferg 1980). Generally, the Danish hunting season comprises the autumn and early winter (Table 1). Reports are legally mandatory, but response rates amongst approximately 170 000 hunters have varied between 56 and 97% with a downward trend over time (1971–2006; Asferg 2008). Due to incomplete reporting, a correction factor is applied to derive the published national annual harvest (Supplementary material Appendix 1). The total national annual harvest for each individual species, used in this paper, is the sum of all the corrected annual hunting harvest totals from the 98 Danish municipalities.

Until 2004, no discrimination was made between common and jack snipe Lymnocryptes minimus, and therefore the ratio of common to jack snipe in the annual harvest of common snipe was estimated using the proportion of common snipe in the Danish wing survey record. The wing surveys are based on material submitted by hunters that voluntarily submit wings from harvested waterbirds identified to species, sex and age class by experts (Christensen and Fox 2014).

Table 1.

Summary of the changes in Danish hunting legislation and general description of hunting for six bird and mammal species. Species codes: GP = grey partridge, CO = coot, CS = common snipe, WP = wood pigeon, RD = roe deer and BH = brown hare.

Since 1976 point counts undertaken during 1 May to 15 June have been used to generate annual population indices of bird and mammal species (i.e. the Danish CBC, Larsen et al. 2011, Supplementary material Appendix 1). At the start of the scheme, only ca 75 routes were covered throughout Denmark, but numbers increased rapidly and have exceeded 300 since 1987. National point count indices were computed using TRIM software designed for analysis of census data (Pannekoek and van Strien 2001). This uses log-linear Poisson regression models to calculate annual indices of population size changes, using a generalised estimating equations (GEE) approach. TRIM offers a method to analyse count data that are zero-inflated and over-dispersed, as well as taking into account serial correlation between counts in consecutive years, coping with data containing missing observations, and the uneven distribution of birds. TRIM also calculates standard errors of indices (SE). A detailed description of the statistical background of TRIM and the formulas for calculation of SE is found in Pannekoek and van Strien (2001). TRIM fits models to generate the point count indices as estimates with SEs for years throughout the time series, with the exception of an anchor year where the index is arbitrarily set to 100 to generate an index of change in relation to this point in the time series. We chose to set the first year (e.g. for snipe 1987) as our anchor baseline year, but for this arbitrary value, we cannot generate a SE. However, because the comparison between the harvest statistics and the count indices require SEs for the point count estimates, we only commenced the comparison of the time series in 1988 for this species.

Comparison of harvest statistics and point count indices

Common snipe, wood pigeon, coot, roe deer, grey partridge and brown hare were selected for analysis as they are not subject to a substantial release of captive-reared animals.

For the purposes of these analyses, we considered that if the maximum SE of a point count index exceeded 30% of the total index range, it was considered too variable for inclusion in the analysis. For species where the full data set for birds (1976–2011) and mammals (1984–2011) did not meet the 30% criterion, the earliest years were excluded one by one until the 30% criterion could be met, but we did not truncate beyond 1987 to ensure a minimum series of 25 years. By choosing a 30% criterion we ensured that less abundant species such as coot and common snipe were not represented in the analysis in years with little data available; for further background information on the point count indices and standard errors see Table 2.

The annual harvest in Denmark for each analysed species was compared against the national point count indices on an arithmetic scale to determine the nature of the relationship between the two variables. We intended to incorporate an autocorrelation component (AR(1)-covariance matrix) in a repeated measures design, given that annual point count data represent a time series. However, a positive final hessian matrix could not be obtained from the AR(1)-models or the models did not converge, most likely because the models were too complex for the data available. A linear mixed model with the annual harvest as dependent variable and the point count index was applied instead. Time (year) was incorporated as a random effect using a variance component covariance structure, which provided stable models. To address indices with large standard errors (SE) we weighted indices by 1/SE², so that observations with large SE were less heavily weighted than those of small SE.

Data were log_e-transformed and a mixed model analysis with the annual harvest as dependent variable undertaken as described above to test if a regression slope (β) of unity existed for the point count index (β not significantly different from 1; t-test of slope) corresponding to a type I response. If β < 1, the regression line would level off (saturation effect) on a back-transformed arithmetic scale, corresponding to a type II response, and if β > 1, it would correspond to an accelerating curve, i.e. the concave part of a type III response.

Table 2.

Summary of the bird count indices, the associated average and maximum standard errors together with a relative measure of variation based on maximum standard errors. Data are provided for the six analysed bird and mammal species. Species codes: GP = grey partridge, CO = coot, CS = common snipe, WP = wood pigeon, RD = roe deer and BH = brown hare.

Unlike in North America, where there are regulated limits on how much game can be harvested by the individual, there are no harvest limits imposed in Denmark. However, changes in hunting legislation (e.g. hunting season length or the banning of hunting coot from motor boats — Table 1) were expected to influence total hunting harvest. For this reason, these changes and their interaction with the point count index were incorporated as covariates in additional analyses with log-transformed data. For example, different periods were treated as categorical variables by coding each period of hunting legislation with a corresponding integer. We used corrected Akaike information criterion AIC_c because of small sample sizes to derive the most parsimonious model and the most likely functional response.

The assumption of linearity was investigated by describing the relationship between the annual harvest and point count indices with Gompertz functions (Winsor 1932) to identify possible non-linear components under the expectations of type II and III functional responses (Beyer et al. 2013). The non-linear relationships were described as: bag = B₁ × exp(-B₂ × exp(-B₃ × index)). From this formula it is evident that if any of the parameters B₁ B₂ or B₃ are zero, a non-linear relationship does not exist between the annual harvest and the point count index. Hence, we used the 95% confidence intervals of the parameters as a proxy to determine non-linear relationships between the annual harvest and the point count index. We here present such results showing the nature of these relationships, i.e. a type II or type III functional response.

There was a high probability that data were serially autocorrelated as both harvest and point count indices represent time series. The two data sets were therefore compared to examine whether they comprised the same autocorrelation structure, as this would support the overall hypothesis that changes in harvest records can be used as a proxy for population changes. Using the approach of Forchhammer and Asferg (2000), we incorporated time (year) as a covariate in the auto-regressive analysis of the annual harvest and the point count index, respectively. Using year as a covariate assumes monotonically linear trends, so we incorporated the term year X year to allow for curvi-linear relationships. Where this changed the AR(x)-estimates very little, we retained the linear models.

The lag 2 partial autocorrelation function was used to estimate the first and second order partial correlation coefficients and these were plotted against each other to derive the fluctuations in the harvest and point count indices (Royama 1992). Initial analysis of Durbin—Watson statistics showed that this was not significantly different from 2 in AR(2)-models (p > 0.05), indicating no delayed density dependence (β₃, β₄ etc.). Pairwise t-tests of the coefficients of the annual harvest and the point count indices for each species were undertaken.

Regression analyses were carried out in SAS, using PROC MIXED and PROC NLIN, while PROC AUTOREG was applied for the autocorrelation analyses (SAS ver. 9.3).

Comparisons of bag statistics and indices

Four of the six species analysed showed significant positive arithmetic relationships between the annual harvest and point count indices (Table 3). The explanatory power of these regression models was reasonably good (r² > 0.50) in common snipe, roe deer and grey partridge. This could be the first indication that the harvest statistics for these species could be used as a measure of changes in populations (i.e. type I functional response).

Log_e-transformed data sets confirmed that the harvest statistics of common snipe showed the best fit of all species to a type I functional response to the point count indices in a direct comparison between the two variables (Table 3, Fig. 1a). The years 2003 and 2005 were excluded from the final analyses as outliers with large residuals (-0.51 and -0.61; Fig. 1a). Values of Cooks D (15.9 and 34.4) were close to or greater than 4/n and |DFFIT| (0.63 and 0.99) were greater than 2sqrt(p/n), where n = sample size and p = number of parameters. In 2003 and 2005 the annual harvest was considerably lower than expected from the point count index values. A post hoc analysis showed that the weather in September and October (the main hunting period) 2003 and 2005 was unusually dry at 70% and 76% of average precipitation respectively (Danish Meteorological Institute) likely reducing common snipe accessibility to hunters across Denmark.

The relationship between harvest statistics and point count indices in wood pigeon also showed a slope that converged to unity in a log-log comparison, but only when changes in hunting legislation were taken into account (Table 3). After 2006, hunting legislation was changed twice (Table 1), leaving too few data points (5) available for effective analysis. Hunting was banned in August during the period 1994–2006, which lowered the harvest size significantly compared to 1977–1993 (t = 3.37, DF = 27, p = 0.002). A corrected hunting harvest was obtained (i.e. as if hunting was still present in August) from a regression model that included the logarithmic annual harvest as dependent variable and the logarithmic point count index and period as covariates (model estimates: 8.7633 + (0.7964 × index) + (0.1945 × period)). The mean logarithmic difference (0.1945) between the mean size of the harvests during the two periods was added to the hunting harvests of 1994–2006. The corrected hunting harvest converged to a slope of unity (Fig. 1b), when compared directly to the point count indices on a logarithmic scale. Following this method, the relationship resembled a type I functional response with enhanced explanatory power (Table 3).

Table 3.

Relationships between annual hunting harvest and annual point count indices (index) indicated by the regression slopes of linear mixed models incorporating periods with differing hunting legislation (period, inserted as a dummy variable) for six bird and mammal species. For each logarithmic model functional response types were determined by the significance of t-tests (p < 0.05) of regression slopes of the indices: if slope = 0, no functional response type; if 0 < slope < 1, type II; if slope = 1, type I; if slope > 1, type III. For all models explanatory power (r ²), the difference in the corrected Akaike information criterion (δAIC_c) to the model with the smallest AIC_c-value and sample size are indicated. Species codes: CS = common snipe, WP = wood pigeon, CO = coot, RD = roe deer, CP = grey partridge and BH = brown hare. * The outliers 2003 and 2005 are not included.

Coot showed two very different relationships between harvest and point count indices (Fig. 1c) which was best explained as before and after the ban of hunting with motorboats in many Danish fjords, where the species occurs in large concentrations (Table 1, 3, Fig. 1c). Prior to the ban, there was evidence of a type I relationship between harvest statistics and indices, while after the ban there were indications of a saturation effect (type II), suggesting that the growth of the annual harvest decreased with increasing point count indices.

Roe deer also showed different patterns before and after 1996 (Table 3, Fig. 1d). In this case, the reduced growth of the hunting harvest with point count indices was associated with a month's extension of the hunting season (Table 1).

Finally, the two farmland species, grey partridge and brown hare, showed the most complex relationships between harvest statistics and count indices. In grey partridge significant non-linear relationships were found (Gompertz function: B₁ = 271 751 [178 490 - 365 011; 95% confidence interval], B₂ = 3.097 [2.482 - 3.713]; B₃ = 0.026 [0.014 - 0.037]; R² = 0.83; n = 34; type III functional response; Fig. 1e). For both grey partridge and brown hare the complex patterns were also reflected in a mixture of response types related to periods with varying hunting legislations (Table 3, Fig. 1f).

Autocorrelation structures

Auto-regressive models were established in order to test whether the harvest statistics and count data showed the same autocorrelation structures for each species using direct density dependence DD with coefficient β₁ and delayed density dependence DDD with coefficient β₂ in AR(2)-models. The AR(2)-coefficients were plotted against each other (Fig. 2) to determine the dynamics of the fluctuations in the data sets and to investigate the within-species correspondence between the auto-regressive parameters. Grey partridge, brown hare and common snipe fell within the lower semi-circle in Fig. 2, indicating these species tended to show cyclic fluctuations in the point count index (see Royama 1992 for further explanation). Distance between points was particularly large in grey partridge (Fig. 2) as there was a significant difference in both DD (t = 3.17, DF = 68, p = 0.002) and DDD (t = 2.66, DF = 68, p = 0.01). A significant difference in DD was also observed in roe deer (t = 3.51, DF = 52, p < 0.001) which indicated that the strength of density dependence in these two species differed between harvest statistics and point count indices. For all other species there were no significant differences in the strength of DD and DDD, suggesting that the autocorrelation structures in the data sets for harvest statistics and point count indices were comparable for these species.

Figure 1.

The relationship between the harvest size and point count indices (compared to a slope of unity) for four avian and two mammal species during periods of differing hunting legislation. For common snipe the outlier observations in 2003 and 2005 were not included.

How reliable are harvest statistics as a proxy for population changes?

Hunting harvest records have been the only data source for studies of population abundance in some species where direct surveys did not exist (Potts et al. 1984, Forchhammer and Asferg 2000, Kerlin et al. 2007). Where managers have no other sources of data with which to evaluate population trends, this situation may be unavoidable, but the present study shows that it is unwise to assume an arithmetic relationship between harvest statistics and relative population abundance.

Figure 2.

Plot of the relationship between first (1 + β₁ and second order delayed density dependence (β₂ coefficient of harvest records and point count indices. Outside the triangle the population will become extinct; inside the triangle the dynamics either reflect point stability or dampened fluctuations; below the semicircle the dynamics are cyclic. Under the semi-circle, the directions of arrows indicate the direction of change in parameters in order to go from short to long periodicity of cycles (Royama 1992, Bjørnstad et al. 1995). Species codes: GP = grey partridge, CO = coot, CS = common snipe, WP = wood pigeon, BH = brown hare and RD = roe deer.

It is useful to consider the harvest by hunters as an analogue to classic predator—prey interactions by incorporating functional response types which are established here via the log-log comparisons. We established two criteria to test whether it could be concluded with reasonable confidence that harvest statistics can be used as a proxy for population changes, namely 1) a slope of unity in a log-log-comparison of harvest statistics and point count indices and 2) comparable autocorrelation structures in the two data sets to establish that the nature of the fluctuations in total harvest and count indices were comparable. Whilst this comparison does not seem important to the findings presented here, we consider this fundamental to the analyses, because it supports the use of 1) above and because it is essential to take into account the autocorrelation in both data sets given that they represent time series with potentially differing autocorrelation structures that could detract from the analysis. Of the six species, only data for common snipe complied with the basic prerequisites 1) and 2) above. It was essential to incorporate changes in hunting legislation into the analyses to improve model fitting between harvest and count data. However, using this information to make general corrections to the harvest data was difficult except in the case of the wood pigeon, because changes in hunting legislation were associated with temporal and spatial changes in hunting effort.

Changes in hunting effort also affected the relationships between the data sets for coot and roe deer in which the overall functional responses resembled type II responses making the harvest statistics for these species unsuitable as a proxy for population size. In the case of roe deer, a type II functional response was predicted given that this species has increased markedly in abundance to the extent that hunters are no longer able to maintain the harvest rate relative to the large population increases of very recent years.

The classic predator—prey type II interaction is associated with predators becoming constrained by the time available to search for and handle prey items when these become very abundant ultimately leading to a decelerating intake rate (Stephens and Krebs 1986). In recreational hunting type II responses may arise if hunters aim for a certain harvest size per season even at high abundance (saturation). The annual harvest in Denmark is remarkably insensitive to changes in season lengths which may indicate that many hunters aim for a relatively fixed harvest size irrespective of the opportunities to kill more prey (Sunde and Asferg 2014). In contrast to North America where harvest limits are applied more directly to individual hunters the opportunity to affect annual harvest by manipulation of hunting season length in Denmark (and in Europe as a whole) is more limited (Williams and Johnson 1995). The precise behavioural mechanisms that generate hunter saturation may be numerous. For instance, once a hunter has shot a certain number of roe deer for their own consumption they may be less inclined to invest more time in hunting, killing and butchering yet more roe deer given restriction on commercial sales of game in Denmark.

We predicted that common snipe would show a type III functional response, because the species has demonstrated a reduction in annual point count indices over the period. However, there has been very little public discussion about the decline of common snipe compared to that of the grey partridge which showed a type III response. For the grey partridge it is well known that hunters have shown restraint because of the widespread awareness of the current population decline. The same was expected for brown hare, although this showed only weak indications of a type III response. Assuming that hunters try to adjust their hunting effort based on assessments of local and regional population density, the complex harvest pattern shown for brown hare may result from the difficulties in assessing changes in population densities effectively. For example, the North American study of white-tailed deer Odocoileus virginiana suggested that the discrepancy between abundance based on perceptions and that in reality would widen as population density declines (Van Deelen and Etter 2003). In addition, the nocturnal behaviour of brown hare makes it difficult to assess changes in its population abundance based on point counts carried out by daytime. Although classic type III predator—prey responses are typically associated with prey switching (Schenk and Bacher 2002), it is unknown to what extent hunters showed restraint in killing, for example, grey partridge or brown hare and to what extent that was associated with switching to other huntable species.

How robust are the comparisons of harvest statistics and point count indices?

Both harvest statistics and point counts are biased indices of true wildlife population size. The magnitude of bias remains unknown. However, for our purposes we assume bias is constant within both indices. We compared point count indices with harvest data, based on the assumption that the former reflected true population fluctuations, but is this assumption correct (Caughley 1980)? A robust count index must retain a constant relationship with overall population abundance over time. Failure to meet this presumption has weakened the arguments for mesopredator suppression in Australia because of serious methodological flaws in generating indices of organism abundance (Allen et al. 2011). In the present study, we attempted to avoid the pitfalls of differential seasonal and habitat abundance and binary measures (presence/absence) by maintaining systematic methods to generate abundance indices based on continuous measures.

Nevertheless, systematic changes in detectability over time could decouple the otherwise constant relationship between population indices and overall abundance. Habitat changes or any other behavioural change may affect detectability, so it is important to consider examples of such bias. The extraordinary difficulties of measuring changes in detectability of birds and mammals mean that few analyses exist in the literature. Nonetheless, those studies that do quantify bias serve as examples to demonstrate the risk of drawing erroneous conclusions, when trying to use changes in hunting harvest as a proxy for population abundance.

Danish farmland has undergone dramatic changes that have reduced open foraging habitats in spring, such as spring sown crops (Fox 2004). These changes may be particularly important for species, which exploit both farmland (open habitats) and woodland (closed habitats) such as roe deer and wood pigeon. We cannot therefore exclude the possibility that the open farmland habitats associated with higher detectability has become less attractive during spring. This would under-record roe deer and wood pigeon during recent point counts, because a greater proportion may remain in closed habitat (woodland), underestimating relative abundance now compared to earlier times. In the initial analysis roe deer and wood pigeon showed indications of a type II and type I functional responses, respectively (Table 3). If the recent point count indices underestimated population size, shifting graph points to the right in recent years in Fig. 1d would conform to a type II response for both roe deer and wood pigeon.

What if the hunting harvest suffers from the same bias as the point count index? For example, autumn farmland is probably more attractive to roe deer than it used to be because the predominance of winter cereals in modern farming ensures more prevalent new green growth than in former times when bare soil and stubbles predominated. In this way roe deer could have become relatively more detectable and/or accessible to hunters in recent years. Moving points from recent years downwards in Fig. 1d also conforms to a type II response.

While the introduction of systematic biases in indices and annual kill of roe deer did not change our conclusions about the functional response, the example of wood pigeon clearly demonstrated the importance of having reliable estimates of population abundance and hunting harvest data which are not inflated by systematic bias in species detectability over time.

Detection bias may also originate from differences in timing; point counts take place in the main breeding season (May/June), whilst harvest statistics are compiled mainly during autumn and early winter. Many bird species occur in high numbers and show clumped distributions so are more conspicuous and potentially accessible to hunters during the autumn and winter than the same species are to ornithologists during their count surveys in May/June. As long as the difference in detection rates amongst hunters and ornithologists does not change systematically over time this would not introduce any differential detection bias. Systematic changes in the annual mortality or reproductive output (e.g. density dependent responses) could create such biases and thereby influence the relationship between harvest statistics and count indices as could permanent changes of migration routes. The hunting harvest of some bird species are known to include individuals that occur in Denmark on migration from elsewhere during the hunting season (e.g. wood pigeon, coot and common snipe). For instance, the Swedish wood pigeon population contributes many birds to Denmark in autumn, but it has a slightly different breeding trend from that in Denmark (Lindström et al. 2011). However, if this were a major source of variation in our study, we would expect many more outliers in the regression analyses and lower explanatory power in the models (Table 3).

Concluding remarks

The relationships between harvest statistics and point count indices in Denmark are complex and varied considerably between species. We therefore conclude that considerable caution should be exercised before total harvest can be considered as an index of population abundance. This is consistent with the conclusions from other recent studies (Ranta et al. 2008, Willebrand et al. 2011), and our own a priori expectations, which predicted that only coot of the studied species would exhibit a type I functional response with which it only partly complied. We therefore recommend detailed investigations of the relationship between harvest records and population abundance, before harvest records are used as an index for population size. It is essential that the detectability of a species to counters or accessibility to hunters does not change systematically over time. Harvest statistics are affected by the behaviour of the hunters themselves, for example as a result of exercising voluntary hunting restraint on species of concern and as a result of saturation effects in the case of very abundant species. Other studies have confirmed that hunting effort is a stronger predictor of harvest statistics than population abundance (Imperio et al. 2010, Willebrand et al. 2011).

Despite these findings, there remain powerful arguments for the continued collection of harvest records and collation of hunting harvest statistics. Denmark remains one of the only countries in Europe that compile such unique data, which are fundamental to the effective adaptive management of sustainable hunting. To fully understand the relationships between harvest statistics and prey population abundance, we recommend undertaking innovative sociological studies to provide further insight into the mechanisms that affect the hunting effort and decision-making of individual hunters in response to changing hunting legislation, hunters' own assessment of game population size, their motivation, choice of prey and effort, all factors known to affect the harvest size (Schmidt et al. 2005, this study).