Browsing Damage on Scots Pine: Direct and Indirect Effects of Landscape Characteristics, Moose and Deer Populations

Abstract

Reducing browsing damages from cervids (Cervidae) on economically valuable tree species is a challenging task in many countries. Apart from cervid abundance, landscape characteristics, such as forest composition, land use, forage availability and climatic conditions, may affect the degree of browsing through both direct and indirect effects. A better understanding of basic mechanisms in this complex system is needed to design efficient and convincing management strategies. Focusing on Sweden as a case, which has been widely studied using regression analyses only, we applied path analysis to test the validity of a model on the indirect and direct links between landscape characteristics, cervid populations, and browsing damages on Scots pine (Pinus sylvestris). Our results validated the tested model in which moose (Alces alces) density and pine availability directly influence browsing damages. Increasing amounts of pine forests, preferred deciduous trees, and young forest had positive direct effects on moose densities and thereby indirectly contributed to increased browsing damages. The density of smaller deer species showed no direct effect on browsing damages on pine. Path analysis corroborated our attempt to disentangle direct and indirect potential causal drivers of browsing damages and shows that the choice of statistical method may alter the understanding of mechanistic driving forces.

1. Introduction

Ecosystem interactions can be illustrated as complex webs with various direct and indirect relationships, which can be difficult to understand and analyze [1,2]. Nevertheless, insights into these complex systems are often needed to design efficient natural resource management strategies that contribute to desired outcomes [2,3]. An example of such complex systems affecting human interests is the dynamic interactions between cervids (Cervidae), forage, and browsing damages to economically important trees. Cervids can cause different types of browsing damage such as top shoot browsing, stem breakage, bark stripping or fraying (hereafter, we will collectively refer to these as browsing damage) which can reduce tree growth and cause stem deformations, thereby negatively affecting timber quality [4,5,6,7].

Several studies have investigated the drivers of browsing damage across Scandinavia (e.g., [8,9,10,11,12]) and other places of the world (e.g., [13,14,15,16,17,18]). However, previous analyses on browsing damages are almost exclusively limited to testing direct relationships without demonstrating the interplay between diverse ecological factors and distinguishing direct from indirect effects. These studies commonly use regression analyses to investigate direct effects or interactions between explanatory variables on the response variable. Testing for such direct relationships may identify correlations and dependence, but it cannot prove causation [19,20]. Furthermore, by neglecting the systemic perspective of ecological mechanisms in regression analyses, research ignores the possibility that a predictor may indirectly influence the response variable via other predictors, so-called mediators [2,21].

In Sweden, both wildlife and forestry are of national importance and separately managed while highly interacting with each other. Several studies that have applied regression analyses to explain browsing damages [11,12,22,23] have triggered a debate on whether increasing pine availability is more efficient in alleviating browsing damages than reducing the moose population [24,25]. This narrative has been reinforced by a recent report claiming that lowering the moose population by about 20% has had little effect on browsing damages on pine [26]. Thus, the present discussion in Sweden regarding browsing damages is focused on basic causal relationships and whether lowering the moose population really would have an effect on browsing damages. Therefore, different analytical approaches than regression analyses might be useful to test patterns of mechanistic relationships from a more systemic and holistic perspective. One such approach is path analysis, as it allows indirect effects to be separated from direct effects by testing a model describing causal relationships [1]. Thus, path analysis expands on regression analysis by adopting the subjective understanding of a system and tests its validity based on assumptions of causality [2].

In this study, we used path analysis to disentangle direct and indirect effects steering browsing damage and to compare our findings to previous results from regression analyses. The dataset used is similar to datasets in previous regression studies. We elaborated and tested a model of basic causal links (see Section 1.2) between landscape characteristics (i.e., forest composition, land use, forage availability, and climatic conditions), cervid densities, and browsing damage levels (Figure 1). This approach provides a better systemic understanding of the ecological processes and mechanisms that determine browsing damages. This may lead to more convincing knowledge support for wildlife management. More specifically, we focused on the ecological chain of events affecting damage to the economically important Scots pine (Pinus sylvestris, hereafter pine) in Sweden.

1.1. Swedish Forest and Ungulate Management

Sweden produces around 5% of the global timber, pulp, and paper production and is one of the world’s largest exporters of these products [27]. At a national level, the forestry sector contributes a vital part of the Swedish economy and has a long-standing history of modifying the landscape. Swedish forests are shaped by silvicultural practices, meaning that forest stands are heavily managed for production purposes, using regeneration, afforestation, fertilization, and even-aged stand management with single dominating tree species to increase yields [28]. In large parts of the country, pine as one of the main conifer tree species (besides Norway spruce, Picea abies) dominates the landscape in form of rather homogeneous production stands. In contrast to Norway spruce, pine is a heavily utilized tree species by cervids in young forest stands [29,30,31]. The consequences of browsing damages on pine are related to high economic losses for the forest industry [32].

Moose (Alces alces) is the largest cervid species in Scandinavia, and its diet contains the biggest share of pine compared to other cervids, especially during winter [31]. Additionally, the Swedish moose population is considered among the densest in the world [33,34]. Moreover, expanding populations of the smaller cervid species, namely roe deer (Capreolus capreolus), red deer (Cervus elaphus), and fallow deer (Dama dama), during the past few decades [35] are of importance for the browsing dynamics in Sweden [36] (hereafter, we will refer to these three species collectively as deer while distinguishing them from moose). The expansion of cervids is also common elsewhere in Europe [37]. To some extent, intensified forest management practices, such as even-aged forest management, contributed to the population increases of cervids in Sweden, where young forests provide high amounts of forage [38,39]. At the same time, large carnivores are unevenly distributed across the country with low densities at the national scale [40,41].

