Probability of Bark Stripping Damage by Red Deer ( Cervus elaphus ) in Austria

Bark stripping by red deer (Cervus elaphus) causes considerable damage to Austrian forests, however, the incidence of bark stripping was never examined from large scale survey data. In this manuscript we present a logistic regression model for bark stripping damage (static model) and a model for recent (5-year period) bark stripping damage to previously undamaged trees (dynamic model) developed from Austrian National Forest Inventory data. Both models showed bark stripping damage to be most frequent in core red deer habitat areas and less frequent in less suitable habitat. Damage was concentrated at elevations of 400–1200 m and in alluvial forests (only static model). Norway spruce (Picea abies), European ash (Fraxinus excelsior), Sweet chestnut (Castanea sativa) and Sorbus spp. had 11–12 times more injuries than all the other species. Red deer preferred the smallest trees with a breast height diameter of 5 cm for bark stripping and damage probability decreased rapidly for trees with a breast height diameter greater than 25 cm. Our static model showed a maximum of bark stripping damage in stands with a mean height of 20 m. In the dynamic model the probability for bark stripping damage decreased with decreasing mean height. Also, in the static model the probability for bark stripping damage increased with increasing spruce proportion and with increasing stand density whereas in the dynamic model the proportion of previous bark stripping damage was a good predictor. Goodness of fit and discrimination of both models were good. In combination with forest growth models, the bark stripping models can be used to predict the risk of damage associated with different forest and habitat management options.


Introduction
Wild ungulates in Europe increased rapidly during the last century (Gill 1990).This has led to an increasing impact on wild and cultivated vegetation.In forest ecosystems it is widely recognised that ungulates can markedly alter a forest's structure and species composition by browsing, fraying trees with antlers and stripping bark.Perhaps bark stripping is the most serious damage to forests because trees affected by bark stripping do not readily recover and a considerable number of trees may be killed (Akashi and Nakashizuka 1999).
After infection there is a delay of several months before fungal species begin to advance from the infection site.Stain development is then very rapid, until it eventually slows and may cease altogether.For Norway spruce the extent of rot is usually 2-3 m after few years, and 4-5 m after 25 years, and can eventually reach 10 m (Vasiliauskas 2001).Consequently, the most valuable lower part of the stems becomes unsuitable for timber within 2-3 years after debarking (Randveer and Heikkilä 1996) and value losses for trees with decay range from 5-15% compared to healthy trees (Roeder 1970, Meng 1978).Apart from these considerable economic losses debarking is also a thread to forest stand stability.There is a risk of breakage at the wound (Randveer and Heikkilä 1996) which can be aggravated by high winds and snow (Fruhmann and Roeder 1981).
Among the ungulate species present in Austria red deer (Cervus elaphus L.), sika deer (Cervus nippon L.), fallow deer (Dama dama L.) and mouflon (Ovis musimon Schreiber) are known to remove bark.In the majority of cases, however, red deer are the damaging agents (Ueckermann 1960, Prien 1997).Bark stripping by red deer in Austria occurs usually in winter (Völk 1997), when the bark is gnawed-off with the aid of the lower incisors and eaten, but sometimes also in spring or summer when the bark can be more easily removed.Red deer typically removes bark from the main stem at a height from 0.8-1.7 m (2.0 m) above ground (Prien 1997).The damage caused by bark-stripping of red deer is a serious problem in Austria.According to the National Forest Inventory 4.5% of the total standing volume, which represents 5.4% of the total standing basal area and 7.9% of the total stem number, are bark stripped (BFW 2005).In certain areas bark stripping has become so common that future stands are placed in jeopardy.
Many hypotheses have been proposed to explain why ungulates strip bark.These include low availability of other food sources, the need for nutrients contained in bark, the need for a constant basal component of fibre and deer ranging behaviour.Despite the number of alternative possibilities most authors have agreed that food shortage is the most frequent cause of bark stripping (Ueda et al. 2002).In cool temperature zones mainly deep snow cover reduces the resource availability and may lead to more bark stripping.The extent of bark stripping damage is known to depend on red deer population density, although other factors such as stand age, tree species, stand density and winter snow cover may modify this relationship (Prien 1997).In addition, bark stripping in one year may increase the probability of bark stripping in the future (Hildebrandt 1959, Trisl et al. 1999).In forest management planning quantitative estimates of the risk of bark stripping damage would be helpful.However, until recently there have been little quantitative studies on bark stripping to explain the variation in bark stripping activity.Further, most existing studies were confined to small areas.
The aim of this analysis was to examine the incidence of bark stripping damage from large scale survey data.Both the accumulated damage and recent bark stripping damage in a five year period was considered.More specifically the following issues were addressed: i) to identify conditions which make trees and forests vulnerable to bark stripping damage, including previous bark stripping, ii) to quantify and compare the effects of individual bark stripping factors.

