Ungulate-adapted forestry shows promise for alleviating pine browsing damage

High densities of ungulates can increase human-wildlife conflicts. Where forestry is an important economy, intensive browsing can lead to browsing damage, resulting in volume losses, poor stand regeneration


Introduction
Food subsidies from agriculture and forestry, selective ungulate harvests, and loss of large carnivores have resulted in an increase in certain wild ungulate species' distribution, abundance, and density (Apollonio et al. 2010).An increase in ungulate densities represents a challenge economically (Putman 1996, Putman et al. 2011) because ungulates are often the main drivers of plant population dynamics, forest structure, and ecosystem processes (Danell et al. 2003, Ross et al. 2016, Speed et al. 2019).Intensive browsing can reduce forest regeneration or cause shifts in plant species composition (Gill 1992, Rooney and Waller 2003, Schütz et al. 2003).As a result, forest owners and forestry planners look for practical, long-term measures to mitigate consequences of intensive browsing.
Measures to mitigate browsing damage include intentional feeding of ungulates away from critical areas such as young forest stands (diversionary feeding; Geisser and Reyer 2004), intentional feeding to increase food availability (supplemental feeding; Milner et al. 2014), exclusionary fences or deterrents, and/or changes to ungulate harvest strategies (Putman et al. 2011).The efficacy of each method is scale and context dependent (Boyce et al. 2017).Further, measures are often difficult to successfully implement because they are costly, ungulates cross jurisdictional boundaries, and researchers often work at fine spatial scales (Hobbs 2003) whereas wildlife management occurs at broader scales (Weisberg and Bugmann 2003), thus creating a mismatch between wildlife movement and management boundaries (Meisingset et al. 2018).
Spatial scale is a critical component to most ecological questions because resource selection occurs at different hierarchical orders, or scales (Johnson 1980, Senft et al. 1987, Boyce 2006) yet studies linking hierarchical scales of foraging to mitigation measures are rare.A foraging ungulate, for example, moves within its geographical range (1st order), establishes a home range (2nd order), and within that home range may make seasonal movements, selecting feeding patches (3rd order), individual plants and which parts of the plant to eat (4th order; Johnson 1980).Similarly, mitigation measures can operate at multiple spatial scales: coarse scales, where seasonal ungulate movements are altered by winter feeding stations (Jones et al. 2014), or fine scales, where resident ungulates are diverted away from high-traffic roads during winter (Milner et al. 2014).However, the scale at which ungulates respond to resources depends on multiple factors including temporal scale, physiology, life history traits, and habitat (Gaillard et al. 2010).
Northern latitudes experience strong seasonality in ungulate food availability and habitat conditions (Dussault et al. 2005), and functional responses to food availability and food quality are scale dependent.Moose (Alces alces), for example, the largest member of the deer family (Cervidae), select for abundant browse irrespective of quality at large spatial scales, but select for higher quality browse at finer scales (van Beest et al. 2010b).However, experimental increases in food availability do not necessarily increase food consumption (Edenius et al. 2014) and large scale selection could constrain available forage at a finer spatial scale (Wilmshurst et al. 1999).
Our study system is in the boreal forests of Scandinavia where forestry is a primary economy, moose have a high recreational and economical value, and moose-forestry conflicts are abundant (Lavsund et al. 2003).Industrialization of commercial forestry in the 1960s, combined with concurrent changes in moose harvest strategies, caused local moose densities to spike in the 1980s and 1990s (e.g., local winter densities of 5-6 moose per km 2 ; Lavsund et al. 2003, Speed et al. 2019).While moose densities have slowly declined over the past twenty years (Speed et al. 2019), as has food availability (Milner et al. 2013), moose continue to negatively affect tree growth and survival by browsing the bark or apical shoot, or by breaking the tree stem.Subsequent tissue damage and changes to growth morphology can reduce the economic value of growing trees and forest stands, resulting in what is commonly termed 'browsing damage' (Hörnberg 2001, Lavsund et al. 2003).In this study, we instead focus on the density of undamaged trees, as it is the density of unaffected trees that result in adequate stand regeneration rather than the damaged trees.We considered a tree 'undamaged' if the tree did not have top-shoot browsing, bark browsing, main stem breakage, or if <60% of a tree's shoots have been browsed (Hårstad 2008).In Scandinavia, browsing damage is primarily applied to Scots pine (Pinus sylvestris) and Norway spruce (Picea abies) because they hold commercial value.Scots pine is a primary moose food in winter and browsing damage most commonly occurs in young Scots pine forests (5-20 years old) during winter when food is limited and where moose congregate at lower elevations (van Beest et al. 2010b).
Economic losses as a result of browsing damage have resulted in conflicts between forest owners, who prioritize timber production, and moose hunters, who harvest on average 196,000 moose annually in Norway and Sweden (for study years 2011/12-2014/15; public data from www.ssb.no and www.älgdata.se).Wildlife managers have used moose population reduction via harvest as the primary strategy to decrease intensive moose browsing.However, since young pine forests are highly selected by moose, moose population reduction does not consistently reduce browsing damage (Reimoser and Gossow 1996).Diversionary or supplemental feeding (typically with silage) are other mitigation strategies, but feeding wild ungulates was recently banned in Norway due to the detection of chronic wasting disease (Stokstad 2017).Managers thus need to be creative in designing alternative, effective, long-term mitigation strategies in Scandinavia.
One mitigation strategy that has been minimally studied is the modification of conventional forestry practices during felling and site preparation to increase available forage.Branches and tree stems <5 cm in diameter from felled trees are left on site because they have low commercial value (Månsson et al. 2010).Conventional logging uses some of the branches as "slash mats" to reduce the compaction of underlying vegetation and soils from the heavy machinery.However, branches are crushed and become inaccessible to moose after snowfall, and shoots no longer have the biting resistance necessary for browsers because they are not attached to a tree (Månsson et al. 2010).A single mature Scots pine in Sweden holds on average 29 kg dry weight of moose forage (Månsson et al. 2010), yet only about 5% of potential forage remains available after the trees are felled and cleaned for hauling.Heikkilä and Härkönen (2000) found that residual Scots pine tree-tops raised above the snow, what we term ungulate-adapted slash piles, were utilized four times more than treetops lying on the ground.Machine operators can thus create ungulate-adapted slash piles with palatable species (e.g., Scots pine, birch; Shipley et al. 1998).This contrasts with traditional slash piles that include all felled species.Despite the increase in food availability, the use of ungulate-adapted slash piles (hereafter, slash piles) have not clearly mitigated browsing damage and require further research (Heikkilä andHärkönen 2000, Edenius et al. 2014).
In addition to slash piles, soil scarification is a common site preparation method in Scandinavia whereby the top organic layer is overturned to expose mineral soil, with the aim to improve seed establishment and increase soil temperature ( Örlander et al. 1996, Béland et al. 2000, Berg et al. 2008).An increase in soil scarification intensity can increase Scots pine seedling establishment (Saursaunet et al. 2018), thus increasing Scots pine density and food availability when trees are within browsing height and before stand thinning ( Örlander et al. 1996).While soil scarification can have many deleterious ecological effects (Atlegrim andSjöberg 1996a, Örlander et al. 1996), previous research in our study area found pine seedling density increased with mineral soil exposure (Saursaunet et al. 2018).Thus, intensive scarification could increase food availability during early-tomid successional stages, creating a forage-rich landscape and reducing browsing damage via resource dilution (Tscharntke et al. 2012).
At each stage of intervention in commercial forestry, actions are taken to optimize timber or pulp production, as described above.We propose 'ungulate-adapted forestry' be an additional step added to this process to optimize ungulate forage production.Here, we test two methods that can be part of 'ungulate-adapted forestry': ungulateadapted slash piles and intensive scarification.Our objective was to examine if ungulate-adapted forestry via slash piles and intensified soil scarification can alter ungulate browsing ecology and forest damage (Fig. 1).We tested if conventional and ungulate-adapted forestry influenced: (1) browsing damage; (2) moose habitat use; (3) bite size; (4) browsing pressure near treatment stands; and (5) whether these changes were similar at different spatial scales.Long-term monitoring plots were placed at varying distances from stands that were logged and scarified, rather than placed directly in logged and scarified stands.Current work (Mathisen et al. unpublished results) addressed the within-stand changes whereas our study investigated responses outside the stands.
Rarely is browsing damage connected to browsing ecology in the literature, yet it could explain many of the mechanisms influencing damage.For example, browsing damage could depend on the abundance of preferred species in the same plot.We expected a diversionary effect of the experiment to lead moose away from young stands, but also a supplementary effect to increase overall food availability, both leading to a decrease in browsing on Scots pine in the studied stands.We focused our analyses on Scots pine since it is a bulk winter food for moose in Scandinavia and has high economic value (Shipley et al. 1998).Our data were collected at the plot scale, which corresponds to the patch selection scale in hierarchical forage selection (3rd order; Herfindal et al. 2015), and at the tree and shoot scale, which corresponds to food selection (4th order; Senft et al. 1987).
Research suggests that pine damage decreases with an increase in pine forage availability (Bergqvist et al. 2014, Herfindal et al. 2015, Pfeffer et al. 2021).We thus hypothesized the density of undamaged pine stems would be higher in the areas close to slash piles (H1).Bite size, which is an index of forage intake (Gordon 2003), can reflect available forage.For example, moose select larger bites as browse density and quality decline, and as distance between patches increases (Vivas and Saether 1987, Shipley and Spalinger 1995, Shipley et al. 1998).Large bites require less handling time per unit biomass consumed, but result in a greater intake of fiber, which increases mastication, rumination, and digestion time.Small bites have less fiber, but require greater handling time per unit biomass consumed (Palo et al. 1992, Shipley 2007).Thus, bite size is a trade-off between food intake and quality.Because ungulate-adapted forestry increases availability of Scots pine (Mathisen et al. unpublished results), we hypothesized bite size to decrease near ungulate-adapted stands due to increased food availability (H2).
We hypothesized ungulate-adapted forestry would decrease browsing pressure (H3) because of the increase in alternative forage via slash piles and intensive scarification (Månsson et al. 2010, Edenius et al. 2014).Further, we assumed that consuming pine shoots from concentrated slash piles would be more efficient than browsing on dispersed trees.We hypothesized ungulate-adapted forestry to increase habitat use (H4) because of the creation of a forage-rich landscape and spill-over effects on surrounding stands, as was found with moose habitat use close to supplemental feeding stations (Gundersen et al. 2004, van Beest et al. 2010a).Based on research from supplemental feeding stations in Norway (van Beest et al. 2010a, Mathisen et al. 2014), we hypothesized stronger effects at smaller spatial scales for all response variables (H5).