The Swedish wildlife management system is to a large extent centered around moose, and hunting is the primary cause of death for the species [42]. Since 2012, moose is managed according to an adaptive co-management system to improve the basis for balancing the number of animals to levels that are acceptable with regard to browsing damage levels and other societal interests. Moose management areas (MMAs, in Swedish Älgförvaltningsområden) represent the focal management unit within the system to initiate management at a larger “ecosystem scale” (see [43] for a detailed description of the moose management system). Boundaries of MMAs are intended to accommodate the migratory behavior of moose by encompassing at least one distinct moose population [44]. Each MMA is managed by a moose management group (in Swedish, Älgförvaltningsgrupp) consisting of representatives for landowners, hunters, and in some places Sami (indigenous people) who collectively formulate adaptive management plans for each MMA. Similar management systems do not exist for the smaller deer species in Sweden, although hunting of red deer also requires a specific license.

Ten years after the introduction of the moose management system, the targets on maximum levels of browsing damages on pine (i.e., maximum 5% annually damaged pine trees; [45]) which have been formally adopted by the Swedish Forest Agency (SFA, in Swedish Skogsstyrelsen) are still not reached in the majority of the MMAs [46]. Although the current system provides more data to support operative management decisions, there still seems to be a need for a coherent interpretation of the data. In this study, we partly used data that originate from the management system and are commonly used as a foundation for decision-making.

1.2. Drivers of Browsing Damage—Theoretical Framework

Previous regression analyses regularly identified cervid populations and in particular moose densities as direct predictors of browsing damage on pine in Scandinavia [8,11,13]. Recent studies showed that more of the variation in browsing damage could be explained when other cervids were included in the analysis [12,18]. However, for Sweden, this may be a consequence of the latitudinal gradient, which is correlated with both pine density and other cervids, specifically roe deer (Capreolus capreolus). Nevertheless, smaller deer species can directly cause damage at lower heights [47], but more importantly, they may influence the feeding behavior of moose due to competition over alternative forage, e.g., ericaceous shrubs (Vaccinium spp.; [36]), and thereby indirectly affect damage levels. Thus, we hypothesize that both moose densities and deer densities will directly affect browsing damage levels (Figure 1, path c) and that interspecific competition over forage resources influences the intake rates and diet composition of moose and deer.

Besides cervid populations, landscape characteristics (i.e., forest composition, land use, forage availability or climatic conditions) have been identified as direct predictors of browsing damage in previous regression analyses [9,10,12,18]. However, here we see a need to apply a more nuanced understanding of the direct and indirect pathways of how landscape characteristics influence browsing damage (Figure 1). Landscape characteristics determine the carrying capacity (e.g., reproduction rate, fitness) for different cervids [48]. Furthermore, landscape characteristics determine habitat selection and foraging behavior, thereby steering densities across spatial and temporal scales [49,50]. For example, landscapes dominated by forests and silviculture are associated with moose [48,51,52], whereas agricultural landscapes and pastures promote the occurrence of deer [53,54]. Thus, we assume that landscape characteristics will directly influence cervid densities and thereby have an indirect effect on browsing damage (Figure 1, path b + c).

At the same time, certain landscape characteristics can have direct effects on browsing damage (Figure 1, path a). For example, forest composition can directly influence the likelihood of browsing damage on certain tree species dependent on their availability. We assume a frequency-dependent relationship on selection, i.e., that selectivity for a preferred tree species will increase if its availability is low at the landscape level [55]. A higher number of available stems will result in a lower relative level of browsing damage due to a dilution effect [22,56]. Another example is climatic conditions, which, besides their indirect effects via cervid densities, can also have direct effects on browsing damage. Snow cover can restrict the mobility of cervids during winter [57,58], causing overexploitation of a local browsing resource due to stationarity, independent of local cervid densities. Thus, snow can restrict the accessibility of forage [58,59,60], directly affecting damage.

In this study, we tested the above-described theoretical relationships between variables according to Figure 1 using path analyses.

2. Materials and Methods

2.1. Study System

Our study system included all MMAs in Sweden (55.6°–68.0° N, Figure 2). For data extractions, we used the boundaries of the MMAs that were in place during 2019 since MMAs undergo continuous administrative changes over time. We first tested a model across all MMAs (n = 144) on a national level because the same moose management system and supporting data are used throughout the country. However, we assumed that ecological mechanisms influencing damage levels might vary in different parts of the country due to differences in landscape and climate properties across Sweden’s large latitudinal gradient. Therefore, in a second step, we divided MMAs along the boreal and nemoboreal vegetation zones [61] into northern and southern MMAs (Figure 2). This border aligns with the general latitudinal differences in tree species regeneration (increasing pine and decreasing spruce regeneration from south to north; [27]), agriculture (decreasing from south to north; [62]), the cervid community (decreasing multi-species cervid communities from south to north; [63]), and climatic conditions (increasing gradient of stronger and snow-richer winters from south to north; [64]). This separation resulted in 68 MMAs belonging to the northern group and 76 MMAs belonging to the southern group (Figure 2).