Data
To develop models for bark stripping damage we used the data of the Austrian National Forest Inventory collected during the periods 1981-1985, 1986-1990 and 1992-1996, which gave us the possibility to control for a large number of confounding variables.The Austrian National Forest Inventory is a systematic permanent inventory.Plots are clustered at the 4 corners of a 200 m × 200 m square.The clusters themselves are distributed over Austria with a distance of 3.89 kilometres.Each year 1100 clusters are assessed so that the total inventory comprises 5500 clusters and 22 000 sample plots (Forstliche Bundesversuchsanstalt 1981, Forstliche Bundesversuchsanstalt 1994).For this study, however, only 5800 plots that are not situated on forest edges or non-forest land, that are undivided by forest roads and that are in a single ownership were used.Permanent plots were established from 1981 to 1985.Plot centres were marked by a hidden iron stake buried underground and plot coordinates were recorded.To relocate the plot centres at the following two forest inventories 5 and 11 years later, a metal detector was used.Red deer habitat suitability and red deer presence was assessed for each cluster (4 sample plots).All other measurements were taken at the individual sample plots with the iron stake as plot centre.Plots size varied depending on the variables measured.
Each cluster of the forest inventory was classified according to red deer habitat suitability and recent presence of red deer into: 1) core habitat, red deer present 2) moderately suitable habitat, red deer present 3) moderately suitable habitat, red deer not present 4) not suitable habitat, red deer present and 5) not suitable habitat, red deer not present.Core habitat areas are defined as areas, that provide habitat of substantial quality year round.There are foraging sites, shelter and hiding cover and rutting stands.Core habitat areas should be managed with the objective of sustaining red deer viability.Core habitat areas are mostly found in the little fragmented landscape of the Austrian Alps.Moderately suitable habitat areas are usually areas adjacent to core habitat.Landscape fragmentation is however higher than in core habitat areas and the land is more intensively managed.Carrying capacity for red deer is accordingly lower.Red deer has spread to these areas because of the high population densities in core areas.Unsuitable habitat areas are highly fragmented areas of intensive agriculture (Malin 1989).Estimates of habitat suitability and recent presence of red deer were provided by regional foresters and are based on data drawn from the number of shot deer and hunter observations.
On each sample plot site variables were assessed using fixed radius plots of 300 m 2 (radius = 9.8 m).Sample trees were selected using two different methods: i) trees with a dbh > 10.4 cm were selected by angle count sampling with a basal area factor of 4 m 2 ha -1 variable radius plot with a different radius for each tree, radius [m] = dbh [cm]/basal area factor [m 2 ha -¹], on average 8 trees per plot), ii) trees with a dbh of 5 to 10.4 cm were assessed on fixed radius plots of 21 m 2 (radius = 2.6 m, on average 1 tree per plot).The following site variables were recorded: elevation to the nearest 100 m, slope to the nearest 10% and aspect to the nearest 45°, slope position, soil type, soil moisture, soil depth and vegetation type.Sample trees were recorded by their polar coordinates and marked by a nail at the base of the tree.Each selected tree was measured for dbh (mm) and height (m) and tree species was recorded.The social position of each tree was assessed using the classification of Kraft (1884).The occurrence of any visible bark stripping damage (fresh and healed wounds) was noted for each tree.Size, location of the wound and wound condition, however, were not assessed (Forstliche Bundesversuchsanstalt 1981, Forstliche Bundesversuchsanstalt 1994).Further, recently bark stripped trees were not distinguished from trees with old wounds at the assessment.Therefore the trees with recent bark stripping damage were defined as trees that did not have a bark stripping damage at the previous inventory, but had a bark stripping damage at the current inventory (Forstliche Bundesversuchsanstalt 1981, Forstliche Bundesversuchsanstalt 1994).
Quadratic mean diameter, Lorey's mean height, basal area of larger trees (Wykoff 1990), the percentile in the basal area distribution, crown competition factor (Krajicek et al. 1961), stand density index (Reineke 1933) and top height (largest 12 m 2 /ha, Pollanschütz 1973) were calculated from the inventory data.Further, the percentage of trees with bark stripping damage and the species proportion for all tree species as percentage of the stand's basal area were calculated.