Study area
Our study area lies between 60.8 • -61.4 • N and 12.2 • -12.7 • E in Innlandet County (Fig. 2; Fig. S1).Elevation ranges from 265 to 750 m above sea level.The area experiences cold, snowy winters (mean January temperature 2011-2018: − 9.3 • C; Norwegian Meteorological Institute) and short, cool summers.Land cover is dominated by boreal forests, which are managed for timber and pulp production.Production forests, which are largely coniferous, typically undergo one precommercial thinning at 10-20 years to remove competing deciduous shrubs and trees.Stands undergo 1-2 thinning events at 40-50 years and 70-80 years to optimize commercial tree density.All time estimates are dependent on site productivity.Natural regeneration from seed trees is most common for pine, whilst spruce are often planted.
In winter, moose typically migrate from summer ranges in highelevation mountainous areas to low-elevation valley bottoms where snow depths are reduced (Sweanor and Sandegren 1988, Bunnefeld et al. 2011, Singh et al. 2012).Roe deer (Capreolus capreolus) and red deer (Cervus elaphus) are present in our study area but occur at low densities.Large carnivores include wolves (Canis lupus), brown bears (Ursus arctos), wolverines (Gulo gulo), and Eurasian lynx (Lynx lynx).Wolves and bears prey on neonate calves in the spring and early summer, and moose are the main prey of wolves throughout the year (Swenson et al. 2007, Sand et al. 2008, Zimmermann et al. 2015).However, annual moose offtake by hunters, which only occurs in the autumn, is 2.4-3.5 times higher than that from wolves, where predation occurs year-round (Zimmermann et al. 2019).