2.2. Dataset

Data collected for this study originated from various sources (see

Table 1. Description of variables used in this study. Sample size representing the number of national moose management areas (MMAs) with available data was 144 for all variables.

Category	Variable	Description	Year	Source
Browsing damage (=dependent variable)	Damaged pine 2	Proportion of annually damaged Scots pine stems	2016–2019 (mean value; see Section 2.2.1)	Moose browsing inventory [46]
Forest composition 1	Pine availability 2	Total number of available Scots pine stems per hectare land area	2016–2019 (mean value; see Section 2.2.2)	Moose browsing inventory [46]
Land use 1	Pine forest	Proportion of pine forest	Raster from 2019	National vegetation cover [65]
	Agricultural land	Proportion of agricultural land	Raster from 2019	National vegetation cover [65]
Forage availability 1	Preferred deciduous 2	Proportion of sampling plots in young forests where rowan, aspen, willow or oak is present	2016–2019 (mean value, see Section 2.2.5)	Moose browsing inventory [46]
	Young forest 2	Proportion of forage producing young forest derived from satellite images	2016–2019 (mean value; see Section 2.2.5)	Forage prognosis [46]
Climate	Snow days	Number of days with snow	Mean annual value for climatic normal period 1961–1990	Climatic raster [66]
Cervid populations	Moose density 2	Simulated number of moose per 1000 ha land area	Hunting seasons 2015/16–2017/18 (mean value; see Section 2.2.3)	Harvest data [67] were simulated in Helge [68]
	Deer index	Sum of shot roe deer, red deer, and fallow deer per 1000 ha land area	Hunting seasons 2015/16–2017/18 (mean value; see Section 2.2.3)	Harvest data [69]

¹ These variable categories are included in the overall term landscape characteristics throughout the text. ² Data are officially available and operatively used within the Swedish moose management.

). Most data are officially available and used operatively within the Swedish moose management by the county administration boards (in Swedish, Länsstyrelserna) and moose management groups.

2.2.1. Browsing Damage

We used the proportion of annually damaged Scots pine stems (i.e., damaged pine; see Table 1) as the response variable. Data on damaged pine are compiled by the SFA in a moose browsing inventory (in Swedish Älgbetesinventering, Äbin), which is an officially used annual monitoring system for browsing damages in young forest regenerations with a mean height of 1–4 m. The method records the number of pine stems with at least one of three types of damage: browsing on the apical leader shoot, damaged bark (e.g., bark stripping or fraying) and broken stem. Each spring, young forest stands are monitored after snowmelt in approximately half of the number of MMAs per county according to a systematic random sampling design (see [12] for details). The method focuses on trees valuable for silviculture, so-called production trees. Therefore, only stems taller than half of the average height of the two tallest production stems are considered per sampling plot. Due to a partly high variation in results from year to year [70], the proportion of damaged pine stems is presented by the SFA as a 3-year average value per MMA to be used for management decisions. We received the raw data from 2016 to 2019 directly from the SFA and calculated the mean proportion of damaged pine stems per MMA across years with available data.

2.2.2. Forest Composition

Based on previous studies (e.g., [12]), we expected pine availability (see Table 1) to have a direct link to the proportion of damaged pine stems (i.e., browsing damage; Figure 3). For each MMA, we calculated pine availability by dividing the estimated total number of pine stems in young forest stands by terrestrial area (ha). Again, we analyzed the mean value per MMA across years with available data from 2016 to 2019.

2.2.3. Cervid Populations

Moose is expected to be the main damage agent on pine when considering the above-mentioned three types of browsing damage and the height span of 1–4 m. Therefore, we expected a direct effect of moose density on damaged pine (Figure 3). Harvest statistics are mandatorily reported by hunters to the county administrative boards [67], whereas moose observations (in Swedish Älgobs) are voluntarily reported (i.e., collected by hunters in autumn; see also [71]). We simulated moose density (i.e., number of moose per 1000 ha terrestrial area) in Helge (Simthinc Version 5.3; see Supplementary Material File S1, [72]), which took harvest statistics and information from moose observations (i.e., proportion of adult males, reproduction as number of calves per female) from the hunting seasons 2015/16 to 2017/18 under consideration. In addition to mortality related to hunting, a fixed additive general mortality rate was assumed equally for all MMAs. The regional differences in mortality due to, for example, the occurrence of predators were not explicitly included in the simulations because of a lack of information at the scale of our study. The simulated density correlates well with harvest statistics (r = 0.795, p < 0.001). We then calculated a mean moose density for these years. We performed this simulation because moose is highly managed in Sweden, meaning that harvest statistics themselves are dependent on the specific management goal of each MMA, which in turn is dependent on recent damage levels. Therefore, we believe that a simulation is closer to the true density index than the annual harvest and comparable between MMAs.

Next to moose, deer (i.e., roe, red, and fallow deer) could potentially directly damage pine trees inventoried by the SFA, or they could, perhaps more likely, indirectly affect damage levels via resource competition with moose. To account for this in our path model, we included a direct effect of the deer species on damaged pine and a residual covariance between moose and deer to represent the relationship between them without indicating a cause–effect association (Figure 3). We calculated a mean deer index (i.e., sum of harvested deer per 1000 ha terrestrial area) representing density purely based on hunting statistics from the seasons 2015/16 to 2017/18. Harvest data of deer are voluntarily reported by hunters and collected in an online database per hunting district (in Swedish Jaktvårdskrets) [69]. To establish harvest data of roe deer, fallow deer, and red deer for each MMA, we used the proportional area overlap of MMAs and hunting districts (see Supplementary Material File S1).

