Supplementary Information for the Paper Titled

Supplementary Information for the paper titled:

Cervid exclusion alters boreal forest properties with little cascading impacts on soils

Anders Lorentzen Kolstad, Gunnar Austrheim, Erling J. Solberg, Aurel M.A. Venete, Sarah J. Woodin, James D.M. Speed

Figure S1. Scatter plot showing the interaction effect of 8 years of ungulate exclusion and annual aboveground tree biomass production (left) and site productivity index (right) on the Canopy Cover Index (CCI) at 12 boreal forest clear-cut sites in central Norway (the box showing the overall exclusioneffect for the right figure). Dots above the dashed line represent sites where ungulate exclusion increased the CCI (measured in full summer). Solid line is the least square fit. The right relationship is shaped by the fact that mean CCI index inside the exclosure was contributing 50% to the specification of the productivity index.

Figure S2. Depth profiles and tables of means of carbon (C) and nitrogen (N) ratios (C:N), concentrations, and stocks, and soil pH. The data are from boreal forest clear-cuts in central Norway, with and without browsing by ungulate herbivores. The litter was not included in the total C, as it was assumed to represent only a transient carbon storage.

Figure S3. Figure showing the interaction effect of 8 years of ungulate exclusion and the Productivity Index on total N stocks (kg m-2). Shaded areas are 95% confidence intervals.

TableS1. Spearman’s rho correlations between plant root simulating (PRS) probes data and tea bag index (TBI) data with burial length, temperature sum (summed daily averages) and precipitation sums (accumulative). The PRS probes were subject to varying burial lengths, whereas the tea bags were buried for approximately similar duration thus making it easier to isolate and interpret the effect of temperature and precipitation on these variables. Asterix refer to P-values and denotes the significance level (* P < 0.05; ** P < 0.01; *** P < 0.001).

PRS probes (N=72)
Burial length / Temperature sum / Precipitation sum
Burial length / 1
Temperature sum / 0.903*** / 1
Precipitation sum / 0.878*** / 0.900*** / 1
NO3 / 0.384*** / 0.513*** / 0.441***
NH4 / -0.295* / -0.432*** / -0.407***
Total N avaliability / -0.226 / -0.242* / -0.257*
Relative nitrification / 0.410*** / 0.556*** / 0.507***
TBI (N=66)
Burial length / Temperature sum / Precipitation sum
Burial length / 1
Temperature sum / 0.013 / 1
Precipitation sum / -0.029 / 0.379*** / 1
S / 0.030 / -0.486*** / -0.276*
k / 0.148 / -0.011 / 0.023

Table S2. Results of linear mixed effects models testing the effect of moose exclusion and site productivity on 26 ecosystem properties. Main effects are reported as results from likelihood ratio tests comparing against a minimum adequate model. Significant parameters are in bold (α=0.05). Arrows indicate the direction of change. Square brackets refer to concentrations.

