Relation of pine crop damage to species-specific density in a multi-ungulate assemblage

Deer management in forest ecosystems requires information on deer densities and impacts to inform culling decisions with a known target density for acceptable damage levels. In multi-ungulate assemblages, managers need knowledge of relative impacts by different species or guilds. In an extensive (195 km2) conifer forest in eastern England, we related Scots pine (Pinus sylvestris) leader damage (% dominant shoots browsed) in 48 restocked stands (1–3 years growth) over multiple years (n = 79 observations) to species-specific annual muntjac (Muntiacus reevesi), roe (Capreolus capreolus), fallow (Dama dama) and “large deer” densities (composite of annual fallow and multi-year mean red deer Cervus elaphus) using generalized mixed effects models. Forest-wide density surface models were calibrated through intensive annual thermal imaging distance transects and local densities around stands resampled within confidence bounds. Models also examined effects of ground vegetation and hare presence (Lepus europaeus). More pine leaders were browsed at higher fallow or large deer densities (22% and 18%, respectively, increased leader damage across inter-quartile range). Leader damage intensity was not influenced by ground vegetation, hares, muntjac (across the range 8.3–41.6 individuals km−2) or roe deer density (1.7–19.4 individuals km−2). To reduce pine crop damage to economically acceptable levels, managers need to reduce fallow deer to a density as low as 0.6 individuals km−2 (CI = 0.06–1.44, which is considered impractical) or reduce combined large deer density to 2.3 (CI = 1.18–3.46) individuals km−2. Reducing muntjac or roe abundance would minimally affect leader damage in this system, but may be important for other tree species. Multi-species deer management requires species-specific understanding of impacts and robust density estimates.