2.2.4. Land Use

Due to better light conditions and the high availability of alternative forage in the field layer [73,74], we considered the proportion of pine forest as a foraging habitat of higher value for moose than spruce-dominated forests. Ericaceous shrubs are important forage for all cervids, and competition seems to influence the foraging behavior of moose by increasing its pine intake [36]. Therefore, we included the variable with a direct link to moose density, which we assume is partly or fully mediating the effect of pine forest on damaged pine (Figure 3).

For deer, pastures and agricultural areas are important foraging habitats [53,54], especially during the vegetation growth period. Therefore, we included a direct effect of agricultural land on deer densities, which could in turn lead to an indirect effect on pine damage. We tested a similar direct link to moose density; however, we expected the opposite effect since agricultural areas represent on the one hand a less utilized moose habitat [75] and on the other hand less forest [65].

We extracted the proportion of pine forest and the proportion of agricultural land per MMA from the national vegetation cover, which was compiled in a 10 m raster by the Swedish Environmental Protection Agency (SEPA, in Swedish Naturvårdsverket) in 2019 [65]. To extract a quantified land use feature (i.e., proportion) for each MMA, raster images were interpreted with image pixel analysis.

2.2.5. Forage Availability

Young forest stands that are established after a mature forest is felled (either by clear-cuts or storm events) provide a high quantity of forage for deer and are therefore a preferred foraging habitat [11,39]. Information on young forest area is presented in an official database by the SFA as part of the reports of the moose browsing inventory. These data are derived from a difference analysis of satellite images, which is a standardized method used by the SFA. It gives information on the annual forested area that has been finally felled (i.e., clear-cut). A model is then used to estimate the aggregated area of forest with an approximated average stand height between 1 and 6 m. The model considers regional differences in site productivity. Data from 2016 to 2019 were directly received for all MMAs from the SFA. We then calculated the mean proportion of young forest area per MMA for this four-year period.

Tree species such as rowan (Sorbus aucuparia), aspen (Populus tremula), willow (Salix spp.) and oak (Quercus robur) are highly preferred and selected deciduous forage by cervids [11,76]. Therefore, we expected that the presence of preferred deciduous species in young forests of MMAs is associated with the proportion of suitable habitat. The presence of preferred deciduous species is also monitored in the moose browsing inventory by the SFA. We used the proportion of sampling plots with the presence of preferred deciduous species and again calculated a mean value for MMAs with available data from 2016–2019.

2.2.6. Climate

In many parts of Sweden, snow cover is prominent for several months. Image analysis was used to extract the proportion of area with different durations of snow cover from a raster map. This map was produced by the Swedish Meteorological and Hydrological Institute (SMHI) representing the defined normal period from 1961 to 1990 by the World Meteorological Organization [66]. We then calculated the number of days with snow cover (i.e., snow days) as a weighted average for each MMA.

2.3. Statistical Analysis

Confirmatory path analysis was applied to test the theoretical model illustrated in Figure 1. Using MMAs (n = 144) as the unit of analysis, we simultaneously estimated a multitude of direct and indirect effects between forest composition, land use, forage availability, climate, cervid populations, and browsing damage. We fitted the same model to the complete dataset (n = 144) to create a national model for all of Sweden and a split dataset to test the model for northern Sweden (n = 68) and southern Sweden (n = 76). All tested effects and evaluated models are illustrated in Figure 3.

Before fitting the model, we inspected the datasets to check if they fulfil the underlying assumptions (i.e., multicollinearity and multivariate normality) for path analysis. The dataset had no missing values, and descriptive statistics for all variables can be found in

Table 2. Summary of variables used in the path analysis, including their sample size (n), mean value (M), standard deviation (SD) and range (Min–Max) within the complete dataset for the national model and the respective samples for the northern and southern Sweden models.

Category	Variable	Unit	All of Sweden (n = 144)				Northern Sweden (n = 68)				Southern Sweden (n = 76)
			M	SD	Min	Max	M	SD	Min	Max	M	SD	Min	Max
Browsing damage	Damaged pine	proportion	0.15	0.07	0.05	0.38	0.13	0.07	0.05	0.38	0.17	0.07	0.05	0.34
Forest composition	Pine availability	1 ha−1	28.60	25.72	0.31	109.02	45.76	26.82	1.80	109.02	13.25	10.65	0.31	54.12
Land use	Pine forest	proportion	0.22	0.10	0.02	0.41	0.26	0.09	0.03	0.41	0.18	0.09	0.02	0.39
	Agricultural land	proportion	0.11	0.12	0.00	0.62	0.05	0.07	0.00	0.33	0.16	0.12	0.02	0.62
Forage availability	Preferred deciduous	proportion	0.49	0.16	0.02	0.83	0.41	0.17	0.02	0.81	0.55	0.13	0.23	0.83
	Young forest	proportion	0.06	0.02	0.00	0.14	0.06	0.02	0.00	0.10	0.05	0.03	0.00	0.14
