Introduction

Increasing deer abundance across Europe and North America has important consequences for humans and ecosystems. Deer contribute to rising incidence of tick-borne encephalitis (Jaenson et al. 2018) and Lyme borreliosis (Gilbert et al. 2012; Kilpatrick et al. 2014). Deer populations influence forest structure, composition, biodiversity, and nutrient cycling (Trdan and Vidrih 2008; Ramirez et al. 2018) and at high densities, can limit forest regeneration (Gill 1992; Côté et al. 2004), with implications for climate change mitigation strategies through land use change, woodland expansion, and increased carbon sequestration (Bastin et al. 2019; UK Climate Change Committee 2020). In the absence of top predator regulation, deer management primarily utilises shooting to control deer abundance and impacts (Department for Environment 2010); however, for ethical wildlife control, it is important to assess whether altering human practices could mitigate problems reducing reliance on culling (Dubois et al. 2017). Immunocontraception remains problematic, while fencing is expensive to install, costly to maintain, can cause access problems, and in some landscapes risks collision mortality of forest grouse (e.g. black grouse Tetrao tetrix, and capercaille Tetrao urugallus, Moss 2001; Baines & Andrew 2003). The local distribution and aggregation of deer and the severity of their impacts can be influenced by habitat and landscape attributes (Hurley et al. 2012; Jarnemo et al. 2014; Royo et al. 2017; Spake et al. 2020) and forest management (Reimoser and Gossow 1996; Vospernik and Reimoser 2008; Reimoser et al. 2009; Redick and Jacobs 2020). If landscape context effects are sufficient, it may be possible to influence local deer activity-abundance away from patches where vulnerable, browsing-sensitive, tree crops are grown, thus reducing the need for lethal control (provided that any impact other than that on vegetation is absent), or conversely, it may be possible to create high deer density areas with few public recreationists, for more effective culling in terms of the effort required.

Although local habitat preferences are well understood for many deer species (Acevedo et al. 2010), landscape context effects are less well studied (though see Saïd & Servanty 2005; Lamberti et al. 2006; Coulon et al. 2008), with the potential influence of recreational activity on deer distribution (Stankowich 2008) often overlooked. Countryside recreation is encouraged for health and human wellbeing, and levels are increasing (Johna et al. 2010); in the UK, through the Countryside and Rights of Way Act 2000 (CROW Act) large areas are now designated as open access, including much of the national forest (Edward and Kazushige 2015) where recreation is balanced with deer control to allow forest regeneration. Sustained high levels of recreational access may constrain deer management through safety concerns over shooting, but if deer are sufficiently disturbance-averse then changing the location of car-parks, access points, and pathways (Mallord et al. 2006) could potentially facilitate deer management. Building on the ‘landscape of fear’ hypothesis (Laundré et al. 2001), there is interest in using perceived predation risk both from stalking on foot and recreationists with dogs, to shift deer distribution away from areas requiring protection (Miller et al. 2001; Cromsigt et al. 2013). Human disturbance can affect deer habitat use and ranging behaviour (Coulon et al. 2008; Bonnot et al. 2013), and deer increase vigilance (Ciuti et al. 2012) and flush away from recreationists (Miller et al. 2001; Taylor and Knight 2003). It is unclear whether individuals in hunted populations consistently perceive recreationists as a predatory threat (Altmann 1958; Behrend and Lubeck 1968; Frid and Dill 2002), as deer may learn to avoid the characteristic behaviour of stalkers but habituate to and tolerate recreational walkers or tourists (Laundre et al. 2010).

