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

**Hisanori Harayama 1,\* , Qingmin Han <sup>2</sup> , Makoto Ishihara <sup>1</sup> , Mitsutoshi Kitao <sup>1</sup> , Akira Uemura <sup>2</sup> , Shozo Sasaki <sup>1</sup> , Takeshi Yamada 1,**† **, Hajime Utsugi <sup>2</sup> and Yutaka Maruyama <sup>3</sup>**


Received: 18 June 2020; Accepted: 24 July 2020; Published: 26 July 2020

**Abstract:** Deer browsing is a major factor causing significant declines in yield in short rotation coppice (SRC) willow, but the resultant yield loss is difficult 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 first 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−<sup>1</sup> yr−<sup>1</sup> . The potential yield in the absence of deer browsing ranged from 2.2 to 7.5 oven dry ton ha−<sup>1</sup> yr−<sup>1</sup> among clones, and was significantly 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 significantly reduce willow biomass yield from the clones we studied, and thus countermeasures to control deer browsing are therefore necessary if sufficient willow biomass yield is to be produced.

**Keywords:** deer browsing; *Salix*; short rotation coppice; woody biomass; yield loss

### **1. Introduction**

Global warming, which is mainly attributable to greenhouse gases such as CO<sup>2</sup> [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–10]; under favorable conditions, it typically yields a biomass of over 10 oven dry ton (odt) ha−<sup>1</sup> yr−<sup>1</sup> , 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−<sup>1</sup> every 2–5 years over a

plantation life span of 20–25 years, with no need for replanting [11,12]. SRC willow can be established on marginal land that is generally unsuitable for food crops [13–15]. In addition to producing biomass for use in the bioenergy and biofuel industries, SRC willow provides positive ecosystem services, such as carbon sequestration in the roots, biodiversity, and water quality [16,17]. However, there are countries or regions where SRC willow has not become established commercially, such as Japan, despite the existence of suitable environmental conditions. This is because of risks and uncertainties concerning its production, management, and marketing, which need to be resolved if the use of SRC willow is to be promoted [18].

One cause of significant biomass loss in SRC willow plantations is the extensive damage due to browsing by large animals, such as deer [19–21]. Smaller SRC willow plantations are more susceptible to deer browsing, which is therefore a greater obstacle to commercial viability, especially in the early stages of SRC willow introduction, when plantations tend to be small [19]. The construction of fences as animal deterrents is effective but expensive [12,20]; in order to evaluate its cost effectiveness, it is therefore important to be able to estimate the yield loss due to deer browsing. In areas where commercial SRC willow has been widely developed, such as in northern Europe, potential yield by region for each willow variety is known; yield loss due to deer browsing can therefore be easily estimated by comparing the browsed actual yield with potential yield in the absence of browsing. In contrast, in areas where SRC willow is underdeveloped and potential yield is not known, estimating yield loss from deer browsing is difficult and requires more extensive investigation by, for example, individually estimating shoot losses based on shoot bite diameters using an allometry equation between shoot diameter and biomass [22].

The use of woody biomass for energy production has progressed in Japan since the Great East Japan Earthquake in 2011 [4,23], but SRC willow is still not commercially cultivated. Hokkaido, the northernmost of Japan's four main islands, is potentially suitable for growing SRC willow because of its flat terrain and cool temperate climate with sufficient rainfall. Additionally, its cold winter temperatures increase the demand for local heat and/or combined heat and power plants, which are typical destinations of harvested SRC willow [19]. In fact, in experimental plantations in Hokkaido, willow yields of more than 10 odt ha−<sup>1</sup> yr−<sup>1</sup> can be obtained with proper cultivation management [24,25]. Meanwhile, Yezo-sika deer (*Cervus nippon yesoensis*) density has increased rapidly in Hokkaido since the late 1980s [26], and the species is currently considered overabundant [27]. Since the late 1990s, Yezo-sika deer damage to agriculture and forestry has been severe [26], with recent economic impacts estimated at approximately 4 billion yen per year [28]. In general, susceptibility to deer damage varies among willow clones and is one of the key properties considered in selecting cultivars in Europe and the United States [29–31]; the selection of appropriate willow cultivars is still under development in Japan, and thus there is limited information on the susceptibility of willow clones to deer damage as well as their potential yield.