Climate	Snow days	count	115.88	41.48	40.00	215.00	150.87	32.17	87.00	215.00	84.57	15.79	40.00	121.00
Cervid populations	Moose density	1000 ha−1	7.21	2.61	0.10	17.15	7.90	2.64	1.61	17.15	6.60	2.45	0.10	12.17
	Deer index	1000 ha−1	7.55	7.75	0.00	53.60	2.00	1.97	0.00	8.51	12.52	7.63	3.98	53.60

. For the complete dataset and the dataset for southern Sweden, Pearson’s product-moment correlations between all variables were below 0.70 (see Supplementary Material File S2, Tables S1 and S3) and variance inflation factors (VIFs) were below 4. Therefore, we deemed the degree of multicollinearity acceptable and applied no remedies. For the northern part of the dataset, correlations between deer density, snow depth, and agricultural land exceeded these thresholds (see Supplementary Material File S2, Table S2). However, as the level of multicollinearity was still moderate (VIF < 6), was only restricted to a sub-section of the model, and did not include the dependent variable (i.e., browsing damage), we decided to fit the hypothesized model to the northern dataset. Nevertheless, as multicollinearity might influence parameter estimates or confound the contribution of independent variables, this model section needs to be interpreted with caution. Mardia’s multivariate skewness coefficient (947.89, p < 0.001) and the multivariate kurtosis coefficient (11.15, p < 0.001) indicated that multivariate normality could not be assumed in the complete dataset, nor in the split datasets. Thus, we decided to use lavaan’s MLR estimator (a robust version of the maximum likelihood estimator; [77]) for our path analysis. MLR calculates robust “Huber–White” standard errors based on the observed information matrix and a robust likelihood ratio test statistic, which is asymptotically equivalent to the Yuan–Bentler T2* test statistic [78]. All indirect effects of land use, forage availability, and climate via cervid populations on browsing damage (see Figure 3) were estimated using the delta method to calculate their standard errors.

Model fit was assessed via a robust version of the chi-square (χ²) results for testing the exact fit hypothesis, the Tucker–Lewis index (TLI, good model fit > 0.95), the root mean square error of approximation (RMSEA, good model fit < 0.08), and the standardized root mean square residual (SRMR, below 0.05 indicating good fit) [20]. All statistical analyses were conducted in R [79] using lavaan [77], a covariance structure analysis package. We treated all variables as continuous and used 0.05 as the cut-off for statistical significance.

3. Results

3.1. National Model

The path analysis showed a good fit of our hypothesized model at the national level with a robust χ²-value of 15.991, (df = 8, p = 0.043) and an SRMR of 0.028 indicating a good absolute model fit. A robust RMSEA of 0.078 and robust TLI of 0.945 confirmed the acceptable fit of our path analysis. Unstandardized estimates, standard errors, and results of significance testing for all estimated parameters of the national model are presented in

Table 3. Unstandardized parameter estimates (est), corresponding standard errors (SE), outcomes of significance testing (Z- and p-value) and standardized estimates (std. est) of the national model. The table contains estimates for the path coefficients, followed by covariance estimates between collaboration variables, and variance estimates for all endogenous variables. Estimates of indirect and total effects have been calculated; standard errors for these have been computed using the delta method. An equivalent table for the northern model can be found in the Supplementary Material File S2, Table S4.

Variable		est	SE	Z-Value	p-Value	Std. Est
Moose density
	Pine forest	5.762	2.391	2.410	0.016	0.213
	Agricultural land	−5.375	1.964	−2.737	0.006	−0.239
	Preferred deciduous	3.906	1.130	3.457	0.001	0.246
	Young forest	39.764	8.753	4.543	0.000	0.369
Deer index
	Agricultural land	4.867	5.118	0.951	0.342	0.073
	Young forest	−45.598	18.579	−2.454	0.014	−0.142
	Snow days	−0.107	0.011	−9.489	0.000	−0.574
Damaged pine
	Pine availability	−0.001	0.000	−3.831	0.000	−0.404
	Pine forest	−0.019	0.078	−0.246	0.806	−0.026
	Moose density	0.007	0.002	3.367	0.001	0.257
	Deer index	0.001	0.001	1.513	0.130	0.094
	Snow days	0.000	0.000	−0.210	0.834	−0.021
Covariance
	Moose density—deer index	−1.296	0.897	−1.444	0.149	−0.113
Variance
	Moose density—moose density	3.874	0.668	5.798	0.000	0.574
	Deer index—deer index	33.689	13.393	2.515	0.012	0.565
	Damaged pine—damaged pine	0.004	0.001	7.367	0.000	0.780
Indirect effects
	Pine forest via moose density	0.041	0.021	1.945	0.052	0.055
	Agricultural land via moose density	−0.038	0.017	−2.265	0.023	−0.061
	Preferred deciduous via moose density	0.028	0.010	2.733	0.006	0.063
	Young forest via moose density	0.282	0.112	2.517	0.012	0.095
	Agricultural land via deer index	0.004	0.005	0.790	0.429	0.007
	Young forest via deer index	−0.040	0.028	−1.442	0.149	−0.013
	Snow days via deer index	0.000	0.000	−1.542	0.123	−0.054
Total effects
	Pine forest	0.022	0.077	0.282	0.778	0.029
	Snow days	0.000	0.000	−0.803	0.422	−0.075

.