Statistical Methods
The data were analysed using logistic regression techniques (Eq.1).Two models, one for the probability of bark stripping damage (static model) and one for the probability of recent (five-year period) damage to a formerly undamaged tree (dynamic model), were developed.Models were weighted by basal area to account for the different sampling techniques in the Austrian National Forest Inventory.Additionally, the dynamic model was weighted by period length.Weighting compensates for differences in the sample and population fraction and creates a representation of the dataset that would have been obtained under simple random sampling (King and Zeng 2003).Inferences about the intercept parameter b 0 (Eq. 1) are not possible without knowledge of the sample fraction (Hosmer and Lemeshow 2000).(E.g.: 152 895 trees were used to calculate the static model, 8363 trees were bark stripped, thus 5.47% of the observed cases were bark stripped.Since trees with a dbh > 10.4 cm were selected by angle count sampling, each tree represented a basal area 4 m 2 ha -1 .Trees with a dbh ≤ 10.4 were selected by fixed radius plots and the represented basal areas ranged from 0.9-4 m 2 ha -1 .The calculation of the basal area proportion of damaged trees yielded 5.39%.A logistic regression model without weighting yields an average damage probability of 0.0547 corresponding to the percentage of bark stripped cases, a weighted logistic regression model (weighted by basal area) yields an average damage probability of 0.0539 corresponding to the basal area proportion of bark stripped trees.) where p = probability for an event e = base of natural logarithm b 0 -b k are model coefficients, ß 0 = intercept, ß 1 -ß k estimated parameters for the kth independent variable x 1 -x k are explanatory variables, independent variables Model parameters were estimated using the maximum likelihood method.Explanatory variables were first selected by listing factors associated with bark stripping from the literature.
We then calculated the cross-classification of the occurrence of damage for different classes of the explanatory variables for the following variables: tree species, dbh, height, height/diameter ratio, stem number per hectare, basal area per hectare, stand density index (Reineke 1933), crown competition factor (Krajicek et al. 1961), dominant height, mean height, quadratic mean diameter, basal area of larger trees (Wykoff 1990), percentile in the basal area distribution, social position (Kraft 1884), species proportion for all tree species, proportion of evergreen trees, elevation, slope, aspect, slope position, red deer habitat suitability classes, alluvial forests and percentage of previous damage (dynamic model).Univariate logistic models were fit to check the assumption of linearity in the logit.Variables were then entered manually into the model.Variables for which the Wald χ 2 -statistic was significant at a = 0.05 were included.Models were tested using the Hosmer-Lemeshow goodness of fit test (Hosmer andLemeshow 1980, Lemeshow andHosmer 1982).The Hosmer-Lemeshow test is a Pearson Chi-square test, which groups the data into 10 classes based on the predicted probability.This test statistic is recommended, if continuous variables are in the model (Hosmer and Lemeshow 2000), because then the minimum of 5 observations per cell required for the calculation of the Pearson or Deviance χ 2 -test is frequently not met.
The final models were selected based on the log-likelihood.For both models pseudo-R 2 values were calculated.Pseudo-R 2 values range from 0 to 1, with 0 indicating that the model is no better than the mean and 1 indicating that it is a perfectly fitting model.The pseudo-R 2 statistic however does not vary linearly with goodness of fit and in contrast to linear regression models values > 0.2 are regarded as good (Litz 2000).
As measure for the model's ability to discriminate, the area under the ROC-curve (Receiver Operating Characteristics curve) was used (Hosmer and Lemeshow 2000).The area under the ROC-curve is 0.5 for models with no discriminative ability, whereas an area of 1.0 indicates perfect discrimination.In practice models with an area under the ROC-curve > 0.8 are considered good models (Hosmer and Lemeshow 2000).Further, to supplement more rigorous methods of assessment of fit, the classification table for the models was calculated.This table is the result of cross-classifying the observed cases of bark stripping, with a dichotomous variable whose values are derived from the estimated logistic probabilities.To obtain the derived dichotomous variable we compared the estimated logistic probabilities to uniformly distributed random numbers.From the classification table we calculated the percentage of correctly classified bark stripped trees, the percentage of correctly classified undamaged trees and the overall percentage of correctly classified trees.This approach of predicting group membership can be useful for model applications, where classification is the stated goal.But it is important to keep in mind that two models can not be compared by measures derived from a classification table since these measures are completely confounded by the distribution of probabilities in the samples upon which they are based (Hosmer and Lemeshow 2000).
The model effects were also expressed as odds ratios, which were calculated from the model parameters (Eq.2).The odds ratio approximates how much the frequency of damaged trees to undamaged trees changes when the variable changes by one unit.
where OR i = odds ratio for the ith explanatory variable b i = estimated parameter for the ith explanatory variable x i2 -x i1 = change of the ith explanatory variable by one unit

