Tracking red deer population size using deterministic cohort analysis

Reliable and cost-effective monitoring tools to track population size over time are of key importance for wildlife management and conservation. Deterministic cohort analysis may be used to this aim, especially in hunted populations, but it requires that all mortality events are recorded and that individual age at death is known exactly. In this study, we investigated the reliability of cohort analysis as a relative index to track over-time variation in red deer (Cervus elaphus) abundance, in the absence of exact information about natural mortality and age. Visual tooth inspection was used to age 18,390 individuals found dead or hunted between 1982 and 2020 within the Trentino sector of the Stelvio National Park and the Val di Sole hunting district (Central Italian Alps). Temporal trend of reconstructed population size was checked using spring spotlight counts as a benchmark, through the Buishand range test and a linear model. Our results showed a significant and positive relationship between reconstructed population size and spring spotlight counts between 1982 and 2013, suggesting that cohort analysis could reliably track red deer population trend up to 7 years in the past. With a relative error of  +  1.1 (SD  =  1.5) years in the estimation of age, and fairly stable hunting pressure, our results support the use of deterministic cohort analysis as a relative index of abundance for monitoring red deer over time, even in the absence of exact information about natural mortality. Under violation of assumptions, however, the performance of deterministic reconstruction should be carefully inspected at the management scale.

