Surveillance and habitat diversity affect European brown hare (Lepus europaeus) density in protected breeding areas

The European brown hare (Lepus europaeus) is an important game species throughout Europe. In Italy, for preventing the introduction of allochthonous strains, the management of brown hare populations has focused on the establishment of small protected areas (ZRCs), appositely managed for disposing of wild-born hares for restocking hunting territories. We investigated the effects of both land cover and surveillance on hare density and habitat preferences in 20 ZRCs, monitored twice per year (pre- and post-breeding periods) between 1997 and 2017. Density, as assessed by spotlight counts, ranged between 2.8 and 47.0 ind/km2 in spring and 5.0 and 68.4 ind/km2 in autumn. Surveillance, percent length of protected boundaries, year of institution and habitat diversity, as assessed by Shannon’s Index, were the main factors affecting hare density. During their foraging activity, hares selected ryegrass, hayfields and lucerne, while avoided maize stubble and ploughed fields and were never recorded in poplar plantations or next to human settlements. While the effects of habitat heterogeneity on hare density have been widely studied, we suggest that the involvement of local stakeholders may be of paramount importance for ensuring effective conservation measures.


Introduction
The European brown hare (Lepus europaeus) is an important small game species (Hacklander and Schai-Braun 2019;Schai-Braun et al. 2019), which, since the 1960s, has suffered such a dramatic decline in western and central Europe (e.g. Broekhuizen 1982;Pielowski 1990;Tapper 1992;Mitchell-Jones et al. 1999) to be listed as "Near Threatened" or "Threatened" in several countries (e.g. Switzerland, Germany and Austria) and in Appendix III of the Berne Convention on the Conservation of European Wildlife and Natural Habitats.
Despite having probably evolved in steppe grasslands, currently the brown hare is generally more abundant in intensively cultivated farmland than in pastures (Vaughan et al. 2003). Since the 1960s, the intensification and mechanization of agriculture have been paralleled by the reduction in hare numbers throughout Europe, pointing at habitat changes as the main cause of hare decline (Hutchings and Harris 1996;Smith et al. 2005). Several studies have related the loss of habitat diversity caused by increased field size to the decrease of food availability, with consequent higher adult mortality and smaller litter size (Frylestam 1980;Smith et al. 2005). In addition, the reduction of permanent cover caused by the removal of hedgerows and grassy field margins may locally exacerbate the impact of predators, mainly the red fox (Vulpes vulpes) (Smith et al. 2004).
Italian hare populations have followed the same trend recorded throughout Europe, showing the greatest reduction in the central and southern peninsula. The decline reached a peak in the 1980s, for the combined effects of agricultural intensification, hunting pressure and European brown hare syndrome (EBHS) (Spagnesi and Trocchi 1993). In an ineffective attempt to halt hare decline, restocking was carried out throughout the peninsula, using hares from source populations in Eastern Europe and locally causing the replacement of native Lepus europaeus meridiei by allochthonous strains, as occurred to a larger scale in France and Denmark (Suchentrunk et al. 2006). Moreover, translocations can spread pathogens and cause non-native gene introgression and mixing of different evolutionary significant units (Champagnon et al. 2012), altering the morphology and behaviour of recipient populations and reducing their ability to adapt to environmental changes.
To prevent the negative effects of restocking with captivebred or non-native hares and enhance natural dispersal, the Italian legislation has promoted the establishment of a network of small protected areas (called "Zones for Restocking and Capture" ZRCs), appositely managed for disposing of wildborn hares for restocking hunting territories through either dispersal or, usually, capture and translocation (Meriggi et al. 2015).
Thanks to ZRCs and habitat improvements supported by EC regulations (Genghini and Capizzi 2005;Canu et al. 2013), currently the negative population trend of the brown hare seems to have stopped or reversed at least in the north of Italy (Trocchi and Riga 2005). Nonetheless, hare numbers are relatively high only within the protected areas, where autumn densities range between 20 and 50 ind/km 2 , while in hunting grounds densities are still much lower (0-5 ind/km 2 ; Trocchi and Riga 2005).
The recovery of hare populations is partially hindered by illegal hunting, a major cause of hare mortality in northern Italy (Meriggi et al. 2001) as elsewhere (e.g. Misiorowska and Wasilewski 2008). The support of local hunting associations and rural communities has been reported to be crucial for assuring surveillance and monitoring of hare populations, and thus the productivity of ZRCs (Meriggi et al. 2001).
In fragmented landscapes, habitat quality of focal sites, such as ZRCs, can have a marked effect on species abundance patterns (Thomas 1994;Thornton et al. 2011), as many species are more influenced by fine-than large-scale features of the environment (Cushman and McGarigal 2004). Hares are expected to cope better in habitat mosaics, providing spatial and temporal heterogeneity in resource (food and cover) availability (Tapper and Barnes 1986;Meriggi and Verri 1990). Accordingly, in northern Italy, hare density was reported to be higher in hilly than intensively cultivated (maize or rice fields) ZRCs (Meriggi and Alieri 1989).
In this study, we investigated the effects of land use on European hare density and habitat selection in 20 ZRCs (province of Lodi, Lombardy region, N Italy), for which data were collected yearly from 1997 to 2017. We expected population density and habitat selection to be affected by both habitat diversity and size of protected areas.