Forty-three percent of the variation in moose density was explained by the occurrence of pine forest (β = 0.21, p = 0.016), preferred deciduous (β = 0.25, p = 0.001), young forest (β = 0.37, p < 0.001), and agricultural land (β = −0.24, p = 0.006). In line with our predictions, the latter had a negative estimate indicating that moose densities are lower in landscapes dominated by agriculture. Surprisingly, agricultural land showed no significant effect on the other deer species (β = 0.07, p = 0.342). Nevertheless, 44% of the variation in the deer index was explained by the negative effects of young forest (β = −0.14, p = 0.014) and snow cover (β = −0.57, p < 0.001); both of these are negatively correlated with agricultural land (see Supplementary Material File S2, Table S1) and reflect the scarcity of these deer species in northern Sweden. The model indicated a negative relationship between moose and deer, but the estimate was non-significant (β = −0.11, p = 0.149).

Our model was able to explain 22% of the variation in pine damage, showing that higher moose densities lead to increased damage levels with a standardized path coefficient of β = 0.26 (p = 0.001). Higher pine availability showed a negative effect on damages (β = −0.40, p < 0.001). Contrary to our prediction, neither the proportion of pine forest in the landscape nor snow days showed a significant direct effect on pine damage (Figure 3a). The calculation of indirect effects showed significant positive effects of preferred deciduous (β = 0.06, p = 0.006) and young forests (β = 0.10, p < 0.012) on pine damage levels via the increase in moose density (Table 3). Accordingly, the indirect effect of agricultural land was negative (β = −0.06, p = 0.023). The indirect effect of pine forest on damaged pine via moose density was just above the threshold for statistical significance (β = 0.06, p = 0.052). Given that deer index had no direct effect on pine damage, none of the remaining indirect effects were significant (Table 3).

3.2. Models for Northern and Southern Sweden

The model for northern Sweden showed a good fit with a robust χ²-value of 8.272 (df = 8, p = 0.407), SRMR of 0.019, robust RMSEA of 0.022 and robust TLI of 0.997. The explained variation of pine damage was 21%, nearly identical to the national model, and both pine availability (β = −0.35, p = 0.013) and moose density (β = 0.27, p = 0.005) remained as the two variables with a significant effect (Figure 3b). Although the standardized estimate for deer index showed a similar value (β = 0.27) to that for moose density, its effect on damaged pine was not significant (p = 0.183). Regarding the direct effects of landscape characteristics on cervid populations, the northern model showed some difference compared to the national results (Figure 3). Thirty-nine percent of the variation in moose density was explained, whereby only preferred deciduous trees showed a significant effect (β = 0.48, p = 0.000), while the effects of pine forest (β = 0.17, p = 0.276), agricultural land (β = −0.10, p = 0.215) and young forest (β = 0.24, p = 0.092) lay above the threshold for statistical significance. Given the high correlations of deer index with snow days and agricultural land, it is not surprising that the northern model showed a drastic increase in the explained variance of deer index (R² = 0.80, Figure 3b). Compared to the national model, young forest was not a significant predictor of deer index in northern Sweden (β = −0.07, p = 0.248), while agricultural land (β = 0.48, p = 0.000) and snow depth (β = −0.49, p = 0.000) showed similarly strong but opposing effects. Out of all tested indirect and direct effects of landscape characteristics on browsing damage, only the presence of preferred deciduous showed a significant indirect effect by increasing moose densities, and thereby the proportion of damaged pine was significant (β = 0.13, p = 0.009, see Supplementary Material File S2, Table S4). All estimates, standard errors and results of significance testing for direct and indirect effects of the northern model can be found in the Supplementary Material (see Supplementary Material File S2, Table S4).

Fitting our hypothesized model to the southern dataset resulted in an unsatisfactory model fit (robust χ² = 29.223, df = 8, p = 0.000; robust SRMR = 0.067; robust RMSEA = 0.162, robust TLI = 0.498). Thus, it is not recommendable to interpret the estimated parameters. The retained R²-values indicated that the model was only able to explain 4% of the variation in deer index in southern Sweden. We see this as an indication of a limited understanding of multi-species communities in southern Sweden and/or potentially shortcoming of current management data in delivering needed information in multi-species systems.

4. Discussion

4.1. General Discussion

Path analysis verified large parts of our model on the dynamics between landscape characteristics, cervid populations and browsing damage levels on Scots pine in Sweden. Our results confirmed that moose density and the availability of pine trees are major direct drivers of browsing damage levels, which is in line with previous regression analyses [9,12,22]. However, path analysis also helped us to disentangle complex effects in a more refined way as compared to ordinary regression models. Thus, the present study contributes to a better systemic understanding through an analysis approach that allows for the testing of indirect pathways.

More specifically, our hypothesis that alternative forage (i.e., preferred deciduous trees) and land use properties (i.e., cover of pine forest, agricultural land and young forest) indirectly affected browsing damage via their effect on moose densities was corroborated by the results. Land use characteristics determine both habitat quality and forage availability for moose and deer, which steer cervid densities [80] and therefore only have an indirect effect on browsing damages.