Results
Coefficients for the models for bark stripping damage are shown in Table 1.In Austria 5.4% of the total standing basal area are damaged by bark stripping and annually recent stem damage is found on 0.13% (0.7% in 5 years) of trees without a previous damage (Table 1).The following variables were significant in both models: red deer habitat suitability class, tree species, dbh, mean height and elevation.Further in the static model spruce proportion, stand density index and a dummy variable for alluvial forests are included.
In the dynamic model the percentage of previous bark stripping damage was important.Norway spruce (Picea abies (L.) Karst.),European ash (Fraxinus excelsior L.), Sweet chestnut (Castanea sativa Mill.) and Sorbus spp.(Sorbus aucuparia L., Sorbus torminalis Crantz, Sorbus aria Crantz, Sorbus domestica L.) were the species that were most severely affected by bark stripping.These species were 11-12 times more frequently bark stripped and recent bark stripping injuries were more frequent than for all other tree species (Fig. 1a, Fig. 1b).Apart from this strong main effect there was, however, no significant interaction of tree species with dbh or mean height.Most bark peeling damage was found in areas classified as core red deer habitat, in these areas also considerably more recent damage occurred.In areas of moderate habitat suitability, where red deer was not present yearround, considerable lower probabilities of bark stripping were observed.The effect of red deer habitat was more pronounced in the static model than in the dynamic model (Fig. 1a, Fig. 1b).
Red deer preferred trees with small breast height diameters for debarking (Fig. 2a, Fig. 2b).A maximum of bark stripping was found in stands with a mean height of 20 m (static model, Fig. 3a).In the dynamic model the probability for bark stripping damage decreased with decreasing mean height (Fig. 3b).
Bark stripping damage was most frequently found at elevations between 400 and 1200 m above sea level (Fig. 4a, Fig. 4b).In the static model the probability for bark stripping damage further increased with increasing stand density (Fig. 3a) and spruce proportion (Fig. 4a), and a higher proportion of damaged trees was found in alluvial forests (Table 1).In the dynamic model, the probability of a new bark stripping damage increased with the percentage of previous bark stripping damage.Only when damage reached 90% did this probability decrease again (Fig. 5).
Goodness of fit and discrimination for both bark stripping models was very good, as was indicated by pseudo-R 2 values > 0.2 and C-statistics > 0.8 for both models (Table 1).The models classified 17% (static) and 5% (dynamic) of damaged trees correctly.This was 3 and 7 times more than without the model, respectively.Of undamaged trees more than 95% were correctly classified, because of their high proportion (Table 1).
An important cause for different vulnerability to bark stripping between species is the process of bark formation.Even though most tree species have been reported to have been bark stripped by red deer, the period of vulnerability varies considerably.Debarking was concentrated to a relatively short period for trees with rapidly thickening bark (oak, Douglas fir, Scots pine, European larch), approximately 5-12 years.In contrast, species with slowly thickening bark, for example Norway spruce and European beech, can be debarked for 15-60 years and consequently these species are more seriously affected (Trisl et al. 1999, Prien 1997).
Another important factor is bark adhesion.Bark of beech is notoriously difficult to detach in winter, which probably explains why Norway spruce in Austria is severely affected whereas little damage is observed on European beech.Bark thickness or bark adhesion alone, however, do not explain the differences between each species.Further differences between species occur due to a different efficiency in wound healing.Wounds of beech, linden and poplar heal more rapidly, and therefore, after a few years, scars are no longer visible.For Norway spruce, wound closure is considerably slower (Molnár and Schmitt 1994).Finally, different damage probabilities between species are influenced by the liking and preferences of red deer (Ueckermann 1960).

