Effects of Bark Stripping on Timber Production and Structure of Norway Spruce Forests in Relation to Climatic Factors

The aim of this study was to assess the effects of bark stripping caused by sika deer (Cervus nippon [Temminck]) on the production and structure of young Norway spruce (Picea abies L. Karst) forest stands (41–43 years). Production parameters, structure, diversity, and the dynamics of radial growth in selected forest stands in relation to climatic conditions were evaluated. Similar to other production parameters, stand volumes showed lower values on research plots heavily damaged by bark stripping (290 m3 ha−1) compared to stands with lower tree stem damages (441 m3 ha−1). A significant decrease in stem volume was recorded for trees with stem circumference damage higher than 1/3 of the stem circumference. In most cases, the trees were damaged between the ages of 10–23 years, specifically the radial growth was significantly lowered in this period. The diameter increment of damaged trees dropped to 64% of the healthy counterparts in this period. Bark stripping damages reached up to 93% of the stem circumference with a mean damage of 31%. Stem rot was found on 62% of damaged trees. In our study area, with respect to the terms of climatic conditions, precipitation had a higher effect on radial growth of the Norway spruce compared to temperature. The main limiting climatic factor of tree growth was the lack of precipitation within a growing season, particularly in June of the current year.


Introduction
Globally, most of the Earth's surface has been changed by human activities.One of the most significant reasons these changes occur is agricultural land use [1][2][3].However, amended agricultural landscapes were successfully colonized, e.g., by increasing wildlife populations of large ungulates which have been documented across many countries in Europe [4][5][6][7][8][9].This increase has been driven not only by changes in land use, but also with the absence of natural predators (large carnivores), which have played a crucial role in the reduction of ungulate populations [10].One of the species with increased population density is sika deer (Cervus nippon [Temminck]), which was first introduced from Japan to Europe after 1860 [11,12].At the moment, the most abundant populations of sika deer within Europe can be found on the British Islands, in Germany, or in the Czech Republic [11,[13][14][15][16][17].
High population density of hoofed game is naturally connected with damages found on agricultural crops [18] or forest stands [19][20][21] within the area.Young trees (of both coniferous Forests 2019, 10, 320 2 of 21 and broadleaved species) are especially affected by browsing [22], as the assimilation shoots and sprouts provide natural aliment for deer [23]-though older trees could also be damaged by bark stripping.Browsed seedlings generally show significant deceleration in growth [24,25], and in extreme cases repeated browsing of young trees could lead to tree mortality [26].The intensity of browsing is also influenced by environmental factors, such as altitude, type of vegetation or the structure of present phytocenosis [22,27].
To date, bark stripping caused by ungulates has not been sufficiently explained.One reason could be the relatively high content of nutrients found in the bark [28].The occurrence of bark stripping is also higher in the areas with intensive human activity, where the deer usually hide under the cover of young forest stands [29].The intensity of bark stripping is not only determined by population density of deer, but also by forest stand characteristics [30,31] including the age of the forest stands [21] or social position of the tree [32].The damages are subsequently reflected in the economic evaluation of damaged wood [33,34], mainly because of frequent infections and spreading of the fungal pathogens in the stem [35][36][37].Nevertheless, economic loss is one of several problems within these stands.Severe infection of multiple trees could cause the disintegration of the entire stand due to reduced resistance to both biotic and abiotic factors [38].If more than 90% of stem circumference is stripped, the whole tree usually dies [27].For a more complex evaluation of bark stripping impacts, the following four aspects are crucial: the length of time over which the bark stripping occurs, the size of the damaged trees or trunks, the frequency of repeated damage, and the recovery of the wounds [34].
Bark stripping in countries with a high population density of sika deer is an important point of interest in forestry because of the aforementioned reasons.In Japan, similar problems were found in local tree species, e.g., Sakhalin fir [Abies sachalinensis (Fr.Schmidt)] Masters or Japanese Larch [Larix kaempferi (Lamb) Carr.] [39], whereas the European populations of sika deer usually cause damage to the Norway spruce [Picea abies (L.) Karst.][40] or the sitka spruce [Picea sitchensis (Bong.)Carr.] found in Scotland [34,41].
One of the largest European populations of sika deer is located in the Czech Republic [13,14].Within the last 20 years, the number of hunted individuals of sika deer increased from 3835 in 1995 to 17,106 in 2018 [42] under the same regulations, but the actual population size is undoubtedly much higher.Damage intensity naturally increases with higher population density, especially in Norway spruce stands across Central and Western Europe [37,43,44].Damaged stands generally cause a significant decrease in timber production, loss of stability, thus lowering the ability to provide their ecological and environmental functions [18,33,[45][46][47][48].The aim of this paper was to evaluate: (1) the timber production, structure, diversity, and interactions of measured parameters in Norway spruce stands exposed to different intensity of bark stripping caused by sika deer; and (2) the effects of bark stripping and stem rot on radial growth of Norway spruce in relation to climatic factors (temperature and precipitation).The study is focused on Norway spruce stands that are heavily influenced by sika deer, and stands without significant deer impact in comparable acidic sites in the western part of Bohemia (Czech Republic).

Study Area
The research was conducted on 16 permanent research plots (PRPs) established in planted monospecific Norway spruce forest stands in the western portion of the Západočeská Upland area (Natural Forest Area 6) in the Czech Republic (Figure 1).PRPs were situated in 4 localities in total (4 PRPs on each locality).Two localities (PRPs 1-8) were not under significant game influence (hunting district Horní Kozolupy, marked "LD plots" as "low game damage" further in text), however, the other two localities (PRPs 9-16) were under strong influence of Sika deer (hunting district Rochlov, marked "HD plots" as "high game damage" further in text).In the hunting district of Rochlov (HD plots), on average, 135 sika deer (Cervus nippon Temminck) individuals and 36 roe deer (Capreolus capreolus plots), on average, 135 sika deer (Cervus nippon Temminck) individuals and 36 roe deer (Capreolus capreolus Linnaeus) individuals were hunted per year between 2013-2017.In the hunting district of Horní Kozolupy (LD plots), these numbers were significantly lower, on average, with 34 sika deer individuals and 26 roe deer individuals hunted per year within the same time period.The studied spruce stands on all PRPs were of similar age (41-43 years), according to the forest management plan.Site quality index and soil type (Cambisol) were identical between all PRPs.The stands were artificially established and protected by an iron fence for 9 to 11 years (depending on locality) after plantation.The altitude of the study locality ranged from 471 to 496 m a.s.l., with the bedrock mostly composed of slate and phyllite.The mean annual temperature in the study area was approximately 7.7 °C [49], with the minimum temperature recorded in January (−1.9 °C), and the maximum recorded in July (17.9 °C).The mean annual precipitation ranged from 500 to 590 mm yr -1 with the lowest precipitation recorded in February (21 mm) and the greatest precipitation in June (67 mm), see Figure 1.The duration of snow-coverage was approximately 40 days, with the average maximum snow depth recorded at 20 cm.The vegetation season lasted ~155 days with a mean precipitation of 360 mm and a mean temperature of 13.9 °C.The study area qualified as a humid continental climate zone, characterized by warm to hot, humid summers and cold to severely cold winters (Dfb) according to the Köppen climate classification [50].According to the detailed Quitt climate classification [51], the area can be classified in the cold climatic region (subregion CH 6).In terms of phytocoenology, all PRPs belong to Querceto-Fagetum acidophilum (Acidic Oak-Beech) and the alliance Luzulo-Fagion sylvaticae Lohmeyer et Tüxen in Tüxen 1954.

Data Collection
In November of 2017, Field-Map technology [52] was used to establish 16 square, individual, 100 m 2 (10 × 10 m) PRPs.The measured characteristics within the tree layer were as follows: position of the tree and its crown projection (minimally at four directions perpendicular to each other) measured The mean annual temperature in the study area was approximately 7.7 • C [49], with the minimum temperature recorded in January (−1.9 • C), and the maximum recorded in July (17.9 • C).The mean annual precipitation ranged from 500 to 590 mm yr -1 with the lowest precipitation recorded in February (21 mm) and the greatest precipitation in June (67 mm), see Figure 1.The duration of snow-coverage was approximately 40 days, with the average maximum snow depth recorded at 20 cm.The vegetation season lasted ~155 days with a mean precipitation of 360 mm and a mean temperature of 13.9 • C. The study area qualified as a humid continental climate zone, characterized by warm to hot, humid summers and cold to severely cold winters (Dfb) according to the Köppen climate classification [50].According to the detailed Quitt climate classification [51], the area can be classified in the cold climatic region (subregion CH 6).In terms of phytocoenology, all PRPs belong to Querceto-Fagetum acidophilum (Acidic Oak-Beech) and the alliance Luzulo-Fagion sylvaticae Lohmeyer et Tüxen in Tüxen 1954.

Data Collection
In November of 2017, Field-Map technology [52] was used to establish 16 square, individual, 100 m 2 (10 × 10 m) PRPs.The measured characteristics within the tree layer were as follows: position of the tree and its crown projection (minimally at four directions perpendicular to each other) measured by Field-Map technology (IFER, Jílové u Prahy, Czech Republic), diameter at breast height (DBH 130 cm ≥ 4 cm; to the nearest 1 mm) measured by Mantax Blue metal caliper (Haglöf, Långsele, Västernorrland, Sweden), and tree and live crown base height (to the nearest 0.1 m) measured by Laser Vertex hypsometer (Haglöf, Långsele, Västernorrland, Sweden).Game damage was evaluated and classified using the methodology of the Institute of Forest Ecosystem Research, Ltd. (IFER, Jílové u Prahy, Czech Republic) and Forest Management Institute (ÚHÚL, Brandýs nad Labem, Czech Republic) using national inventory measurements [53,54].Bark stripping was measured along the girth of the tree stem using a tape (to the nearest 1 mm) at DBH and subsequently divided into three categories: healthy trees (damaged circumference ≤ 1/8 of the stem girth, 134 trees in total), low game damage (damaged circumference > 1/8 and ≤ 1/3 of stem girth, 62 trees in total), and high game damage (damaged circumference > 1/3 of stem girth, 54 trees in total) [53,55].
In 2017, increment core samples of 16 spruce trees on each of the 4 localities (4 cores on each PRP) were taken at DBH using a Pressler auger (Haglöf, Långsele, Västernorrland, Sweden).Thirty-two core samples were taken perpendicularly to the tree stem axis on the site damaged by bark stripping (HD plots, hunting district Rochlov) and another 32 samples were taken in an opposite, undamaged site (LD plots, hunting district Horní Kozolupy).Trees were randomly chosen by the RNG function (Excel, Microsoft).Tree-ring widths were measured to the nearest 0.01 mm with an Olympus (Tokyo, Japan) stereo microscope on a LINTAB measurement table (RINNTECH, Heidelberg, Germany) and recorded by the TSAPWIN software (RESISTOGRAPH, Heidelberg, Germany).Throughout the dendrochronological analysis, the occurrence of bark stripping damage (year) and stem rot (the percentage of diameter affected by rot) were visually recorded by significant color changes (fungi, black coloring in stripping damage, etc.), separations of ring widths, occurrences of defects (resin, rind galls, etc.), and changes in the mechanical structure of each sample.

Data Analyses
Based on the measured dendrometric data, stand characteristics including stand volume [56], stand density index (SDI) [57], crown closure [58], and crown projection area were computed for the tree layer.Height curves were constructed using the Näslund height-diameter function [59] (see Supplementary File S1-List of indices), and also the relationship between crown size and DBH was evaluated via linear regression.
To assess stand diversity, the following indices were computed using Sibyla 5 software [60]: Arten-Profile index as an indicator of vertical structure (using ratio of the basal area of tree species in stand layers) [61], diameter and height differentiation as indicators of structure differentiation (using ratio between tree diameters, heights and distances) [62], crown differentiation index (using minimum and maximum crown diameter), and total diversity index as indicator of complex stand diversity (using tree species diversity, diversity of vertical structure, diversity of tree spatial distribution, and diversity of crown differentiation) [63].The PointPro 2 software (Zahradník & Puš, CULS, Prague, Czech Republic) was used for the computation of horizontal structure index.The spatial distribution of trees on PRPs was evaluated according to the Clark-Evans aggregation index (using real distances between two nearest trees and plot parameters) [64], such as horizontal structure of damaged trees.The criteria of diversity indices are shown in Table 1, and a detailed description is given by Pretzsch [65].Detailed information about the computation of selected indices is provided in Supplementary File S1-List of indices.Tree-ring increment series were individually cross dated (removal of errors connected with the occurrence of missing tree rings) using t-tests, and subsequently, visually checked according to the Yamaguchi [66] method.Particular curves from PRPs were detrended in a standard way and mean tree-ring series were created from collected data in the R software [67].Annual tree-ring width was standardized by tree age and represented by tree ring width index.Package dplR [68] was used for creating smoothing spline for a different time window using built-in method to detrend, filter, and visualize the tree-ring series data [69].
Standard age detrended chronology was interleaved and plotted with a variety of splines (8 and 32-year splines) for removing short-term fluctuations.The analysis of negative pointer years was conducted according to Schweingruber et al. [70].For each tree, the pointer year was identified as an extremely narrow tree ring that did not reach 40% of the average increment from the previous four years.A negative year was defined as the occurrence of a strong increment reduction in at least 20% of the trees on their respective plot.For modeling the diameter increments in relation to climate characteristics (monthly precipitation and temperature), the DENDROCLIM2002 software [71] was used for statistical relationships between climate and tree radial growth using correlation and response functions.The software used bootstrapped confidence intervals to estimate the significance of both correlation and response function coefficients.In this case, we used correlation and response functions for monthly precipitations and temperatures from May of the previous year to September of the current year in range with our climatic data.Climate data  were obtained approximately 30 km away from the meteorological station site in Kralovice (468 m a.s.l.; WGS84 49 • 58 51 N, 13 • 29 06 E).
The relationship between production parameters (DBH, stand volume, height) of trees with different types of bark stripping damage (healthy trees, low game damage and high game damage; 250 trees in total) was tested using the Kruskal-Wallis test with subsequent multiple comparisons (according to Siegel & Castellan [72]).The assumptions of ANOVA (normality) was not met in all cases.For modeling the stem volume of each particular tree (based on DBH and stem circumference damage), a generalized linear model with gamma distribution was used.For evaluating the differences in stem volume predicted by the model between LD and HD plots, separate models for LD and HD plots were constructed, and their differences were evaluated.Parameters for the model were chosen to be easily measurable in the field (DBH and circumference damage).An inverse value was used as a link function.All computations related to the Kruskal-Wallis test and the generalized linear model were performed in R software [67].Three-dimensional plots were made using Gnuplot 5.2.The unconstrained principal component analysis (PCA) was performed in the CANOCO 5 program [73] to analyze the relationships between production parameters, diversity, structural characteristics, game damage, and similarity of all 16 PRPs.Data were centered and standardized during the analysis.The results of the PCA analysis were visualized in the form of an ordination diagram.Situational maps were created in the ArcGIS 10 software (ESRI).

Production of Spruce Stands
In 2017, the number of trees per hectare was comparable (difference of 6.1%) in PRPs between both hunting districts, specifically ranging from 1100 to 2200 trees in PRPs established on LD plots, and 1200 to 2000 trees in PRPs established on HD plots.However, a greater difference (19.0%) was observed in the stand density index (SDI from 0.81 to 0.99 on LD plots and from 0.49 to 0.97 on HD plots) with both averaging higher indicators of stand density in favor of HD plots.Crown closure ranged from 77.3-93.9% in PRPs and showed higher variability in HD plots (by 8.0%), such as in SDI (by 62.9%).Differences were also found in production characteristics.A mean DBH of 20.03 cm (±3.42 SD) was observed in LD plots, and a mean DBH of 17.49 cm (±1.37 SD) was found in HD plots.A similar trend was found with respect to mean height: 20.07 m in LD plots and 17.53 m in HD plots.The mean total stand volume on LD plots was 440.5 m 3 ha −1 (±112 SD) and 289.9 m 3 ha −1 (±73 SD) on HD plots.The total stand basal area ranged between 41.0 and 64.7 m 2 ha -1 within LD plots and 23.7-51.0m 2 ha -1 within HD plots.The dominant height (HDOM) was comparable in all localities (1LD-4LD: 26.5 m, 4LD-8LD: 26.3 m, 9HD-12HD: 22.0 m, and 13HD-16HD: 24.2 m), however, in HD plots, where individual tree tops had an occurrence of bark damage and rotting, HDOM was lower.Detailed information is presented in Table 2. Näslund height-diameter functions were computed to evaluate the differences in height to diameter ratio between selected damage levels, see Figure 2. The coefficient of determination was observed to be 0.61 with respect to healthy trees and 0.54 with respect to trees with high game damage (values are also depicted in the plot).Trees with small damage showed a low coefficient of determination.Computed parameters for healthy trees: a = 1.00, and b = 0.18; trees with small damage: a = 0.46, and b = 0.22; trees with high game damage: a = 1.09, and b = 0.18).The relationship between crown diameter and DBH was assessed by linear regression (Figure 3).Observed functions altogether with coefficients of determination are presented in the respective figure.Crown diameter of healthy trees appeared to follow a linear trend with stable variation around the regression line and a relatively high coefficient of determination (R 2 = 0.61) when related to DBH, whilst trees with greater game damage showed a rather high variation at lower DBH resulting in a relatively low coefficient of determination.The relationship between crown diameter and DBH was assessed by linear regression (Figure 3).Observed functions altogether with coefficients of determination are presented in the respective figure.Crown diameter of healthy trees appeared to follow a linear trend with stable variation around the regression line and a relatively high coefficient of determination (R 2 = 0.61) when related to DBH, whilst trees with greater game damage showed a rather high variation at lower DBH resulting in a relatively low coefficient of determination.