Our main objectives were to estimate: (1) the amount of yield loss due to deer browsing in willow clones with unknown potential yields, using a simple method counting the number of browsed and non-browsed stems; and (2) the yield of each clone in the absence of deer browsing—that is, potential yield.

#### **2. Materials and Methods**

#### *2.1. Study Site*

The study was conducted at an experimental SRC willow plantation in Shimokawa town, Hokkaido prefecture, northern Japan (44◦250 N, 142◦420 E, 249 m a.s.l., Figure 1a). This plantation was established to modify and improve SRC willow management procedures [24] and to select high yield clones. The parent material of soil in the study site was quaternary alluvial sediments, with rounded and sub-rounded fragments covering the surface [32]. The study site was 56 m × 60 m, almost level, and was plowed and harrowed in late October 2012. In May 2013, one-year old

cuttings (length, approximately 20 cm) of *Salix pet-susu* Kimura and *S. sachalinensis* F.Schmidt clones (*n* = 7 and 5, respectively) were planted by hand in a double-row system at a density of 20,000 cuttings ha−<sup>1</sup> ; the distance between and within double rows was 1.5 and 0.5 m, respectively, and plant spacing within rows was 0.5 m (Figure 1b). The cuttings were collected from experimental plantations in Shimokawa town and Sapporo city (Figure 1a) a few days before planting and refrigerated in plastic bags. Since willow clones appropriate for SRC were still under selection and no commercial cuttings existed in Japan, there was no guarantee the cuttings used in this study were high-yielding. All planting rows were sheeted with a plastic mulch (width, 1.5 m, and thickness, 0.021 mm; Sunshat, C.I. TAKIRON Corporation, Tokyo, Japan) to suppress weed growth. However, we did not conduct additional weeding between adjacent plastic mulches, except for half of the S.I.27, S.I.67, P.I.81, and P.I.82 clones, resulting in lower yields due to weed overgrowth in the first harvest rotation [24]. A total of 28 double rows (approximately 60 m long) were planted, each comprising plants of a single clone; the rows were randomly arranged (Figure 1c). Almost all the willow plants were cut back in November 2013 to promote sprouting, except for half of the S.I.27, S.I.67, P.I.81, and P.I.82 clones, so that the impact of the cut back procedure on biomass yield could be evaluated [24]. All plants were harvested in November 2016 using chain and hand saws, and then coppiced in spring 2017. The study site was enclosed by an electric fence, except in the snow period (November to April) to avoid destruction of the fence by snow coverage, usually over 1 m. *Forests* **2020**, *11*, x FOR PEER REVIEW 3 of 12 respectively) were planted by hand in a double-row system at a density of 20,000 cuttings ha−1; the distance between and within double rows was 1.5 and 0.5 m, respectively, and plant spacing within rows was 0.5 m (Figure 1b). The cuttings were collected from experimental plantations in Shimokawa town and Sapporo city (Figure 1a) a few days before planting and refrigerated in plastic bags. Since willow clones appropriate for SRC were still under selection and no commercial cuttings existed in Japan, there was no guarantee the cuttings used in this study were high-yielding. All planting rows were sheeted with a plastic mulch (width, 1.5 m, and thickness, 0.021 mm; Sunshat, C.I. TAKIRON Corporation, Tokyo, Japan) to suppress weed growth. However, we did not conduct additional weeding between adjacent plastic mulches, except for half of the S.I.27, S.I.67, P.I.81, and P.I.82 clones, resulting in lower yields due to weed overgrowth in the first harvest rotation [24]. A total of 28 double rows (approximately 60 m long) were planted, each comprising plants of a single clone; the rows were randomly arranged (Figure 1c). Almost all the willow plants were cut back in November 2013 to promote sprouting, except for half of the S.I.27, S.I.67, P.I.81, and P.I.82 clones, so that the impact of the cut back procedure on biomass yield could be evaluated [24]. All plants were harvested in November 2016 using chain and hand saws, and then coppiced in spring 2017. The study site was enclosed by an electric fence, except in the snow period (November to April) to avoid destruction of the fence by snow coverage, usually over 1 m.