Forestry activities
Our study had three forestry activity levels; no forestry actions, conventional forestry, and ungulate-adapted forestry.Conventional stands had conventional logging (no slash piles) and low intensity soil scarification (see Section 2.2.2).Ungulate-adapted stands had ungulateadapted logging (ungulate-adapted slash piles) and high intensity soil scarification (see Section 2.2.2).

Logging
Conventional logging does not make residual forest materials available for moose.In ungulate-adapted stands, harvester operators created slash piles from discarded Scots pine branches and tree-tops (stem diameters <5 cm) (Fig. 1).From November to March between 2012 and 2015, slash piles were created during felling.Slash piles varied in size and frequency by stand, particularly with available forage.In our study area slash piles doubled the amount of Scots pine biomass available, on average, compared to conventional logging (Mathisen et al. unpublished results).Harvester time spent in conventional and ungulate-adapted logging stands did not differ (Mathisen et al. unpublished results).Ungulate-adapted logging occurred in 46 stands (268 ha) while stands (204 ha) received conventional logging (Table S1).While we aimed for an equal number of conventional and ungulate-adapted stands, previous research has shown that many factors influenced when, where, and why small private Norwegian landowners harvested timber (Bashir et al. 2020) and these may not align with forest research goals.
Logging occurred at varying distances to plot centers (Fig. 2).To assess how spatial scale influenced our response variables, we created 250 m, 500 m, and 1000-m radius buffers centered on each monitoring plot in ArcMap (Environmental Systems Research Institute 2011).We limited our distance to 1000 m to reflect the maximum distance at which feeding stations have an effect on browsing and moose density in Norway (Mathisen et al. 2014).For each buffer size we created three unique logging variables: (1) a 4-level factor with 4 treatment levels (absence of logging, presence of conventional logging, presence of ungulate-adapted logging, and presence of both conventional and ungulate-adapted logging), (2) the area (m 2 ) with conventional logging; and (3) the area (m 2 ) with ungulate-adapted logging.All variables refer to logging that occurred the winter prior to spring data collection.See Table 1 for details.

Scarification
Scarification occurred during the summer from 2011 to 2014, one to two years after logging.Ungulate-adapted scarification occurred in stands (306 ha) while 47 stands received conventional scarification (390 ha; Table S1).Conventional scarification practices in Norway typically expose 13-20% mineral soil to improve seedling regeneration (Øyen 2002).We classified 0-20% mineral soil exposure as low intensity (conventional) and >20% as high intensity (ungulate-adapted) (Øyen 2002).Using the same buffer sizes as for logging, we created six scarification variables for each buffer size: variables (1) and ( 2) were 5-level factor variables indicating no scarification or the presence of 1-4 year old (time since scarification) conventional and ungulate-adapted scarification; variables (3) and ( 4) represented the area with ungulate-adapted scarification 1-2 and 3-4 years old in each buffer; and variables ( 5) and ( 6) represented the area with conventional scarification 1-2 and 3-4 years old in each buffer (Table 1).For the area scarified variables, we grouped stands by age because we would expect younger stands to have lower browsing biomass relative to older scarification stands.

Pellet counts
We counted ungulate pellet groups (Neff 1968) in 100-m 2 circular plots during late spring each year from 2012 to 2015.Pellet groups represent habitat use, or the time animal(s) spent in a plot (Månsson et al. 2007).Although moose were the dominant browser, we counted pellets for all ungulates present.We identified ungulate species by morphological characteristics of the pellets (Spitzer et al. 2019).To include a pellet group in our counts, >50% of the pellets from an individual group had to fall within the plot.Only piles with ≥20 pellets for moose and red deer, and ≥10 pellets for roe deer were counted.We visually distinguished between current winter and old (prior to winter) pellets.Typically, winter pellets were brown, in pellet form, and positioned on top of leaf litter and forest debris, while summer pellets were often in patty form, had leaf litter on top, or had mold or fungus growth (Zimmermann et al. 2015).We counted only current winter pellets, which corresponded to the winter browsing period.Pellets were removed from the plot each spring to avoid double counting the following year, except for the year prior to the start of our study.