Tracking variations in population abundance is fundamental for wildlife management and conservation (Lawton 1993;Yoccoz et al. 2001). To this aim, absolute or relative estimates of population size are needed, though the trade-off between their costs and benefits must be considered (Msoffe et al. 2007). For example, absolute estimates of population size can be obtained through statistical methods such as capture-mark-recapture (Schwarz and Seber 1999) or distance sampling (Buckland et al. 2015), which, however, may be expensive and time-consuming. Less expensive relative abundances indices (RAIs) can be obtained using simpler methodologies, for example ground counts (e.g., block counts: Corlatti et al. 2015), camera traps (Palmer et al. 2018), track counts (Fryxell et al. 1988) or ecological indicators (Morellet et al. 2007). The use of RAIs, however, requires validation against some known standard (Morellet et al. 2007).
Ungulates, and cervids in particular, are of utmost importance from the scientific, economic and conservation standpoint for consumptive and non-consumptive reasons (Mysterud et al. 2007;Putman et al. 2011). The use of techniques that allow for robust monitoring of their numerical variations over time is, therefore, crucial (Putman et al. 2011). For example, estimates of absolute densities can be obtained using mark resight (Gould et al. 2005), distance sampling (Focardi et al. 2002) and pellet-group counts (Fattorini et al. 2011). In open areas, ground counts may also be employed (Putman et al. 2011), though direct counts often have imperfect detection probabilities; therefore, they tend to underestimate population size and are best seen as RAIs 1 3 (cf. Corlatti et al. 2016). An efficient method to track relative abundance of cervids over time may be through spring spotlight counts, at least in mountainous environments (Garel et al. 2010;Corlatti et al. 2016). Alternatively, the use of individual age-at-harvest data (catch-at-age data) in cohort population analysis can also be employed to track variations in population abundance (Fry 1949;Pope 1972;Fryxell et al. 1988).
A cohort is defined as a set of individuals sharing the same characteristics at the same time and includes individuals of the same age (Manson and Wolfinger 2001). Cohort analysis is a retrospective analysis of abundance (Solberg et al. 1999) and assumes that individuals belonging to a given cohort that die either from natural mortality or human causes (typically harvesting) are recovered. The age at death can then be used to reconstruct population size, as the number of individuals of a specific age will correspond to the sum of the number of individuals from that age and of the individuals in subsequent years (Ueno et al. 2009). The reconstruction of population size through cohort analysis was originally developed in fishery (Fry 1949), but it has been extended to other taxa, including ungulates such as red deer Cervus elaphus (Lowe 1969;Mysterud et al. 2007), white-tailed deer Odocoilenus virginianus (Fryxell et al. 1991) and chamois Rupicapra rupicapra (Reiner et al. 2020). In its simplest form, this deterministic method assumes that the population is closed, all mortality events are recorded and that, for each dead animal, exact age is known. When these basic assumptions are met, cohort analysis returns absolute abundance. Deterministic approaches, however, may be very sensitive to violation of assumptions (Millspaugh et al. 2009) and in the absence of full knowledge on the number of individuals that die of natural causes, estimates of natural mortality can be obtained to assess true abundance (Fryxell et al. 1988;Solberg et al. 1999). Furthermore, in the absence of exact knowledge on age, age classes can be used (Skalski et al. 2012). Uncertainty on natural mortality and exact age, therefore, requires to step from a deterministic approach to more robust statistical approaches (Skalski et al. 2005;Millspaugh et al. 2009) which, however, increase modeling complexity.
Given that simple methodological approaches are most desirable for management purposes, in this study we take advantage of a 38-year time span and 18,390 individuals of red deer harvested and found dead in the Italian Alps to assess the robustness of deterministic cohort analysis to track population size over time, in the absence of exact information about natural mortality and age estimation. To this aim, we reconstructed the population abundance of red deer from 1982 to 2020 based on visual age assessment and without adjusting for unknown natural mortality; under these circumstances, deterministic cohort analysis reduces to a relative abundance index and requires validation against some known standard applied to the same population over the same time span. In particular, spring spotlight counts were used as a benchmark (cf. Corlatti et al. 2016) to estimate if, and for how many years, population abundance can be reliably tracked using deterministic cohort analysis under violation of basic assumptions.
The study area was located in the north-western part of the Autonomous Province of Trento, Central Italian Alps (10.78613 E, 46.32888 N; cf. Fig. 1 in Bonardi et al. 2017). The boundaries of the area include the Trentino sector of the Stelvio National Park (28% of the area) and the hunting district of Val di Sole (72% of the area). The area extends for about 640 km 2 , from 500 to 3700 m a.s.l., and encompasses the distribution range of a red deer population, whose movements were monitored using radio tracking (Bonardi et al. 2017). The climate is alpine, with mean temperature ranging between − 0.4° in winter and 13.2° in summer and yearly mean precipitation of 900 mm. The lower part of the area is forested, while above the tree line, the habitat is mainly composed of alpine and subalpine meadows (Bonardi et al. 2017). At the beginning of the nineteenth century, red deer was considered extinct in the southern part of the Alps, mainly due to hunting pressure. Red deer naturally recolonized these areas, and rapidly increased their number; the current density in the study area is about 8 individuals/km 2 in summer, with seasonal peaks of about 20 individuals/ km 2 in winter when individuals concentrate locally (Bonardi et al. 2017). In addition to red deer, three other ungulate species are present within the study area: roe deer Capreolus capreolus, chamois and ibex Capra ibex. Large predators are rare, with stable presence of wolf Canis lupus only in recent years, and brown bear Ursus arctos. Hunting is not allowed in the protected portion of the area (i.e., within the Stelvio National Park) but it is allowed in Val di Sole. Within the Val di Sole hunting district, hunting takes place from September to December. The hunting quota is fairly stable and never exceeds 35% of the red deer population size estimated by spring spotlight counts, and it is biased towards younger animals (1 year of age), while sex ratio is unbiased (1:1).
To reconstruct red deer population size with cohort analysis, we used three different types of data sources: natural mortality, hunting and car accident data collected between 1982 and 2020. Information was collected by hunters during the hunting season and by forest rangers of the Autonomous Province of Trento and of the Stelvio National Park during ordinary work and surveys in the field (cf. Bonardi et al. 2017). For each animal, age was assessed by inspecting tooth eruption and wear; the accuracy of age determination was assessed by counting cementum layers on the first inner molar of a subsample of 237 individuals (121 males and 116 females) (De Marinis 2018). Each individual feeds a global dataset, updated annually, which includes animal ID, date of finding or shooting, location, causes of death, sex, age and age class. In some cases, for example when only the skeleton without the skull was found, assessing the age of the individual was difficult. The age of unknown individuals (10% of the entire dataset) was assigned based on the proportion of dead individuals of a given age in the same year. The retrospective reconstruction of population size followed Mayle et al. (1999), excluding calves from the reconstructed data within each year to ensure comparability with spring counts (see below).
To assess the goodness of variation in red deer population size reconstructed with cohort analysis, the number of red deer counted using spring spotlight counts was used as a benchmark (cf. Corlatti et al. 2016). Each year between 1982 and 2020, in April-May, at least three surveys were conducted on separate days. During each survey, from 11.00 pm to 3.00 am, three operators drove along the same pre-defined routes in the lower part of the study area, where deer clump during spring nights. In addition to the driver, one operator was in charge of spotlighting, while a third noted down the number and composition of the animal groups. For analysis, only the survey with the greatest number of animals was used (cf. Bonardi et al. 2017). As raw counts usually underestimate true population size, data were adjusted for an underestimate of 0.35 to return absolute estimates (the coefficient of underestimation was determined by comparing spotlight counts with mark-resight estimates, cf. Bonardi et al. 2017). Furthermore, to disentangle process and observation errors, adjusted counts were filtered using a statespace model (cf. Bonardi et al. 2017) and a Kalman filter approach (Dennis et al. 2006) with the package 'MARSS' (Holmes et al. 2012(Holmes et al. , 2013. At the time of spring spotlighting, calves are not yet born; therefore, counts only included individuals of 1+ years of age. To compare the two methods, cohort analysis vs. spotlight counts, we used the Buishand range test (Buishand 1984) in the R-package 'trend' (Pohlert 2020) with R v. 4.0.4 (R Core Team 2020) in RStudio v. 1.3.1056 (RStudio Team 2020). This non-parametric test allows to compare the variation of the mean of a time series, where the null hypothesis of no variation, i.e., Δ = 0 (homogeneous time series) is tested against the alternative hypothesis Δ ≠ 0 (i.e., there is a shift in the mean of the time series). As applied to our study, the null hypothesis assumes that the annual difference between spring spotlight counts and cohort analysis remains constant over time ( Δ = 0 ). Therefore, the test would return the year when the series diverge ( Δ ≠ 0 ). To this aim, we calculated the rescaled adjusted partial sum (S k ), i.e., the deviation from the mean, as: where k refers to the number of observations of a time series x 1 , x 2 , … x n with mean x , and x i refers to the difference between filtered spring spotlight counts and number of individuals from cohort analysis. When S k ≅ 0 , the series are homogeneous without any divergence (Jaiswal et al. 2015). The shift can be assessed by computing the rescaled adjusted range (Rb) as follows: Rb = max (Sk)−min (Sk) x . The ratio between Rb and √ n can be compared using pre-defined values given by Buishand (1984) and the p value is estimated via Monte Carlo simulations using 2000 replicates (Pohlert 2020). The value returned by the test defines when filtered spring spotlight counts and number of individuals from deterministic reconstruction diverge, thereby indicating for how many years the population abundance obtained through cohort analysis is reliable, as compared to spotlight counts. From the change point onward, the data obtained through reconstruction are assumed to no longer reliably reflect population size. To further support the results of the Buishand test, we compared the two monitoring methods with a linear regression between reconstruction data (response variable) and filtered spring spotlight data (explanatory variable), including only the years before the point of divergence. Preliminary data exploration showed that the ordinary least squares (OLS) assumption of constant variance in the residuals was violated, thus leading to heteroskedasticity issues. Therefore, we decided to fit a linear model using a weighted least squares (WLS) approach with weights = 1 (fitted values) 2 , where the fitted values were obtained from a regression of absolute residuals vs. fitted values of the corresponding OLS model (cf. Cohen et al. 2003). We also compared the linear WLS model with a quadratic WLS model to test for saturation effect (cf. Garel et al. 2010) using the AICc (Akaike Information Criterion corrected for small samples) (Hurvich and Tsai 1989).
Between 1982 and 2019, a total of 18,390 (yearly min = 115 in 2019, yearly max = 1841 in 2000) individuals were recorded for use in cohort analysis. Among these, 13,760 died from hunting, 3832 died from natural mortality (mostly due to severe winters; Bonardi et al. 2017) and 798 died from accidents with vehicles. Visual inspection of tooth wear and eruption overestimated the exact age with a mean error of 1.1 years (SD ± 1.5). Figures 1S-2S-3S in the supplementary file show details of the sex-and age structure of the animals used for cohort analysis. Yearly adjusted abundance estimates using spring spotlight counts were very similar to filtered data, and varied between 491 and 2977 (Fig. 1).
The Buishand range test suggested the year 2013 (i.e., the 31st year of monitoring) as the change point (Rb = 2.107; p < 0.001) (Fig. 2). The AICc values of the two WLS models were similar (linear WLS model: AICc = 402.35; quadratic WLS model: AICc = 403.13). With similar predictive power, the simplest model should be selected (Burnham and Anderson 2002). Thus, we selected the linear WLS model as the best model, which showed a significant positive relationship ( = 0.81;p < 0.001) between the outcomes of the two estimation methods and explained much of the variance in the reconstructed data R 2 = 0.96 (Fig. 3). The residual diagnostics of the selected model did not show major violations of the main assumptions of linear models.
To obtain reliable results with cohort analysis, rules of thumb suggest that data should include an age where only 1% of a cohort is still alive in the population (Solberg et al. 1999). Therefore, some authors used a cut-off value to estimate population size, which for red deer corresponds to 7 years in males and 14 years in females (Solberg et al. 1999;Mysterud et al. 2007). Consistently, and without using any pre-defined cut-off value, our data showed that deterministic cohort analysis can be successfully used to track red deer population abundance up to 7 years in the past, despite the violation of some crucial assumptions.
In our study, it was not possible to detect all individuals that died in a given year; thus, the reconstruction of population size through cohort analysis did not return the 'true' size of the population. While statistical methodologies can be used to estimate natural mortality, thus returning estimates of absolute abundance (Skalski et al. 2005), these approaches are more computationally demanding than simple, deterministic cohort analysis. Our study suggests that when population variation, rather than absolute population size, is of interest, deterministic cohort analysis can still return a reliable index of relative abundance Red deer abundance from spring spotlight counts Red deer abundance from cohort analysis Fig. 3 Relationship between abundance from reconstruction and abundance from spring spotlight counts for the red deer population in Val di Sole and in the Stelvio National Park between 1982 and 2013. Gray shaded area represents 95% confidence interval even in absence of exact data on natural mortality. Furthermore, while in other ungulate taxa, such as bovids, age determination is easy and exact, thereby allowing for a reliable population reconstruction (cf. Reiner et al. 2020), age assessment in cervids is problematic. In our study, age estimation through tooth eruption and wear proved relatively accurate. In addition, the red deer is a long-lived animal, as individuals may live up to 20 years and have a generation length of 14 years (Pacifici et al. 2013) and small errors in age estimation of long-lived animals possibly have negligible consequences on population reconstruction. Under these circumstances, our results support the use of a simple, deterministic population reconstruction approach as a RAI to investigate temporal variation in population abundance, despite the violation of basic methological assumptions.
This information is particularly important for wildlife managers, as simple cohort analysis may, for example, allow to retrospectively reconstruct the sex-and agespecific vital rates of a population, thereby helping the decision-making process without the need for complex modeling. Great care, however, is needed when interpreting our results: deterministic approaches may be very sensitive to violation of assumptions and scale of analysis (Millspaugh et al. 2009), and their results may depend on hunting regimes (cf. Reiner et al. 2020) as well as on the proportion of animals that died of natural causes. Furthermore, the robustness of our deterministic cohort analysis may owe to relatively high accuracy in age determination and long lifespan of the target population. The same approach may not be robust if age estimation has lower accuracy or the species has a shorter lifespan. More generally, the performance of deterministic cohort analysis under violation of basic assumptions should always be validated at the local scale against some known standard, before it is used for management purposes.