Previously it has been suggested that damage levels on pine are more closely linked to the availability of pine than to moose density [12,22]. The introduction of a mediator (i.e., moose density) between landscape characteristics and browsing damage increased the statistical inference for moose density associated with browsing damages, indicating that the conclusions drawn are influenced by the setup of the statistical model. Path analysis has the advantage of elaborating more on the causal mechanisms between variables, which could be of more value for decision support in management. This is perhaps of specific importance when trying to build and test models for data that partly highly correlate with each other and are collected at large spatial scales and across different areas as in our study. In the case of Sweden, we tried to capture important north–south gradients via our variable selection and by splitting our model since these gradients influence, for example, stand productivity, tree species composition, snow cover, and deer species abundances.

Contrary to previous research [12], we did not find a direct effect of deer index on browsing damage. We argue that the separation of direct and indirect effects with path analysis, as performed in this study, made it possible to denote with higher accuracy the role of moose in the complex interplay of herbivores, habitat, forage, and damages. The covariance between the smaller deer species and moose was not significant. This could be explained by a moose population that is more evenly distributed across MMAs as compared to deer, which exist in higher densities towards the south. Although these species are taxonomically close, they differ in terms of forage selection, such that moose and roe deer include a higher proportion of browse in their diets as compared to red deer and fallow deer [31]. Based on this, the lack of a clear relationship between moose and deer in our analyses was expected, as were the contrasting effects of landscape characteristics on these animal groups. We find it likely that the direct contribution of deer to damaged pines in the stand height interval analyzed in this study (1–4 m) is marginal, and they do not compete with moose over the same habitats. Nevertheless, because some important forage species, such as Vaccinium shrubs, are shared between the four cervid species, deer may play a moderating role in the relationship between moose and pine damage [36].

Although preferred deciduous trees (i.e., rowan, aspen, willow or oak) are highly selected as forage items by moose, these trees are usually of relatively low abundance (approximately 10% of total tree abundance) within commercially used young forest stands [46,76] and therefore only of marginal importance when considering the overall intake of food by moose [31]. Nevertheless, quality aspects such as body weight seem to be influenced by the composition of available forage species [81]. Our model highlighted that preferred deciduous had a positive effect on moose density, which indirectly affected browsing damage on Scots pine. An indirect reason for this could be that the occurrence of preferred deciduous tree species also may be linked to more nutritious areas in general, therefore offering browse of better quality [75,82,83]. If this is the case, we expect that general forest productivity is associated with both moose density and damage. However, we cannot conclude that forest stands with a higher presence of preferred deciduous trees are associated with a higher risk of browsing damage on pine, although other studies indicate this [84].

Neither pine forest, representing a foraging habitat with abundant understory vegetation, e.g., ericaceous shrubs, nor snow days, representing winter harshness, had a direct effect on browsing damage, which was contrary to our assumptions. With testing a direct effect of snow days on browsing damage, we assumed that on average colder and longer winters in MMAs would affect the browsing behavior of moose due to potential limited movement [85,86] as well as making the field layer less available for foraging [60]. We expected that snow days would lead to higher browsing damage levels independent of respective cervid densities. Instead, snow days had an effect on deer index indicating that the distribution of the smaller deer species is limited by stronger winters. Regarding the proportion of pine forest in the landscape, we expected that this habitat with a higher abundance of favorable forage as compared to the darker spruce forest types [87] would reduce the browsing pressure on young forest stands. However, we did not find such an effect, which could imply that our variable might not sufficiently represent forage availability in the landscape or that most of the variation is already included in the variable pine availability.

4.2. Comparison of National and Regional Models

When we split the dataset into two parts, our attempt to fit the same model to the northern and southern regions of the country separately showed that the strength of direct and indirect drivers of browsing damage varies across the country. This may be related to variation within variables along national gradients and how variables differ between the national and regional scales. For example, the strength of the effects of landscape characteristics on cervid populations differed considerably between the northern model and the national model. While all forage and land use variables significantly influenced moose densities in the national model, only the amount of preferred deciduous trees showed a significant effect in the north. This might be explained by a lower variation in land use, more specifically pine forest, within the northern region than at the national scale and therefore less influence within the model. At the same time, this means that the lower occurrence of agricultural areas in the north becomes of higher relevance for deer that are known to utilize both forest and agricultural land where available [53,54]. Besides these differences in the role of landscape characteristics on cervid populations, the northern model showed similar results to our national model concerning browsing damage, with moose densities and pine availability being the only direct drivers.

Our theoretical model did not fit the southern part of the dataset. Apart from a lower sample size and less variation across MMAs when dividing Sweden into two parts, there may be dynamics and processes of a more complex nature in southern Sweden due to a higher diversity in land use and a more pronounced multi-species cervid system with considerably higher densities of all deer species. For example, the abundance of deer species is not to the same degree explained by landscape characteristics that follow a latitudinal gradient, as it is at the national scale. Furthermore, there is a general east–west gradient of pine forest in southern Sweden, which is negatively associated with site productivity and annual precipitation. Therefore, we find it likely that a finer resolution or better quality of the data may be needed to be able to map the detailed processes and improve model fit in southern Sweden.

4.3. Management Implications

The degree of browsing damages within forests is dependent on a multitude of factors [88] where food availability and density of cervids play a central role. Although there are several known measures that can be undertaken to alleviate browsing damages [89], national targets are still far from being reached in Sweden [45]. The conflicting interest between silviculture and game hunting creates deep-rooted positioning on how to resolve the conflict. In Sweden, however, it is not only an issue that involves forest owners and hunters. Central authorities also have an interest as silviculture is of high economic and strategic importance for sustainability while game hunting is of high importance for recreational hunting, cultural heritage, and multiple uses of the forest ecosystem.