Introduction
Increasing forest extent is key to meeting carbon targets and limiting global climate change (Griscom et al. 2017;Bastin et al. 2019; UK Committee on Climate Change 2019; Lewis et al. 2019). However, across much of temperate Europe, North America, Japan and New Zealand, high deer abundance threatens successful tree establishment (Takatsuki 2009;Ramirez et al. 2018;Lesser et al. 2019;Latham et al. 2020) as well as profoundly altering forest structure and ecosystem function (Gill 1992a;Côté et al. 2004;Carline et al. Communicated by Christian Ammer. 1 3 2005). Controlling such impacts in the absence of natural predators requires deer to be managed (Crete 1999;Hothorn and Müller 2010); ensuring such management is effective and justifying it to a wider public requires an appropriate evidence base (Reimoser 2003). Deer browsing is affected by a range of factors (Spake et al. 2020), of which local deer density is usually the most important. Deer densities can be driven by management regimes, habitat and human disturbance (see examples in Hemami et al. 2005;Bonnot et al. 2013;Coppes et al. 2017); however, manipulating local densities through habitat management or recreational access would only subtly affect local deer numbers (Zini et al. 2021). Identifying species-specific deer density thresholds below which forest crop damage is either acceptable, or where marginal costs of culling at such low numbers outweigh benefits of further reducing damage, supports deer managers by defining a sustainable density (Welch et al. 1992(Welch et al. , 1995 allowing appropriate cull targets to be set. A complicating issue is that deer species have been introduced beyond their natural range (Dolman and Wäber 2008;Carden et al. 2011;Champagnon et al. 2012). Evidencebased management of natural and modified multi-species ungulate assemblages requires understanding of which species cause the greatest impacts and of the species-specific relation of damage to density.
Although exclosure experiments, which contrast controls exposed to deer and experimental plots lacking deer, can usefully demonstrate the dramatic nature of vegetation impacts (see, for example, Simončič et al. 2019), they cannot identify what target density of deer is sustainable or desirable. To achieve this, local impact intensity should be examined across a range of known local deer densities, ideally for free-ranging deer at landscape scale not in small or poorly replicated exclosures. Heinze et al. (2011) related browsing damage on tree saplings to local densities of free-ranging fallow (Dama dama) and roe deer (Capreolus capreolus), but with estimates of plot-specific density estimates having wide uncertainty through faecal pellet counts (Gill et al. 1997;Campbell et al. 2004;Hemami et al. 2007). In England, combined impacts of muntjac (Muntiacuus reevesi), roe, fallow, red (Cervus elaphus) and sika (C. nippon) deer on tree regeneration were examined across multiple lowland coniferous and deciduous woodland blocks, but without separating species-specific contributions . There remains a need for studies that evaluate species-specific forest crop impacts across a range of deer densities using reliable methods of density estimation.
We examined impacts of three deer species: muntjac, roe and fallow deer, in 48 forest stands re-sampled in multiple years, across an extensive forest landscape in eastern England, to predict species-specific densities for successful regeneration of planted Scots pine Pinus sylvestris. Annual thermal imaging distance sampling, which reliably quantifies deer densities (Gill et al. 1997;Morelle et al. 2012;Wäber and Dolman 2015), allowed us to relate local impacts to species-specific annual densities of each species extracted from fine-grained density surface models (though for wideranging, low-density red deer, an aggregate multi-year density was used-see details in Methods). We tested an explicit a priori expectation that damage intensity would be more strongly predicted by the density of fallow, or the combined density of large deer (annual fallow and multi-year mean red density), rather than either of the two smaller concentrate selectors, roe and muntjac. Crop impacts may be modified by the local availability and relative palatability of other forage and ground-layer vegetation (Welch et al. 1991;Gill 1992a;Augustine and Jordan 1998;Palmer and Truscott 2003;Harmer et al. 2010); therefore, we also tested whether damage intensity was reduced with greater availability of preferred browse and forage within stands, as well as the presence of European (brown) hares (Lepus europaeus; hereafter "hare").

Study system
The study was conducted in Thetford Forest, an extensive (195 km 2 ) lowland conifer-dominated forest landscape in eastern England ( Fig. 1) managed by rotational clearfelling and replanting of even-aged stands (mean stand area = 8.1 ha, SD = 9.1). Four deer species occur (Wäber et al. 2013): muntjac and roe are the most abundant (see Fig. 2 and Table 1), both are relatively sedentary (Chapman et al. 1993), and their density varies within and between 12 forest "blocks" (Hemami et al. 2005;Wäber et al. 2013), which vary from contiguous core to outlying areas with greater access to the surrounding arable crops and grassheath. Fallow deer occur at lower densities and are concentrated in the southern forest blocks, although having a larger home range than the two smaller species (Borkowski and Pudełko 2007) fallow deer tend to be hefted (Apollonio et al. 2003), further contributing to local variance in density. For each of these three species, it was possible to map local densities annually and examine their independent contribution to impact levels. In contrast, in Thetford forest red deer occur in small groups (mean group size across years = 7.1, SD = 5) at a low density across the landscape (mean = 1.4 individuals km −2 , SE = 0.26) and are wide ranging (Reinecke et al. 2014) resulting in low spatial variance (see Fig. 1S Supplementary Information, Appendix A). Furthermore, red deer numbers are slightly higher in the southern part of the forest, making the red deer density variable correlated to fallow density. For this reason, stand-level crop impacts could not be analysed in relation to local red deer density. However, 1 3 we combined the local density of red deer and fallow deer (with greater spatial variance) in an aggregate measure of "large deer" density, to assess whether this better predicted damage than the influence of fallow deer alone. Deer are considered to be crepuscular animals (Náhlik et al. 2009;Sándor et al. 2011;Krop-Benesch et al. 2013; Ensing et al.  Mori et al. 2020;Zanni et al. 2021), and our surveys, carried out after dusk, measured the spatial distribution of the deer when active, thereby reflecting the distribution of their foraging pressure. The forest is set in a wider mosaic of arable, small woodlands, heathland and grassland (some extensive) also occupied by resident deer. Fallow, roe and muntjac density did not differ significantly between thermal imaging distance transect surveys of two neighbouring estates (measuring around 160 km 2 ; with grasslands, grass-heath, arable and scattered forests) and that measured in the forest (unpublished; data cannot be shared for private land), while the local density of roe and muntjac deer within the forest is not related to local access to external arable habitats (Zini et al. 2021). Fallow and red deer are relatively mobile and feed in open habitats at night; however, both utilize young growth and grasslands within the forest landscape with no systematic pattern of groups "commuting" between forest interior and the outside landscape. Importantly, any local influence on densities of dispersal (settlement and demography) or of individual daily foraging movements is included in the landscape density surface maps that are calibrated with spatially explicit observations of deer during their nocturnal active period in winter when crop impacts can occur. However, our measures of functional density (nocturnal numbers active in the forest) of larger deer could underestimate total numbers present diurnally in the forest, in which case crop damage might be mitigated before fallow density is pushed fully to zero.
The forest landscape is dominated by Corsican pine (Pinus nigra var. maritima, 59% of forest area), Scots pine (15%) and other conifers (5%); and limited deciduous tree crops (10%). Despite dominating second rotation planting during the last decades of the twentieth century, replanting of Corsican pine ceased from 2006 due to increased prevalence of red-banded needle blight (Dothistroma septosporum), most likely owing to milder and wetter recent winters . Subsequent replanting has used a variety of tree species, with sufficient restocked stands of Scots pine available for impacts to be related to speciesspecific deer density.

Crop impact data
We assessed leader (dominant shoot) damage as this commonly causes the greatest effect on forest crop harvest values (Motta 1996). Young conifers are particularly vulnerable while deer are still able to reach the leader (Welch et al. 1992), with greatest browsing damage in winter (Miller et al. 1982;Welch et al. 1992).
We collected impact data in 2018, 2019 and 2020, from restocked Scots pine stands aged 1-3 years, as Scots pine is most vulnerable to leader browsing within 3 years of planting (Holloway 1967). In 2018, 29 suitably aged stands were sampled, these were resurveyed in subsequent years unless they exceeded 3 years age, and additional newly planted stands were added (sampling 23 in 2019 of which 17 had been sampled in 2018, and 27 in 2020 of which 12 were sampled in 2019), providing 79 observations across 48 unique stands. Between March and April, we recorded winter browsing that removed leader buds from the previous growing season and recorded major fraying damage sufficient to cause leader death. In each replicate stand, a transect was walked on the longest diagonal, assessing the nearest tree every five steps, sampling at least 60 saplings per stand.
It was not possible to consistently distinguish deer browsing damage from possible contributions from hares or rabbits, Oryctolagus cuniculus; however, most damage is likely attributable to relatively abundant deer; nevertheless we controlled for the local presence of hares in analyses (see "Deer density and hare presence" section), but rabbits were excluded as they are unlikely to browse pine (see Supplementary Information, Appendix B). Other small mammals were considered unlikely to have contributed to leader damage as, although bank vole Myodes glareolus may consume coniferous buds (see review in Lyly et al. 2014), they did not consume Pinus sylvestris foliage in cafeteria experiments (Hjältén et al. 2004), the number of branches or height of Pinus sylvestris seedlings were not reduced by exposure to both bank and field vole Microtus agrestis in the field (Lyly et al. 2014) and small mammals most commonly damage bark rather than browsing pine foliage (Gill 1992b).
In order to measure alternative sources of food, we recorded the relative frequency of ground cover classes (noted at each transect point at the tip of the surveyor's shoe) classified as: bramble Rubus fruticosa agg., other shrubs (mainly silver birch, Betula pendula), grass, herbaceous vegetation (hereafter "herbs") or non-forage (e.g. bare ground, brash, or bracken Pteridium aquilinum litter), with the proportionate contribution of each class considered in analyses. Herbs and grasses may divert browsing impacts in early winter before dying back; in contrast, bramble is semievergreen, offers forage all year round and is an important predictor of roe and muntjac winter distribution (Hemami et al. 2004(Hemami et al. , 2005. As calcareous soils support greater plant species richness (Eycott et al. 2006) and more palatable grass species (Rodwell 1992), we calculated the percentage of calcareous soil in each restock stand (extracted from Corbett 1973) as a proxy for palatability of forage. Fig. 2 The distribution within the forest landscape of the intensity of Scots pine leader damage and densities of fallow deer, combined large (fallow and red) deer, roe and muntjac (predicted at 100 × 100 m grid from DSMs including latitude and longitude) shown separately for each of three study years ◂

Deer density and hare presence
Annual species-specific winter deer density surfaces were generated in 2018, 2019 and 2020 from distance sampling data, obtained in late winter immediately before sampling tree crops. Intensive nocturnal thermal imaging transects (see Supplementary Information, Appendix C) were conducted by thermal imaging (Pulsar Helion XP50) between January and March at the end of the cull season when vegetation concealment is minimal. Methodology followed Wäber and Dolman (2015), but with distances to detected deer measured using laser rangefinder binoculars (Leica Geovid 10 × 42). Each one-sided transect section was surveyed twice in both 2018 and 2019 (mean annual effort: 567 km, or 3.02 km km −2 across the forest), and once in 2020 owing to resource limitations (mean effort 1.5 km km −2 ). Across the 3 years, the mean number of observations (single or groups of deer) was 1193 (SD = 415) for muntjac, 312 (SD = 110) for roe, 73 (SD = 12) for fallow and 30 (SD = 18) for red deer.
Species-specific detection functions that account for the reduced probability of detection of individuals at greater distances (Buckland et al. 2001) were fitted pooling data from the 3 years (using the R package "Distance", Miller 2016), excluding 5% of most distal observations (following Buckland et al. 2001) and controlling for three structural visibility classes a priori defined by crop age (see Zini et al. 2021). The best-fitting key detection function was selected and adjustments for group size (to account for lower probability of detecting smaller groups at further distances) evaluated, by comparing model AIC (see Supplementary Information, Appendix D for details). We then fitted the selected detection functions separately for each survey year to obtain annual density estimates for muntjac, roe and fallow deer for each of 1162 surveyed composite transect "segments" (mean length = 462 m, SD = 297; for details on segment demarcation see Zini et al. 2021). These were used to calibrate landscape-wide DSMs, using a GAM (package "dsm" Marshall et al. 2017) that related abundance per segment to geographical coordinates, using penalized thin-plate regression splines to produce the best-fitting smoothed surface at a broad scale, greater than that of fine-scale stand structure. Thus the DSM implicitly accounts for a range of factors affecting local deer density, including local habitat, neighbouring deer management regimes, and ease of access to surrounding arable crops. The GAM was interpolated to predict densities for the entire study area as 100 × 100 m rasters using packages "raster" and "sp" (Pebesma and Bivand 2005;Bivand et al. 2013;Hijmans 2017), allowing resampling for the stands in which crop impacts were measured. Uncertainty around deer density estimates comes from both the detection function and the GAM. Assuming the detection function and GAM are independent, we calculated estimates of variance of predicted deer density per raster cell following Miller et al., (2013). In Supplementary Information, Appendix D (Table 4S) we report the yearly, species-specific coefficient of variation (CV) associated with the detection function, the CV associated with the GAM and the total CV (comprising both detection function and GAM CV). Red deer are a wide-ranging species (Reinecke et al. 2014), and the limited number of observations did not allow annual density surfaces to be modelled. Therefore, red deer observations from 2018 to 2020 were pooled and modelled with a single detection function (again incorporating the visibility class detectability covariate) to fit a single GAM, generating a multiyear average DSM (Supplementary Information, Appendix A). The proportion of deviance explained by DSMs was 0.13 for muntjac, 0.14 for roe, 0.40 for fallow (multi-annual means) and 0.70 for red deer (single pooled DSM).
For each sampled stand in each year, species-specific deer densities (and pooled "large deer" density) were resampled from the annual density surfaces using packages "rgeos" and "raster" (Hijmans 2017; Bivand and Rundel 2018) within a buffer around the centroid of each stand at a scale relevant to the species winter home range. For muntjac and roe deer, buffer radii of 250 m (18 ha) 500 m (75 ha), respectively, were based on the mean of sex-specific home ranges measured in Thetford Forest using radio-telemetry in November-December and January-February (Chapman et al. 1993). For fallow deer (and for large deer), a buffer radius of 850 m (219 ha) was based on home ranges of radio-collared male and female fallow deer in southern Poland (Borkowski and Pudełko 2007). Fallow and large deer densities were square-rooted to avoid leverage during modelling. During thermal imaging sampling, hare and rabbit were also recorded when seen but without distance measurement owing to time constraints; hare were considered in analysis on the basis of presence/absence. Hares are long-lived sedentary animals (Hulbert et al. 1996;Macdonald 2009) and their home ranges (11.8 ha in agricultural areas in Austria, 22.2 ha in Scotland, Hulbert et al. 1996;Schai-Braun and Hackländer 2014) are two to three times the area of a stand. For modelling, we considered presence of hares (as a dummy variable, zero = absent, 1 = present) if they were recorded (either during nocturnal distance sampling or during impact assessment fieldwork) from the focal stand sampled for tree impacts, or a directly adjacent stand, in any of the 3 years.

Analysis
We examined the percentage of sampled leaders recorded as damaged, hereafter "damage intensity" and incorporated stand identity to control for lack of spatial independence of stands resampled more than one year. We examined two model sets: one including fallow deer density and the second including large deer density, as correlation between these (r = 0.9) precluded combining them within the same models (Freckleton 2002). All remaining candidate explanatory variables were non-colinear (r ≤ 0.7, Freckleton 2002). Given the relatively large number of candidate variables (Table 1), we conducted analysis of each model set in three phases. In each phase, model-averaged variable coefficients and CIs were averaged across 1000 iterative models calibrated by resampling (independently for each stand-year) the 95% CI of species-specific deer density (derived from the maps of the coefficient of variation, see "Deer density and hare presence" section), to incorporate uncertainty in deer local deer estimates. We included a random effect of stand ID (n = 48) in all models to account for pseudoreplication of stands sampled in multiple years.
First, we related damage intensity (as % per stand) to all combinations of species-specific deer densities, presence of hares and tree crop age using generalized mixed effects models with random intercept and negative binomial error distribution, using MMI (Burnham and Anderson 2002) in package "MuMIn" (Bartoń 2018). Candidate variables were considered supported if included in the 95% confidence subset of models and their model-averaged coefficient did not span zero, following (Boughey et al. 2011); any lacking support in this initial MMI was excluded from subsequent analyses. Second, for each set (fallow or large deer), we built a series of further models each comprising the base herbivore model (incorporating species variables supported at phase one) and one of the vegetation modifiers. For each model set, we then subjected the provisionally supported variables (base herbivore model plus any supported vegetation modifier) to MMI. Last, for each model set, a final model comprising the supported herbivore and vegetation variables was used to predict the threshold deer density above which pine leader damage exceeded 20% intensity per annuum (that has been proposed as economically acceptable, Prien and Müller 2010), and the damage intensity expected at deer densities of 1, 5 or 10 individuals km −2 . The explanatory power of the two final model sets was compared in terms of their conditional (proportion of total variance explained through both fixed and random effects) and marginal (proportion of total variance explained through fixed effects) R 2 . All analyses were performed in R (R core Team 2018).

Deer density
Deer density in buffers around the crop impact stands was highly variable: ranging 8.2-42.0 individuals km −2 for muntjac, 1.8-19.3 km −2 for roe, 0-28.1 for fallow km −2 , 0.9-2.0 km −2 for red deer and 0.9-29.9 km −2 for large deer (see Table 1)-but with large deer overwhelmingly comprising fallow. Across the crop impact stands, muntjac were the most evenly distributed deer species, while spatial variance in density was greatest for fallow deer (Table 1). Relative spatial distribution of deer was broadly consistent between years, the mean of the correlation coefficient between density in stand buffers in 2018 versus 2019 and 2019 versus 2020 was 0.77 for muntjac, 0.92 for roe, 0.96 for fallow deer and 0.75 for large deer.

Models of crop damage
The mean leader damage per stand per annuum was 32.1% (SD = 25.9). The spatial distribution of leader damage intensity and species-specific density of fallow, large, roe deer and muntjac are shown in Fig. 2. MMI showed that the intensity of leader damage was greater in restocked stands with younger trees and with greater fallow deer density in the surrounding buffer, or with greater density of large deer (Figs. 3, 4). In contrast, no support was found for any effects of muntjac density, roe deer density or presence of hares in either model set (Fig. 3); these herbivores were therefore excluded from subsequent modelling.
No support was found for an effect of any ground vegetation variable, in the subsequent series of models that included the supported variables from either of the herbivore base models (fallow deer or large deer density; tree crop age) and each vegetation variable in turn (Fig. 3). Final models therefore comprised effects of tree crop age and the density of fallow deer plus or density of large deer; these provided similar explanatory power (for fallow deer conditional R 2 = 0.37, marginal R 2 = 0.65, for large deer 0.35 and 0.65, respectively; delta AIC for large deer relative to fallow = 1.8).

Discussion
Scots pine leader damage was primarily attributed to fallow deer or the combined effect of large deer (primarily comprising fallow), with interquartile effect sizes, + 22% and 18%, respectively, but we found no evidence of any effect of introduced muntjac or native roe deer. Although fallow deer are bulk feeders with their diet mostly comprising grass and graminoids, conifer leaves comprised up to 20% of winter rumen content in southern England (Jackson 1980) and pine foliage comprised 6% or 10% of winter diet in southern or northern Poland, respectively (Borkowski and Obidziński 2003). In addition to pine being a nontrivial component of fallow deer winter diet, their herding behaviour (Apollonio et al. 1998) means many individuals can Fig. 3 Relation of Scots pine leader damage to candidate explanatory variables of herbivore abundance and ground vegetation composition. Shown are model-averaged coefficients (filled symbol: supported variables in red, variables lacking support in black; bars represent 95% CI) from 1000 MMI iterations, each randomly resampling deer densities across the CI of density surfaces generated from intensive thermal imaging distance sampling. Model coefficients of herbivore species and crop age MMI are shown in (a and c). Variables retained in the herbivore base model (fallow deer or large deer density; tree crop age) were also supported in MMI analysis of models incorporating each ground vegetation variable in turn, therefore, for each of presentation only the ground vegetation coefficients are reported (b and d). All models controlled for the random effect of forest stand identity (n = 48). (Color figure online) simultaneously aggregate in a restocked stand, increasing the potential for substantial impact. Red deer have a similar herding behaviour and dietary strategy (Clutton-Brock et al. 1982) and have been shown to browse on Scots pine (Palmer and Truscott 2003;Pépin et al. 2006;Krojerová-Prokešová et al. 2010), but their overall density in Thetford Forest is low. However, models still predicted moderate residual damage even at very low density (e.g. 1 km −2 ) of fallow, or of combined large deer; this suggests additional impacts not being captured by our analysis, potentially due to other deer species or hares, or early-winter presence of fallow in areas where they were not detected by late-winter thermal imaging survey.
Pine leader damage intensity was not affected by muntjac or roe deer density in this landscape. For roe, this is consistent with evidence from both highland Scotland, where winter browsing incidence on Scots pine was unrelated to roe density (Palmer and Truscott 2003), and north-eastern Germany, where browsing on pine trees also was not related to roe density, across a similar range as observed in the current study (approximately 1-15 individuals km −2 ) (Heinze et al. 2011). A previous study in north-east England, however, showed a positive association between conifer leader browsing and an index of roe deer occupancy (Ward et al. 2008). Roe deer have high nutrient demands (Gordon and Illius 1994) and fibrous, tannin-rich conifer foliage does not represent optimal browse (Gonzalez-Hernandez and Silva-Pando 1999). Conifer needles in late winter or spring make up to 40% of roe deer diet in colder climates than those of our study area (de Jong et al. 1995), but other studies carried on in areas with comparable climatic condition show conifer only forms a negligible part of roe deer diet (Tixier et al. 1997;Duncan et al. 1998;Storms et al. 2008).
The lack of any relation of leader damage to muntjac density is notable, as elsewhere pine has been found to be browsed by this species. For muntjac in south-east England, rumen analysis showed Scots pine was mostly browsed in late winter when alternative forage was scarce (33% of individuals culled in March and April had Scots pine in their rumen, declining to 8% in May and June; Chapman and Tutchener 2017). Given muntjac's small size (for which we can assume a similar dietary strategy to the roe deer; Bodmer 1990) and the evidence showing Scots pine is only browsed when food is scarce, it is plausible that muntjac could also have fed on Scots pine but not enough to be detected by our analysis. The lower spatial variance in muntjac density may Fig. 4 Intensity of Scots pine leader damage related to the density of fallow deer, density of large deer and tree age (separately for each model set). Shown are predictions (continuous line) from a model including supported variables only (fallow or large deer density and crop age), incorporating the random effect of forest stand identity (n = 48). Shading represents (95% CI) have reduced our power to detect any additive impact if this was subtle relative to that of the larger fallow.
Leader damage was also unrelated to presence of hares though these were only sparsely distributed weakening our ability to detect any contribution. In the Scottish Highlands, browsing impact of mountain hare on saplings (mostly birch) was only detectable at high densities (45 hares km −2 ), while below one hare km −2 no effect was found (Miller et al. 1982). Although hares are capable of browsing pine, grasses comprise the majority of their diet (Schai-Braun et al. 2015) and their impact on Scots pine is likely negligible where alternative forage is widely available.

Effects of ground vegetation composition
After accounting for the density of fallow or large deer, the intensity of leader damage was not modified by ground vegetation composition, availability of alternative browse, or percentage of calcareous soil, contrary to the "neighbour contrast defence" hypothesis (Bergvall et al. 2006). In contrast, in highland Scotland, damage to Scots pine by red deer was less when high quality graminoid forage was available but greater when ground vegetation predominantly comprised heather (Palmer and Truscott 2003). It is possible that in the current study ground vegetation composition may not have influenced crop damage due to the generally low availability of bramble in restocked stands (that could shield Scots pine and offer alternative forage), in contrast to local bramble dominance in pre-thicket and mature stands. But fallow, being generalist feeders (Obidziński et al. 2013) may not be diverted from browsing on Scots pine saplings. In contrast, more specialist browsers such as roe and white-tailed deer are more likely to be diverted by the presence of palatable species in the understorey (Tuomi and Augner 1993;Hjältén et al. 1993;Augustine and Jordan 1998;Ward et al. 2008).

Deer density thresholds to alleviate crop damage
To reduce leader damage intensity on restocked Scots pine to only 20% per annuum would require control effort supressing fallow density towards zero (0.6 individuals km −2 ) or the combined large deer density to a low level (2.3 individuals km −2 ). Tolerable levels of damage depend on crop value relative to management costs; damage can be tolerated at higher levels (> 55%) if further reducing impacts (through culling) is expensive (Ward et al. 2004). Density thresholds for successful regeneration in this landscape may not apply to some other systems. Successful regeneration of lowland mixed woodland in England exposed to varying mixtures of sika, red, fallow, roe deer and muntjac occurred with a combined overall density of 10 individuals km −2 , higher than in our study . However, for lower productivity upland or moorland areas, impact density thresholds are around 4-5 large deer or 25 roe deer km −2 for native woodland regeneration (Putman et al. 2011), while for regeneration of Scots pine in highland Scotland, red deer density should not exceed 3-4 km −2 (Scott et al. 2000), similar to the large deer threshold of our study. For deciduous forest in north-western Pennsylvania regeneration occurred at a density of 7-9 white-tailed deer (Odocoileus virginianus) km −2 (Horsley et al. 2003), for balsam fir (Abies balsamea) in Canadian boreal forest with white-tailed deer 15 km −2 (Hidding et al. 2012). However, all of these studies consider natural regeneration; in contrast browsing predisposition may be greater with clear-felling and replanting (Reimoser and Gossow 1996;Partl et al. 2002), particularly as herding species can aggregate and feed for a prolonged period in a large restocked stand offering concentrated food with reduced travel cost. In contrast, higher fallow deer density may potentially be tolerated in forest management by small-gap creation and natural regeneration (Peterken 1996).

Conclusions
In the forest landscape we investigated, cultural preferences have led to an emphasis on controlling muntjac, that is regarded as un undesirable and damaging invasive species. However, our evidence instead shows large numbers of fallow deer are causing most crop damage. Fallow have increased to high abundance since colonising the study area in the late nineteenth century (Wäber 2010), partly due to a tolerance due to their aesthetic appeal and economic value for recreational stalking, and have similarly increased to high density in many other parts of the UK. Achieving the low density required to supress impacts to an acceptable level would be a considerable challenge for deer management. Fallow and red deer are wider ranging than roe or muntjac and, if heavily culled, can move to areas where they perceive a lower risk of predation (Tolon et al. 2009;Proffitt et al. 2009). Although fallow and red deer management must operate at a landscape-scale, requiring coordination across landowners (POST 2009), reductions to a negligible level may still be impracticable. In addition, a reduction in fallow deer numbers may lead to more roe deer establishing home ranges in the forest due to reduced inter-specific competition (see: Ferretti et al. 2011). Greater roe deer densities may then impact more substantially on the forest crop than at present, and some management of roe or muntjac should, therefore, also take place despite focussing most of cull efforts on fallow deer. In the forest landscape considered here, forest managers should also consider alternative approaches to mitigation, including planting tree species more resilient to deer damage (such as larch-Larix spp.-and sycamore-Acer pseudoplatanus; Perks et al. 2005;Hein et al. 2009), greater reliance on natural regeneration instead of planting, or use of fencing in areas where crop impacts are unavoidable and crop value offsets this expense. More generally, our findings emphasize the crucial importance of taking an evidence-based approach when managing deer in a multispecies system, and of understanding which deer species are responsible for impacts.