Study area
The Po-Venetian plain is the largest Italian plain (ca. 46,000 km 2 ), intensively managed for the production of rice, maize, wheat, sugar beet, fruit and cattle husbandry (Ajassa et al. 1997). In its central part (Lombardy region), 77% of the plain is farming land, while only 3.65% is covered by residual woodland. Throughout the study period, the maize was by far the most widespread crop (42-54% of cultivated land), and the percent cover of maize fields still increased (on average + 0.6% per year; Balestrieri et al. 2019).
In the study period , the climate was continental to moderate sub-littoral, with a mean ± SD annual temperature of 14.1 ± 0.7°C (min monthly mean: − 5.4°C, in February 2009; max monthly mean: 35.4°C in August 2003) and mean ± SD annual rainfall of 867.9 ± 245.8 mm (calculated from daily values obtained from the historical archives of the Regional Agency for Environmental Protection).

Methods
From spring 1997 to spring 2017, hare density was assessed twice a year, during both the pre-breeding (March-April) and post-breeding periods (October-November), for a total of 532 surveys (on average 26.6 ± 10.6 surveys per ZRC; Tables 1 and 2). Monitoring was carried out using spotlight counts, a widely used and effective method for assessing hare population densities (Langbein et al. 1999). Hares were searched for from a moving car (speed: 5-10 km/h), by two observers, starting each survey ca. 2 h after sunset. Dirt roads crossing cultivated land were selected as to survey the largest possible portion of each ZRC (70-90%, depending on size) while avoiding double counts, and kept constant throughout the study period. Both sides of the road were scanned by a handle lamp (100 W; Frylestam 1981; Barnes and Tapper 1985). Both the number of observed hares and actual size of lighted sub-areas were marked on aerial orthophoto maps, recording the type of habitat where each hare was first observed. The total area surveyed during any session was calculated by summing all lighted sub-areas. Variation in hare detectability may result in the non-random distribution of the hares missed during spotlight counts (Schai-Braun et al. 2019). We assumed this bias to be negligible, since structured habitat types such as copses comprised, on average (± SE), only 0.6% ± 0.1% of our study areas. Methods and staff were consistent throughout the study period.