Influence of Tree Size
Susceptibility to bark stripping was strongly size dependent.The smallest trees with a dbh of 5 cm (minimum dbh at assessment) were most severely affected and damage probabilities decreased rapidly for trees with a dbh greater than 25 cm.Damage can begin at a dbh of 2-4 cm (Prien 1997), when the main stem becomes rigid and accessible.As dbh increases damage probabilities peak.Rinner (1996) reported highest damage percentages at a dbh of 4-11 cm.Comparably, Faber and Thorson (1996) found a mean diameter of 4 cm for damaged trees, when investigating moose bark stripping damage.Akashi and Terazawa (2005) also reported decreasing damage probabilities with increasing dbh for sika deer.Other researchers suggest rather higher peak diameters of 8-16 cm (Kamsker 1979) or 11-15 cm (Greutter 1985).These results might, however, be biased by the minimum diameter recorded in the study.
The smallest trees are preferred by red deer because of their thin bark with a favourable water content (Prien 1997).In general, bark decreases in digestibility as it ages, thus possibly explaining the lower attractiveness of larger trees; nutrient level and fibre content are other factors (Gill 1992a).Norway spruce is vulnerable to bark stripping at an age of 10-60 years (Prien 1997, Trisl et al. 1999), it reaches the most vulnerable diameters on average at an age of 25-35 years (Trisl et al. 1999), tough the most vulnerable age may vary considerably from site to site (Trisl et al. 1999, Welch et al. 1987).On sites with good tree growth the period of vulnerability is reached early, on poor sites, especially at high elevations, the age of vulnerability is reached later and lasts longer (Trisl et al. 1999).
Some authors suggest that percentages of bark stripping depend on the social position of the tree, there are, however, differing views whether suppressed or dominant trees are preferably used (Prien 1997).Welch et al. (1987) observed that below a mean diameter of 6.3 cm bark stripped trees were significantly larger than the plot mean (predominant and dominant trees), in a range of mean diameter of 6.3-12.7 cm stripped trees had average diameters (dominant and codominant trees) and in stands with a mean diameter above 12.7 cm, bark stripped trees were significantly smaller (suppressed trees).This supports the hypothesis that tree diameter is important for bark stripping, whereas the social position plays no or a minor role.
In the static model most damage was found at a mean height of 20 m, approximately at the end of the stem exclusion phase.This indicates that the damage accumulates over a long time period, before it can be reduced by thinning.Similarly, also other workers observed most bark stripping damage in pole stands (Lichtenwanger 1983).In contrast the probability of a recent damage decreased with mean height.But it has to be kept in mind that the dynamic model adjusts for previous damage percent, and that the damage percent is still increasing for small tree heights.