Transformation / # Obs. / # Sites / Moose
Exclusion / Productivity
Index / Interaction
effect / ICC
C stocks 0-10 cm / Log / 30 / 15 / χ2(1)=1.8564
P = 0.1730 / χ2(1)=0.0471
P = 0.8283 / Site: 0.4643
C stocks 10-20 cm / Sqrt / 26 a / 13 / χ2(1)= 4.4350
P = 0.0351 / Site: <0.001
C stocksTotal / Log / 30 c / 15 / χ2(1)=2.5964
P = 0.1071 / χ2(1)=1.7571
P = 0.1850 / Site: 0.3399
[Nitrogen] litter / 30 / 15 / χ2(1)=0.1408
P = 0.7074 / χ2(1)=3.759
P = 0.0525 / (↑) / Site: 0.5319
N stockslitter / 30 / 15 / χ2(1)=0.5804
P = 0.4462 / χ2(1)=6.2769
P = 0.0122 / ↑ / Site: 0.2270
N stocks 0-10 cm / 30 / 15 / χ2(1)=1.2519
P = 0.2632 / χ2(1)=1.1406
P = 0.2855 / Site: 0.5360
N stocks 10-20 cm / 26 a / 13 / χ2(1)= 5.0861
P = 0.0241 / Site: 0.0252
N stocks total / Log / 30 / 15 / χ2(1)= 6.4679
P = 0.0110 / Site: 0.4661
C:Nlitter / 30 / 15 / χ2(1)=0.6291
P = 0.4277 / χ2(1)=5.1033
P = 0.0239 / ↓ / Site: 0.7887
pH litter / 30 / 15 / χ2(1)=1.0834
P = 0.2979 / χ2(1)=4.0724
P = 0.0436 / ↑ / Site: 0.6290
pH 0-10 cm / 30 / 15 / χ2(1)=1.2261
P = 0.2682 / χ2(1)=0.0170
P = 0.6455 / Site: 0.4252
pH 10-20 cm / 28 c / 14 / χ2(1)=1.8842
P = 0.1699 / χ2(1)=0.0727
P = 0.7874 / Site: 0.6191
aTwo sites removed due to lacking soil depth. bOne location only went to 10cm depth, but was still included. cOne sites removed due to lacking soil depth ICC = intraclass correlation coefficient.

Table S3. Table with mean values (SE) of ecosystem properties analysed in relation to 8 years of ungulate exclusion from boreal forest clear-cut sites in central Norway. N=12 for all except C:N, C and N stocks, pH and bulk density where N=15.

Variable / Unit / Open plots / Exclosure
S / - / 0.375 (0.011) / 0.365 (0.017)
k / - / 0.010 (0.001) / 0.010 (0.001)
NO3 availability / µg 10cm-2burial length-1 / 4.2903 (2.583) / 2.829 (1.342)
NH4availability / µg 10cm-2burial length-1 / 16.611 (2.986) / 15.528 (2.057)
Total ION availability / µg 10cm-2burial length-1 / 20.918 (4.064) / 18.326 (2.698)
Relative nitrification / - / 0.123 (0.047) / 0.164 (0.061)
Average soil temperature a / Degrees Celsius / 11.621 (0.141) / 11.065 (0.136)
C:N (organic soil) / - / 28.768 (1.735) / 30.778(1.735)
pH (organic soil) / - / 4.091 (0.080) / 4.101 (0.087)
Total soil C stocks / kg m-2 / 8.296 (0.713) / 9.137 (0.498)
Total soil N stocks / kg m-2 / 0.323 (0.040) / 0.339 (0.019)
Bulk density (organic soil) / kg m-3 / 132.738 (16.259) / 119.208 (13.171)
D:C biomass ratio / - / 0.639 ( 0.154) / 3.260 (1.327)
Canopy Cover index / percent / 14.117 (2.950) / 32.735 (5.015)
a5cm depth, between 15June-25July. S = litter stabilization factor; k = decomposition rate; ION = inorganic nitrogen; D:C = deciduous : coniferous.

Extended methods (EM)

Biomass models

A dataset containing heights (height classes of 50 cm; precise heights only for year 2016), and diameter at ground level (DGL; 2016 only) recorded in permanent subplots since 2009 (see Speed and others (2013) for details) were used as input for locally calibrated biomass models to estimate above-ground biomass of individual trees for each year of the experiment. To calibrate the biomass model, we harvested some trees of birch (Betulapubescens), rowan (Sorbusaucuparia), pine, and spruce, from outside the experiment, but in close proximity to four of the sites (site numbers 5, 6, 8, 12; Table 1), using a stratified random approach to cover the entire relevant size range for each species. Biomass functions were not calculated for species rarely encountered within the experiment: The biomass of Salix capreaand Betula pendula was estimated using the biomass function for Betulapubescens, and similarly, Juniperuscommunisbiomass was estimated using the function originally fitted for spruce. No rowan above 79 cm could be found anywhere outside the exclosures, so we sampled a sufficient number of trees from inside the fences (close to the fences) to enable us to fit separate biomass models for trees inside and outside the exclosures. We also sampled one birch and two pine individuals from inside the fences. The vertical height of the live shoot apex and the DGL were recorded before the trees were harvested at ground level and dried to constant weight at 60ºC and then weighed. Biomass models were obtained for each species by starting with a saturated model including cubic, quadratic and linear terms for both height and DGL, and performing backwards elimination by removing the least significant terms first, and ending with the model with the highest adjusted R2-value (Table EM1).

