Assessing red deer hunting management in the Iberian Peninsula: the importance of longitudinal studies

Study area

The two study populations are located in Castilla-La Mancha, South Central Spain. The Mediterranean woodlands and scrublands in the area constitute a continuous distribution of independently managed private or public hunting estates, including natural reserves. One of the estates is a private property devoted to game hunting (LM; 38°55′N, 0°36′E; 600–850 m a.s.l.) and the other is a public hunting preserve (LQ; 39°45′N, 4°15′E; 600–1,100 m a.s.l.). Since 1989, the LM deer population has received supplementary forage (protein-rich pellets) all the year-round, and in summer, alfalfa (Medicago sativa) has been supplied on the ground. The deer population in the LQ estate did not receive supplementary forage. They were selected as being representative of the more frequent typologies of big-game hunting estates in Mediterranean Spain. Table 1 summarises the main characteristics of the study areas.

The Mediterranean habitat is characterised by savannah-like landscape called ‘dehesas’, with evergreen oak, Quercus ilex, woodlands, and scrublands (dominated by Cystus spp., Erica spp., Pistacia spp., Phyllirea spp. and Rosmarinus spp.), scattered pastures and small areas of crops.

Estimation of red deer density and population growth rates

Between 1989 and 2015, the annual harvest rates and deer population sizes were estimated at both hunting estates. The deer density was estimated in the LQ estate by spotlight censuses from a 4 × 4 vehicle that were performed during August–September. Each transect was an average of 20.3 km ± 2.34 (S.E.) in length, and the speed of the vehicle was <15 km/h. The distance between the observer and the deer was measured with a Leica LRF 1200 Scan telemeter (Solms, Germany) (range 15–1,100 m; precision ±1 m/±0.1%). The annual abundance of the deer was estimated using the distance sampling method (Buckland et al., 2004) with Distance 6.0 software (Thomas et al., 2010). Half-normal, uniform and hazard rate distributions for the detection function were fitted against the data using cosine, Hermite polynomial and simple polynomial adjustment terms, which were fitted sequentially. The selection of the best model and adjustment term were based on Akaike’s Information Criterion (AIC). More details on the analytical framework can be found in Acevedo et al. (2008).

The size of the deer population at the LM estate (Table 1) was estimated annually by using repeated direct counts at feeding sites at the end of July and during the rutting season (total = 8 days/year), that is, September and October (see Rodriguez-Hidalgo et al., 2010). In the summer, most deer congregate around feeding sites and are counted at dusk. The estimates obtained from direct counts are equivalent to those obtained using the same approach used in LQ, that is spotlight counts and distance sampling. We evaluated the consistence between distance sampling and direct counts at feeding sites in three different years in LM (2003, 2006 and 2020). The results showed a high equivalence in the estimates obtained; 0.32 (CV 16.8%)–0.36 deer/km² (distance sampling—direct counts), 0.41 (CV 24.1%)–0.42 deer/km² and 0.41 (CV 20.4%)–0.48 deer/km² for 2003, 2006 and 2020, respectively.

We obtained the sex- and age-specific structure of the deer populations during 1989–2015 period using the number of animals recorded for each sex and age during the transects (LQ) and direct counts (LM). We registered the number of calves (<1 year old), yearlings (between 1 and 2 years) and adults (>3 years) (Landete-Castillejos et al., 2004). We combined the numbers for both sexes and defined the sum of yearlings plus adults as the total population size in year t (N_≥1,t). We defined the number of yearlings as the number of recruits at the end of the year. We calculated the adult sex/ratio (females/males), the birth rate (calves/females), the survival rate (yearlings/calves from the previous year)*100, the recruitment rate in year t: (N₁,_t+1/N_≥1,_t) (Ueno, Kaji & Saitoh, 2010) and the exponential population growth in year t (ln(N_t+1/N_t)) (Hebblewhite, Pletscher & Paquet, 2002); where N_t represents the population of a given year and N_t+1 is the population of the next year to characterise the population dynamics. Henceforth, we refer to the exponential population growth rate as the population annual growth rate. We specify that our study was based on the logistic population growth model (which is not linear). However, as abovementioned, ln(N_t+1/N_t) (i.e. Exponential growth rate or intrinsic rate of increase) is a linear function of N (Hebblewhite, Pletscher & Paquet, 2002). The intrinsic rate of increase, as defined in the exponential equation, is not a constant number at all but rather is itself a function of the density of the population, and it can be assumed a linear approximation.