Stand Diversity and Structure
The prevailing spatial pattern (horizontal structure) of trees was random with a tendency of regular distribution on all PRPs, and a significantly regular structure was observed on three plots (7 LD, 14 HD, and 15 HD; see Table 3).The vertical structure (according to A index) was low to moderately-varied on LD plots (A = 0.103-0.494)and moderately to strongly-varied on HD plots (A = 0.218-0.627).Diameter differentiation was low on both variants of plots (TMd = 0.108-0.225),with respect to height differentiation (TMh = 0.063-0.170).Higher structural diversity was generally observed on HD plots compared to LD plots.Total stand diversity, B, (including species diversity, crown differentiation, vertical and horizontal structure) identified monotonous structures on all PRPs, while the highest total diversity was observed on PRP 11 (high game damage plot, B = 3.041) and the lowest on PRP 3 (low game damage plot, B = 1.952).For detailed information see Table 3.

Stand Diversity and Structure
The prevailing spatial pattern (horizontal structure) of trees was random with a tendency of regular distribution on all PRPs, and a significantly regular structure was observed on three plots (7 LD, 14 HD, and 15 HD; see Table 3).The vertical structure (according to A index) was low to moderately-varied on LD plots (A = 0.103-0.494)and moderately to strongly-varied on HD plots (A = 0.218-0.627).Diameter differentiation was low on both variants of plots (TMd = 0.108-0.225),with respect to height differentiation (TMh = 0.063-0.170).Higher structural diversity was generally observed on HD plots compared to LD plots.Total stand diversity, B, (including species diversity, crown differentiation, vertical and horizontal structure) identified monotonous structures on all PRPs, while the highest total diversity was observed on PRP 11 (high game damage plot, B = 3.041) and the lowest on PRP 3 (low game damage plot, B = 1.952).For detailed information see Table 3.