Table EM1. Best-fitting models to explain sapling biomass based on height and diameter at ground level (DGL). Only rowan had sufficient samples of un-browsed trees (from inside exclosures) to allow two models. N = sample size (of which are un-browsed trees in parentheses).

N / Height range (cm) / DGL range (mm) / Equation / Multiple R2
Birch / 28 (1) / 18 - 298 / 3.73 - 37.00 / y = 0.078843DGL2 + 0.009197HGT2 / 0.9458
Pine / 26 (2) / 23 - 256 / 4.00 - 91.43 / y = 0.325839DGL2 + 0.0007434DGL3 / 0.9856
Spruce / 24 (0) / 25 - 296 / 4.58 – 63.58 / y = 0.020293HGT2 + 0.006092 DGL3 / 0.9671
B Rowan / 19 / 23 - 79 / 2.83 - 24.22 / y = 0.006664HGT2 + 0.082983DGL2 / 0.8123
UB Rowan / 12 / 51 - 405 / 6.00 - 41.95 / y = 0.0053962HGT2 / 0.9841
Birch = Betulapubescens; Pine = Pinussylvestris; Spruce = Piceaabies; Rowan = Sorbusaucuparia; B = browsed; UB = un-browsed.

Intercepts were supressed (set to zero) in all models. Model parameters with negative coefficients were removed to ensure positive biomass estimates, as long as the reduction in model fit as a result of their removal was negligible.

The models were used to estimate standing biomass inside and outside the exclosures in 2016 for deciduous and coniferous species separately. Trees above 6 meters height were excluded from the analysis because they had potentially large effects on the results due to measurement error and because they are remnants from before the clear-cutting and therefore not interesting in terms of the impact of recent browsing release.

Productivity Index

To characterise the sites in terms of productivity potential, a possibly important covariate, we estimated the mean annual increment in biomass since the start of the experiment for inside the exclosures only. Because DGL was only available for one of the years in the dataset, we also constructed biomass models for un-browsed trees of birch, rowan and pine, using only height as the predictor (Table EM2).For spruce, the model was built on samples from both inside and outside the fences as the growth forms did not differ (Figure EM1).

Figure EM1. Relationship between diameter at ground level (DGL) and height for trees measured in Trøndelag, central Norway, in 2016. Diverging lines (least square fit) indicate different growth forms inside and outside the exclosures, so that for all species except spruce, trees inside herbivore exclosures were taller for the same diameter at base compared to trees from outside the exclosure.

For un-browsed rowan we had sufficient samples from inside the fence (n= 12) to fit a model separately (Table EM1). For un-browsed pine and birch, we opted to refrain from cutting trees inside the permanent exclosures and instead fitted the biomass models using a ‘back-fitting’ procedure explained next. First, we used the original biomass models (Table EM1) to estimate the 2016 standing biomass inside the exclosure (Figure EM2-A). Second, using the measured height and the estimated biomass for 2016, we back-fitted another set of biomass models, using the same model selection protocol as described for the original biomass models (Table EM2; Figure EM3), and then estimated the biomass inside the exclosures for all years before 2016.By doing this we inescapably lose track of the associated residual errors of the models, but because the R2 of the ‘best-fitting models’ were so high (0.9458 and 0.9856 for birch and pine, respectively), and because visual checks confirmed that the few data points originating from un-browsed individuals did not appear as outliers (Figure EM4), we believe the ‘back-fitting’ approach was justified, and indeed the best compromise. Finally, we calculated the estimated mean annual increment in tree above-ground biomass inside the exclosure (standardised to per area unit).