Data analysis: theoretical maximum sustainable yields harvest (MSYt)

Ecological carrying capacity can be approached by the number of individuals that the environment ‘can support’ in a given area (Chapman & Byron, 2018). Therefore, in this study the ecological carrying capacity, represented by the K of the logistic equation, is useful to study the balance between forage resources, negative impacts and population productivity. The theoretical Maximum Sustainable Yields harvest (MSY) is the maximum number of deer that can be harvested from the population on a continual basis (McCullough, 1992), in this case, based on the logistic population growth model. Theoretically, and to establish a reference practical point for internal comparisons in our study populations, we considered that the population level that produces the MSY is approximately 50% of K (Jensen, 1996), which has been referred to as the Economic Carrying Capacity (Caughley, 1979), or I-Carrying Capacity (i.e. K/2 or population abundance at MSY, MacNab, 1985; McCullough, 1992; Buckland et al., 1996). When the growth rate (as density-dependence relationship) is modelled by logistic equation, r₀ is the maximal rate of population growth from low density and K is indicated where the regression cuts the X axis (i.e. exponential population growth is 0). Under the logistic model, a negative linear relationship between the growth rate and density (including confidence intervals) is expected, which is used as a simple test of statistical density dependence (Wolda & Dennis, 1993; Reynolds & Freckleton, 2005; Lebreton & Gimenez, 2013).

A linear model was used to evaluate the effects of density on the population dynamics and calculate the related parameters. Beyond analysing the complete set of potential factors determining the growth rate (the authors, in preparation) our model included as covariables the age and sex structure of the harvest. Age and sex structure are critical factors regarding the demographic parameters affected by population density or supplementary feeding, and they are quite relevant in determining culling strategies (Trenkel, 2001). For this purpose, age class was simplified to individuals <1 year/(calves) and >1 year (broadly, adults). These extractions, separately by sex and ages, are expressed in terms of proportion of individual extracted respect to available numbers. In order to determine the potential impacts of factors depending on the management regime (supplement vs unsupplemented), the population density and their respective interactions with other explanatory factors were included in the model.

The fitted linear regression on the exponential growth rate on the density helped us to calculate the theoretical K parameter, that is the size of the population for which the population growth rate is 0. The residual population accommodating MSY was also calculated as a proxy to a sustainable extraction of the populations (because it meets the requirements of maximising hunting while reducing the impacts of overabundance). According to the logistic model, at low densities, the population growth rate is the maximum, but the product (rN = net recruitment (number of net individuals that is added to the population per unit of time), r is referred to as the intrinsic rate of increase) is low. As density increases, the growth rate decreases and N increases, the net recruitment increases to a maximum K/2. This point, the MSY, is the theoretical maximum number of individuals that can be extracted from the population on the long-term. Since this theoretical parameter was considered a proxy for purposes of comparison between the populations in our illustrative scenario, following Stephens & Sutherland (1999) we were conservative in calculating and establishing a safety margin to avoid overexploitation (data for deer in the literature approaches 60% of K) (Clutton-Brock et al., 2002; Kaji et al., 2010). The calculation of such parameters allows us to compare real population extraction rates and those based on the assumption of the MSY value over years in two sets of contrasting data under different management regimes.

Characterization of red deer population dynamics in the contrasted study populations