Diameter Growth of Norway Spruce
A comparison of the average tree-ring curves for all 16 PRPs shows their high rate of agreement among them (t-tests ≥ 5.1).This consistency allowed us to develop a regional standard chronology for the spruce stands in the Zapadoceska Upland region.The regional standard tree-ring chronology indicated a difference in the dynamics of radial growth, especially on HD plots (Figure 4).Until 1992, radial increments of spruce on damaged plots was considerably lower that on LD plots.For example, in 1983 the mean radial growth reached only 28.8% of normal (healthy) diameter increments.Since the impact of game attack subsided in 1992, radial growth of both variants (healthy vs trees with high damage) was relatively balanced.

Diameter Growth of Norway Spruce
A comparison of the average tree-ring curves for all 16 PRPs shows their high rate of agreement among them (t-tests ≥ 5.1).This consistency allowed us to develop a regional standard chronology for the spruce stands in the Zapadoceska Upland region.The regional standard tree-ring chronology indicated a difference in the dynamics of radial growth, especially on HD plots (Figure 4).Until 1992, radial increments of spruce on damaged plots was considerably lower that on LD plots.For example, in 1983 the mean radial growth reached only 28.8% of normal (healthy) diameter increments.Since the impact of game attack subsided in 1992, radial growth of both variants (healthy vs trees with high damage) was relatively balanced.
The low radial increments that occurred in 1998 and 2015 qualified them as negative pointer years.The first half of 1998 had the warmest recorded temperatures (8.2 °C, long-term mean 6.1 °C) since 1961 (the beginning of meteorological measurements).In 2015, low increments were caused by extremely high temperatures (annual temperature 9.9 °C, long-term mean 7.9 °C), with low amounts of precipitation (for example in February 2 mm, long-term mean 21 mm).Moreover, PRPs established on HD plots were also found to be a negative pointer year in 2003.This year was the driest year recorded from 1961 through 2017 (304 mm, long-term mean 488 mm), especially with respect to the vegetation season (196 mm, long-term mean 325 mm).Regarding the effects of precipitation and temperatures, significant (p < 0.05) correlations were observed in 5 months in cases of trees with large damage (5 significant values) compared to healthy since 1961 (the beginning of meteorological measurements).In 2015, low increments were caused by extremely high temperatures (annual temperature 9.9 • C, long-term mean 7.9 • C), with low amounts of precipitation (for example in February 2 mm, long-term mean 21 mm).Moreover, PRPs established on HD plots were also found to be a negative pointer year in 2003.This year was the driest year recorded from 1961 through 2017 (304 mm, long-term mean 488 mm), especially with respect to the vegetation season (196 mm, long-term mean 325 mm).
Regarding the effects of precipitation and temperatures, significant (p < 0.05) correlations were observed in 5 months in cases of trees with large damage (5 significant values) compared to healthy trees (1 significant correlation).Diameter increments recorded between 1983 and 2017 positively correlated with the precipitation in June of the previous year (r = 0.30) and the current year (r = 0.40), while it was negatively correlated with precipitation in January of the current year (r = −0.38).Temperature had significant positive effects on radial growth in July of the previous year (r = 0.37) and in February of the current year (r = 0.29) for trees with high game damage.However, with respect to healthy trees, there was only one positive correlation of radial growth found with precipitation, specifically in April of the current year (r = 0.29).

