Time series analysis plays an important role in the detection of mechanisms that drive population fluctuations. However, long time series are rare, with ungulate data sets usually not exceeding 50 years. In this article, we describe a long-term data set of population density indices of five ungulate species obtained from the analysis of bag records collected in the Castelporziano Preserve, Rome, Italy. Hunting statistics are often used as proxies for population density; however, in the case of long time series for large mammals, there are no comparative studies to assess the validity of such data. We evaluated the ungulate time series, using two different approaches: we 1) compared hunting statistics with independent animal counts, and 2) assessed whether or not habitat composition of the drive areas was representative of habitat availability in the whole estate. Regressions between bag data and animal counts gave significant results only for three species, whereas bag data corrected for hunted area were significantly correlated to animal counts for all five species. The results suggest that use of bag records not corrected for hunting effort and without any previous validation could lead to misleading estimates of abundance indices. Finally, our analysis showed that density indices of the five species were not significantly affected by the selection of habitats where hunting drives were organised. Our data set may contribute to the understanding of ungulate ecology in the Mediterranean environment.
Time series analyses play an important role in the detection of mechanisms leading animal populations to fluctuate in space and time (Royama 1992, Turchin 2003).
For ungulate populations, long-term time series are needed for at least two reasons. First, available data sets are often represented by population indices rather than by absolute measures of population size. Since these indices are prone to observation errors that enhance noise, it may be difficult to make inferences about patterns of population change with short data sets (Solow & Steele 1990, Dennis et al. 2006). Second, if we are investigating second-order dynamics (caused by feedbacks with a time lag or delayed density-dependent mechanisms), we need a data set longer than a complete cycle (Turchin 2003). Ungulates can exhibit cycles with periods of several decades. For example, Peterson et al. (1984) and McLaren & Peterson (1994) showed a 16-34 year cycle for a moose Alces alces population in North America, and Ogutu & Owen-Smith (2005) found cyclic variability with half-periods of 10-18 years in 12 African ungulate species. Finally, longer cycles (30-50 years) may be predicted by theoretical models involving tri-trophic (vegetation-deer-predator) interactions (Turchin 2003).
However, long-term time series are seldom available, with most data sets usually not exceeding 50 years (see for example Brook & Bradshaw 2006). Only bag records are available for longer periods, such as the data set described by Jones (1914) concerning deer dynamics in the 18th and 19th century in North America or the Greenland hunting records of caribou Rangifer tarandus available since 1900 (Meldgaard 1986, Forchhammer et al. 2002).
Hunting statistics are often used as proxies for population density (Grøtan et al. 2005, Mysterud & Østbye 2006), but attention must be paid to variations of effort (e.g. the number of hunters), and to the relationship between bag data and actual animal abundances. Strictly speaking, an index of abundance, such as any harvest index, should be linearly correlated to actual abundances in order to yield unbiased estimates of demographic parameters (such as population trends; Gibbs 2000). Some attempts to validate harvest indices using independent sources of information have been made for several species of deer, with data covering short time-spans (Grøtan et al. 2005, Mysterud et al. 2007). To our knowledge, there are no comparative studies using long-term time series for large mammals to assess the validity of hunting statistics.
In our paper, we describe a long-term data set of population density indices for five ungulate species obtained from the hunting records collected in the Castelporziano Preserve in Rome, Italy, during the period 1878-1986. The studied guild included two endemic Italian subspecies, i.e. the Italian roe deer Capreolus capreolus italicus (Festa 1925), the Maremma wild boar Sus scrofa majori (de Beaux & Festa 1927), the reintroduced red deer Cervus elaphus, and the two introduced species fallow deer Dama dama and the Asian antelope nilgai Boselaphus tragocamelus, which is now extinct in the area. The ecology of this ungulate community in the Castelporziano Preserve is well known since it has been studied as part of a project of global monitoring of the area, which is quite important for conservation of the biodiversity of Mediterranean habitats (Scarascia Mugnozza 2001). Other features of the Castelporziano Preserve make it an ideal site to obtain long-term population records for ungulates. First, the ungulate management system has remained virtually unchanged for more than one century, based on strictly controlled hunting parties whose records have been conserved with few losses. Second, an independent set of ungulate count data is available for a relatively long period (i.e. 37-51 years depending on the species), which could be used as an independent source of information to validate bag records. Finally, environmental changes (e.g. vegetation cover), which represent important covariates of animal population dynamics, are quite well documented since the late 19th century.
The aim of our paper is to evaluate the collected time series and to document the presence and intensity of sampling biases in the bag data. We used two different approaches:
We compared hunting statistics (total number of killed animals and density indices) with independent animal counts.
Since the different species are known to prefer specific habitat types, a non-random habitat use by hunters could result in over- or underestimation of population indices. To estimate possible sources of bias, we compared the yearly habitat use by hunters (i.e. the proportion of each habitat sampled in hunting drives) with the availability of each habitat type in the whole estate.
Material and methods
Study area
The Castelporziano Preserve, situated about 20 km southeast of Rome, Italy, is a completely fenced area of 48 km2, excluding the area of Capocotta (Fig. 1). The climate is Mediterranean, with a period of summer drought during May-August and maximum rainfall in November and December (Blasi 1994).
Historical background
In 1872, the estate was purchased by the Italian government and managed primarily as a Royal hunting reserve, with some animal husbandry and agricultural activity. At that time three species of ungulates were present: roe deer, fallow deer and wild boar. In 1889, red deer and nilgai were introduced to Castelporziano. During the First World War (1915-1918), royal hunts were suppressed and few drives were performed by rangers. During the Second World War (1940-1945), hunting was very intense to satisfy the demand of meat for charities. In 1943, the estate was occupied by the German army for about six months, and for the first three of these months we lack detailed records of bags. An estimated total bag of 1,000 head was given later by the head keeper. During the German military occupation nilgai and red deer became extinct.
In 1948, Castelporziano became a hunting reserve for the President of the Republic, but the wildlife management remained largely unchanged. After the many wildfires that ravaged the area during the war, intense forestation took place to recover the habitat. Red deer were reintroduced again in 1958, but this new population was not hunted and hence is not included in our investigation.
No large predator was present in the area during our study period.
Vegetation
The plant communities currently present in Castelporziano were described by Pignatti et al. (2001). The holly oak Quercus ilex grove represents the typical meso-mediterranean evergreen wood of the study area (27%), characterised by a rich understory. On recent dunes close to the sandy shore, this habitat assumes the aspect of a typical Mediterranean maquis. The deciduous oak wood (34%) consists of Q. cerris and Q. frainetto, in most cases associated with a high coverage undergrowth of Oriental hornbeam Carpinus orientalis. The few open areas in Castelporziano (8.3%) are largely characterised by arid pastures, often degraded by overgrazing. Commercial stands (21%) consist of pure or mixed (with holly oak) coetaneous woods of domestic pine Pinus pinea or cork oak Q. suber plantations. These habitats are invaded by the typical meso-mediterranean undergrowth. A few Eucalyptus sp. and Populus sp. plantations are also present.
We used the present land cover of the Preserve and, in particular, the habitat classification of Grignetti et al. (1997) as a template for interpretation of past time data. The first datum available is a cadastral map (dated 1867) annexed to the bill of sale of the Preserve to the Italian government. This is a freehand watercolour map with an accompanying table reporting the surface areas of two vegetation classes (wooded and open areas). For later periods, we used six aerial photos taken at about 10-year intervals from 1930 to 1980 (aerial photographs refer to sheets 149 II NE (Acilia) and 149 II SE (Castelporziano), series 25V (1:25,000) of IGM maps). Each of the six aerial photos were interpreted with the aid of a stereoscope to reconstruct the land cover (at forest parcel scale) of Castelporziano through time (Petit & Lambii 2002, Velázquez et al. 2001; see Fig. 1).
We used seven categories: 1) mixed oak wood, 2) cork oak wood - both natural and reforested, 3) maquis, 4) pastures and cultivated fields, 5) clearings and open areas - including uncultivated fields and scantily wooded areas, 6) pine plantations and 7) ‘other artificial woods’. Buildings and wetlands were not considered. In our analysis, we pooled mixed oak woods and maquis (collectively referred to as ‘natural woods’), while pastures and cultivated fields were grouped with clearings and uncultivated fields to form the ‘open areas’ class.
Data collection
Our data sources are the records conserved in the national archives of Rome (i.e. Central State Archive with records from 1878 to 1945 and Quirinale Library and Castelporziano Preserve Archive with records from 1946 to 1986).
Our basic sampling unit is the number of animals of each species killed during each hunting drive. Drives started in 1878 and usually 10-12 drives were organised each year. Hunting was carried out from November of year t to March of year t+1, and each hunting season is referred to using index t. The area to be driven was surrounded by mobile curtains to reduce the possibility of animals escaping. Hunters arranged in a single line shot all the animals driven forward by dogs and personnel. Figures on the proportion of animals that escaped shooting were unavailable, but gamekeepers reported cases when unusual numbers of animals fled. The participants in a hunting drive (e.g. guests of the King or of the President, or gamekeepers) were usually numerous (up to 60), as were the personnel driving the animals (up to 30). Unfortunately, the number of participating hunters and gamekeepers was not recorded in all occasions. Most drives were documented by sketches made by rangers, indicating the location of the drive, and the position of the curtains and hunters. Most of the areas were used year after year (median time between subsequent use of the same drive was 395 days), but rarely more than once per hunting season. The boundaries of 92 drives (corresponding to 1,099 of 1,430 hunting reports available in the archives) were determined and their surface area computed (average: 65.2 ± 4 ha). In the case of repeated hunting in a single plot during two or more consecutive days (N = 37), we pooled all killed animals. Drives that were judged by the gamekeepers as unsuccessfully carried out (e.g. because of bad weather or a small number of hunters; N = 12) were discarded. In total, we used 1,087 records (1878-1986) for the analyses. Unfortunately, information about sex, age and weight of the shot animals was not collected.
Using the habitat cover maps, we could compute the habitat composition of each drive, and which land cover map was used to establish the habitat composition for each hunting period is shown in Table 1.
A small number of animals were shot out of the hunting drives, especially by the King himself. The total number of specimens harvested every year, both during drives and free hunting (bag records) was carefully reported in the registers, even during the World Wars.
Rough counts of the ungulate populations were available for some periods. There were two types of surveys: 1) counts of all ungulate species performed by 5-8 pairs of gamekeepers from vantage points during October-November from 1906 to 1942, and 2) wild boar counts carried out during summer (July-August) at supplemental feeding sites from 1892 to 1919. Ungulate counts were carried out in one day, from several vantage points placed in open areas. Gamekeepers surveyed almost the same areas every year, however, the exact number of vantage points (for ungulate counts) and feeding sites (for wild boar counts) used is unknown. Both methods are still used as cheap monitoring systems (starting in 1988), and recent analyses showed a good correlation with population estimates based on Distance Sampling and Capture-Mark-Recapture (Focardi et al. 2002, 2005, 2008). For wild boar, we used the two types of count data separately in the analysis.
Data analysis
We calculated density indices for each species as the number of animals killed during the hunting season per km2 of driven area. To avoid zero data when applying natural logarithms, one individual of each species was added to the total number of killed animals. We calculated confidence intervals according to the derivation of Burnham et al. (1987) for log-normal distributions.
We performed linear regressions between raw bag records (total number of killed animals in each hunting season) and absolute numbers of counted animals (on all vantage points or feeding sites). Furthermore, we carried out linear regressions between density indices and count data, after transforming both data sets into natural logarithms.
To determine gamekeepers' habitat selection of drive area, we used the Ivlev electivity index (Krebs 1989). Electivity (E) for habitat I is
where ri is the proportion of the ith habitat type driven and pi is the proportion of occurrence of the ith habitat type in the environment. E ranges from −1.0 to +1.0, with positive values indicating preference, negative values indicating avoidance, and values close to 0 indicating no selection. We estimated the availability of each habitat type from its proportion in the pertinent land cover map (see Table 1). We computed the use by hunters as the proportion of each habitat type inside the driven plots used in that season. Hunting densities of the five species were regressed against Ivlev electivity indices for natural woods (Enw) and open areas (Eoa). To check for possible bias in the regression results due to collinearity between the two independent variables, we computed the Variance Inflation Factor (VIF) according to Neter et al. (1990).
To control the variation in population density of the five species while investigating the effect of habitat selection by hunters on density indices, we also performed a multiple regression between hunting densities and both counted animals and electivity indices for natural woods and open areas (for the years from which count data were available).
With time series data, the residuals derived by an ordinary regression are often correlated over time, hence the regression assumptions are violated. In particular, positive autocorrelation of the errors generally tends to make the estimate of the error variance too small, so that confidence intervals are too narrow and true null hypotheses are rejected with a higher probability than the one stated by the selected significance level. Contrary to this, negative autocorrelation of the errors generally tends to make the estimate of the error variance too large, so that confidence intervals are too wide and the power of the tests is reduced. We performed a Durbin-Watson test to detect positive/negative autocorrelation of residuals. When the test gave significant results (P < 0.01), we used the Yule-Walker method to correct for autocorrelation (Gallant & Goebel 1976).
We analysed data using SAS 9.1 (SAS Institute Inc., Cary, North Carolina, USA) and performed manipulation of geographical data using ArcView 3.2 (Environmental Systems Research Institute, Redlands, California, USA).
Results
Time series
We reconstructed the time series of hunting densities of the five ungulate species (Fig. 2). We could not obtain reliable indices of population sizes for 9 of the 109 years (8.3%), and data were missing mainly during the two World Wars. Fallow deer (see Fig. 2A) reached the highest densities among all species (about 100-300 individuals/km2) during 1915-1940, but the population rapidly dropped off during the Second World War. After this war, the fallow deer population increased again, but its abundance appears to be quite variable. Roe deer and wild boar (see Fig. 2B,C) were more abundant (3-30 individuals/km2 and 13-133 individuals/km2, respectively) in the 19th century up to 1915. After a period of decline between the two World Wars, the two populations recovered (especially wild boar), but never reached previous densities. Finally, red deer and nilgai (see Fig. 2D,E) always occurred in relatively low densities (up to 12 individuals/km2 for both).
Comparison with count data
The total numbers of killed and counted animals were statistically correlated only for fallow deer, roe deer and wild boar (Table 2). We found positive autocorrelation of residuals for three species (fallow deer, red deer and nilgai). However, the results of the regressions were not influenced by the correction. The (ordinary) regression coefficient of harvest as a function of counts (see Table 2) was significantly > 1 for fallow deer (one-sample t-test: t = 3.81, P = 0.0006, N = 35). Since the y-axis intercept was significantly negative (-1598.3 ± 291.36) for this species, the total harvest of fallow deer was not a fixed percentage of counts. In other words, harvest was more intense when counts were higher. On the contrary, for the other species, the regression coefficient (see Table 2) was < 1 (roe deer: t = −409.9, P < 0.0001, N = 35; wild boar autumn counts: t = −25.1, P < 0.0001, N = 35) and the y-axis intercept was near 0 (roe deer: 1.7 ± 2.05; wild boar: −19.5 ± 39.6).
Density indices obtained from bag records were statistically correlated to counts for all species (Table 3, Fig. 3), but for fallow deer, roe deer and nilgai, regression errors were positively autocorrelated. Again, the regression coefficient for fallow deer was higher than for the other species.
Electivity index
Castelporziano underwent significant habitat modifications over time (Fig. 4). Natural woods and clearings/uncultivated lands constantly decreased (from 71% in 1930 to 58% in 1954, and from 11% in 1930 to 5% in 1959, respectively) in favour of pastures and cultivated fields. From the late 1950s, however, many open areas were reforested with pines and cork oaks, covering up to 13% and 10% of the Preserve in 1980, respectively. Other artificial woods were less relevant (< 2% of the total area in 1980). Natural woods recovered slightly in the last years of the investigated period, covering about 60% of the study area in 1980.
Also the use of each habitat type by gamekeepers varied among years (see Fig. 4). Natural woods was the most used category (70.2% ± 4.4 per year during the first period; 47.7% ± 1.9 and 25.6% ± 1.6 for mixed oak wood and maquis during all other periods). Together, cultivated and uncultivated fields represented 14% ± 1.6 of the total hunted area, pine plantations and cork oak woods had an average proportion of 7.3% ± 1.7 and 9.3% ± 1.3, respectively, and other artificial woods represented only 2% of the total hunted area.
The use by hunters of natural woods and open areas was not always proportional to their availability. The mean values of the electivity indices (Enw = 0.05 ± 0.018; Eoa = −0.42 ± 0.03) suggested, on average, a use of natural woodlands proportional to their availability and avoidance of open areas.
All multiple regressions of hunting densities as a function of Enw and Eoa presented a highly significant and positive autocorrelation of errors. Yule-Walker estimates showed that none of the ungulate species densities was dependent on the annual electivity index values, except for a negative effect of open areas index on wild boar population (Table 4). However, the variance of wild boar density accounted for by electivity was < 5%. For all regressions, the estimated variance inflation factor values were at maximum 1.6, well below the value of 10, indicating the presence of multicollinearity problems.
When including count data in the model, multiple regressions displayed statistically significant effects, but the contribution of the electivity indices was limited (Table 5); only for roe deer we obtained a negative effect of open areas. For wild boar, only summer counts were considered, since they showed the best agreement with the density indices of the species (see Table 3).
Discussion
The availability of reliable estimates of yearly hunting effort, coupled with independent count data, allowed us to evaluate the density indices of five ungulate species, obtained from time series of hunting data, and to explore possible sources of bias.
As it often happens with time series of population densities, most of the regression analyses showed a positive autocorrelation of residuals. This has to be considered, since with ordinary regressions we run the risk of rejecting too often a null hypothesis. Results of our regressions did not change significantly after autocorrelation correction, most likely because of the relatively large sample size.
Raw bag records and count data were closely associated only for the most abundant species (fallow deer, roe deer and wild boar). This probably depended on the hunting strategy adopted in the Preserve. Harvests were planned before the hunting season on the basis of count outcomes, but usually only counts of the species with the highest densities were considered. Fallow deer were particularly pursued by hunters and the proportion of killed animals rose as the number of counted animals increased. In contrast, roe deer and wild boar harvests represented a fixed percentage of counts. For the two species with lower densities, red deer and nilgai, the numbers of killed animals did not reflect actual population variations. We showed that we introduced systematic biases, if we did not correct the bag for the hunting effort.
Hunting records, not corrected for hunting effort, are one of the primary sources of data in vertebrate population ecology (see for example Brook & Bradshaw 2006), but care must be taken when using them without any previous validation with independent information. We agree with Pettorelli et al. (2007) that bag record statistics cannot be used as surrogates of density based on the assumption that they are unbiased; the assumption must first be tested.
The high correlations between counts and hunting densities suggest that harvest data, corrected for hunted area, can be used as a reliable proxy for population abundance. However, the relationship between the two variables is not simple. We obtained large R2 values when log-counts and log-density indices were compared, indicating a non-linear relationship between counts and bags. Similar findings have been described by Cattadori et al. (2003) for red grouse Lagopus lagopus. This effect could depend on a bias in count data (difficulty in detecting animals at low densities) or in hunting (inability at high density to shoot all the animals crossing the hunters' line). Hunting density is an underestimate of actual density, since some of the animals can escape from shooting, and this underestimation could increase when density is higher. Hunted area can be considered an index of hunting effort where the main hunting technique is the drive. In other cases, other measures of hunting effort (such as the number of hunters or the number of hunting days) could be useful for adjusting bag data, depending on the management practices. Given the peculiarity of the study area, more studies are requested in other situations (non-fenced areas, different habitats and hunting techniques) in order to establish whether bag data corrected for hunting effort can be generally considered a good index of population density for ungulates.
When bag data are used as a proxy for animal density, bias can also arise from non-representative habitat sampling in hunting drives, particularly when the studied species have strong preferences for different habitats. In our study area, open areas appeared to be undersampled by hunters, and this may have biased the population estimates. The density indices that we calculated from bag records can therefore be regarded as minimum densities in the wooded areas of the Preserve. We presume that avoidance of open areas was partly due to the need to maximise hunting yield, since drives were performed only during daylight when ungulates move and rest inside wooded areas. In this regard, avoidance of open fields by hunters could be partly compensated for by ungulate behaviour, but exact figures to estimate this bias are not available. However, we did not detect any effect of changes in hunters' electivity on hunting densities (except for a negative effect of the Ivlev index for open areas on wild boar density index), suggesting that the ‘sampling strategy’ of gamekeepers was unlikely to produce strong biases in population indices. Even when we accounted for variation in population densities by including count data in the multiple regressions, we did not find any effect of hunters' habitat selection. The ‘fair’ sampling by gamekeepers was related to the fact that almost all the preserve was used for hunting and that gamekeepers rotated hunting among the different zones.
Our data set has the potential to contribute to the understanding of ungulate ecology in the Mediterranean area, which is poorly known compared with the ecology of populations living in temperate or boreal environments or in open ranges in the tropics. The presence of wild ungulates in Mediterranean ecosystems is in fact quite limited, so that there are few possibilities to study the population dynamics of these communities (but see Gaillard et al. 2003, Toïgo et al. 2006). Long-term studies are required to understand factors affecting ungulate populations, and long-term hunting data may be an important starting point to investigate their dynamics.
Acknowledgments
Acknowledgements - we would like to thank the Preserve of Castelporziano, in particular the Director Dr. Demichelis for providing support and accommodation, A. Tinelli and the Osservatorio per gli Ecosistemi Mediterranei for finding additional material and aerial photographs, and the gamekeepers for useful information. We are also grateful to the employees of the archives for providing consultation.
References
Figure 1.
The Castelporziano Preserve and the Capocotta area. The land cover obtained from the aerial photo taken in 1980 is shown only for the Castelporziano area.
Figure 2.
Time series of hunting densities of the five species fallow deer (A), roe deer (B), wild boar (C), red deer (D) and nilgai (E). Population densities (individuals/km2) are transformed into natural logarithms and dashed lines indicate 95% confidence intervals.
Figure 3.
Log-linear relationship between density indices (D) and counted animals (C) of the five species fallow deer (A), roe deer (B), wild boar (C autumn and D summer), red deer (E) and nilgai (F).
Figure 4.
Availability (A) and average use by gamekeepers for drive hunts (U) of each vegetation type during each period.
Table 1.
Periods of the time series and the corresponding land cover maps.
Table 2.
Analysis of the relationship between the total number of killed animals and the number of the animals counted expressed by number of observations (N), coefficient of determination (R2), probability value (P) and regression coefficient (for both ordinary regressions and regressions with autocorrelation correction). In the Durbin-Watson test the Durbin-Watson statistic (DW), the probability values for testing positive autocorrelation (P < DW) and probability value for testing negative autocorrelation (P > DW) are reported; when P < 0.01, the regression with autocorrelation correction (Yule-Walker method) was performed. The number of observations refer to the number of years for which both kinds of data were available. For wild boar, results of regressions with both autumn and summer counts are given.
Table 3.
Analysis of the relationship between the log-transformed hunting densities and the log-transformed number of the animals counted, expressed by number of observations (N), coefficient of determination (R2), probability value (P) and regression coefficient (for both ordinary regressions and regressions with autocorrelation correction). In the Durbin-Watson test the Durbin-Watson statistic (DW), the probability values for testing positive autocorrelation (P < DW), and probability value for testing negative autocorrelation (P > DW) are reported; when P < 0.01, the regression with autocorrelation correction (Yule-Walker method) was performed. The number of observations refer to the number of years for which both kinds of data were available. For wild boar, results of regressions with both autumn and summer counts are given.
Table 4.
Results of multiple regressions with autocorrelation correction (Yule-Walker method): log-transformed hunting densities vs Ivlev index for natural woods and open areas. Probability values of each independent variable (Enw and Eoa) are also given.
Table 5.
Results of multiple regressions: log-transformed hunting densities vs log-transformed counted animals (c), Ivlev index for natural woods (Enw) and open areas (Eoa). For each species either the ordinary regression (when residuals are not autocorrelated) or the Yule-Walker regression (when residuals of ordinary regression are positively autocorrelated) is shown. Probability values of each independent variable are also given. For wild boar, only summer counts are considered.
