Estimation of Yield Loss Due to Deer Browsing in a Short Rotation Coppice Willow Plantation in Northern Japan

: Deer browsing is a major factor causing signiﬁcant declines in yield in short rotation coppice (SRC) willow, but the resultant yield loss is di ﬃ cult to estimate because it requires extensive investigation, especially when the standard yield is unknown. We investigated a simple method for estimating yield loss due to deer browsing. We enclosed an experimental SRC willow plantation in Hokkaido, northern Japan, planted with 12 clones, with an electric fence; deer browsing did, however, occur in the ﬁrst summer of the second harvest cycle. We counted the number of sprouting stems and deer-browsed stems per plant and, after three years, the yield of each clone was analyzed using a generalized linear model with the above two parameters for the numbers of stems as explanatory variables. The model explained the yield of 11 out of the 12 clones, and estimated that browsing of a single stem per plant could reduce yield to 80%. Losses due to deer browsing were estimated to be as much as 6.0 oven dry ton ha − 1 yr − 1 . The potential yield in the absence of deer browsing ranged from 2.2 to 7.5 oven dry ton ha − 1 yr − 1 among clones, and was signiﬁcantly positively correlated with the estimated yield loss due to deer browsing. Our results suggest that a generalized linear model can be used to estimate the yield loss due to deer browsing from a simple survey, and deer browsing could signiﬁcantly reduce willow biomass yield from the clones we studied, and thus countermeasures to control deer browsing are therefore necessary if su ﬃ cient willow biomass yield is to be produced.


Introduction
Global warming, which is mainly attributable to greenhouse gases such as CO 2 [1], has recently led to increased efforts to utilize woody biomass energy from managed forests, which are treated as carbon neutral, in Europe and the United States [2,3], as well as in Japan [4]. Short rotation coppice (SRC) willow is the most successful woody biomass energy crop in the cool temperate regions of Europe and North America [5][6][7][8][9][10]; under favorable conditions, it typically yields a biomass of over 10 oven dry ton (odt) ha −1 yr −1 , with a harvest rotation interval of 2-5 years. This capacity is associated with the coppicing ability of willow, which produces multiple stems and vigorous post-harvesting regrowth from the stump, resulting in repeated biomass yields of 20-50 odt ha −1 every 2-5 years over a Figure 1. Study site and experimental design. Location of the study site in Shimokawa town, Japan (a), planting design with double-row system (b), and arrangement of willow clones in the 28 rows (c). Clone names beginning with "P" (white rows) and "S" (gray rows) in panel (c) represent Salix petsusu and S. sachalinensis, respectively.

Deer Browsing Damage to Sprouting Willow Stems
In August 2017, the first year of sprouting, we found deer browsing damage to willow stems following a problem with the electric fence. Similar willow damage during growing season was previously reported to be due to browsing damage by Yezo-sika deer in a study using sensor cameras at the same study area [33]. In November 2017, we investigated the impact on the plants. An initial visual inspection suggested that browsing damage differed significantly among double rows containing different clones, and that damage to individual plants in a single-clone double-row was generally the same throughout the 60-m length (Figure 2a). Therefore, in each of the 28 double rows, the total numbers of stems and sprouting stems browsed by deer were counted for each plant Figure 1. Study site and experimental design. Location of the study site in Shimokawa town, Japan (a), planting design with double-row system (b), and arrangement of willow clones in the 28 rows (c). Clone names beginning with "P" (white rows) and "S" (gray rows) in panel (c) represent Salix pet-susu and S. sachalinensis, respectively.