Browsing surveys
In 2012 and 2015, we assessed browsing from the same plot centers as those where we conducted pellet counts but used 50-m 2 plots.We identified the tree species and counted the number of browsed and unbrowsed shoots from the previous growing season.A shoot was defined as live, woody growth ≥1 cm long.We registered browsing from only the current winter season (i.e., "fresh") where browsed shoots were still wet with resin and the wood had not died nor become grey (Ball and Dahlgren 2002).We counted browsed and unbrowsed shoots between 30 cm and 3 m above ground.The lower height represents average snow depth, below which trees are not available for browsing during winter, and the upper height represents the maximum browsing height for moose (Nichols et al. 2015).We counted browsed and unbrowsed shoots on up to 10 trees per species, working from the plot center to the edge in a spiral pattern.On each browsed tree, we used digital and manual calipers to measure twig diameter (to the nearest 0.1 mm) just below the bite for up to five bite diameters.When >5 browsed shoots were present, which was rare, technicians closed their eyes and grabbed shoots to measure.We also assessed browsing damage but only in young forests (cutting class 2; Table 2).See Table 2 for a summary of field measurements.

Statistical analysis
We used four response variables: density of undamaged Scots pine, bite diameters, browsing pressure, and moose habitat use (Fig. 1, Table 2).Browsing response variables were for Scots pine only.Undamaged pine were assessed only in young forests whereas the other three response variables represent forests of all age classes.

Model fitting
Model formulation differed by response variable.We defined a priori models using hierarchical regression models in a Bayesian framework (Gelman et al. 2013a).We fit presence of treatment (logging and scarification) and area of treatment in different models, resulting in four models per response variable per buffer size (48 models total).All logging and scarification models included additional habitat variables which were standard across each response variable (Tables 3 and 4).We fit logging and scarification models separately due to different expected temporal responses.With logging, we would expect slash piles to offer food only during the winter of logging and thus included data from all

Table 1
Forestry variables used for modeling.Each variable was extracted using 250, 500, and 1000-m buffers, which were centered on each plot.Scarification age refers to time since scarification (in years).Area (m 2 ) within buffers containing 1-2-year-old conventional scarification stands Area conventional scarification 3-4 yr Area (m 2 ) within buffers containing 3-4-year-old Area ungulate-adapted scarification 1-2 yr Area (m 2 ) within buffers containing 1-2-year-old ungulate-adapted scarification stands Area ungulate-adapted scarification 3-4 yr Area (m 2 ) within buffers containing 1-2-year-old conventional scarification stands years.In contrast, scarification offers increasing amounts of food with age.We therefore used browsing data only for 2015 for scarification models, since we would not expect sufficient regrowth to occur shortly after scarification.All continuous predictor variables were mean-centered and scaled to one standard deviation.We did not include correlated (Pearson r > | 0.7|) variables in the same models (Dormann et al. 2013).We included spline-based smoothers (k = 5) on all continuous variables, as we expected non-linear responses.
We built hierarchical models (Tables 3 and 4) by including different spatial scales (shoot, tree, plot; Table 2) in different sub-models (Szewczyk and McCain 2019).For example, we included variables measured at the tree-level as 'fixed effects' (population level) and variables measured at the plot level as 'random effects' (group level).For model fitting, we used weakly informative (default) priors on all parameters (Appendix III) and randomly generated initial values.We fit all models using package brms (Bürkner 2017(Bürkner , 2018) ) which uses the Stan programming language (Stan Development Team 2018).We ran 4 chains with 2000 iterations with 1000 warmup each, which resulted in 4000 posterior samples.We checked parameter convergence by visual inspection of the chains and with the Gelman-Rubin diagnostic (Gelman et al. 2013a).We evaluated model fit of the top model with posterior predictive checks.
Our interests with model fitting were two-fold: first, we wanted to compare conditional effects of logging and scarification across buffer sizes to identify scale-specific responses.We thus interpret all models.Second, we did model selection to identify the scale explaining the most variation for each response variable.We did separate model selection for logging and scarification models since the datasets were different.We determined the most parsimonious model using Watanabe-Akaike information criterion (WAIC) weights.WAIC weights (w i ) can be interpreted similar to the more familiar Akaike information criterion (AIC) as the relative support for the model, given the data (Gelman et al. 2013a).WAIC is appropriate for Bayesian approaches as it averages over the posterior distribution rather than conditioning on a point estimate (Gelman et al. 2013b).All models were fit in program R version 3.6.1 (R Core Team 2018).