In an extensive study across a large heterogeneous landscape, we quantified the local abundance at scales of 49.2 ha and 37.2 ha (SD = 11.5, SD = 9.6 for roe deer and muntjac respectively, representing aggregation or settlement, not sub-populations) of re-introduced roe deer (Capreolus capreolus, Ward 2005) and introduced invasive non-native muntjac (Muntjacus reevesi; DEFRA 2018), by extensive thermal imaging distance sampling with robust sample sizes. Local abundance was related to fine-grained habitat suitability using spatially explicit modelling that accounts for spatial intercorrelation (Miller et al. 2013). We tested habitat and landscape context features known or likely to influence local abundances of roe deer and or muntjac (Sinclair 1997; Putman and Moore 1998), and whether local deer distribution was affected by levels of recreational activity, using a spatially explicit extensively calibrated recreational path-network model. Roe deer model controlled for muntjac density and vice versa muntjac models controlled for roe deer density. We a priori considered that roe would respond more strongly to recreation, as they are considered sensitive to disturbance (Hewison et al. 2001; Coulon et al. 2008; Bonnot et al. 2013), while muntjac appear tolerant of human contact at least in their introduced UK range. Last, we examined the effectiveness of landscape, habitat, and recreation as tools for vegetation impact mitigation, contrasting the spatial variation in local deer density accounting for the fine-grained variables against a null model ignoring these variables.

Methods

Study area

The study was conducted in a 187-km2 conifer-dominated, lowland plantation forest landscape in eastern England, Thetford Forest, characterised by a semi-continental climate and sandy soils varying from acidic podsols to calcareous redzinas. Four deer species are resident: roe deer, muntjac, red deer (Cervus elaphus), and fallow deer (Dama dama) (Zini et al. 2019); these are managed by a professional ranger team to mitigate browsing impacts, particularly on vulnerable tree crop species. The forest landscape spans 12 ‘blocks’ that vary in predominant soil type and configuration — from contiguous blocks in ‘core’ forest-dominated landscape to outlying fragments providing easier access to surrounding arable fields and grasslands (Supplementary material, Appendix A). Configuration and soil type are independent (Zini et al. 2019) providing landscape-scale replication. Planted in the 1930–1950s, plantations are managed to maturity, clear-felled, and replanted (Eycott et al. 2006). With 63% of the forest in the second crop rotation, the felling-replanting cycle has been slowed to sustain the local forestry industry and important open-habitat species (Lin et al. 2007; Wright et al. 2009; Sharps et al. 2015), creating a fine-grained mosaic within each block that comprises even-aged stands (mean area = 4.5 ha, SD = 6.5; Zini et al. 2019) permeated by an extensive trackway network. In addition to timber, carbon, and biodiversity, recreational amenity is an important social benefit of multi-criteria forestry (Oh et al. 2017) and most of the forest (75%) is designated open access being visited by approximately 1 million people annually (1064 visitors km−2 year−1; Armour-Chelu, Brooke & Riley 2014).

Survey method

We used extensive nocturnal thermal imaging distance sampling conducted in late winter (January to March 2018) when low temperatures enhance thermal contrast and detection is improved by leaf-drop of deciduous under-canopy and collapse of herbaceous vegetation and bracken Pteridium aquilinum (Wäber and Dolman 2015). One-sided (from passenger-window) vehicle-based distance-transect surveys were conducted on the trackway network, using a thermal imager (Pulsar Helion XP50), while driving at ≤ 16 km/h with headlights off, between 20.00 and 02.00 h. Each transect was driven twice totalling 567 km; overall sampling intensity (3.02 km/km2) was similar in each of the 12 blocks (mean = 3.3 km/km2, SD = 0.8; see Supplementary material, Appendix A for transect routes). Each deer detected was briefly illuminated by a roof-mounted spotlight (‘lamped’) allowing its distance (distancesighting) to be measured using Leica Geovid 10 × 42 laser rangefinder binoculars. Lamped animals remained stationary or slowly walked away; this lack of urgent escape behaviour is consistent with a lack of night hunting (professional rangers do not operate at night, and poaching gangs were not detected during extensive multi-annual fieldwork). The angle of the observed animal relative to the surveyed transect line (anglesighting) was estimated to 5° using an angle board and perpendicular distance (from the transect line) was calculated as.