Habitat diversity was assessed by Shannon's index
To assess the degree of surveillance, we classified ZRCs into three types of patrolling frequency/intensity (1: patrolling by land owners and local hunters; 2: patrolling by officers of hunting associations; 3: no patrolling). We also considered the total area and year of institution of each ZRC, assuming that Table 1 Hare densities (ind/km 2 ) in spring, as recorded in 20 protected areas (ZRCs; numbers correspond to those in Fig. 1 larger and older protected areas host better preserved hare populations. Considering that, during the hunting period, dogs can be used to flush hares out of ZRCs and into hunting grounds, we calculated the proportion of the total perimeter of each ZRC delimited by watercourses, canals, paved roads or railway tracks ("Prot bound" = length of protected boundaries/total perimeter). We assumed that these landscape elements may hinder hunting dogs from penetrating into protected areas. For each ZRC, the relationship between pre-and postbreeding densities was tested by Spearman rank correlation (rho), while the effect of land use on post-breeding hare density was tested by random forest, using a subset of the data including only ZRCs and years for which both densities and habitat covers were available (N = 101). Random forest is a machine learning, non-parametric method which is highly suitable for analysing compositional data (Brűckner and Heethoff 2017) and performs well also with longitudinal settings (Adler et al. 2011;Konerman et al. 2015).
A random forest model is made up of hundreds of unpruned classification and regression trees, each trained by selecting a random bootstrap subset 'Xi' (i = bootstrap training iteration of the database X, ranging from 1 to t) and a random set of predictor variables. Predictor variables are evaluated by how much they decrease node impurity or how often they make successful predictions in the forest of classification and regression trees (Breiman 2001). Random forest regression (RFR) is of particular interest for identifying non-linear relationships amongst both continuous and categorical variables without processing (no need to rescale or normalise the inputs), thus allowing the analysis of variables that are difficult to be defined using other traditional statistical methods (Cutler et al. 2007;Siroky 2009;Vincenzi et al. 2011). We assessed variable importance by the percent Gini increase of mean square error in nodes (%IncMSE) that use randomly permuted values of a predictor in the model, and increase in node purity (IncNodePurity), which is calculated based on the reduction in sum of squared errors whenever a variable is chosen to split. The larger the increase in the node purity values, the greater the importance of the predictor (Breiman 2001). A permutation test was included to provide a significance level for each predictor, with α = 0.05 as significance threshold value. The RFR model was applied 1000 times, and the 95 percentile of the ordered distribution of node impurity values was taken to assess the significance level of each individual variable (Balestrieri et al. 2013). Partial dependence Table 2 Hare densities (ind/km 2 ) in autumn, as recorded in 20 protected areas (ZRCs; numbers correspond to those in Fig. 1)  plots were used to show how the RF model predictions were influenced by major predictors (Cutler et al. 2007).
To assess hare (N = 1653) habitat preferences during their foraging activity, the chi-squared (χ 2 ) goodness-of-fit was used to test for non-random habitat use by comparing the observed frequencies of use of the different habitat types with the expected ones (White and Garrot 1990). Expected frequencies were calculated based on the percent cover of each habitat and compared to those observed in the same years in which habitat cover was recorded. To determine whether a habitat was selected or avoided, Bonferroni simultaneous confidence intervals for the proportion of use of each habitat were calculated (White and Garrot 1990).
RFR was also used to test the relationship between mean hare densities and year of institution, area, type of surveillance and percent length of protected boundaries of ZRCs. Mean hare density was compared among the three classes of increasing surveillance by Kruskal-Wallis test.

Results
Throughout the study period hare density ranged between 2.8 and 47.0 ind/km 2 in the pre-breeding period and 5.0 and 68.4 ind/km 2 in the post-breeding period (Tables 1 and 2). Considering the 11 ZRCs for which the longest series of density estimates were available (nos. 1, 3-8, 12, 14-16 in Tables 1 and 2), post-breeding densities were significantly correlated with those of both the previous and following springs (rho = 0.44-0.87, P < 0.01 for all correlations), with only one exception (no. 16; P = 0.1).
The RFR model also included the percent cover of ryegrass (%IncMSE: P = 0.13; IncNodePurity: P = 0.20) and ploughed fields (P = 0.08 and 0.37, respectively), as top variables. For the first variable, hare density increased with increasing values of cover, while for ploughed fields optimal values ranged between 10 and 30% (Fig. 3).
During their foraging activity, hares selected ryegrass, hayfields, lucerne and fallow land, while tended to avoid maize stubble and ploughed fields and were never recorded in soy fields, poplar plantations, woods and next to urban areas or infrastructures (Fig. 4).

Discussion
Lacking the assessment of hare detectability, our survey method probably did not allow to obtain actual estimates of European hare density; nonetheless, our main aim was to find the relationship between hare abundance and land use types and we are confident that the rank order of estimates for the 20 study areas was not affected by the chosen technique (Langbein et al. 1999).
Translocations may have represented a further source of bias, whether hare removal affected the density of the source population. In agreement with the Ministry's guidelines (Trocchi and Riga 2005), translocations were never carried out when density was lower than 20 ind./km 2 , occurred routinely in two ZRCs (nos. 6 and 14 in Fig. 1), on alternate years in one (no. 3) and only occasionally in a further three (nos. 1, 8, and18); in all cases, the number of captured hares never exceeded 20% of the post-breeding population. As translocations precede a period of high mortality (winter), these low rates of harvest should have no effects on population dynamics (Boyce et al. 1999;Lebreton 2005).
Consistently with previous studies (Tapper and Barnes 1986;Meriggi and Verri 1990;Vaughan et al. 2003;Smith et al. 2004Smith et al. , 2005, crop diversity enhanced hare numbers, as hares require a mosaic of vegetation structures for both sheltering and facing seasonal or man-driven variation in food availability. The availability of some ploughed fields may benefit hare density by providing open areas which allow oversight of the surroundings and the prompt sighting of mammal predators (Hewson 1977;Slamečka 1991;Pépin Fig. 2 Variable importance, as assessed by both the percent increase of mean square error and increase in node purity, in the Random Forest Regression model; significant variables (permutation test) are shown by red bars (Prot_bound: percent length of protected boundaries) Fig. 3 Partial dependence plots for the four most important variables in the Random Forest Regression model (see Fig. 2). The Y-axis (D, density) value is determined by the average of all the possible model predictions with the dataset when the value of the objective predictor is X. Small ticks on the X-axis indicate deciles of the variables and Angibault 2007). In agricultural areas, hares, particularly males, select ploughed fields as daily resting sites (Meriggi et al. 2015), possibly because they are able to camouflage themselves into the furrows.
Accordingly, in the study area foraging hares selected ryegrass, together with hayfields and lucerne fields, which provide good quality food, positively affecting female body condition and reproductive success (Wincentz 2009). Consistently with previous studies, hares avoided both ploughed fields (Meriggi and Alieri 1989;Vidus-Rosin et al. 2009) and maize (Pavliska et al. 2018), which does not provide good forage (Katona et al. 2010).
As expected, hares were never recorded next to rural buildings, where disturbance caused by men and domestic animals (dogs and cats) can be expected to be the highest, and in tree stands, supporting their preference for open areas.
While the role played by land use patterns in shaping the densities of hare populations has been widely studied, the effects of surveillance against poaching, i.e. the illegal taking of hares in protected areas, has received little attention, probably because protected areas for hare reproduction are a prerogative of the Italian game management system. We pointed out that lack of protection can largely affect hare density and thus the success of ZRCs. Poaching occurs mainly during the hunting season, and mortality is more likely to be additive rather than compensatory for ZRCs hosting low-density populations (Bartmann et al. 1992). It is widely recognized that local farmers and hunters can play a major role in the conservation and improvement of habitat quality for game species, with positive effects for biodiversity (Oldfield et al. 2003;Reid et al. 2010). Although differences were not statistically significant, our results suggest that the involvement of stakeholders may ensure the most effective protection against the incursions of outsiders, who are perceived to threaten the interests of the local community (Bell and Morse 2007).
As expected, older ZRCs tended to host larger hare populations, suggesting that two decades of effective protection may be necessary to set up viable populations. The positive effect of barriers at the boundaries of the protected areas may depend on the reduction in hunting pressure at the ZRChunting ground interface during the hunting season; alternatively, barriers, at least roads with high traffic densities and wide canals (Underhill and Angold 2000), may limit hare dispersion and enhance philopatry.

Conclusions
Habitat restoration aims to support declining species as well as societal demands for ecosystem goods and services (Hanberry et al. 2015). In this regard, being a traditional game species throughout Europe, the brown hare can play a pivotal role. Fig. 4 Habitat selection by brown hares during their foraging activity in the pre-breeding season. Observed and expected frequencies (F) were compared by the chi-squared test with Bonferroni's confidence intervals for the proportion of use (P < 0.001) Through the analysis of a large dataset, we confirmed that habitat heterogeneity is a major factor affecting hare density in intensively cultivated land. Our results support those obtained by landscape simulation models, which suggest that any management measure increasing farmland diversity has the potential to increase hare abundance (Langhammer et al. 2017).
In contrast, in the central Po-plain, in the last decade, the cover of widespread maize and rice fields has increased, to the detriment of barley and winter wheat fields, together with the use of herbicides and insecticides (data from the archives of the national institute for statistics and Lombardy region; Balestrieri et al. 2019). Increasing cultivation of intensive crops, and associated loss of rotational set-aside, are expected to have a negative impact on European hares (Gevers et al. 2011).
In the national authorities' intention, the dispersion of wildborn hares from ZRCs should ensure the natural restocking of surrounding hunting areas, averting the release of non-native strains or captive-reared individuals. To achieve this goal, it is paramount to ensure greater public awareness and the involvement of local stakeholders in the management of both protected areas and hunting territories.
Acknowledgements A. Meriggi and an anonymous referee improved the first draft of the manuscript with their comments and suggestions. S. J. Lovejoy kindly provided language editing.
Funding information Open access funding provided by Università degli Studi di Pavia within the CRUI-CARE Agreement.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.