Bark Stripping Damage and Stem Rot
A large difference was found between the bark stripping damage on forest stands in studied hunting districts.On LD plots, 16.6% of trees were damaged by game with a mean damage of 3.1% of stem circumference.The opposite situation was observed on HD plots, where a ratio of 92.7% with a mean damage of 31.2% of damaged stem girth was observed.Overall, stem rot was found in 62.1% of trees that were damaged by bark striping.Fresh bark stripping was found only in 7.3% of damaged stems.Stem rot was not detected in unharmed trees.According to dendrochronology, damages by bark stripping occurred on trees between 10-23 years old.Stem rot reached up to 17.6% (±14.7 SD) of the tree diameter in trees damaged by bark stripping.The analysis of spatial patterns of damaged trees (by Clark-Evans index) showed higher aggregations compared to spatial distributions of all trees on the PRPs.Significant differences were observed on all PRPs with respect to mean tree height (Kruskal-Wallis test, chi-squared = 55.2, p < 0.001), DBH (Kruskal-Wallis test, chi-squared = 42.7,p < 0.001), as well as in tree volume (Kruskal-Wallis test, chi-squared = 50.1,p < 0.001, see Figure 5) when comparing healthy trees (circumference damage ≤ 1/8 of stem girth) to those with less damage (circumference damage ≥ 1/8 and ≤ 1/3 of stem girth), and those with more damage (circumference damage ≥ 1/3 of stem girth).Mean DBH was significantly higher in healthy trees (19.6 cm) compared to trees with high damage (15.3 cm, p ≤ 0.001), with similar results found in mean height (19.8 m vs. 16.9 m), and tree volume (0.27 m 3 vs.0.14 m 3 ).
A generalized linear model (for more information see Material and Methods, Data Analyses section) was presented in the following form (simplified notation of model from R software): For results of the model, see Table 4.The goodness-of-fit was evaluated using the ratio of residual to null deviance (both values are present in respective table).
well as in tree volume (Kruskal-Wallis test, chi-squared = 50.1,p < 0.001, see Figure 5) when comparing healthy trees (circumference damage ≤ 1/8 of stem girth) to those with less damage (circumference damage ≥ 1/8 and ≤ 1/3 of stem girth), and those with more damage (circumference damage ≥ 1/3 of stem girth).Mean DBH was significantly higher in healthy trees (19.6 cm) compared to trees with high damage (15.3 cm, p ≤ 0.001), with similar results found in mean height (19.8 m vs. 16.9 m), and tree volume (0.27 m 3 vs.0.14 m 3 ).For a graphical representation of the model see Figure 6A.A rapid decrease in stem volume with respect to the presented surface plots, was illustrated when approximately one third of the stem circumference (measured at breast height, ca.130 cm) was damaged by bark stripping.The decrease was naturally more striking with increasing stem damage percentages.Healthy trees could, according to the model, grow up to 0.5 m 3 (or even more) of the stem volume for DBH at 30 cm, while these values could fluctuate to 50% or lower for seriously damaged trees.A similar decline was observed in trees with lower DBHs.
The differences between models for HD and LD plots are generally on very low level (Figure 6B) when applied to trees with DBH up to 25 cm.The overall model presented in Table 4 could be used for described conditions.the number of trees and spatial pattern of trees.Vertical diversity was observed to positively correlate with diameter and height differentiation.Total diversity negatively correlated with canopy closure and mean DBH.In terms of production, stand volume positively correlated with basal area, stocking, mean height, and DBH, while these parameters negatively correlated with bark stripping damage.The influence of game (LD plots vs. HD plots) showed low significance to mutual relationships among diversity indices compared to greater differences among production parameters.The LD plots with higher production parameters of timber occupy the left part of the diagram, while the high structural diversity and share of damage were typical for PRPs on HD plots.represented the number of trees and spatial pattern of trees.Vertical diversity was observed to positively correlate with diameter and height differentiation.Total diversity negatively correlated with canopy closure and mean DBH.In terms of production, stand volume positively correlated with basal area, stocking, mean height, and DBH, while these parameters negatively correlated with bark stripping damage.The influence of game (LD plots vs. HD plots) showed low significance to mutual relationships among diversity indices compared to greater differences among production parameters.The LD plots with higher production parameters of timber occupy the left part of the diagram, while the high structural diversity and share of damage were typical for PRPs on HD plots.