Deer Browsing Damage to Sprouting Willow Stems
In August 2017, the first year of sprouting, we found deer browsing damage to willow stems following a problem with the electric fence. Similar willow damage during growing season was previously reported to be due to browsing damage by Yezo-sika deer in a study using sensor cameras at the same study area [33]. In November 2017, we investigated the impact on the plants. An initial visual inspection suggested that browsing damage differed significantly among double rows containing different clones, and that damage to individual plants in a single-clone double-row was generally Forests 2020, 11, 809 4 of 12 the same throughout the 60-m length ( Figure 2a). Therefore, in each of the 28 double rows, the total numbers of stems and sprouting stems browsed by deer were counted for each plant positioned 5 m from the edge, that is, typically 20 plants per double row, with the exception of some rows that included dead plants (we did not find any plants that had died due to deer damage at the time of the browsing investigation). The number of plants investigated for each clone ranged from 34 to 108. The percentage of browsed stems per plant (%) was calculated by dividing the number of browsed stems by the total number of stems per plant. With the exception of August 2017 and the snow period, electric fencing protection from deer was functional, and there was no evidence of extensive deer browsing in other periods during the study.  Planting rows of a clone with minimal browsing (white inverted triangles) adjacent to rows of a severely browsed clone (black inverted triangles) (a). Willow stems broken by Yezo-sika deer while browsing willow stem tips (b).

Willow Yield Measurement
In late October 2019, three years after the first harvest, we harvested all the willow plants, except for those that were sampled for fresh mass to dry mass conversion. We did this using a sugarcane harvester (UT-120K, Uotani-tekko Inc., Nara, Japan), which cuts willows into billets approximately 20 cm long and harvests them into a harvesting bag. We harvested each 60-m single-clone double row into a single bag, which was then immediately weighed using a digital crane scale (resolution: 0.5 kg; 1ACBP-K, Shuzui Scales Co., Ltd., Aichi, Japan). The fresh yield (fresh t ha −1 yr −1 ) was determined by subtracting the weight of the harvesting bag (11.5 kg) and dividing by the planted area in a double row (0.011 ha) and the rotation period of three years. To obtain the oven-dry mass yield from the fresh yield in the bags, we harvested nine typical willow plants using a saw before harvesting with the sugarcane harvester and weighed each plant (1.63 fresh kg per plant on average, 14.63 fresh kg in total). We then transported the sample plants to the laboratory, dried them at 70 °C to constant mass, and reweighed them. The oven-dry yield of each clone (odt ha −1 yr −1 ) was calculated by multiplying the fresh yield mean of each of the 12 clones by the mean of the ratio of dry to fresh weight of the sample plants.

Data Analyses
All statistical analyses were performed using R version 3.6.2 [34]. We estimated willow yield and yield loss due to deer browsing from generalized linear models (GLMs) with a gamma distribution and log link function. The actual browsed yield of each clone was analyzed by a GLM in which the objective variable was the mean dry yield for each clone, and the explanatory variables were, for each clone, the mean number of stems and mean number of browsed stems per plant. We also analyzed the same variables with GLMs with a gamma distribution and inverse link function and a Gaussian distribution and identity link function. We compared the three models using Akaike's information criterion (AIC) (Table S1); the GLM with a gamma distribution and log link function had the smallest AIC, and so was adopted as the willow yield estimation model. We determined that a clone (P.I.82) was an outlier from the residual plot, the quantile-quantile plot, and scale-location plot using the plot function in R ( Figure S1). We reanalyzed the GLM for the dataset excluding the outlier clone, and estimated the GLM coefficients. The coefficient of determinations for each of the GLMs (r 2 GLM) [35] for all 12 clones and for the 11 clones excluding the outlier was calculated using the req function of the "rsq" package in R. The GLM excluding an outlier clone was used to simulate willow

Willow Yield Measurement
In late October 2019, three years after the first harvest, we harvested all the willow plants, except for those that were sampled for fresh mass to dry mass conversion. We did this using a sugarcane harvester (UT-120K, Uotani-tekko Inc., Nara, Japan), which cuts willows into billets approximately 20 cm long and harvests them into a harvesting bag. We harvested each 60-m single-clone double row into a single bag, which was then immediately weighed using a digital crane scale (resolution: 0.5 kg; 1ACBP-K, Shuzui Scales Co., Ltd., Aichi, Japan). The fresh yield (fresh t ha −1 yr −1 ) was determined by subtracting the weight of the harvesting bag (11.5 kg) and dividing by the planted area in a double row (0.011 ha) and the rotation period of three years. To obtain the oven-dry mass yield from the fresh yield in the bags, we harvested nine typical willow plants using a saw before harvesting with the sugarcane harvester and weighed each plant (1.63 fresh kg per plant on average, 14.63 fresh kg in total). We then transported the sample plants to the laboratory, dried them at 70 • C to constant mass, and reweighed them. The oven-dry yield of each clone (odt ha −1 yr −1 ) was calculated by multiplying the fresh yield mean of each of the 12 clones by the mean of the ratio of dry to fresh weight of the sample plants.