Density of undamaged pines
The number of undamaged pine trees per 50 m 2 was calculated as: where a i is the number of undamaged pine assessed in plot i (field protocols capped this at 10), b i is the number of pine assessed (capped at 10), and c i is the total count of pine.This measurement is used to assess national forest regeneration regulations (Regulation of Sustainable Forestry: https://lovdata.no/dokument/SF/forskrift/2006-06-07-593).
We excluded trees >10 m to be consistent with the definition of a young forest (Table 2).
We used a gamma distribution y i Gamma(n i , α) with a log link where n i are the combination of predictor variables for each plot-level observation i and α is the shape parameter.We added 1e-5 to our response variable, as zeros were present (Zuur et al. 2009).We included fixed effect variables collected at the plot level: moose pellet counts, available shoots summed, birch density, and RAW density (Tables 3 and 4).Year was included as a fixed effect parameter because it had only two levels (2012,2015) and would be difficult to estimate variance if included as a random effect.Pine density was excluded as a predictor variable as it was included in the response.To account for the nested sampling design, our grouping structure was an intercept of quadrat nested within site (Tables 3 and 4).

Bite diameter
We used a gaussian distribution gaussian(μ f , σ) with an identity link where μ f is the linear combination of predictors for the tree-level observation f andσis standard deviation.A gaussian is appropriate to use on strictly positive values when the tail of the distribution has a low likelihood of overlapping zero, which is the case for bite diameters (Zuur et al. 2009).For fixed effect variables collected at the tree level, we included year, tree height, number of browsed shoots, number of unbrowsed shoots, damage presence, and accumulated browsing (Tables 3 and 4).Moose and pine density were collected at the plot scale and were used as slope terms, with the plot as the intercept.Because bite sizes are a trade-off between food intake and quality, we expected bite

Table 2
Field measurements used for modeling.The plot area was 50 m 2 for browsing variables and 100 m 2 for pellets.All variables were measured in each plot, except undamaged pine, which was assessed only in young forests (cutting class two).Response variables are indicated with a 'Y'.sizes to increase where habitat use was high and decrease with increasing pine density.We included quadrat and tree ID separately as intercepts.We included a unique tree identifier as we have repeated measures per tree (Tables 3 and 4).

Browsing pressure
We restricted analyses to Scots pine <10 m to be consistent with the density of undamaged pine analysis (n removed = 70; n total = 8067) and excluded trees where the number of available shoots was zero (i.e., no browse available; n = 2735 removed).We used a beta regression y f beta(μ f , φ)with an identity link forφ and a logit link for u f where f is the tree-level observation.We transformed the response variable to exclude zero and one: where n is the number of observations (Smithson andVerkuilen 2006, Cribari-Neto andZeileis 2010).We initially modeled browsing pressure with a binomial distribution, but the extreme tail led to poor model convergence and fit.We used a nested grouping structure to account for the design and repeated measures of multiple trees per plot.We included moose pellet counts and pine density as random slopes.Candidate models are presented in Tables 3 and 4.

Table 3
Candidate logging models for the four response variables density of undamaged pines, bite diameters, browsing pressure, and moose habitat use.Each model was run with 250-m, 500-m, and 1000-m buffer data.Subscripts refer to the level at which data were collected (i = plot, site = j, k = quadrat, f = tree).For random effect variables, values on the left of the | are used as slope terms and values on the right are intercept terms.Nested random effect terms are represented by /.Variables are defined in Tables 1 and 2.

Moose habitat use
We used the number of pellet groups as the response variable, which represents the time animal(s) spent in a plot.Månsson et al. (2007) identified pellet counts as an unbiased estimator of habitat use.We used a Poisson distribution y i Poisson(λ i ) with a log link, for each plot-level observation i.We included five plot-level population variables: year, pine density, birch density, and RAW density (Tables 3, 4).We lacked tree density data in years 2013 and 2014, when only pellet counts were conducted.Rather than excluding 2013 and 2014, which represented 45% of the full dataset, we imputed missing tree densities separately for each species (pine, silver birch, downy birch, rowan, aspen, willow spp.).Tree density in Scandinavia does not change substantially between years unless forestry activities occur (Hedwall et al. 2019).We used two complementary datasets to identify if a stand had forestry activities: (1) field data set (see Table 2) and ( 2) spatial dataset (i.e., boundaries of forestry activities in our study area; see Fig. 2).We used multiple imputation as a robust means to impute missing data (Sterne et al., 2009;White et al., 2011).Multiple imputation creates several imputed data sets based on other variables in the dataset.We created 10 different datasets (i.e., multiple imputation) and fit models separately to each dataset.This practice does have pitfalls (Rubin, 1996) but is becoming standard in the field of medicine, for example (Azur et al., 2011).We used multiple imputation by chained equations using the random forest algorithm from package mice (van Buuren and Groothuis-Oudshoorn 2011).We compared distributions of the univariate and multivariate datasets for each variable to evaluate prediction accuracy.Each model was fitted to the 10 datasets separately and posteriors were pooled across models.We evaluated sub-model convergence via r-hat values.

Density of undamaged pine
We registered browsing damage in 424 young forest plots (Table S2, Fig. S3).There were no ungulate-adapted scarification stands older than two years within the 250-m buffers.From the hierarchical models, the density of undamaged pine was highest near ungulate-adapted logging and where both conventional and ungulate-adapted logging occurred, regardless of buffer size (Figs. 3 and 4).The area logged had minimal effect on the density of undamaged pine (Fig. S4).The most parsimonious logging model was logging presence within 1000 m (w i = 0.82).
Plots near conventional and ungulate-adapted scarification stands showed contrasting relationships with the density of undamaged pine: plots close to one and four-year-old conventional scarification stands had the highest and lowest density of undamaged stems, though we note there was high uncertainty (Fig. 3).The area scarified had minimal effects on the density of undamaged pine (Fig. S5).The most parsimonious scarification model was scarification presence within 250 m (w i = 1).From the top model, an increase in moose habitat use increased the density of undamaged pine from 13.9 to 21.2 undamaged pine per 50 m 2 when moose pellet group counts changed from 0.4 to 4.1 (Fig. S6).All logging and scarification models suffered from relatively poor model fit for large values of the response (Fig. S7).
For scarification analyses, bite diameter data were restricted to 2015 (n = 1111 bites).From the data, maximum bite diameters were highest where scarification did not occur regardless of scarification type (Fig. S9).From the hierarchical models, bite diameters did not decrease as the scarified stand aged as expected (Fig. S9).The area scarified had minimal effects on bite diameters (Fig. S12).The most parsimonious scarification model was scarification presence, but weights were split among buffer sizes (500 m: w i = 0.30; 1000 m: w i = 0.30; 250 m: w i = 0.20).

Browsing pressure
We assessed browsing pressure on 5252 Scots pine (Table S2, Fig. S13).From the hierarchical models, the area with conventional logging within 250 m had the highest browsing pressure at intermediate stand area treated (4.8 ha).The area with ungulate-adapted logging had a weak negative effect on Scots pine browsing pressure (Fig. S14).Parameter uncertainty increased as the buffer size decreased.The presence of logging had the strongest effect on browsing pressure at 250 m, where areas near ungulate-adapted logging stands had 27% lower browsing pressure than conventional logged stands (Fig. 5).The most parsimonious model was the area logged within 250 m (w i = 0.95).
The area scarified had no apparent effects on browsing pressure, except for young ungulate-adapted scarified stands (1-2 years old) within 250 and 500 m where browsing pressure spiked at intermediate stand area (Fig. S15).The presence of scarification stands had the greatest effect at 250 m where browsing pressure was lowest at age one (0.06, 90% CI = 0.04-0.10)and highest at age two (0.24, 90% CI = 0.20-0.29)(Fig. 5).The most parsimonious model was the area scarified within 250 m (w i = 1).

Moose habitat use
We registered moose pellet groups in 3630 plots (Table S2).Mean pellet groups were higher where logging occurred (logging = 0.43, SD = 1.09; no logging = 0.29, SD = 0.77) within 250 m (Fig. S16).From the hierarchical models, habitat use within 250 m was 67% lower near ungulate-adapted logged stands relative to conventional stands (Figs. 4  and 6).The area logged had minimal effect on habitat use (Fig. S17).The most parsimonious model was logging presence within 250 m (w i = 0.87).From the top model, predicted habitat use was 1.8 times as high in young (cutting class two) and mature forests (cutting class five) relative to clear cuts (cutting class one), and habitat use peaked with an optimal pine density (Fig. S18).
For scarification models, habitat use declined as the conventional scarification stand aged whereas habitat use increased as ungulate- adapted scarification stand aged (Fig. 6).This effect was most pronounced at 1000 m.Moose habitat use decreased with the area of ungulate-adapted scarification, with strong non-linearities for plots near conventional scarification (Fig. 7).The top model was scarification presence within 250 m (w i = 1).

Discussion
We tested whether ungulate-adapted forestry changed browsing damage, bite diameters, browsing pressure, and moose habitat use across three spatial scales.For logging, the density of undamaged pine increased near ungulate-adapted stands (supporting H1), bite diameters increased where logging did not occur (not supporting H2), and browsing pressure (supporting H3) and habitat use (supporting H4) decreased near ungulate-adapted logged stands.For scarification, results were more equivocal, but habitat use decreased over time near conventionally scarified stands but increased near ungulate-adapted stands (not supporting H4).Conditional effects of logging and scarification across response variables were most pronounced at 250 m (H5).
As expected, the density of undamaged pines was highest near ungulate-adapted logging, and where both conventional and ungulateadapted logging occurred, regardless of buffer size (H1).The density of undamaged pine in plots was 1.4-1.7 times higher near ungulateadapted relative to conventional logging stands across buffer sizes (Fig. 3).The mechanism for the damage decreases near ungulateadapted stands could be explained by the concurrent decrease in browsing pressure and habitat use (i.e., pellet counts).A previous standlevel analysis found that available forage biomass doubled and mean biomass removed was higher in our ungulate-adapted logged stands (available biomass: 98 kg per ha) relative to conventional logged stands (available biomass 53 kg per ha; Mathisen et al. unpublished results).This would suggest that moose were able to maintain intake rates of Scots pine from slash piles, while concurrently reducing time spent feeding in the surroundings (H4).
Our results show promise for ungulate-adapted logging as a measure to increase the density of undamaged pine stems.Most of our study area occurs in forest productivity zones (site index) F6-8 and F11-14 (Astrup et al. 2019).According to the National Norwegian Regulation of Sustainable Forestry, F11-14 areas should have a minimum of 1900 pines per hectare, but 2280 to 4560 are recommended.F6-8 areas should have a minimum of 950 pines per hectare, but 1520 to 2470 are recommended (https://lovdata.no/dokument/SF/forskrift/2006-06-07-593).From our hierarchal models, conventional logging within 250 m reduced densities of undamaged pine (2254 pine per ha, 90% CI = 718-9394) outside the recommended range.However, ungulateadapted logging increased densities (3882 pine per ha, 90% CI = 2048-7582) to within the recommended tree density.These results show strong support for ungulate-adapted logging as a mitigation measure against moose browsing damage.
While we found that the presence of logging treatments had a positive effect on the density of undamaged pine, the area treated rarely affected undamaged pine densities (Fig. S4).It is possible that we did not treat large enough areas, as has been suggested by other studies where extensive feeding is necessary to see any effects (Putman and Staines 2004).We suggest however that more research is required.For example, we expect forage saturation would occur (e.g., type II functional response) assuming many new individual moose do not move into the study area.Future research should focus on treating a range of stand sizes, particularly above our maximum ungulate-adapted logged stand size of 57 ha, to identify optimal stand sizes for specific contexts and scales.
We found browsing pressure was 27% lower near ungulate-adapted logged stands relative to conventionally logged stands (H3; Fig. 5).This was as expected due to the doubling in forage biomass via slash piles (Mathisen et al. unpublished results), increased browsing efficiency from slash piles versus dispersed trees, and potentially reduced plant secondary metabolites of slash contents because most shoots developed above moose browsing height.However, chemical responses to browsing remains poorly understood; for example, Burney and Jacobs (2012) found only one tree species (Thuja plicata) of a multi-species study that increased terpene production in response to simulated browsing.For unbrowsed trees, secondary metabolites differed by tree height and species: Nordengren et al. (2003) found secondary metabolites from field-collected tree samples (152-727 cm in height) increased with height for willows but decreased for birch.Our generalized understanding of tree chemical defenses in response to browsing remains limited.Despite the support of our hypothesis that browsing pressure decreased near ungulate-adapted stands, we should expect a two-fold decrease in browsing pressure relative to conventional stands if moose foraging matched the doubling in forage biomass.Our results did not support this.This mismatch could be because food availability is not the limiting factor in browsing damage.Our hierarchical models indicated food availability had a strong positive effect on the density of undamaged stems only when food availability was low.When food was highly abundant, there was no apparent effect on undamaged pine.
Another explanation for the mismatch is that slash pile contents are of lower quality than shoots available within browsing height, not higher quality as we suggested earlier.If this is the case, digestion of slash pile contents would take longer, thus limiting intake rates (Belovsky 1984).Moose may also require complementary diets (Felton et al. 2016(Felton et al. , 2020; but see Hjeljord and Histøl 1999), which additional pine browse would not facilitate.The complementary diet approach contrasts with typical supplementary or diversionary feeding designs, where the goal is often to maximize energy intake or assumes that food is the limiting agent, without regard to balanced diets.For example, moose that eat carrots and potatoes in supplemental feed may be fiber deficient.This may increase their propensity to browse tree bark, which is high in fiber (Felton et al. 2020).We suggest that diet mixing and alternative forage in supplemental feed requires further research.
Our study design contrasts with many typical supplemental and diversionary feeding studies.First, slash piles and scarification stands resulted in dispersed resources, rather than few point locations where only the most dominant individuals can feed (Ozoga 1972, Putman andStaines 2004).Dispersed resources should also decrease time spent at the feeding site, as moose could avoid intense competition for resources.This would also reduce risk of disease transmission since individuals can avoid overcrowding (Mysterud et al. 2019).This is pertinent in Norway, where chronic wasting disease was recently detected in moose (Stokstad 2017).Second, the 'feed' in our study is a natural part of a moose's winter diet so the potential for individuals to suffer from pH imbalances or insufficient fiber are reduced (Mysterud 2010).This also reduces the potential of affecting behavioral traits or 'natural' selection of fed individuals or populations (Mysterud 2010).
As expected, we found logging and scarification effects on browsing pressure, and logging effects on moose habitat use, were strongest at the smallest spatial scale (250 m; H5).This corresponds to third-order patch selection whereby moose adjust movements within their winter home range to feed or rest near recently modified stands.Similar feeding patterns have been found for white-tailed deer (Odocoileus virginianus) (Ozoga and Verme 1970) and migratory moose in Sweden (Sahlsten et al. 2010).This makes sense both from a movement and energy maximization standpoint.On average, moose move very little while on their winter range (on average 2 km per day in winter in Scandinavia; Pfeffer et al. 2018).We would thus not expect moose to make longdistance winter movements 'in search' of our treatment areas; rather, they would adjust patch selection from within their seasonal home range.From an energy maximization standpoint, moose could feed in logged or scarified stands, which are on average quite small in our study (133 ha for conventional and ungulate-adapted logging and scarification stands) and still be close to resting sites in mature forest stands where this is protective cover.Based on our results, we cannot identify at what distances moose are influenced by our treatments and therefore cannot suggest at which distances ungulate-adapted stands should be placed from young forests.To help answer this question, we suggest future studies use concurrent GPS-collar data to evaluate multi-scale responses to ungulate-adapted forestry.
One unexpected result was the difference in peak habitat use by time since scarification: habitat use near conventional stands increased over time whereas habitat use near ungulate-adapted scarified stands decreased (Fig. 6).One explanation for this could be that both conventional and ungulate-adapted stands attract moose, but because ungulate-adapted stands have more forage, moose spend more time in the scarified stands.Because there is less forage in conventionally scarified stands, moose instead forage more in the surrounding scarified stands.Indeed, pioneering trees such as birch, which are attractive browse for moose, dominate regrowth in Scandinavian boreal forests (Wam et al. 2010).Another explanation for the decrease in habitat use over time near ungulate-adapted stands could reflect nutrient loss, which is facilitated by intensive mineral soil exposure from ungulateadapted scarification.Nutrient loss could result in a slower regrowth period (Bergquist andÖrlander 1998, Knudsen 2014).This is supported by an analysis in our study area : Saursaunet et al. (2018) found current annual growth of Scots pine and downy birch declined as soil scarification intensity increased.Slower regrowth in ungulate-adapted stands could influence not only the biomass available but feeding preferences by moose: previous research has shown that moose browsing increases as the plant reaches moose chest height (Bobrowski et al. 2015), so older stands (e.g., four versus one-year-old stands) may offer not only more abundant browse, but that which has lower handling time.
Despite the possible benefits of soil scarification to increasing ungulate food availability, it has extensive negative ecological effects: it facilitates soil carbon release, intensive site preparation may increase nutrient loss and decrease long-term site productivity, and is detrimental to understory species like bilberry, which are an important food source for herbivores (Atlegrim and Sjöberg 1996b, Jandl et al. 2007, Bergstedt et al. 2008, Maillard et al. 2010).Thus, the possible benefits of increased food availability of Scots pine must be weighed against the many detrimental effects for intensive soil scarification to be justified.For our study, we had more ambiguous signals from scarification effects on increasing the density of undamaged pine.As such, we recommend ungulate-adapted forestry should focus more on creating slash piles versus extensive implementation of intensive scarification.Regardless of the benefits of feeding as discussed above, there are certainly risks associated with supplemental and diversionary feeding.Previous research has shown that feeding can change foraging patterns (van Beest et al. 2010a), restrict movement (Guillet et al. 1996), change the amount of time spent on seasonal ranges (Jones et al. 2014), increase the risk of disease transmission (Sorensen et al. 2013), and increase population productivity (Milner et al. 2014).A study reviewed by Putman and Staines (2004) found that the only variable correlated with red deer browsing damage was the presence of supplementary feeding: aggregations of deer around feeding stations produced high local densities, resulting in a significant increase in forest damage.Similarly, Mathisen et al (2014) found browsing pressure increased over time, likely due to an increase in carrying capacity from supplemental feeding.Feeding programs thus often require a simultaneous increase in ungulate harvest with concurrent population and forage monitoring.Feeding is also associated with a poor mismatch with the timing of migration and plant phenology (Jones et al. 2014): movements of migrating elk that did not use winter feed grounds in Wyoming, USA closely matched spring greenup.In contrast, fed elk stayed longer on stop-over locations resulting in a poor mismatch of green-up and later arrival to summer ranges.Fed elk also departed summer ranges early, resulting in fed elk spending nearly a month less on summer ranges than unfed elk.Migratory moose in Scandinavia (Singh et al. 2012) could display similar patterns, and the extended period on winter and transitional ranges could thus intensify browsing on the natural forage before winter feeding starts, thus counteracting the intended effects of supplemental feeding.
There are of course challenges in doing large-scale, long-term forestry experiments such as coordinating among different land tenures, having protocols followed at all levels of operation, and changing timber prices which make meeting forestry research goals unrealistic (Bashir et al. 2020).For example, we attribute much of our model uncertainties to few non-zero data points, meaning less logging and scarification occurred than we expected.In total, only 4.72 km 2 was logged during our study.Other items which researchers may consider in future research include how sawmills purchase timber (e.g., continuous supply, so forest owners are motivated to log both in the summer and not just the winter) as well as reduced costs with logging larger, continuous stands rather than more smaller, dispersed stands, which we had hoped for in our study.Despite these shortcomings, our data suggest that larger logged stands could increase density of undamaged pines.While beyond the range of our data, perhaps even larger ungulate-adapted logging stands would illicit a stronger positive effect on undamaged pine.While large-scale experiments are difficult, they are important since we will experience large scale habitat and wildlife range shifts with climate change, which will influence wildlife ecology, agriculture, and forestry.

Conclusions
Our results suggest that ungulate-adapted forestry can reduce browsing damage, but more work is needed to determine how the area logged can produce detectable effects on browsing damage and at what spatial scale we can see moose movement is affected.We found that the intensive scarification can reduce browsing damage as the stand ages, but this comes at a cost as soil scarification can have strong negative effects.We suggest that supplementary feeding should be followed by careful population and forage monitoring.Provided that feeding is dispersed, natural forage from adapted forestry could be a better alternative to silage feeding.Future research should focus on whether feeding ungulates with mixed forage (e.g., deciduous and coniferous) could better account for a complementary diet.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Overview of the experimental study design.The four response variables (bite diameter, browsing pressure, density of undamaged pine, moose habitat use) are pictured at center.The conventional and ungulate-adapted logging and scarification treatments, which occurred at varying distances from plot centers, are featured at left and right.Response variables are described in Sections 2.4.2-2.4.5 and in Table 2. Illustration by Heidi Loosen (loosenstudio.net).

Fig. 2 .
Fig. 2. Map of forestry activities from to 2015 at one of our three study sites (Gravberget).Scarification occurred from 2011 to 2014 and logging occurred from 2012 to 2015.Scarification occurred one to two years after logging.Only winter logging stands are pictured here.Inset map shows our study area (white rectangle) in southern Norway.Each site contained 20-22 quadrats.Each quadrat contained plots (grey circle).At each plot we measured the density of undamaged pine, browsing pressure, bite diameters, and moose pellet groups.See Fig. S2 for a map of all study sites.

Fig. 3 .
Fig. 3. Conditional effects of the presence of logging (top row) and scarification (middle and bottom rows) in 1000 m, 500 m, and 250-m buffers on the density of undamaged pine per 50 m 2 .Logging data were from Norway in 2012 and 2015 (n = 424).Scarification data were from 2015 (n = 177).Bars represent 90% credible intervals.Scarification age represents time (in years) since scarification.Note the different y-axis limits for the middle row, left panel and lack of scarified stands >2 years of age in the lower right panel.In the top row of panels, 'conv' stands for conventional and 'ung' stands for ungulate-adapted logging.

Fig. 4 .
Fig. 4. Posterior probability distributions for presence of conventional and ungulate-adapted logging across four response variables.Posteriors are from 250-m buffer models.Highest density intervals (HDI) are drawn at 90 (orange) and 100% (red).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) A.E. Loosen et al.

Fig. 6 .
Fig. 6.Conditional effects of the presence of logging (top row) and scarification age (middle and bottom rows) on moose pellet counts.Data were collected from 2012 to 2015.Error bars represent 90% credible intervals.Scarification age represents time (in years) since scarification.Note the different y-axis limits for the upper right and center panels.In the top row, 'conv' stands for conventional and 'ung' stands for ungulate-adapted logging.

Table 4
Candidate scarification models for response variables density of undamaged pines, bite diameters, browsing pressure, and moose habitat use.Each model was run with 250-m, 500-m, and 1000-m buffer data.Subscripts refer to level at which data were collected (i = plot, site = j, k = quadrat, f = tree).For random effect variables, values on the left of the | are used as slope terms and values on the right are intercept terms.Nested random effect terms are represented by /.Variables are defined in Table1.