**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-***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.

#### *2.2. Deer Browsing Damage to Sprouting Willow Stems*

*susu* and *S. sachalinensis*, respectively.

*2.2. 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 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

containing different clones, and that damage to individual plants in a single-clone double-row was

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. *Forests* **2020**, *11*, x FOR PEER REVIEW 4 of 12 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.

**Figure 2.** Willow stem tips browsed by Yezo-sika deer in the experimental short rotation coppice. 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**). **Figure 2.** Willow stem tips browsed by Yezo-sika deer in the experimental short rotation coppice. 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**).

#### *2.3. Willow Yield Measurement 2.3. 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 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−<sup>1</sup> yr−<sup>1</sup> ) 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−<sup>1</sup> yr−<sup>1</sup> ) 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.

#### weight of the sample plants. *2.4. Data Analyses*

(*r*

*2.4. 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 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,

using the *plot* function in R (Figure S1). We reanalyzed the GLM for the dataset excluding the outlier

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. *Forests* **2020**, *11*, x FOR PEER REVIEW 5 of 12 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

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. 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.

#### **3. Results 3. Results**

#### *3.1. Deer Browsing and Willow Yield across Clones 3.1. Deer Browsing and Willow Yield across Clones*

*3.2. Estimation of Yield Loss Caused by Deer Browsing*

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.81 to 8.1 for S.I.27, across clones (Figure 3a). The clonal means of the number of stems browsed by deer per plant ranged from 0.61 for S.T.27 to 8.1 for S.I.27 across clones (Figure 3b). The clonal means of the percentage of browsed stems per plant ranged from 11% for P.C.51 to 99% for S.I.27, across clones (Figure 3c). The actual yield of each clone ranged from 1.05 odt ha−<sup>1</sup> yr−<sup>1</sup> for clone S.T.27 to 4.26 odt ha−<sup>1</sup> yr−<sup>1</sup> for P.I.82 (Figure 3). There was no significant relationship observed between actual yield and the number of stems per plant, number of browsed stems per plant, and percentage of browsed stems per plant across clones (Figure 3). There were no significant interspecific differences in the number of stems per plant (Mann–Whitney U test; *p*=0.53), browsed stems per plant (*p* = 0.53), percentage of browsed stems per plant (*p* = 0.34), and actual yield (t-test; *p* = 0.08). 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.81 to 8.1 for S.I.27, across clones (Figure 3a). The clonal means of the number of stems browsed by deer per plant ranged from 0.61 for S.T.27 to 8.1 for S.I.27 across clones (Figure 3b). The clonal means of the percentage of browsed stems per plant ranged from 11% for P.C.51 to 99% for S.I.27, across clones (Figure 3c). The actual yield of each clone ranged from 1.05 odt ha−1 yr−1 for clone S.T.27 to 4.26 odt ha−1 yr−1 for P.I.82 (Figure 3). There was no significant relationship observed between actual yield and the number of stems per plant, number of browsed stems per plant, and percentage of browsed stems per plant across clones (Figure 3). There were no significant interspecific differences in the number of stems per plant (Mann–Whitney U test; *p* = 0.53), browsed stems per plant (*p* = 0.53), percentage of browsed stems per plant (*p* = 0.34), and actual yield (t-test; *p* = 0.08).

**Figure 3.** Relationships of the means of number of stems per plant (**a**), number of browsed stems per plant (**b**), and percentage of browsed stems per plant (**c**) to actual willow yield across 12 willow clones. Clone names, coefficients of determination (*r* 2 ), and level of significance (*p*) of the standardized major axis regression are shown in each panel. Circles represent *Salix pet-susu* clones and triangles *S. sachalinensis* clones; odt, oven dry ton. Error bars represent standard errors of the mean (*n* = 34–108). **Figure 3.** Relationships of the means of number of stems per plant (**a**), number of browsed stems per plant (**b**), and percentage of browsed stems per plant (**c**) to actual willow yield across 12 willow clones. Clone names, coefficients of determination (*r* 2 ), and level of significance (*p*) of the standardized major axis regression are shown in each panel. Circles represent *Salix pet-susu* clones and triangles *S. sachalinensis* clones; odt, oven dry ton. Error bars represent standard errors of the mean (*n* = 34–108).

A simple GLM incorporating the number of stems and number of browsed stems could

2GLM increased to

### *3.2. 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−<sup>1</sup> yr−<sup>1</sup> 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−<sup>1</sup> yr−<sup>1</sup> for SRC willow, yield losses caused by deer browsing were estimated at approximately 4, 6, 7, and 8 odt ha−<sup>1</sup> yr−<sup>1</sup> 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.


<sup>1</sup> Residual deviance: 0.891 on 9 degrees of freedom, *r* 2 GLM = 0.43; <sup>2</sup> Residual deviance: 0.619 on 8 degrees of freedom, *r* 2 GLM = 0.74.

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−<sup>1</sup> yr−<sup>1</sup> , assuming an absence of deer browsing, from 1.1–4.3 odt ha−<sup>1</sup> yr−<sup>1</sup> 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−<sup>1</sup> yr−<sup>1</sup> for S.T.27 to 6.0 odt ha−<sup>1</sup> yr−<sup>1</sup> 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

the estimate. SE is the standard error.

the *r*

(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−<sup>1</sup> yr−<sup>1</sup> , which was almost the same as the actual yield with browsing (Figure 5a). Number of browsed stems per plant −0.199 (0.82) 0.048 −4.186 0.003 <sup>1</sup> Residual deviance: 0.891 on 9 degrees of freedom, *r* 2GLM = 0.43; <sup>2</sup> Residual deviance: 0.619 on 8 degrees of freedom, *r* 2GLM = 0.74.

**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

GLM for 12 clones (including an outlier) <sup>1</sup>

Number of stems per plant 0.411 (1.51) 0.118 3.488 0.007 Number of browsed stems per plant −0.223 (0.80) 0.055 −4.037 0.003 GLM for 9 clones (excluding an outlier) <sup>2</sup>

Number of stems per plant 0.373 (1.45) 0.100 3.724 0.006

**Factor Estimate SE** *t* **Value** *p* **Value**

Intercept −1.113 (0.33) 0.616 −1.808 0.104

Intercept −1.009 (0.36) 0.519 −1.943 0.088

*Forests* **2020**, *11*, x FOR PEER REVIEW 6 of 12

0.74, indicating substantial predictive power (Table 1). In addition, for the model excluding an outlier,

25%, 50%, 75%, and 100% of the stems per plant were browsed (Figure 4c).

<sup>2</sup> 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

**Figure 4.** Relationship between actual and estimated willow yield with reference to deer browsing, using a generalized linear model (GLM) (**a**), and GLM simulations for estimated yield with browsing (**b**) and estimated yield loss due to browsing (**c**), assuming 0%, 25%, 50%, 75%, and 100% of browsed stems per plant. In panel (**a**), circles and triangles represent *Salix pet-susu* and *S. sachalinensi* clones, respectively; and an open symbol represents the clone that were determined to be an outlier based on the GLM residuals. Clone names are shown. The coefficient of determination (*r* 2 ), level of significance (*p*), and regression line of the standardized major axis analysis for 11 clones excluding an outlier clone **Figure 4.** Relationship between actual and estimated willow yield with reference to deer browsing, using a generalized linear model (GLM) (**a**), and GLM simulations for estimated yield with browsing (**b**) and estimated yield loss due to browsing (**c**), assuming 0%, 25%, 50%, 75%, and 100% of browsed stems per plant. In panel (**a**), circles and triangles represent *Salix pet-susu* and *S. sachalinensi* clones, respectively; and an open symbol represents the clone that were determined to be an outlier based on the GLM residuals. Clone names are shown. The coefficient of determination (*r* 2 ), level of significance (*p*), and regression line of the standardized major axis analysis for 11 clones excluding an outlier clone are also shown in panel (**a**). The dotted line represents a 1:1 relationship. Loess regressions were applied in panels (**b**) and (**c**). odt, oven dry ton. 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).

significantly correlated with actual yield across clones. The yield loss due to deer browsing was

**Figure 5.** Relationships between actual and estimated yield without deer browsing (**a**), between estimated yield without browsing and estimated yield loss due to browsing (**b**), and the percentage of browsed stems per plant and estimated yield loss due to browsing (**c**) across willow clones. An open symbol represents the clone that was determined to be an outlier based on the GLM residuals. Clone names are shown. The coefficient of determination (*r* 2 ), level of significance (*p*), and significant regression line of the standardized major axis analyses for the 11 clones excluding an outlier clone are **Figure 5.** Relationships between actual and estimated yield without deer browsing (**a**), between estimated yield without browsing and estimated yield loss due to browsing (**b**), and the percentage of browsed stems per plant and estimated yield loss due to browsing (**c**) across willow clones. An open symbol represents the clone that was determined to be an outlier based on the GLM residuals. Clone names are shown. The coefficient of determination (*r* 2 ), level of significance (*p*), and significant regression line of the standardized major axis analyses for the 11 clones excluding an outlier clone are shown. The dotted line in panel (**a**) represents a 1:1 relationship. odt, oven dry ton.

#### shown. The dotted line in panel (**a**) represents a 1:1 relationship. odt, oven dry ton. **4. Discussion**

#### **4. Discussion** *4.1. Yield Loss Caused by Deer Browsing*

*4.1. 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 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

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

11 clones during three years of rotation. As a result, it was possible to predict the yield only by the

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.

#### *4.2. 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−<sup>1</sup> yr−<sup>1</sup> in the first harvest rotation (average of approximately 8 odt ha−<sup>1</sup> yr−<sup>1</sup> ) 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−<sup>1</sup> yr−<sup>1</sup> 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–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

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−<sup>1</sup> yr−<sup>1</sup> ), 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.