Table EM2. Biomass models for un-browsed saplings of birch (Betulapubescens), pine (Pinussylvestris), rowan (Sorbusaucuparia) and spruce (Piceaabies) constricted to only height as predictor variable.

Equation / Multiple R2
Birch / Y = 0.170274HGT + 0.010018HGT2 / 0.9937 a
Pine / Y = 0.0149667HGT2 / 0.8841 a
Rowan / y = 0.0053962HGT2 / 0.9841
Spruce / y = 0.038068HGT2 / 0.9481
aModel was fit on estimated (not measured) values for biomass

Initial analyses revealed that annual biomass production did not fully reflect site productivity potentials. This is based on comparing the ranking of sites in terms of annual biomass production against our own perception of site productivity after multiple visits to all sites, and by observing that the mixed effects models produced predictable outliers. To overcome this, values of average annual biomass increments and the CCI for the exclosed treatment (Figure EM2-C) were standardised by dividing by the maximum value held by any site, before being added together to make up a compound productivity index PI (Figure EM2-D; Table 1). By using two metrics together, we were able to soften the randomness produced by the finite sampling area of a single metric. Visual inspection confirmed a sensible ranking of the sites in terms of PI compared to impressions of perceived productivity potential from multiple field visits. The compound productivity index ranked the study sites in a way that was correlated with organic soil bulk density (Spearman’s rho rs=0.62), C:N ratio of the organic soil (rs=-0.69), decomposition rate (rs=0.46), pH of the organic soil (rs=0.42), and litter stabilization factor S (rs=0.48), all in a predictable manner, indicating its appropriateness as an explanatory variable.

Figure EM2. Standing biomass in year 2016 (A), average annual biomass increments (B), Canopy Cover Index (C), and Productivity Index (D) for the exclosed plots at all 15 sites in Trøndelag, central Norway. The site numbers refer to Table 1 and Figure 1.

Figure EM3. Back-fitted biomass models for un-browsed pine (Pinussylvestris; left) and birch (Betulapubescens; right).

Figure EM4. Best-fitting biomass models (from Table EM1) for pine (Pinussylvestris; left) and birch (Betulapubescens; right) with samples from un-browsed trees in black (for pine) or blue colour (for birch).

Canopy Cover

To estimate canopy cover we directly analysed photos in the field using the Gap Light Analysis Mobile Application (Tichý 2015) and a normal smartphone with an inbuilt camera (Samsung XCover 3, Android v. 5.1.1). A photo was taken, and a modified Canopy Cover Index (CCI) was calculated, for each of the four sides of 10 randomly placed 50 × 50 cm permanent vegetation quadrats per plot (total 40 pictures), and these were later averaged to produce a single value and CCI estimate for each plot. The smartphone was positioned horizontal by using a spirit level and the camera lens was placed 50 cm from the ground (the height of the vegetation frame). Tall herbs and grasses from the field layer were bent aside before the photo was taken. Shrubs (including Rubusidaeus) and trees were left as they were, but single leaves or branches were removed if they were closer than 2 cm from the lens. Each photo captured a 38º view of the hemisphere, and photos were excluded if they captured part of the fence (very few, <1%). The smartphone application was set to analyse all pixels, and in the blue colour spectrum only. Photos were taken in all the non-thinned sites during peak foliage, between 20th June – 24th July 2016.

Tea Bag Index