The densities in LM population increased since 1989 and ranged around 40 ind/km² with considerable peaks (reaching above 50 ind/km² in 1995, Fig. 1). LQ population ranged consistently between 30 and 40 ind/km² overall, until 2009 with a consistent declining trend at the end of the study period, although with considerable peaks (e.g. high 2009), reaching the lowest level in 2014 (slightly above 20 ind/km²). With the exception of the first 4 years, the density was consistently higher in LM than in LQ population (Table 1; Fig. 1). Overall, densities were lower in LQ population, and population changes in both populations did not show any regular and/or cyclical pattern but marked ‘saw-tooth-like’ curves.

The exponential growth rate achieved values close to 0 in both populations (LQ = 0 y LM = 0.05); however, annual variations were significant, with positive and negative values (Table 1). LM population presented higher birth, recruitment survival rates, and proportions of calves and yearlings; while the proportions of females and adults were higher in LQ population (Table 1). In addition, the proportion of the population that was annually harvested and the harvest density were higher in LM population, with relevant annual variations in both study areas (Table 2; Fig. 1). Both populations presented sex ratios balanced slightly favourably to females.

Regarding harvesting, the extraction rates followed irregular patterns (see Fig. 1), alternating between periods of high and low hunting. The proportion of harvested tended to increase in LQ population during the last period of study resulting in a reduction of density (significant negative relationship between the annual proportion of harvest and density along the study period in LQ population, Beta = −0.35, F = 12.9, p < 0.01). There was an association between harvest density and population density in LM population, especially during 1989–1999 (p < 0.05, T = 2.21); however, low extraction rates during 2000-2007 led to an irregular, but increasing, density pattern.

Harvesting parameters in relation to the MSYt value

The linear model evaluating the relationships between the exponential annual population growth rates and densities calculated for both populations for the period 1989–2015 evidenced a decrease in population growth at high densities in both situations, whether artificial feeding was provided or not (F = 30, p < 0.001) (Table 3), while age and sex class specific extraction rates, as well as the interactions between area and the other parameters (including density) were not significant.

The linear regressions fitted to relate the exponential growth rate and population density in each population (Fig. 2) evidenced that both populations presented significant linear adjustments (LQ: F = 13.15, p = 0.0017; LM: F = 33.62, p = 0.001; R² = 0.57 and 0.39, respectively) with similar slopes 0.72 Lm (CI [0.47–0.96]) and 0.65 LQ (CI [0.26–0.98]). In LM population, the test of statistical density dependence established K at 39.45 deer/km², while for LQ population it was 32.64 deer/km². These values are similar to the average densities occurring during the study period: 37.9 and 32.9 deer/km² for LM and LQ, respectively. Therefore, the range of I-carrying capacity assuming 50 of K value was between 19.73 deer/km² and 16.32 deer/km² for LM and LQ populations, respectively.

Figure 3 plots the harvested density and the proportion harvested for both populations, indicating as a reference their values at MSY (for K_50% according to our previous estimations ofItheoretical values). Table 2 shows the total number and Table 4 the proportion harvested in each population. The total number of animals hunted was slightly higher in relation to what was expected under an MSY assumption for both populations. In the case of the extraction rates, the real observed numbers were below the percentage of sustainable extraction for the whole study period, except at 2 years in the case of LM and at 4 years in LQ (Table 2).

Overall, an average of an extra 36 deer per year were hunted (0.52 deer annually hunted/km²) from LQ population compared with a sustainable MSY situation (345 individuals, 5.03 hunted/km²). This represents 10.4% more hunting bag compared to MSY. On the other hand, an extra 21.7 deer were hunted per year (2.5 hunted/km²) from LM population compared with a sustained MSY situation (58 individual, 6.68 hunted/km²), which represents 37.4% more hunting bag compared to (artificial) MSY. During the total study period (27 years) 976 and 586 fewer deer would have been hunted under an MSY situation for LQ and LM respectively (see Table 2). In terms of average extraction rates, a 13.3% and 10.21% lower proportion of the population would have been extracted compared to the MSY situation in LQ and LM populations, respectively (Table 4). In other words, the extraction rate would be 1.75 and 1.43 times higher under an MSY situation, for LQ and LM populations, respectively.

Density dependence and management