Data Analyses
All statistical analyses were performed using R version 3.6.2 [34]. We estimated willow yield and yield loss due to deer browsing from generalized linear models (GLMs) with a gamma distribution and log link function. The actual browsed yield of each clone was analyzed by a GLM in which the objective variable was the mean dry yield for each clone, and the explanatory variables were, for each clone, the mean number of stems and mean number of browsed stems per plant. We also analyzed the same variables with GLMs with a gamma distribution and inverse link function and a Gaussian distribution and identity link function. We compared the three models using Akaike's information criterion (AIC) (Table S1); the GLM with a gamma distribution and log link function had the smallest AIC, and so was adopted as the willow yield estimation model. We determined that a clone (P.I.82) was an outlier from the residual plot, the quantile-quantile plot, and scale-location plot using the plot function in R ( Figure S1). We reanalyzed the GLM for the dataset excluding the outlier clone, and estimated the GLM coefficients. The coefficient of determinations for each of the GLMs (r 2 GLM ) [35] for all 12 clones and for the 11 clones excluding the outlier was calculated using the req function of the "rsq" package in R. The GLM excluding an outlier clone was used to simulate willow yield and yield loss when the number of stems per plant varied from 1 to 10, and the percentage of browsed stems per plant varied from 0% to 100%. In addition, we estimated the non-browsed yield by converting the number of browsed stems in the GLM to zero for each clone in the experimental plot. Yield loss by browsing was estimated by subtracting the estimated browsed yield from the estimated unbrowsed yield, that is, 0% of browsed stems per plant.
Standard major axis (SMA) regressions were performed to analyze bivariate relationships between willow yield and other variables across clones using the "smart" package in R. Interspecific differences between S. pet-susu and S. sachalinensis in the number of stems and number of browsed stems and the percentage of browsed stems per plant were analyzed by a Mann-Whitney U test using the wilcox.test function in R, and the differences in actual yield, estimated non-browsed yield, and estimated yield loss from browsing were analyzed by a Student t-test using the t.test function in R. The level of statistical significance was set at p = 0.05.

Deer Browsing and Willow Yield across Clones
Yezo-sika deer typically browsed the tips of willow stems ( Figure 2a). In addition, deer broke some stems in order to access the tips (Figure 2b). The number of sprouting stems, the extent of deer browsing, and actual willow yield varied largely among clones. Clonal means for the number of stems per plant ranged from 4.8 for P.I.  yield and yield loss when the number of stems per plant varied from 1 to 10, and the percentage of browsed stems per plant varied from 0% to 100%. In addition, we estimated the non-browsed yield by converting the number of browsed stems in the GLM to zero for each clone in the experimental plot. Yield loss by browsing was estimated by subtracting the estimated browsed yield from the estimated unbrowsed yield, that is, 0% of browsed stems per plant. Standard major axis (SMA) regressions were performed to analyze bivariate relationships between willow yield and other variables across clones using the "smart" package in R. Interspecific differences between S. pet-susu and S. sachalinensis in the number of stems and number of browsed stems and the percentage of browsed stems per plant were analyzed by a Mann-Whitney U test using the wilcox.test function in R, and the differences in actual yield, estimated non-browsed yield, and estimated yield loss from browsing were analyzed by a Student t-test using the t.test function in R. The level of statistical significance was set at p = 0.05.