Relative measures of organic matter decomposition was evaluated using the Tea Bag Index (TBI; Keuskamp and others 2013). Pairs of Lipton green tea and Lipton rooibos tea were buried 8 cm deep and retrieved after c. 90 days (mean 92; range 85-97). Initial air-dry weight and final oven-dry weight (70ºC, 48 hours) of tea bags were obtained with 1 µg precision, and a correction factor was found for each tea type to correct for the difference between air-dry and oven-dry weight (factors used in calculating TBI are given in Table EM3). The TBI method relies on using two types of tea with very different decomposition rates (Keuskamp and others 2013). The set burial time is necessary to secure the assumption that the easily decomposable green tea has reached an approximately stable weight (hydrolysable biomass fraction is decomposed and only non-hydrolysable compounds remain) whilst the relatively recalcitrant rooibos tea is still in the exponential decomposition stage. By knowing the fraction of hydrolysable biomass of each standardized tea type, one is able to calculate both a stabilisation factor S indicating the completeness of decomposition (conceptually similar to the ‘limit value’ (Berg and McClaugherty 2008)), as well as the decomposition rate k, without the need for a time series and repeated visits to field sites.

Table EM3. Parameters used to calculate the tea bag index of decomposition (Keuskamp and others 2013).

Replicates used in estimate
Hydrolysablefraction green tea / 0.842 / Taken from the referenced article
Hydrolysablefractionrooibos tea / 0.552 / Taken from the referenced article
Weight (g) of nylon tea bag / 0.120319 / 18
Weight (g) of cord / 0.032975 / 4
Weight (g) of paper label / 0.093805556 / 18
Correction factor for green tea / 0.997 / 2
Correction factor for rooibos tea / 0.981 / 2

PRS probes

Four cation-anion pairs were deployed within each subplot and later combined before lab analysis, yielding 3 samples per experimental plot and 72 samples total (thinned sites were excluded). Plastic cylinders (Polyvinyl chloride (PVC) or polypropylene (PP)) of 10 cm diameter were driven into the soil around each probe pair in order to exclude roots which otherwise would compete with the probes for the same ions (Figure EM5). After retrieval, the probes were cleaned in the field using distilled water and a scrubbing brush. Composite samples were stored in moist conditions in zip-lock plastic bags and shipped to Western AG Innovations Inc (Saskatoon, Canada) where NO3-N and NH4-N was determined colorimetrically using a flow injection analysis system. The units of N availability is reported as µg 10 cm-2 burial length-1 and the burial length varied between 49 to 92 days (identical burial length for each exclosure-open plot pair). The length of burial was assumed long enough for the probes to have archived equilibrium with the soil solution, making the estimates less dependent on diffusion processes and more dependent on nutrient release.

Figure EM5. Picture of two PRS probes inside a plastic root exclusion cylinder.

Soil core sampling and analysis

The three subsamples from each plot were randomly chosen out of nine potential sampling points arranged as a grid around the centre of each experimental plot. If the soil was too shallow to get a full soil core, we tried to get a new core in close proximity. A mineral increment sample was only combined with the composite sample if it was more than 5 cm long/deep and the depth was noted down to adjust the sample volume accordingly. The three subsamples from each experimental plot were bulked into one composite sample, keeping the soil layers separate.

The soil samples were stored in individual plastic zip-lock bags and kept cool until they could be transferred to paper bags and dried at 60 °C for a minimum of 48 hours. In the lab, the litter samples were cleaned of green moss, woody material >5 mm diameter, and live biomass, before they were homogenized using a stone mortar and passed through a 2 mm soil sieve. Stones and roots that were caught in the soil sieve were removed from the sample, washed, dried, and weighed, and the volume was estimated using the water-displacement method. Due to low quantities of root biomass, root volume was measured on all roots combined (by submerging the roots using a stone of known volume) which was then divided by total root weight, and multiplied with the specific root volume of each sample to give the volume. Soil and litter samples were homogenized into fine dust using a steel ball mill (Retch, Haag, Germany) with 30 revolutions min-1 for 3 min and analysed for total C and N using a dry combustion type elemental analyser (mod.NA2500; CE instruments Ltd, Wigan, UK). Note that carbonates are not generally present in these types of soils.

Supplementary Information for the Paper Titled