$${Distance}_{perpendicular}=\cos({angle}_{sighting})\;.\;{distance}_{sighting}$$
(1)

Perpendicular distance was reduced (near-truncated) by the mean verge width of the trackway class (forest roads = 8 m, SD = 6; fire routes = 3.8 m, SD = 2.1; forest rides = 0.60 m, SD = 0.2) following Wäber and Dolman (2015). Stands were categorised into three visibility classes (open, dense, mature) depending on their age and whether they had been thinned (see Supplementary material, Appendix B for details), following Zini et al. (2019). Data were recorded from 2414 single-sided transect ‘sections’ defined by the adjacent surveyed stand (mean length = 225 m, SD = 128), but these were not suitable for density surface modelling due to zero inflation and overdispersion (with potential for biased parameter estimates; Tu 2006). Consequently, sections of the same visibility class lying within 600 m of each other (along either side of the transect centre-line) were combined for analysis as 1162 composite ‘segments’ (mean length = 462 m, SD = 297), irrespective of whether they were contiguous or interrupted by sections of another visibility class (see Supplementary material, Appendix B for further details). Thus, segments may partially overlap, or be nested within, others that differed in detectability class but sampled similar neighbourhood characteristics; their lack of spatial independence was accounted for by spatially explicit analysis (see below).

Candidate variables

Local abundance of muntjac and roe deer were then related to the local intensity of recreational activity, roe deer density (only in the muntjac model), muntjac density (only in the roe deer model), and four environmental variables quantifying habitat and landscape context (Table 1), informed by the ecology of the two deer species modelled, and extracted from buffers around each transect-segment (hereafter ‘localities’; for details of buffer radii, see below).

Table 1 Candidate variables examined in density-surface models of local roe and muntjac deer abundance, showing for each: the mean, standard deviation (SD), and coefficient of variation (CV, SD: mean) across n = 1162 transect segments sampled for local deer abundance
Full size table

Both muntjac and roe are small- to medium-sized deer (10–18 kg and 17–23 kg respectively, Macdonald 2009) that, as concentrate selectors (Tixier et al. 1997), require more nutritious forage relative to grazers and bulk feeders. This, combined with a requirement for thermal shelter and security (Mysterud and Østbye 1995), affects the habitat selection preferences of both species (Hemami et al. 2005). In conifer-dominated plantations, density of roe deer is greater in young or mature forest stands, while density of muntjac is greatest in mature stands; both have low densities in closed-canopy (pole-stage) stands (Hemami et al. 2004) in which ground vegetation is suppressed by shade (Eycott et al. 2006). We therefore extracted the proportion of forest habitat in each locality comprising young (aged 0–5 years since planting; Table 1) and mature (> 45 years, following Hemami et al. 2004) stands from the Forestry England (FE) forest management database (mean polygon area = 2.9 ha, SD = 3.0) and, for surrounding areas, from the Centre for Ecology and Hydrology Land Cover Map 2015, hereafter LCM2015 (Rowland et al. 2017). Forest extracted from LCM2015 was classified as mature forest after ground-truthing. LCM2015 land cover is classified from composite summer–winter satellite images (pixel resolution 20–30 m), spatially referenced to digital cartography (Ordnance Survey MasterMap topographic layer) with a minimum polygon resolution > 0.5 ha.

Arable fields provide high-quality food relative to the forest (Putman and Moore 1998), such that fragmented forest in arable landscapes may enhance habitat suitability for deer, and population growth (Sinclair 1997). Adult female roe deer in forest stands closer to arable have greater body mass (Zini et al. 2019). Cover of arable lands was extracted from the 2018 Crop Map of England (CROME2018) maps of crop types from supervised classification of Sentinel-1 Radar and Sentinel-2 Optical Satellite images. As both deer species avoided entering outdoor pig units (confirmed through extensive thermal imaging survey), these were mapped manually using satellite data derived from Sentinel-2 (Copernicus Sentinel data 2018) and excluded; remaining CROME2018 crop types were pooled as ‘arable’.