The average deer densities of the study populations might be considered high in the context of the values described in the international scientific literature, where large range of values exist for this plastic species even locally (Apollonio, Andersen & Putman, 2010), for example 3 deer/km² in Sweden (Jarnemo, 2008), 3.7 deer/km² in Norway (Mysterud et al., 2007), 8.1 deer/km² in Scotland (Pérez-Barbería, Hooper & Gordon, 2013) or 9.6 deer/km² in the Czech Republic (Prokesová, Barancekova & HoMolka, 2006). Although, Buckland et al. (1996) showed densities between 20 and 34 deer km² in Scottish red deer populations. It should be remembered that Central European literature generally considers that much lower densities than those described here (e.g. densities above 4 deer/km²) are not compatible with other activities of agricultural and forestry production and may have adverse effects on the conservation of the natural environment (Gill, 1992). Similar high ranges of densities have been found in other studies in Mediterranean ecosystems and areas with a temperate climate, for example up to 40 deer/km² in Portugal (Silva et al., 2014), >30 deer/km² in other areas of Spain (Perea, Girardello & San Miguel, 2014) or 26 deer/km² in Italy (Casula et al., 2009), 20–34 deer/km² in Scotland (Buckland et al., 1996). These values indicate that the carrying capacity of red deer in Mediterranean and temperate oceanic environments may present relatively high values, compared to northern or continental latitudes. Climatic conditions and the quality of the food available (e.g. Real Bioclimatic Index, Martínez-Jaúregui et al., 2009) are adverse in the summer in Southern Europe due to drought (with a second minimum in winter), showing significant interannual variability. However, the vegetation productivity curve in our study latitudes varies little over the year due to the presence of evergreen trees and shrub species (Gazol et al., 2018), which contrasts with the Northern and Central latitudes of Europe, where hard winter conditions limit winter survival, determining lower average carrying capacities (Martínez-Jaúregui et al., 2009). However, NDVI varies considerably among plant species, droughts, and the location of the studied forests (Gazol et al., 2018). These authors also indicated than NVDI data are less sensitive metrics of forest resilience to drought than other metrics as TRWi (absolutely dated ring-width indices).

The relative rate of increase per individual decreased linearly with population density in both study populations. However, equilibrium population density, K, is increased by supplemental feeding. This happened even when natural conditions or habitat quality apparently was lower in the supplemented population (Milner et al., 2014). López-Montoya, Moro & Azorit (2017) detected that density-dependence may have a stronger effect on growth than hunting or climatic conditions in populations which are closer to carrying capacity (these aspects are being addressed for these two populations in a separate article). A high population density experienced early in life may cause a reduction in subsequent development, thus affecting antlers size over time (Torres-Porras, Carranza & Perez-Gonzalez, 2009; Pérez et al., 2011), which may be detrimental in terms of the objectives of hunting. Miranda et al. (2015) showed that supplemental feeding provided to ungulates in LM led to increased browsing on plant species whose nutritional composition complemented that of the supplement provided. In our study populations, Acevedo et al. (2008) reported high levels of browsing in both study sites (0.95 and 0.86 browsing indexes for highly palatable species, and 0.85 and 0.83 browsing indexes for species with reduced palatability, in LQ and LM populations, respectively), which are high according to previous studies (Morellet et al., 2001). In addition, aggregation at feeders may increase disease transmission (Barasona et al., 2017). Vicente et al. (2007) showed that supplemental feeding increases the body condition at the cost of increased host contact rates, which implied an ecological trade-off between acquiring resources and the concurrent exposure to mycobacterias (the causative agent of animal tuberculosis, in this case). Understanding these long-term complex links represents an important contribution to the planning and prioritising of population management actions.

We also verify that the neither the population growth rate nor its negative relationship with population density in both study populations was affected or mediated by the rates of harvesting different age and sex classes. This maybe result of the hunting strategies implemented in the study sites, where not specific plans intended population reduction by targeting specific sex by age classes.

Harvest rates and proactive sustainable management