Stand, Site and Habitat Factors
A higher incidence of bark stripping was found in areas of high stand density, probably because the bark remained thinner for longer.Also in dense stands senescence and death of lower branches occurs earlier (Prien 1997).Additionally the dense canopy cover in evergreen stands efficiently intercepts snow fall and thus provides ideal shelter, when snow depth is high and limits the range of deer movement in winter.Pichler (1984) showed for a forest enterprise with serious bark stripping damage, that stands with recent incidence of damage had considerably higher stem numbers than stands were no bark stripping occurred.Further in Slovenia (Debeljak et al. 2001) and Austria red deer seems to select young stands with a high spruce proportion for winter habitat.
Even though the environment of subalpine meadows is exceptionally attractive during summer (Völk 1997, Homolka andMatouš 1999), red deer usually migrates to lower elevations during winter because the climate is milder.Regions at very low elevations are also avoided because of a high pressure from settlements and agriculture, which limits red deer activity (Debeljak et al. 2001).This is reflected in both models, where most damage is found at an elevation of 800-1400 m.Analogue patterns, even though absolute elevation may vary between different regions, are frequently found (e.g.Lichtenwanger 1983, Debeljak et al. 2001).In contrast, Rinner (1996 ) found the most severely damaged stands at elevations ranging from 1500-1800 m.
The results indicate distinct differences in damage percentages for different red deer habitat classes.This suggests that habitat classification coincides with population density, resulting in higher damage percentages where red deer is more frequent.A relatively high red deer population density is undoubtedly a contributing factor.Strong direct relationships between damage and ungulate population are usually reported, where population levels change considerably in a short time period.Motta (1997), for example, found an increase in bark stripping damage the first years after the introduction of red deer in the Italian Alps.The association between ungulate density and bark stripping is, however, often masked and confounded when populations fluctuate less.Then often only weak correlations (Völk 1997, Jalkanen 2001) or no correlations between ungulate density and damage are found.Methodological difficulties in animal census in forest ecosystems usually aggravate the problem (Motta 1997).
A higher incidence of bark stripping damage was found in alluvial forests.In areas of intensive agriculture at lower elevations, only alluvial forests provide suitable habitat for red deer, because of their spatial connectivity.Accordingly, most damage is also concentrated there.
In the dynamic model, strong correlations were found between previous and recent stem damage.This may reflect the habit of red deer to return to the same winter feeding areas, thus increasing the local risk.Bergqvist et al. (2001) also reported strong correlations between recent and previous damage for bark stripping by moose.Previous damage percent might also be an indicator for predisposing factors that could not be considered in the model, because the respective variables were not available from the forest inventory data.Such factors may include the distance to winter feeding areas, edge density, predation or tourism.
Winter feeding concentrates deer and results in damage to nearby trees (Kamsker 1979), approximately at a circumference of 100 m (Gossow 1988).Also, in areas adjacent to forest edges, more damage is observed (Pichler 1984).If deer are disturbed by human activities (hunting, tourism, wood harvesting) it has to "wait" excessively long under cover, before foraging and consequently causes a lot of damage (Gossow 1988).
Adverse weather (long cold periods, snow height and dryness) also cause higher percentages of bark stripping damage (Prien 1997) and in winter deer might concentrate in areas of favourable radiation, wind and snow (Gossow 1988).Aspect however did not qualify as predictor in the two models.Deer sometimes also selects moderately steep areas which are undisturbed by humans but avoids very steep areas (Lichtenwanger 1983, Gossow 1988).In this study, however, we were unable to demonstrate any effect of slope.

Model Application
For predictions of bark stripping probabilities both, the static and the dynamic model could be used, depending on the available information on previous damage percent, which is necessary for the dynamic model.Since the models are based on large-scale inventory data, predicted probabilities of bark stripping damage are within a realistic range.Predictions may not be exact for each single tree, but they may be on average for a larger area.Predictions of damage can be made for different areas of red deer habitat suitability and red deer presence or for changes in the structure and management of forests.The models could also be linked with models of growth and yield or models of habitat suitability and risk levels for different forest and habitat management options could be estimated.For economic evaluations both models could be combined with models, that estimate the expected loss in tree volume or tree quality caused by bark stripping.
In general both models suggest that damage at a given deer density can be reduced by silvicultural management.The distinct differences between species suggest that planting of not susceptible species in areas with critical deer densities might be one way to reduce damage.The preference of red deer for dense stands indicates that thinning might be another way to reduce bark stripping damage.Thinning would, however, have to start very early, because deer prefer trees with small breast height diameters.

Fig. 1 .Fig. 2 .Fig. 3 .
Fig. 1.Odds ratio (OR) in the static (a) and the dynamic (b) model for bark stripping for different red deer habitat suitability classes (core habitat, red deer present; moderately suitable habitat, red deer present; moderately suitable habitat, red deer not present; unsuitable habitat, red deer present; unsuitable habitat, red deer not present) and for susceptible and not susceptible species.Susceptible species are: Norway spruce, European ash, Sweet chestnut and Sorbus spp.RD = Red deer

Fig. 4 .
Fig. 4. Probability (Prob) of bark stripping in relation to elevation in the static model (a) and the dynamic model (b).For the static model the relationship is plotted for different spruce proportions.

Table 1 .
Coefficients and standard errors of significant variables (a = 0.05), model statistics and the percent of trees correctly classified by the model (probabilities were compared to random numbers) for the model for bark stripping damage (bark stripping static) and for the model for recent bark stripping damage to previously undamaged trees in a five-year period (bark stripping dynamic).