Deer Browsing and Willow Yield across Clones
Yezo-sika deer typically browsed the tips of willow stems (Figure 2a). In addition, deer broke some stems in order to access the tips (

Estimation of Yield Loss Caused by Deer Browsing
A simple GLM incorporating the number of stems and number of browsed stems could accurately estimate actual willow yields with deer browsing (Table 1, Figure 4a). In the GLM calculated for all 12 clones including an outlier, r 2 GLM was 0.43, indicating a moderate predictive power (Table 1). When we removed the outlier clone from the analysis (Figure S1), r 2 GLM increased to

Estimation of Yield Loss Caused by Deer Browsing
A simple GLM incorporating the number of stems and number of browsed stems could accurately estimate actual willow yields with deer browsing (Table 1, Figure 4a). In the GLM calculated for all 12 clones including an outlier, r 2 GLM was 0.43, indicating a moderate predictive power (Table 1). When we removed the outlier clone from the analysis ( Figure S1), r 2 GLM increased to 0.74, indicating substantial predictive power (Table 1). In addition, for the model excluding an outlier, the r 2 value in the SMA regression between the actual and estimated yield was 0.74 (p < 0.001), and the regression line almost overlapped the 1:1 relationship line (Figure 4a). Yield was underestimated for the outlier clone P.I.82. Clone P.I.82, which had the highest actual yield among the willow clones studied, was estimated by the GLM to have a predicted yield of more than 1 odt ha −1 yr −1 lower than actual yield. The numbers of stems and browsed stems per plant had a significantly positive (approximately 1.5-fold) and negative (approximately 0.8-fold) effect on yield for each clone in both models, including and excluding an outlier clone ( Table 1). The GLM simulation excluding an outlier clone predicted that yields increased with increasing number of stems per plant, but decreased significantly with the percentage of browsed stems per plant (Figure 4b). The effect of browsing damage on the yield was large; for example, even if 25% of the number of stems per plant were browsed, yield decreased to at least 61% of the non-browsed yield (Figure 4b). The GLM simulation also predicted that yield loss by deer browsing increases with the estimated yield and percentage of browsed stems per plant. For example, with a standard target yield of 10 odt ha −1 yr −1 for SRC willow, yield losses caused by deer browsing were estimated at approximately 4, 6, 7, and 8 odt ha −1 yr −1 if 25%, 50%, 75%, and 100% of the stems per plant were browsed (Figure 4c). Table 1. A summary of the generalized linear models (GLMs) of estimated willow yield under deer browsing conditions for each clone. A gamma distribution with a log-link function was used for the models. One model includes all 12 clones studied, and another model comprised 11 clones, excluding an outlier. The models incorporate the number of stems and browsed stems per plant as explanatory variables. The numbers in parentheses in the "Estimate" column represent the exponential value of the estimate. SE is the standard error. When the GLM was applied to actual data from the experimental plantation, willow yield was estimated to increase to 2.2-7.5 odt ha −1 yr −1 , assuming an absence of deer browsing, from 1.1-4.3 odt ha −1 yr −1 with deer browsing (Figure 5a). The estimated yield without deer browsing was not significantly correlated with actual yield across clones. The yield loss due to deer browsing was estimated to range from 0.1 odt ha −1 yr −1 for S.T.27 to 6.0 odt ha −1 yr −1 for S.I.27, and was significantly correlated with estimated yield under non-browsing conditions across clones (Figure 5b) and the percentage of browsed stems per plant (Figure 5c). There was no significant interspecific difference in the estimated yield without browsing (p = 0.27) and estimated yield loss due to browsing (p = 0.15). In P.I.82, with the maximum actual yield among the clones studied, although an average of 1.7 stems (i.e., 30% of the average 6.1 sprouting stems per individual) were subject to browsing, the potential yield without browsing was estimated to be approximately 4 odt ha −1 yr −1 , which was almost the same as the actual yield with browsing ( Figure 5a).

Factor
Estimate SE t Value p Value GLM for 12 clones (including an outlier)    When the GLM was applied to actual data from the experimental plantation, willow yield was estimated to increase to 2.2-7.5 odt ha −1 yr −1 , assuming an absence of deer browsing, from 1.1-4.3 odt ha −1 yr −1 with deer browsing (Figure 5a). The estimated yield without deer browsing was not significantly correlated with actual yield across clones. The yield loss due to deer browsing was estimated to range from 0.1 odt ha −1 yr −1 for S.T.27 to 6.0 odt ha −1 yr −1 for S.I.27, and was significantly correlated with estimated yield under non-browsing conditions across clones (Figure 5b) and the percentage of browsed stems per plant (Figure 5c). There was no significant interspecific difference in the estimated yield without browsing (p = 0.27) and estimated yield loss due to browsing (p = 0.15). In P.I.82, with the maximum actual yield among the clones studied, although an average of 1.7 stems (i.e., 30% of the average 6.1 sprouting stems per individual) were subject to browsing, the potential yield without browsing was estimated to be approximately 4 odt ha −1 yr −1 , which was almost the same as the actual yield with browsing ( Figure 5a).

Yield Loss Caused by Deer Browsing
Although the actual willow yield with browsing could not be estimated by a single stem parameter (Figure 3), we were able to estimate the yield for 11 out of 12 clones using a simple GLM that incorporated both the number of sprouting stems and browsed stems as explanatory variables (Figure 4a). The high predictive accuracy of this model, despite the lack of stem diameter and/or height data, suggests that growth in each browsed or unbrowsed stem was comparable among the 11 clones during three years of rotation. As a result, it was possible to predict the yield only by the data on the numbers of stems. These results suggest that when an SRC willow plantation is damaged by deer browsing, the biomass loss can be easily estimated using a GLM by counting the number of browsed and unbrowsed stems per plant and measuring the yield in places where the degree of

Yield Loss Caused by Deer Browsing
Although the actual willow yield with browsing could not be estimated by a single stem parameter (Figure 3), we were able to estimate the yield for 11 out of 12 clones using a simple GLM that incorporated both the number of sprouting stems and browsed stems as explanatory variables (Figure 4a). The high predictive accuracy of this model, despite the lack of stem diameter and/or height data, suggests that growth in each browsed or unbrowsed stem was comparable among the 11 clones during three years Forests 2020, 11, 809 8 of 12 of rotation. As a result, it was possible to predict the yield only by the data on the numbers of stems. These results suggest that when an SRC willow plantation is damaged by deer browsing, the biomass loss can be easily estimated using a GLM by counting the number of browsed and unbrowsed stems per plant and measuring the yield in places where the degree of browsing damage is different in the plantation. This approach will be particularly useful in areas where SRC willow is in the developmental stage and standard yields are unknown. In addition, the fact that the number of browsed stems in the first year of rotation significantly explained yield in the final year of rotation in the GLM (Table 1) would indicate that browsing damage, which occurred only once in the first year of rotation, affected biomass increase during the three-year harvest intervals and was not offset during that period. Although the influence of deer on the SRC willow has sometimes been underestimated [21], these results indicate that deer-browsing control can be very important for the success of the SRC willow plantation. On the other hand, yield of P.I.82, which had the highest actual yield among the 12 clones studied, was greatly underestimated by the GLM (Figure 4a). This could possibly result from greater biomass per stem and/or a smaller reduction in biomass per browsed stem in this clone than in the other clones. The stem of P.I.82 might be taller than that of the other clones at the time of browsing damage and thus be less affected by browsing damage per stem, although we did not measure the height.
The model simulation showed that a single browsing event (which occurred during the first summer of our study period) resulted in a significant yield loss over three years of a harvest cycle; the potential for willow yield to be halved if 25% of the number of stems per plant that were browsed by deer was also shown (Figures 4 and 5). This high percentage of yield loss, even when the percentage of browsed stems per plant was small, is attributable not only to direct biomass loss due to deer browsing, but also to the indirect negative effects of browsing damage on biomass increase. In a conifer crop tree, Pseudotsuga menziesii, deer browsing and weed overgrowth inhibited the crop tree height growth, whereas deer browsing on broadleaf competitors resulted in a slightly higher tree height in comparison with that in the plot which excluded deer under high weed suppression by herbicides, indicating potential interaction between competing vegetation and browsing [36]. Because willow growth is very susceptible to weeds [24,37], and willow had a higher preference for deer than weeds at the present study site [33], the reason for the indirect negative effect of deer browsing on willow biomass in this study could be that deer browsing damage delayed the closure of the willow canopy, resulting in weed overgrowth, thus suppressing the growth of the willow susceptible to weeds [24,38]. Because our results are based on only one browsing event that occurred on a small plantation, future studies need to examine how the results of the simulation vary depending on various conditions and situations of SRC willow plantation, such as the extent of weed growth, plantation size, browsing frequency, browsing timing, etc.

Potential Yield of Studied Willow Clones
Only two of the 12 clones studied have a known potential yield in the absence of browsing damage in the study area. The two clones, S.I.27 and S.I.67, yielded approximately 5-11 odt ha −1 yr −1 in the first harvest rotation (average of approximately 8 odt ha −1 yr −1 ) with no deer browsing at the same plantation in this study [24]. We estimated the potential yield to be 7.5 and 7.3 odt ha −1 yr −1 for S.I.27 and S.I.67, respectively, in the second rotation using the GLM (Figure 5a). SRC willow yields are generally greater in the second than in the first rotation [39][40][41], but, in our study, more yield data from the second rotation for these clones were included for areas not weeded between mulches and therefore with potentially reduced yields compared with the first rotation [24]; therefore, the potential GLM yield estimates for these two clones would be reasonable in comparison to the first rotation yield. It might therefore be reasonable to assume that the estimated potential yields by the GLM of the other nine clones other than the outlier are also probably close to their actual potential yield in the absence of browsing, although we were not able to compare the estimated potential yield with the first rotation yield due to the lack of the first rotation data. Future verification with actual potential yield that was Forests 2020, 11, 809 9 of 12 reliably protected against deer browsing will allow us to clarify the accuracy of the estimated potential yield from the model in this study.
In this study, clones with the highest estimated yield tended to show results of greater yield loss due to deer browsing (Figure 5b). In general, there is a trade-off between plant growth and defense against herbivory because of limited resources [42,43]. Therefore, it would be difficult and time consuming to search for and breed super willow clones that achieve both high-yielding and low damage from deer browsing; though susceptibility to browsing animals has already been used as one of the traits for selecting willow commercial variety in Europe [29]. Compared with Europe and North America, hunting is less popular in Japan, where hunters are aging and decreasing in number [44]. Therefore, it will be potentially difficult to use deer hunting as a major management tool for maintaining deer browsing at low levels in SRC willow plantations in Hokkaido, Japan. The potential yields of most clones were estimated to be low (<5 odt ha −1 yr −1 ), even considering that the overgrowth of weeds among mulches could have negatively affected willow yields. If commercial SRC willow plantations are to be promoted in Hokkaido, Japan, the selection of high-yielding clones that are less susceptible to deer browsing, together with the development of inexpensive measures protecting against deer entry to plantations, will therefore be desirable. Further research about various protective measures against deer browsing, including the introduction of physical barriers such as fencing and/or repellents, establishment of large plantations, and deer population management in cooperation with local hunters [20,21,45,46], should lead to a higher willow yield by reducing yield loss due to deer browsing and improved profit forecasts at SRC willow plantations for the studied willow clones in the study site area.

Conclusions
This study demonstrated that the GLM analysis, which incorporates two easily measured parameters-the number of stems per plant and the number of browsed stems by deer-can estimate the biomass loss due to deer browsing in an SRC willow plantation, which cannot be estimated by a single parameter. The GLM also can be used to predict yield in the absence of browsing damage, although it needs further validation in the future. In the 12 clones studied, there was a trade-off between potential biomass yield and susceptibility to browsing damage, suggesting that the importance of preventing deer browsing to obtain adequate biomass yield in an SRC willow plantation.