The high average densities (close to theoretical K values) only yielded 10.4% more hunting bag compared to the MSY in LQ population (this population can be considered as a reference for the natural carrying capacity in our study area). However, the disproportional increase in deer density to achieve such a limited hunting bag increase may in turn lead to negative impacts (see above). As for LM population, the hunting strategy yielded 30.4% more hunting bag compared to MSY due to very high densities (which peaked at up to 50 deer/km²) reached in a sub-optimal habitat for red deer. Therefore, this population cannot be considered a reference of natural circumstances. Theoretically, in LQ population, a density of 16.3 deer/km² may allow extraction rates of up to 30.3% (MSY); however, the actual extraction rate (17.5%) is well below the MSY threshold. Previous studies (such as DeCalesta (2017) or Mihalik (2016)) showed similar situations in wild ungulate populations where the extraction rate is below the sustainable rate under high density situations, and where an increased effort in the harvest rate reduced deer density, thus reducing impacts on ecosystems, while improving deer health and hunter satisfaction (Ueno, Kaji & Saitoh, 2010), although this depends on hunter’s preferences (e.g. trophy, shoots, meat; see Heberlein, 2002).

Our harvest rates compare with the range described by previous authors (Kiffner et al., 2008; Saint-Andrieux et al., 2009; Apollonio, Andersen & Putman, 2010; Martínez & Martín, 2017). They respond to local situations in space, time and management (Milner et al., 2006); however, we evidenced that the significant changes in harvest rates were a consequence of reactive decision taking. Buckland et al. (1996) suggest that the female (hind) cull should be restricted to non-lactating adult (yeld) and immature hinds, reducing hind numbers so that there are just two hinds of age >1 for every stag of age >6 in the population, provided that the population is maintained well below the carrying capacity (roughly estimated at 1,400 deer for LQ; Millas-Prendergast & Martín-Fernández, 2015), which approaches the situation we described for LQ population (K/2 population was estimated at about 1,268 deer). Proactive adaptive hunting management can maintain population levels within an optimum range, maximising the balance between forage resources and population productivity (Pérez-González et al., 2010; Rodriguez-Hidalgo et al., 2010; Herruzo & Martínez-Jauregui, 2013; Herruzo et al., 2016); that is when the population is at the optimal level, growth is high, and deer will be large (McCullough, 1979). For example Santín-Durán et al. (2004) showed a positive relationship between the intensity of parasite infections and population density of red deer in LQ, and Gómez et al. (2001), who showed the impact of red deer on vegetation, also recommended limiting the deer population in LQ. In this population, previous studies found that density was the most important factor explaining the variation in conception dates, with greater densities causing later conception dates (Peláez et al., 2017), as well as having a negative effect on antler quality (Peláez et al., 2018). In the same way, Santos et al. (2018) found a negative density-dependent response of deer body reserves when considering the combined effects of population density and dietary quality, especially in populations without a supplementary feeding regime (as in LQ). Perea, Girardello & San Miguel (2014) showed that some plant species were highly preferred and were not to be found in the areas where the deer were present in great densities (>30 ind/km²) for a long period of time (for example LM).

Priorities for managing hunted red deer populations in Mediterranean areas

Density-mediated factors determined the population growth in these two situations of high red deer densities in a Mediterranean habitat, even when artificial feeding was provided. Our research has shown that, in these particular situations, an increase in the extraction rates (as is already happening in LQ) is needed.

Management implications

The extraction rates reported here can be considered as indicative (rather than reference) values for the two prevalent hunting management regimes in the south-central Iberian Peninsula. Hunting plans must be based adaptively on population monitoring, thereby proactively preventing resource damage. For these purposes, a simple initial approach can be to fix sustainable hunting rates to mitigate the impacts of deer browsing. We must promote awareness of the problem among policy makers and managers, but also develop and apply available population monitoring tools, to be applied in proactive management plans, with active surveillance and the early detection of the consequences of overabundance. Adaptive management involves assessing the effectiveness of management actions in relation to the plan’s goals and objectives (Lyons et al., 2008).