Although our study confirmed that the principal driving factors behind browsing damage levels on young pine were the density of moose and availability of pine, many of the additional management variables included in this study are certainly relevant for a better understanding of the pathways leading to high browsing damage. Our findings call for a more refined systemic perspective on browsing damage as compared to previous studies, which assessed direct relationships only [9,12,90]. Several landscape characteristics, of which some are available for decision support within the Swedish moose management system and some are not, can be attributed to indirect effects on browsing damage levels mediated via moose density. This emphasizes the need to include landscape data in moose management for effective mitigation of browsing damage. In fact, the current Swedish management system considers the entire terrestrial area of an MMA when planning harvest goals. However, the proportional share of the habitat (e.g., forest) of the animals (e.g., moose) within an MMA is rarely considered. A better knowledge of how landscape factors translate into suitable habitats and ultimately the reproduction potential for moose is a way to include ecology into management in a more refined way.

Our study showed that the understanding of the causal relationship between moose density and browsing damages is improved when viewing landscape characteristics as explaining factors for moose density and therefore only indirectly influencing browsing damages. Landscape characteristics should be actively considered in wildlife management as these characteristics affect the carrying capacity of cervid populations [48]. For example, a change in land use may not directly affect browsing damage, but it may alter the population size and fitness of cervid populations, which in turn may affect damage levels. At the same time, landscape changes might increase browsing damages independent of cervid densities when they enhance the predisposition of certain forest stands to browsing or lead to unsuitable spatial distribution of cervid populations [89].

A key factor that is difficult to disentangle is the impact of moose management on the model interactions. Management has a strong impact on the number of cervids and will respond to forest damages, to some degree. As the moose population and browsing damages vary considerably across MMAs, one can draw the conclusion that there are other, perhaps more influencing, factors that managers adhere to. Quality indices in the moose population, including, e.g., reproduction rate and slaughter weight, are examples of such factors. The number of large carnivores present is another important factor, which also influences effective reproduction (i.e., calf survival). If reproduction rate is of more importance for moose managers than browsing damages, one would expect a stronger link between browse availability, predator density, and moose density. Path analysis can be used to disentangle such a relationship.

4.4. Limitations of the Study and Future Research Needs

In comparison to regression analyses, path analysis tests previously defined cause-and-effect relationships between variables, and thus rests on basic assumptions concerning causality, the included variables and tested models. Strictly speaking, causation requires a sufficient degree of association (i.e., covariance) between two variables, that one variable occurs before the other (i.e., that one variable is clearly the outcome of the other) and that no other reasonable causes for the outcome are present [20]. Within our study, a sufficient degree of association was given within the variable set (see Supplementary Material File S2, Tables S1–S3). When it comes to the time order between variables, we see it as reasonable to assume that landscape characteristics are more slowly changing variables compared to yearly estimates of cervid densities or browsing damages. Thus, we see landscape variables as the pre-existing conditions that shape population development. Concerning the aspect of other reasonable causes, we followed common procedure and used existing theory, previous research and the collective experience of the author group to decide which variables should be included in a model to explain browsing damages. Limitations to our model and this study are partly data availability and the use of a confirmatory modeling approach, which excludes the testing of alternative models. Drivers of browsing damage are multifactorial and complex. While higher cervid densities are one of the central factors that can lead to increased browsing pressure, other density-independent aspects such as increased disturbances or nutritional imbalances could also increase browsing levels. At the same time, forest stands could have an enhanced predisposition for browsing due to spatial patterns in the landscape or unsuitable forest conditions [89]. Unfortunately, we were not able to include all of these aspects within our model. We hope that our study can stimulate a scientific discussion and encourage other researchers to test other, perhaps more refined, path models explaining browsing damages. These could include time-lags, such as the recovery time in a forest that has been heavily utilized by cervids, or include data collected through experiments to better match the source and scale of analyzed variables.

This study is, to the best of our knowledge, the first attempt of using path analysis with publicly available national data used in management to explain browsing damages across areas. Despite its limitations, this study has validity and importance to the managers and will add to the public debate on Swedish browsing damages by focusing on the basic causal links.

5. Conclusions

Path analysis can to some extent disentangle complex and interrelated factors by verifying a theoretical model of causal dependencies [20]. Our model showed that there are different direct and indirect pathways that are influencing cervid populations and browsing damages on economically valuable tree species. This allows for multiple potential management strategies to decrease browsing damages. Regulating moose density can offer an effective tool to directly control browsing damage levels, especially in a short time frame. Managing landscape characteristics (except for snow), via, e.g., large-scale changes in silviculture, will eventually also have indirect impacts on browsing damage since cervid densities (e.g., moose within the current Swedish management system) will respond to changes in landscape properties. Forest management acts on a longer time frame than hunting, and there are no collective objectives to adhere to among forest owners. However, in collaboration with forest companies and landowner associations, the SFA advocates increasing the amount of pine forest in southern Sweden. Such advice may fall short if moose populations are managed without consideration of such transitions. Overall, our study shows that even if the management of landscape characteristics and cervid populations currently are detached from each other in the Swedish management regime, they are ecologically intricately intertwined, which needs to be acknowledged for more effective management of browsing damages. Integrating a multi-species perspective (i.e., all four deer species including moose) within game management would be necessary to successfully alleviate browsing damages in forest stands.