#### **5. 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.

**Supplementary Materials:** The followings are available online at http://www.mdpi.com/1999-4907/11/8/809/s1, Table S1. Akaike's Information Criterions (AICs) of the generalized linear models predict willow yield with deer browsing. Figure S1: Residual plot, quantile-quantile plot, and scale-location plot of the generalized linear model used to determine an outlier clone.

**Author Contributions:** Conceptualization, H.H.; methodology, H.H. and S.S.; formal analysis, H.H., S.S., and T.Y.; investigation, H.H., Q.H., M.I., M.K., S.S., T.Y., H.U., and Y.M.; data curation, H.H. and T.Y.; writing—original draft preparation, H.H.; writing—review and editing, H.H., Q.H., M.I., M.K., A.U., S.S., T.Y., H.U., and Y.M.; project administration, H.H., Q.H., A.U., and H.U.; funding acquisition, H.H., Q.H., and H.U. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Forestry and Forest Products Research Institute (FFPRI) grants.

**Acknowledgments:** We thank our collaborators from Shimokawa town and the Shimokawa forest association, especially Yuji Takahashi and Toshio Yamamoto, for assistance with the establishment of this trial, maintenance and harvesting, and also Oji Holdings Cor. for offering willow cuttings, and Uotani-tekko Inc. for assistance with the harvest. We also thank Haruka Yamamoto, Eriko Ito, Tomomasa Amano, Naoyuki Furuya, Toshimitsu Nagasawa, Megumi Yasumoto, Masayoshi Takahasi, and Kiyohiko Fujimoto for their assistance in data collection.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*