Within forest habitats, calcareous soils support mineral-rich and palatable forage (important during gestation and lactation, see Brown 1990) and greater plant species richness (Eycott et al. 2006; Dolman et al. 2012), while acidic soils are often dominated by less-palatable Deschampsia flexuosa (Scurfield 1954), or bracken (Marrs and Watt 2006). High-resolution mapping of forest soils (Corbett 1973) comprising 18 soil series and Cranfield NATMAP data (coarser: 9 series) for surrounding areas were reconciled to a common classification, following Zini et al. (2019). For each locality, we calculated the percentage of the combined forest and grassland area overlaying calcareous soil, excluding arable lands that provide nutrition irrespective of underlying soils that are modified by marling and fertiliser application.

We extracted the proportion of grassland from each locality, combining permanent open habitat within the forest (from the FE management database, mostly grassland but including areas of bracken and scattered scrub and, for the surrounding landscape, pooling LCM2015 classes of improved grassland, rough grassland, neutral grassland, calcareous grassland, acidic grassland, and heather (as ground-truthing showed finer-scale ecological-resolution was unreliable, Zini et al. 2019).

The local intensity of recreational activity in each locality was extracted from an empirical GIS-based Network Analysis considering path element class (forestry road, fire route, ride), path network-distance from the lowest impedance access point (considering behavioural impedance accounting for path classes, as a weighted distance, following Oh & Jeong 2007), car parking capacity at this access point, and population density in concentric buffers around the access point, controlling for hour of day and day of week (Hornigold 2016). Calibration data comprised 1536 observations (single recreational disturbance events) collected during 1713 one-hour survey visits (controlling for time of day, day of week, and season), across 180 locations (intersection nodes in the pathway network) over five years (138 points surveyed in 2007, 180 points in 2008 and 2009, and 26 in 2013–2014, for further details of methodology and recreational modelling see Supplementary material, Appendix C). The distribution of access points and housing have been stable from this earlier survey to the collection of deer data, so relative distribution of recreational intensity is expected to be unchanged. The model related the observed frequency of dog walkers for each footpath element, with mean marginal and conditional R2 from tenfold cross validation of 0.13 and 0.41, respectively, was interpolated across the forest-wide trackway network to predict the mean frequency of disturbance events (DE) for each footpath element.

An index of the intensity of recreational activity across the footpath elements i within each segment-buffered locality j, was calculated as:

$$\sum\nolimits_{i=1}^n\frac{pde\;i\ast length\;i}{areaj}$$
(2)

where pdei is the predicted mean frequency of DE (per hour) on element i, lengthi is the length (m) of footpath element i intersected by the locality, and area j is the area (m2) of the locality.

Species-specific deer densities were resampled from 100-m resolution density surface model (DSM) rasters, generated from distance sampling data that incorporated geographical coordinates only (see below for more details on DSM models).

Statistical analysis

We obtained a large sample size with 1590 muntjac groups and 400 roe groups observed. We built separate DSMs for muntjac and roe deer as a two stage process, following (Miller et al. 2013). First, a detection function, accounting for the decline in detection probability with increasing perpendicular distance from the transect line, was fitted to the observations using package “Distance” (Miller 2016), accounting for segment detectability class (categorical covariate, three levels), after right-truncating data excluding the 5% most distal observations, following Buckland et al. (2001). The best-fitting function (uniform; half-normal; or hazard-rate), adjustments (hermite; cosine; polynomial), detectability covariate, and whether group size was related to perpendicular distance (accounting for lower probability of detecting distal smaller groups) were selected by comparing model AIC (see Supplementary material, Appendix B). The selected model structure provided detection probabilities per segment, from which density per segment was calculated (accounting for effective strip width). Then, a generalised additive model (GAM) was fitted, relating density per segment to geographical coordinates of latitude and longitude using a penalised thin plate regression spline accounting for spatial autocorrelation giving a smoothed surface best-fitting the data, simultaneously considering the set of a priori candidate variables at each locality as linear predictors with Gaussian error, using package “dsm” (Marshall et al. 2017).

Variables were extracted from a buffer radius around each locality, informed both by home-range size reported in the literature and statistical analysis of our dataset. For roe deer, winter home-ranges measuring 75 ha were reported from radio-telemetry in our study area (Chapman et al. 1993), corresponding to a circular home-range of 500-m radius. A 500-m radius represented a coarser grained scale than that of individual stands (mean width = 355 m, SD = 60); thus, we inspected a series of full DSMs (including all candidate variables, geographical coordinates as smoothing parameters, and the random effect of forest block, n = 12) across a range of buffer radii (100–650 m at 50 m increment step, where 650 m was the maximum roe deer home-range reported see Appendix D). This showed a clear AIC minima at 300 m that was selected for subsequent modelling. For muntjac, the reported winter home-range from radio-telemetry of 18 ha (Chapman et al. 1993) corresponds to a circular home-range of 250-m radius. Inspecting AIC values across DSMs built with incremental buffers ranging from 100 to 700 m (at 50 m increments, maximum roe deer home range reported, see Appendix D) showed an inflexion at 250 m that was selected for subsequent modelling (for more detail see Supplementary material, Appendix D). Analyses were conducted in R software (R core Team 2018). Variables were non-collinear across all candidate buffer widths (r ≤ 0.7; Freckleton 2002).

Whether candidate variables contribute to explaining local deer density, was assessed using multi-model inference in an information theoretic framework (Burnham and Anderson 2002), applied across the full set of alternative models comprising all possible combinations of candidate variables at the selected buffer radius (geographical coordinates were included in all models as smoothing parameters), using ‘MuMIn’ package (Bartoń 2018). A variable was considered supported if included in the 95% model confidence set (the sub-set of candidate models for which cumulative Akaike weights amount to 95% of the total) with a model-averaged CI that did not span zero (Boughey et al. 2011). We report model-averaged coefficients and quantify effect sizes as the magnitude of change in local deer density for an increase in the candidate variable from its lower, to its upper, distribution quartile (hereafter ‘inter-quartile range’, IQ). As a measure of explanatory power, we report the adjusted R2 of a model constructed comprising geographical coordinates and supported variables, as well R2 of the model comprising geographical coordinates only — where thin plate regression splines will model spatial variation in deer density encompassing both the unmodelled effect of candidate variables as well as other variance due to unexplained effects. We expect that a model specifying the candidate variables directly would have a greater R2 than a model using geographical coordinates only if these variables greatly affect deer density. Lastly, we predicted muntjac and roe deer densities using the supported variables and geographical coordinates in a raster of 100 × 100 m and compared it to the prediction of the simpler model, comprising geographical coordinates only. To compare the change in the spatial prediction of deer numbers generated by the supported variables and the model including geographical coordinates only, we calculated the correlation coefficient between the two raster predictions.

Results

Roe deer

For roe deer, a hazard-rate detection function incorporating segment-specific detectability class best explained the decline in detectability with distance. A model that incorporated a decline in observed group size (mean = 1.7, SD = 0.8) at greater detection distance received similar support (ΔAIC = 0.14, Supplementary material, Appendix B) and was therefore excluded by parsimony, so that no group size adjustment was made in subsequent density estimation. Effective strip (half) width averaged across the three detectability classes was 90 m (SE = 4 m). Of the 1162 segments, 894 (77%) had zero roe deer detected, while 268 (23%) had one or more deer detections. Roe density across the forest (estimated for 100 × 100 m raster cells) ranged from a 0.9 to 55.4 deer/km2 (Fig. 1). Averaging DSMs across combinations of candidate variables showed that roe deer density was lower in localities (segment-buffers of mean length = 475 m, SD = 298, radius = 300 m, mean area = 49.2 ha, SD = 11.5) containing a greater intensity of recreational activity (Figs. 2–3A), or a greater cover of grassland (Figs. 2–3B) and was higher in localities with a greater cover of young stands (Figs. 2–3C). Roe deer density was not affected by the extent of mature stands, arable lands, calcareous soil or muntjac density (Fig. 2). Adjusted R2 of the DSM including supported candidate variables was 0.17, compared to 0.13 for the model including geographical coordinates only (which nevertheless captures some effect of local variables through the geographic model). The roe deer density raster predicted by the DSM comprising geographical coordinates and supported variables was positively related to that predicted using geographical coordinates only (R2 = 0.74).

Fig. 1
figure 1

Muntjac and roe deer density (individuals /km2) predicted for a 100-m × 100-m raster, by models including supported candidate variables and geographical coordinates or by models including geographical coordinates only. All models included random effects of forest block (n = 12)

Full size image
Fig. 2
figure 2

Models of local roe deer and muntjac density showing effects of candidate habitat and landscape context variables. For each variable, model-averaged coefficients following multi-model inference (red filled symbols non-supported variables, black filled symbols supported variables) with 95% Confidence Intervals (bars) are shown, calculated across the 95% model confidence set. All models included the random effect of forest block (n = 12)

Full size image
Fig. 3
figure 3

Relation of local roe deer and muntjac density to supported habitat and landscape variables. Densities were predicted from density surface models including geographical coordinates and those variables supported by multi-model inference. Shown are local density of roe (a–c) and muntjac (d–g) deer in relation to: recreational intensity (number of disturbance events per hour; a, d): sum of disturbance events per hour across trackways per buffer weighted by buffer area; percentage of cover of grasslands (b, e); young stands (c); mature stands (f); and calcareous soil (g) within localities. Deer densities were predicted across the full range of each predictor variable

Full size image

Muntjac deer

For muntjac, a hazard-rate detection function, incorporating segment-specific detectability class, best explained the decline in detectability with distance. Greater group size (mean = 1.16, SD = 0.4) increased the probability of detection at greater distance (ΔAIC =  + 10.75 on removal of group-size term, Supplementary material, Appendix B); an adjustment for group size was therefore made in subsequent density estimation. Effective strip (half) width averaged across detectability classes was 75 m (SE = 3 m). Of the 1162 segments, 481 (41%) had zero deer detected, while 681 (59%) had one or more deer detections. Muntjac density (at 100 × 100 m) ranged from 4.7 to 77.5 deer/km2 (Fig. 1). Averaging DSMs across combinations of candidate variables showed that muntjac density was greater in localities (segment-buffers of mean length = 475 m, SD = 298, radius = 300 m, mean area = 37.2 ha, SD = 9.6) containing greater intensity of recreational activity (Figs. 2–3D) and with greater cover of mature stands (Figs. 2–3F) or of calcareous soil (Figs. 2–3G), but was lower in localities with a greater cover of grassland (Figs. 2–3E). Muntjac density was not affected by the extent of young stands, arable lands, or muntjac density (Fig. 2). Adjusted R2 of the DSM including supported candidate variables was 0.33, compared to 0.30 for the model including geographical coordinates. The muntjac density raster predicted by the DSM comprising geographical coordinates and the supported variables was positively related to that predicted using the DSM comprising geographical coordinates only (R2 = 0.81).

Discussion

In this extensive landscape-scale study, with intensive replication of local habitat across multiple localities, roe deer had lower abundance in areas with highest levels of human recreational activity, while the reverse was true for muntjac. Although ecologically interesting, the effects on local deer density of recreation and other habitat variables (see below) were relatively small and thus unsuitable to be used as tools for deer impact mitigation.

Local habitat and density effects

The two species differed in their responses to forest habitat. Roe were more abundant in localities with greater percentage of young stands (inter-quartile, IQ, effect size, + 1.3 deer/km2), while muntjac, as previously reported (Hemami et al. 2004, 2005), were more abundant in localities with greater percentage of mature stands (IQ effect size, + 2.4 deer/km2). Younger stands have previously been shown to be important for roe deer populations, with increased recruitment after extensive (75%) conifer re-planting of a woodland in Dorset, and subsequent population declines as the forest matured (Gill et al. 1996). Earlier research in our study landscape similarly found roe deer density was higher in young stands, and muntjac density was higher in mature stands (Hemami et al. 2004). In a forest managed by clear felling and replanting, young stands offer the greatest amount of ground forage (Smolko et al. 2018), but mature stands (with open canopy) offer both ground vegetation and cover (Hemami et al. 2005). Bramble (Rubus fruticosa agg.) is widely exploited by deer (Holt et al. 2010; Wäber 2010) including roe and muntjac (Chapman et al. 1985) that both occur in higher densities where more bramble is available (Hemami et al. 2004, 2005); bramble is restricted to early and late forest growth stages by shading (Battles et al. 2001; Hemami et al. 2004). Differences in habitat use between these two species may be further explained by physiology and predator escape mode: while roe deer use open habitats to forage (Putman 1986) muntjac have shorter legs and more curved backs, adaptations that do not favour sustained running as an escape response but facilitate sheltering in dense low cover. The association of the smaller sub-tropical muntjac with dense shrub may be further explained by a greater requirement for thermal cover, relative to the temperate-boreal distribution of roe deer (Dolman and Wäber 2008).

Muntjac abundance was higher in localities with more calcareous soil (IQ effect size, 7.1 deer/km2) that have more diverse forage (Eycott et al., 2006). Muntjac are smaller-bodied than roe, potentially placing them at the extreme of the browser-grazer continuum (Bodmer 1990) compared to roe, further increasing the need to select nutritious food (sensu Gordon & Illius 1994). Roe deer and muntjac abundance were both lower in localities with a greater percentage of grassland (IQ effect size, − 2.9 and − 2.6 deer/km2, respectively). This may be a direct consequence of habitat suitability, as both species have a concentrate-selector diet in contrast to larger grazing deer such as fallow and red deer (Gordon and Illius 1994; Obidziński et al. 2013). It is also possible that grassland extent may be a proxy for higher local abundance of these larger and potentially competitively dominant (Dolman and Wäber 2008) species that congregate around and graze within grassland areas (Thirgood 1995; Trdan and Vidrih 2008). Despite the two species using similar resources and the likelihood of these species to compete for resources (see Hemami et al. 2005), muntjac density was not related to roe deer density and roe deer density was not related to muntjac density.

Recreational intensity

Roe deer local abundance was negatively affected by the local frequency of dog walkers (IQ effect size, − 1.1 deer/km2). According to the ‘human-caused predation risk’ hypothesis, animals perceive disturbance from recreational activities as a predation risk (Frid and Dill 2002). This can ultimately result in avoidance of disturbed areas and increased energetic costs of antipredator behaviour reducing effective habitat quality (Gill et al. 2001). Animals may trade-off the degree of perceived risk against the resource value of a patch (Charnov 1976). For example, a study of flight responses in southwestern France showed roe deer were more vigilant when foraging in more open landscapes or when close to human settlements (Bonnot et al. 2017); in the same study area, GPS data from collared roe deer show this species uses open habitats (where forage is abundant) more during the night while using forests more during the day, likely to hide from humans; use of open habitats during the day also decreases with proximity to roads and human settlements (Bonnot et al. 2013). Similarly, avoidance of buildings and roads have been shown to be the most important factors in roe deer habitat use (Coulon et al. 2008).

Across three study sites in northern and one in southwestern France, roe deer occurred more frequently in areas with less human infrastructure (Hewison et al. 2001). Red deer and fallow deer have also been shown to be affected by recreation (Boer et al. 2004; Sibbald et al. 2011; Westekemper et al. 2018) and in Germany, red deer used zones with lower recreational activity (Coppes et al. 2017); however, differences in recreational access were applied across a bigger scale than the one examined here.

Local abundance of muntjac was positively related to the local frequency of dog walkers (IQ effect size, + 1.5 deer/km2); their smaller size makes it easier for them to hide in dense shrubs and avoid recreationists. Rangers in Thetford Forest avoid culling deer in areas popular with recreationists over safety concern; animals more tolerant of recreation can therefore reduce the risk of being culled by establishing their home-range in these parts of the forest. While comparison of flight distances in different environments is not an accurate method for inferring behavioural differences it is notable that in a high-throughput ‘safari-park’, in southern England, muntjac flight distances were only 8 m in non-public areas, and 11 m in public areas (Champion 1997), significantly less than those reported for roe deer in the Netherlands (64 m when walking down wind, 41 m upwind, 44 m in calm wind (Boer et al. 2004). This short flight distance, combined with the difficulty of culling deer in areas of high recreational pressure, suggests the trade-off of occupying areas with high number of recreationists is positive for muntjac. Numerous records report muntjac in urban gardens (Harris et al. 1994), suggesting muntjac to be tolerant of proximity to humans. When comparing results of surveys of muntjac and roe deer numbers in urban environments and rural areas muntjac seems to outnumber roe deer in urban environments. Data from the BTO UK-wide Breeding Bird Survey (covering 3295 and 3941 1 km2 grid cells in 2006 and 2017 respectively) provide observations of some mammal species including both deer, this shows that both roe and muntjac increased in rural areas by 81% and 117% respectively between 2006 and 2017 (Harris et al. 2018).

Models incorporating geographical coordinates and the set of supported candidate variables performed better than simpler models including solely geographical coordinates; however, the increase in variance explained was small, and the majority of the variance is to be attributed to unexplained spatial variation. Densities of fallow and red deer were not accounted for in our models and are likely to account for some of the unexaplained variation in local deer density. It is possible that the influence of variables reported here, relating to local deer distribution (aggregative response, Sutherland 1983), may not affect underlying demographic processes (driven by fertility or neo-natal survival) and therefore overall numbers across the landscape. Regarding forest management, this study showed that felling and restocking tree crops may increase local roe deer abundance (as they tend to congregate more in young stands) while decreasing local muntjac abundance. As recreation appears to negatively affect roe deer local distribution, it may be possible to manage recreation (changing the spatial distribution of access points and paths accessibility) and change local density of roe deer to create ‘safe refuges’ favouring greater local abundance. Culling in these areas would then be simultaneously safer and more cost-effective (in terms of the effort required to achieve a given harvest, Hatter 2001). Managing recreation however would only slightly improve forest management if roe was the only deer species present, or the species having the greatest impact. This solution is not suitable for managing muntjac, for which the opposite response to recreational intensity was found. Managing recreation at a fine-scale does not appear to provide a useful tool for muntjac and roe deer management.

Conclusions

Our results show that forest management affects local deer density, but the effects found here were not strong enough to justify adapting forest management to mitigate impacts and would not allow the removal of lethal control by shooting. The variance explained by the models including the variables supported as well as the effect sizes of each variable were relatively low, showing little effect on local deer aggregation. Our conclusion is that culling deer continues to offer the most effective tool for local density and impact reduction. An ethical and effective culling program should start from an assessment of deer numbers to build a demographic model; culling should then be coordinated and take place at a relevant spatial scale for each culled species (Fattorini et al. 2020) and evaluation of impacts and population size will also be necessary to monitor its effectiveness (Dubois et al. 2017).