Discussion
Under European conditions, bark stripping is a common damage caused to forest stands, especially by red deer (Cervus elaphus Linnaeus) or sika deer [37,[70][71][72].In the last decades, the intensity of damage caused to forest stands has been on the rise [21,36,37], due to increasing abundance of the aforementioned species [37,43,74].Extensive damages to forest stands lead to significant loss in economic revenues [33], especially in regions where the Norway spruce is the main economic source with respect to forest management [75,76].The monetary value of stands damaged by bark stripping ranges from 70-95% of the hypothetical value of healthy stands [77].A positive relationship between bark stripping and population density of red deer was documented, e.g., throughout Europe [37,78,79], with respect to sika deer, with higher damage intensity found in Austria [74] or in native extension area (central Japan) [80].
In this study, the differences in growth parameters between stands in locations heavily affected by sika deer and forest stands with significantly lower sika deer activity were statistically supported.As a result of reduced growth of trees on forest sites with large observed damage was their overall stand volume of only 289.9 m 3 ha −1 (± 73 SD) compared to healthy stands with mean stand volume of 440.5 m 3 ha −1 (± 112 SD).The decrease in stand volume of damaged stands could be explained by their long-term damaging.The spreading of sika deer in studied territory happened after WWII, after destruction of game enclosures nearby where sika deer were reared [13,81].The negative effects of sika deer browsing, or bark stripping, on forest stands have also been described by other authors

Discussion
Under European conditions, bark stripping is a common damage caused to forest stands, especially by red deer (Cervus elaphus Linnaeus) or sika deer [37,[70][71][72].In the last decades, the intensity of damage caused to forest stands has been on the rise [21,36,37], due to increasing abundance of the aforementioned species [37,43,74].Extensive damages to forest stands lead to significant loss in economic revenues [33], especially in regions where the Norway spruce is the main economic source with respect to forest management [75,76].The monetary value of stands damaged by bark stripping ranges from 70-95% of the hypothetical value of healthy stands [77].A positive relationship between bark stripping and population density of red deer was documented, e.g., throughout Europe [37,78,79], with respect to sika deer, with higher damage intensity found in Austria [74] or in native extension area (central Japan) [80].
In this study, the differences in growth parameters between stands in locations heavily affected by sika deer and forest stands with significantly lower sika deer activity were statistically supported.As a result of reduced growth of trees on forest sites with large observed damage was their overall stand volume of only 289.9 m 3 ha −1 (± 73 SD) compared to healthy stands with mean stand volume of 440.5 m 3 ha −1 (± 112 SD).The decrease in stand volume of damaged stands could be explained by their long-term damaging.The spreading of sika deer in studied territory happened after WWII, after destruction of game enclosures nearby where sika deer were reared [13,81].The negative effects of sika deer browsing, or bark stripping, on forest stands have also been described by other authors [13,29,46,82].Detailed information comes mostly from Japan where bark stripping caused by sika deer affects population dynamics and species composition of native cool-temperate mixed forest.The large scale of bark stripping (more than 70% of circumference of the stems) leads to a higher tree mortality of damaged trees [83].Debarked trees are also susceptible to wind damage which results in wind-induced falling of coniferous trees through subsequent decay [80].However, to date, the production of forest stands has not been satisfactorily described.
Not only production, but also stand structure and diversity were influenced by game [84].In our study, higher structural diversity was observed on plots damaged by bark stripping compared to undamaged plots due to higher differentiation of slow growing damaged trees.In HD plots, higher crown differentiations were observed, as well as higher variability in crown closure and SDI between PRPs.Deer preference for trees with smaller DBHs was major contributing factor of higher canopy and stocking variability.Subsequently, the mortality of these heavily damaged trees by bark stripping accelerated the increase in canopy-gap ratio [21,80,85].The creation of larger openings increased mosaicity of forest stands on HD plots, which created unsuitable even-aged forest management planning because of lower timber production.Similar to the vertical structure, the horizontal structure of forest stands was largely influenced not only by silvicultural practices, but also by ungulates [86,87].As another pointer of health status and wood production, the relationship between crown diameter and DBH was analyzed as the structure of tree crown had a major effect on stemwood production [88].Heavily damaged trees show high variation in crown diameter at lower DBHs compared to healthy trees, pointing out the irregularities or deformations of crowns of heavily damaged trees.
Månsson and Jarnemo [36] documented, that stems of Norway spruce damaged by red deer had thinner bark and lower branching in the bottom portion of the stem compared to healthy trees.Most commonly, the stands with a mean DBH of ~8-16 cm were damaged by red deer [21], although Welch et al. [32] showed a preferred DBH of 20 cm with respect to sitka spruce.Such values of DBH correspond to values recorded in our study: the actual mean DBH in LD plots was 20.0 cm with ±3.4 SD, and in HD plots the actual mean DBH was 17.5 cm with ±1.4 SD.Regarding age, the most damaged stands were between 18 and 38 years of age [89] with detailed information about critical stand age for bark stripping mentioned by Koltzenburg [90] and Čermák et al. [91], who estimated this critical age to range between 15-30 years old.
Furthermore, the scent of disintegrated bark layers has been found to stimulate red deer to eat it [34].Repeated bark stripping occurs quite regularly on Norway spruce, as Büchsenmeister [92] documented, repeated bark stripping or different types of stem damage occur on 38% of previously damaged trees every year.
Previously published studies have primarily focused on describing bark stripping intensity and stem circumference damage [32,33,41].However, it is also crucial to understand how bark stripping impacts wood quality, wood color [93] and fungal infection, which can cause decomposition of wood mass and stem deformations [18,38,[93][94][95].Wood from previously damaged trees are often infected by Stereum sanguinolentum (Alb.& Schwein.)Fr. in the first year after bark damage occurrences.Only 50 cm 2 of debarked area on the stem is sufficient for penetration of fungal infection [38].The probability of developing infection is higher in cases of wider damage or when trees are damaged at early ages [38,93,95,96].The most suitable temperatures for infection development in Stereum sanguinolentum, range from −8.3 to 5.0 • C [93,97], which correspond to the mean temperatures in the Czech Republic during the non-vegetation season (1.7 • C according to Czech Hydrometeorological Institute (CHMI) [49]).Moreover, deer cause stem damages mainly in winter at time of high snow cover [98], when food availability is low [99], facilitating the possibility of infection.
In terms of climatic factors, the radial growth of spruce in lowland areas was influenced more by precipitation compared to temperature.Similarly, studies from other locations in the Czech Republic [100,101], Germany [102], Poland [103], and Slovakia [104] showed that temperature is the limiting factor of radial growth in mountainous regions, respectively the lack of precipitation in lowlands.Specifically, weather in June and July significantly impacts annual radial growth response [105,106], and is consistent with results reported in this study.Our data also revealed notable effects of climatic conditions in January, which could be due to frost.Temperatures in June and July positively affect growth of trees at higher altitudes [107,108].The limiting effects of low temperatures are more significant, not only at high-altitude sites, but also at northern latitudes, while the importance of precipitation generally increases in southern latitudes [109].
Radial growth can be significantly affected by a variety of other factors, including air pollution, fungal pathogens, and bark stripping [101,108,[110][111][112].Comparing radial growth of healthy trees and trees damaged by bark stripping in our study, annual ring width of sensitive damaged spruce was significantly impacted by climatic factors compared to trees without damages.In our case of study of healthy trees, there was only one significant value observed (positive correlation with precipitation in April), probably due to the location of the study area at an intermediate elevation of the spruce range, similar to sites evaluated in Germany [102].Damaged trees had reduced transpiration flow and were more stressed by abiotic factors compared to healthy trees without bark stripping.Growth of undamaged trees was positively correlated with precipitation in April only, likely due to this period coinciding with the beginning of growing season, therefore exhibiting the greatest growth dynamics and suitable climatic conditions.Similarly, in Switzerland, precipitation had a larger effect on the radial growth of damaged trees compared to health-dominant and health-suppressed trees [105].Tree-ring analyses from Finland showed a strong positive correlation with June precipitation in damaged stands compared to a low observed correlation in healthy stands [113].On our study site, June precipitation correlated with diameter increment and was much weaker in healthy stands compared to heavily damaged stands.
The sensitivity of Norway spruce forests damaged by bark striping is further enhanced by ongoing global climate change, including Czech Republic [106,[114][115][116]. Furthermore, the frequency and duration of extreme climate events (droughts, temperature backlashes, and windstorms) has recently increased [117,118].In combination with global warming and climatic stress, secondary damages to forest stands are often caused by bark beetles and fungal pathogens [119,120].Decreased forest vitality makes stands more susceptible also to windthrow damage [119,121].Specifically in the Czech Republic, the probability of drought between April and June increased by 50% since 1961 [122].Starting 1961, the warmest six years for the last 20 years were observed with an average temperature increase of 0.04 • C per year (the largest increase was recorded for July [49]).Climate change has a significant negative effect on timber production of spruce forests and causes severe economic losses [120,123,124], especially in vulnerable spruce stands previously damaged by bark striping.

Conclusions
This study has estimated significant differences in growth characteristics of Norway spruce with different levels of stem circumference damage caused by sika deer.The total stand volume on sites with lower damage intensity was significantly higher compared to forest stands that were heavily debarked by sika deer.Damaged trees were also more susceptible to fluctuations in precipitation and temperature, while radial increments of healthy trees were relatively stable.
Therefore, it is possible to conclude that damaged Norway spruce stands show higher vulnerability to extreme climatic conditions, which could reduce their stability, resulting in reduced resistance to biotic and abiotic threats.Given the higher vulnerability to climatic extremes highlighted in this study and the increased climatic extremes predicted with climate change, the impacts of bark stripping on forest stands is likely to be greater in the future.Damaged stands are not only vulnerable to wind, frost or heavy snow damage, but also to drought, fungal infections or insects, such as bark beetles.Due to the numerous threats of damaged stands (especially in the context of global climatic changes), it is necessary to pay attention to hunting management practices so game damage does not exceed an ecologically acceptable limit.Preserving ecologically and economically effective management tools of local forests require a significant reduction in the density of local sika

Figure 1 .
Figure 1.Localization of permanent research plots established in Norway spruce stands in the hunting districts of Horní Kozolupy (LD plots) and Rochlov (HD plots) in the Západočeská Upland area and monthly climatic characteristics (mean temperature and sum of precipitation between 1961 and 2017).Wider areas of interest are marked by (▲) symbols.The grey layer in the map symbolizes forest areas in the Czech Republic.

Figure 1 .
Figure 1.Localization of permanent research plots established in Norway spruce stands in the hunting districts of Horní Kozolupy (LD plots) and Rochlov (HD plots) in the Západočeská Upland area and monthly climatic characteristics (mean temperature and sum of precipitation between 1961 and 2017).Wider areas of interest are marked by ( ) symbols.The grey layer in the map symbolizes forest areas in the Czech Republic.

Figure 2 .
Figure 2. The relationship between diameter at breast height (DBH) and tree height of healthy spruce trees (parameter a = 1.00, and b = 0.18), trees with small damage (parameter a = 0.46, and b = 0.22), and trees with high game damage (parameter a = 1.09, and b = 0.18).

Figure 2 .
Figure 2. The relationship between diameter at breast height (DBH) and tree height of healthy spruce trees (parameter a = 1.00, and b = 0.18), trees with small damage (parameter a = 0.46, and b = 0.22), and trees with high game damage (parameter a = 1.09, and b = 0.18).

Figure 3 .
Figure 3.The relationship between crown diameter and DBH for healthy trees and trees with high damage.The bold lines represent linear regression functions (presented in relevant color).

Figure 3 .
Figure 3.The relationship between crown diameter and DBH for healthy trees and trees with high damage.The bold lines represent linear regression functions (presented in relevant color).

Figure 4 .
Figure 4. Standardized annual tree-ring increments of Norway spruce after removing the age trend (ring width index) and relevant smoothing splines (8-and 32-year spline).

Figure 4 .
Figure 4. Standardized annual tree-ring increments of Norway spruce after removing the age trend (ring width index) and relevant smoothing splines (8-and 32-year spline).The low radial increments that occurred in 1998 and 2015 qualified them as negative pointer years.The first half of 1998 had the warmest recorded temperatures (8.2 • C, long-term mean 6.1 • C)

Figure 5 .
Figure 5. Testing the differences in mean DBH, mean height, and mean stem volume between selected levels of stem girth damage.The Kruskal-Wallis test with relevant multiple comparisons was used; selected level of significance was α = 0.05.Indices above bars depict statistically significant differences between variants (different letters indicate significantly different variants).Error bars depict 95% confidence interval.A generalized linear model (for more information see Material and Methods, Data Analyses section) was presented in the following form (simplified notation of model from R software): stem volume ~ DBH + Circumference damage

Figure 5 .
Figure 5. Testing the differences in mean DBH, mean height, and mean stem volume between selected levels of stem girth damage.The Kruskal-Wallis test with relevant multiple comparisons was used; selected level of significance was α = 0.05.Indices above bars depict statistically significant differences between variants (different letters indicate significantly different variants).Error bars depict 95% confidence interval.

ForestsForests 21 Figure 6 .
Figure 6.Graphical representation of the generalized linear model for stem volume prediction based on DBH and stem circumference damage (A) and the differences between separate models for LD and HD plots (B).

3. 5 .
Relationships among Production, Diversity, and Game Damage Results of the PCA analysis are presented in the ordination diagram in Figure 7.The first ordination axis explained 43.2% of data variability, the first two axes combined represented 61.1% of variability, and the first four axes all together explained 82.7% of data variability.The first x-axis represented stand volume and basal area with damages caused by game.The second y-axis

Figure 6 .
Figure 6.Graphical representation of the generalized linear model for stem volume prediction based on DBH and stem circumference damage (A) and the differences between separate models for LD and HD plots (B).

3. 5 .
Relationships among Production, Diversity, and Game Damage Results of the PCA analysis are presented in the ordination diagram in Figure 7.The first ordination axis explained 43.2% of data variability, the first two axes combined represented 61.1% of variability, and the first four axes all together explained 82.7% of data variability.The first x-axis represented stand volume and basal area with damages caused by game.The second y-axis represented Forests 2019, 10, 320 13 of 21

Figure 7 .
Figure 7. Ordination diagram showing results of Principal Component Analysis (PCA) for the relationships among production parameters of stand (mean DBH, mean Height, Basal area, Volume, Stocking SDI, Canopy Crown projection area), diversity (A Arten-profile index, R aggregation index, TMd diameter differentiation index; TMh height differentiation index, B total diversity index), and game damage (Damage).Plot symbols: • plot with low game damage (LD plot), and ▼ plot with high game damage (HD plot).Numbers: 1-16-plot ID, and I.-IV.-locality.

Figure 7 .
Figure 7. Ordination diagram showing results of Principal Component Analysis (PCA) for the relationships among production parameters of stand (mean DBH, mean Height, Basal area, Volume, Stocking SDI, Canopy Crown projection area), diversity (A Arten-profile index, R aggregation index, TMd diameter differentiation index; TMh height differentiation index, B total diversity index), and game damage (Damage).Plot symbols: • plot with low game damage (LD plot), and plot with high game damage (HD plot).Numbers: 1-16-plot ID, and I.-IV.-locality.

Table 1 .
Overview of indices describing the biodiversity and their common interpretation.

Table 2 .
Stand characteristics of the tree layer on permanent research plots 1-16 in 2017 (LD-plots with low game damage; and HD-plots with high game damage).

Table 3 .
Diversity indices of tree layers on permanent research plots 1-16 in 2017 (LD-plots with low game damage; HD-plots with high game damage).

Table 3 .
Diversity indices of tree layers on permanent research plots 1-16 in 2017 (LD-plots with low game damage; HD-plots with high game damage).
Notes: R: aggregation index, A: Arten-Profile index, TM d : diameter differentiation index, TM h : height differentiation index, K: crown differentiation index, B: stand diversity index; * statistically significant (α = 0.05) result for horizontal structure (R index).

Table 4 .
Results of used generalized linear model.Provided p-values were related to the null hypotheses with the respective parameter equal to zero.