Effects of intensive biomass harvesting on forest soils in the Nordic countries and the UK: A meta-analysis

The use of biomass from forest harvesting residues or stumps for bioenergy has been increasing in the northern European region in the last decade. The present analysis is a regional review from Nordic and UK coniferous forests, focusing on the effects of whole-tree harvesting (WTH) or whole-tree thinning (WTT) and of WTH followed by stump removal (WTH + S) on the forest floor and mineral soil, and includes a wider array of chemistry data than other existing meta-analyses. All intensified treaments led to significant decreases of soil organic carbon (SOC) stock and total N stock in the forest floor (FF), but relative responses compared with stem-only harvesting were less consistent in the topsoil (TS) and no effects were detected in the subsoil (SS). Exchangeable P was reduced in the FF and TS both after WTT and WTH, but significant changes in exchangeable Ca, K, Mg and Zn depended on soil layer and treatment. WTH significantly lowered pH and base saturation (BS) in the FF, but without apparent changes in cation exchange capacity (CEC). The only significant WTH-effects in the SS were reductions in CEC and BS. Spruceand pine-dominated stands had comparable negative relative responses in the FF for most elements measured except Mg and for pH. Relative responses to intensified harvesting scaled positively with growing season temperature and precipitation for most variables, most strongly in FF, less in the TS, but almost never in the SS, but were negative for P and Al. The greater reduction in FF and TS for soil organic carbon after intensive harvesting decreased with time and meta-regression models predicted an average duration of 20–30 years, while many other chemical parameters generally showed linear effects for 30–45 years after intensified harvesting. Exchangeable acidity (EA), BS and pH all showed the reversed effect with time, i.e. an initial increase and then gradual decrease over 24–45 years. The subsoil never showed a significant temporal effect. Our results generally support greater reductions in nutrient concentrations, SOC and total N in forest soil after WTH compared with SOH in northern temperate and boreal forest ecosystems.


Introduction
In the northern European region, there is considerable interest in using biomass from forest harvesting residues or stumps for bioenergy. In conventional timber harvesting (stem-only harvesting including stemonly thinning, SOH), branches, tops and stumps are left in the forests.
Removal of these biomass compartments for bioenergy in whole-tree harvesting at the final felling (WTH), whole-tree thinning (WTT) and WTH followed by stump removal (WTH + S) may have consequences for the functioning of the forest ecosystem. As a large part of the nutrients in trees are in the foliage, twigs and branches, removing these tree components will reduce nutrient supply to the soil, which might in the long term increase the risk for nutrient imbalance and reduced forest production (Raulund-Rasmussen et al., 2008;Helmisaari et al., 2011;Tveite and Hanssen, 2013). This may lead to reduced future biomass carbon (C) sequestration and stocks in living forest biomass. In addition, there is a risk for reduction of deadwood and soil organic carbon (SOC) stocks when inputs of organic matter to soils decrease with removal of tree biomass during WTH that was previously left to decay on site. There has also been concern that stump harvesting might lead to reductions in SOC stocks, due to the missing input from stumps as well as by an enhanced decay following the soil disturbance (Persson, 2017).
Forest floor C stocks have been reported to respond mostly, but not always, with a larger decrease after intensification of biomass harvesting (Clarke et al., 2015). Several single site field studies have examined these questions and contrasting results have been found, likely due to variation in site-specific factors. With reference to the two most common tree species in northern Europe, Scots pine (Pinus sylvestris L.) is often planted on shallower and more nutrient-poor soils compared to Norway spruce (Picea abies (L.) Karst.) Thus, it may often not be possible to determine if an effect of WTH compared to SOH depends on either tree species or soil properties (or both).
A worldwide meta-analysis by Achat et al. (2015a) suggests that WTH is likely to lead to reduced contents of soil nutrients and slower growth compared with SOH. A further study (Achat et al., 2015b) found that WTH unlike SOH led to loss of SOC not just in the forest floor but also in the mineral soil when compared to no-harvest control. The metaanalyses on harvesting effects on SOC in temperate forest soils by Nave et al. (2010) and at global scale by Wan et al. (2018) suggest, on the other hand, that at least for SOC the effect in the mineral soil depends on soil type and texture. Similarly, in a recent global meta-analysis Wan et al. (2018) identified soil texture as an important co-variable for the SOC response. Hume et al. (2018) in their global meta-analysis found that SOC and N tended to decrease over time since harvest, while P did not, resulting in a higher forest floor N/P ratio.
Recently, there have been several meta-analyses published on the effects of intensified biomass harvesting in which the scale is worldwide (Achat et al., 2015a,b;James and Harrison, 2016;Hume et al., 2018;Wan et al., 2018). However, results obtained on a worldwide scale may not necessarily apply in the same way to all regions, as there are clear regional differences in factors such as climate and soil types. We suggest that there is still a need for more knowledge on the factors that control the observed differences at regional level, and of how variation in these factors affects long-term site productivity and SOC stocks. Northern Europe (here comprising the Nordic and Baltic countries and northern Great Britain, but not including Russia) is a highly relevant region due to the widespread use of residual forest biomass for energy, and a strong focus on intensive harvesting methods and their consequences in the past years in these countries (e.g. Clarke et al., 2015). Although many of the studies used by Achat et al. (2015a,b), Hume et al. (2018) and Wan et al. (2018) were from this region, it remains unclear to what extent results from their worldwide studies are applicable to northern Europe specifically.
We compiled available data from northern European field experiments on the effects of intensive forest biomass harvesting on SOC, nutrients and further soil chemical properties, including a few studies that were not used in the aformentioned meta-analyses. Additionally, this enabled us to constrain factors such as tree species, soil type and climate to a more limited set of categories. A minor part of the data for this meta-analysis was found in 'grey' literature (published in a national report) or was as yet unpublished, and thus not previously included in meta-analyses. Some authors of published studies also sent auxillary data from their studies upon request. The inclusion of all available data meeting the qualification criteria for inclusion in this study ('grey' as well as peer-reviewed literature) is important for meta-analyses to avoid flawed conclusions as a result of publication bias (Thornton and Lee, 2000).
The main hypotheses of the study were: (1) Intensive biomass harvesting methods used in northern European forests reduce soil nutrient contents; (2) Intensive biomass harvesting decreases forest SOC stocks which reduces the climate mitigation potential of production forests; (3) The most intensive harvesting methods in terms of the amounts of removed nutrients and C (clear-cut WTH and WTH + S compared to SOH) lead to greater reductions in soil nutrient and SOC stocks than WTT that is performed as part of forest management during stand development and where only part (<40%) of standing trees is typically removed; (4) Soils under Norway spruce (Picea abies (L.) Karst.) and Scots pine (Pinus sylvestris L.) have different sensitivity to the intensity of biomass harvesting; (5) Reductions in soil nutrients and SOC contents after intensive biomass harvesting diminish with time elapsed since harvest.

Data compilation
The data search was done consulting Web of Science, Google Scholar and websites from relevant research institutions, as well as by contacting corresponding authors directly. The following set of criteria had to be fulfilled in order to include a study in the meta-database: 1) sites had to be located in the Nordic-Baltic countries or northern Great Britain, i.e. at similar latitude; 2) studies should include pairs of comparable plots on the same site consisting of stem-only harvest (SOH) or stem-only thinning (also termed SOH) control plots and intensively harvested plots; 3) the intensively harvested plot should include at least one of three types of harvesting methods: whole-tree thinning (WTT), whole-tree clear-cut harvest (WTH) or whole-tree harvest followed by stump removal (WTH + S); 4) the study should provide measured data on one or several soil nutrients or carbon, forest floor thickness and depth of the mineral soil layers sampled, soil type, harvested tree species and time between harvesting and soil sampling. Additional recorded information included geographic (latitude, longitude and altitude) and climatic information: mean annual air temperature (MAT), mean annual precipitation (MAP) and mean air temperature and precipitation of four months representing the main growing season (T May-Aug and P May-Aug , respectively). The climate data were retrieved from the WorldClim Global Climate Database using the climatological normal period 1961-1990 and 30 arcseconds (~1 km) resolution (Hijmans et al., 2005). Available information on number of replicates and standard deviations, standard error estimates or P-values of the differences in measured soil nutrients, the nutrient extraction method and site preparation method was also collected. Data from intensified harvesting experiments in Estonia, Latvia and Lithuania did not meet all the above criteria. As a result, this study focuses on the following countries: Denmark, Finland, Norway, Sweden and the United Kingdom (Fig. 1). All the information in the meta-database was collected between 2012 and the spring of 2015. Site information is summarised in Table 1.

Data structure
Altogether our search yielded data from a total of 26 studies ( Table 1) that included 48 sites with paired plots, of which 13, 8, 2, 2 and 1 studies were in Sweden, Finland, Norway, the UK and Denmark, respectively (some studies included more than one site). The studies contained 34 to 203 (average 99) observations (k) for various properties from different soil layers. All selected studies made use of controlled experimental designs or paired plots to study the effects of intensified harvest on the various soil properties.
The meta-analysis was performed separately within three soil layers: 1) forest floors (FF): the organic horizon(s) consisting of undecomposed (L) or partly decomposed (F/H) organic material on top of the mineral soil; 2) topsoil (TS): generally, the mineral soil layer from 0 to 10 cm; 3) subsoil (SS): mineral soil layer >10 or 15 cm and down to 20 or 30 cm. Soil here refers to the fine-grained fraction, from which large organic or inorganic particles had been excluded, for example by sieving. The number of paired observations for each soil property in each layer is shown in Table 2. The maximum sampling depths and layering of the mineral soil differed amongst the studies and thus only some studies included subsoil layers. More than 90% of the data in the topsoil provided information to a depth of 10 cm; however, Olsson et al. (1996b) reported data for only one pooled mineral soil layer of 0-20 cm, which was included as a topsoil, and Egnell et al. (2015) and Vesterdal et al. (2002) reported topsoil layers down to 15 cm and subsoil below that depth. For the subsoil, approximately 90% of the observations provided data to a maximum depth of between 20 and 30 cm. The study of Egnell et al. (2015), however, provided four observations to a maximum of 70 cm which were included as subsoil, but since this study focused on the harvesting treatment WTH + S, these data were only included to calculate effect sizes from harvest treatment type (see Results section). Both concentration units and stock units were included in the database, and both types of unit were used in the analysis. Exceptions were (1) for SOC and total N in the forest floor where only stock units were used and (2) to find out whether there was a systematic difference between concentration units and stock units, where the data were divided up between the two unit types. For studies that reported results both in stock and concentration units, only the stock units were included. Consequently, for most elements stock data dominated the database (Table 2).
Only in a few cases were the variance and standard deviation (SD) estimated from P-values based on a standard normal distribution (Wald test). If the P-values were given as n.s. then a P-value of 0.05 was assumed and if they were given as < 0.05 or < 0.001, a value of 0.02 and 0.001 respectively was assumed and the same SD was then given to both treatments. When data for the forest floor (FF) were provided in the form of separate horizons (e.g. litter (L), fermentation (F) and humus (H) horizons), FF stocks were calculated by adding up the stocks of the respective layers to obtain a single value, whereas weighted averages of concentrations were calculated based on the thickness of each layer. If SD values were provided for the separate layers within the FF, the overall SD for the FF was based on the SD of the dominant layer. Stock estimates in different studies were based on measurements of soil bulk density (fine-grained fraction relative to total volume) rather than pedotransfer functions. SOC and total N (TN) were mostly determined using dry combustion, although in some studies loss on ignition was used for SOC (Wall, 2008;Walmsley et al., 2009) and the Kjeldahl method for TN (Wall, 2008). Soil pH was in almost all studies (except Vesterdal et al., 2002) determined in water. The extractants and digestion methods for P, Al, Ca, Mg, Na, K, Mn, and Zn used in the different studies varied, as did the method of chemical analysis. To avoid the risk of misinterpreting differences in results that were due to use of different extraction/digestion methods, we categorized the majority of these based on the strength of the extractant/acid used: NH 4 -based extractants, BaCl 2 , and strong acids. Available P was mostly determined in the same extracts as the exchangeable cations and the same procedure for categorization was followed, although the Sibbesen and Bray methods were used in addition in one study (Nykvist 2000). Because of the small number of studies within other categories, only the category of NH 4 -based extractants (giving information on "exchangeable" concentrations) was used in the meta-analysis for base cations, Al, Mn, Zn and P, as well as for exchangeable acidity (EA), base saturation (BS) and cation exchange capacity (CEC).
Most studies were carried out in pure stands of either Norway spruce or Scots pine, the most common tree species in the focal area. Thus, the studies that included data on mixed coniferous stands (k ≤ 36 pairs, Table 1) were excluded from the evaluation of tree species effects.
The minimum number of observations to be included in the calculation of an effect size was set to four (k ≥ 4), provided from at least two different studies and two different study sites. This prerequsite ensured that the calculated response ratios were based on a minimum degree of representativity.

Imputation of missing values
Multiple imputation was used to derive standard deviations for studies that did not report this (37% of all observations across all response variables). First, coefficients of variation were calculated for each observation with reported standard deviation, dividing the standard deviation by the reported mean for control and treatment, respectively. Missing coefficients of variation were then imputed from these by random sampling with replacement, either among all (when continuous predictor variables were tested) or a sub-set of observations within the same category level (when categorical predictor variables were tested). Finally, each imputed value was converted back to standard deviation by multiplying with the reported mean of the imputed observation, allowing all observations to be included in the following meta-analysis. This procedure was repeated 500 times, and final parameter estimates were obtained as the average across runs (Wiebe et al., 2006).

Effect size
Log-transformed response ratios were used as effect sizes to quantify the effect of harvesting intensity on each response variable (soil prop- Table 1 Information on sites and treatments in the meta-database. Where the same sites have been used in more than one study, these studies are grouped together and as far as possible listed consecutively in order of increasing time since harvesting.   17 (11) 16 (12) 16 (12) erty), having ln(RR) = ln(X I /X C ), where X I and X C are the mean values of the variables from the intensive (I) and control (C) harvest treatments, respectively. Control treatments were defined as stem-only harvest, including stem-only thinning (SOH). A negative ln(RR) indicated a negative effect of intensive harvesting on a given response variable compared to the control, whereas a positive ln(RR) indicated a positive effect. Assuming that harvest treatments were independent, the variance of ln(RR) was calculated as Hedges et al., 1999), where SD I and SD C are the standard deviations of the mean response values, and n I and n C their sample sizes.

Meta-analysis
For each response variable, overall meta-estimates were obtained from random effects meta-analysis models in which the calculated effect sizes were weighted with the inverse of their respective variances, giving more precise results a larger influence on the overall result (Hedges and Olkin, 1985). Furthermore, a range of mixed effects meta-regression models was used to evaluate (based on regression coefficients) the change in effect size against different a priori selected predictor variables (moderators, see Viechtbauer, 2010, for details). Where appropriate, multiple meta-regression was used to compare the influence of continuous predictor variables among categorical groups. Predictor variables tested included time (years since last harvest), climate variables (growing season and annual precipitation and temperature), longitude and latitude. Uncertainty in the regression coefficients was quantified, using 95% confidence intervals. The level of residual variation when fitting each meta-regression model was used to estimate the degree to which the predictor variables could explain differences among studies, applying restricted maximum-likelihood (REML) estimation of a Q statistic (Viechtbauer, 2007). In order to provide meta-estimates more clearly, ln(RR) was back-transformed (as e ln(RR) ) and expressed as percentage difference at intensive biomass harvesting compared to the less intensive alternatives.

Modelling of non-independence
One of the central assumptions in the statistical framework of metaanalysis is that effect sizes are mutually independent. Conflicting with this assumption, common sampling designs in forest studies introduce correlation because effect sizes may be spatially clustered, e.g. within chronosequences, paired-plot designs, and repeated sampling designs. Some of the studies in our database resampled the plots over a period following harvest or in other ways included multiple observations (sublayers) within the same soil layer, thus providing correlated measurements. To account for this interdependence, we used a modified statistical framework, placing a covariance term in the appropriate diagonal entries of the covariance matrix and using duration between measurements as a scaled linear distance (Lajeunesse, 2011). As a result, multiple (repeated) effect sizes from the same study were downweighted, reducing their influence on the overall result.
All analyses were run in R, version 3.2.1, using matrix notation as implemented in the metaphor Package (Viechtbauer, 2010). R code is available upon request.

Effects of intensive harvesting on SOC and total N, exchangeable elements, pH, base saturation, cation exchange capacity and exchangeable acidity
The mean response ratios for SOC stocks differed significantly from zero for the three intensive harvesting methods (WTT, WTH and WTH + S) (Fig. 2). In the forest floor, a significantly lower level in mean C stocks was observed with all the intensive harvesting treatments, with the highest mean deviation of − 12% for WTH + S (normal scale, as for the subsequently given values). In the topsoil, mean SOC stocks or concentrations were lower by − 21% for WTH; on the other hand, no significant effect of WTH + S was observed on the topsoil SOC. For WTT there was a small but significant positive difference (5%) in the mean SOC in the topsoil. No significant effects of intensive harvesting were found in the subsoil SOC.
For soil TN stocks, significantly lower levels were found in the forest floor after WTH and WTT compared with SOH, amounting to a mean of − 10% and − 4%, respectively (Fig. 2). No significant effects were observed in TN stocks or concentrations in the topsoil or the subsoil following WTH, while, as for SOC, there was a modest significant positive difference (mean of 7%) after WTT in the topsoil. The WTH + S treatment had too few observations to evaluate the response in TN.
In addition to SOC and TN, we found significant effects of WTH and WTT on exchangeable soil nutrients compared with SOH (Fig. 3). In the forest floor, significantly lower mean levels of exchangeable concentrations were observed for P, K, Ca, and Mg for WTT, whereas P, Ca and Zn were significantly lower for WTH relative to SOH. The largest relative negative difference was observed for P after WTT (mean of − 30%). Al was the only element that was significantly increased in the forest floor layer following WTH relative to SOH.
In the topsoil, P was at a significantly lower level after WTH and WTT relative to SOH and this effect was of similar magnitude in both treatments (mean of − 16% and − 20%, respectively; Fig. 3). Other elements that showed significantly lower levels in the topsoil for WTH relative to SOH were K, Ca, and Mn. Al increased significantly in both the forest floor and the topsoil after WTH. While Ca concentrations were at lower Fig. 2. Log-transformed response ratios, ln(RR), of soil organic carbon (SOC) and soil total nitrogen (TN) in the forest floor, topsoil, and subsoil layers following whole tree thinning (WTT), whole tree harvesting (WTH), and WTH with stump removal (WTH + S) as compared to conventional stem-only thinning or harvesting (SOH). Enough observations for the WTH + S treatment were only available for SOC. Significance levels are indicated as ns (P ≥ 0.10), (ns) (P = 0.05-0.10), * (P < 0.05), ** (P < 0.01), *** (P < 0.001) and number of observations included is shown in parentheses. levels in the topsoil after WTH (mean of − 13%), they were greater by a mean of 24% for WTT. On the other hand, a significantly lower level of Mg was observed in the topsoil for WTT.
Finally, for the subsoil, no significant differences were detected for any of the intensive harvesting treatments. As was the case for TN, there were not enough observations for testing the effects of the WTH + S treatment, as well as for Zn, Al and Mn for WTT and for P in the subsoil for all treatments.
The only significant response in pH was found in the forest floor, where pH was significantly lower for WTH compared to SOH (mean of − 2%, Fig. 4). This coincided with a significant increase in EA (mean of 16%) and a significant reduction in BS (mean of − 9%) in forest floors after WTH. CEC did not change significantly as a response to WTH. In the topsoil, both EA and BS significantly decreased (mean of − 21% and − 19%, respectively) following WTH compared with SOH; however, a similar reduction was not observed for CEC (Fig. 4). Lastly, in the subsoil, significant negative log-transformed response ratios were observed for BS (mean of − 10%) and CEC (mean of − 18%).

Differences in intensive harvesting response ratios in forest floors calculated from element stocks versus concentrations
To check if there was systematic bias in observed effects for different studies that expressed elements as concentrations compared to stocks, we compared response ratios for the forest floor following WTT or WTH, in cases where k ≥ 6 for both units (Fig. 5). Different studies based on stocks or concentrations showed similar mean responses for Ca, Zn, Mn and Al. Of these elements, Zn and Mn were dominated by concentration Fig. 3. Log-transformed response ratios, ln(RR), of exchangeable nutrients and Al in forest floor, topsoil and subsoil layers following whole-tree thinning (WTT) and whole-tree harvesting at final felling (WTH) as compared to conventional stem-only thinning or harvesting (SOH). Significance levels are indicated as ns (P ≥ 0.05), * (P < 0.05), ** (P < 0.01), *** (P < 0.001) and number of observations included is shown in parentheses. For P, too few observations were available for analysing in the subsoil. For Zn, Al and Mn, response ratios were only calculated for WTH, as there were too few observations to calculate the effect of the WTT treatment. Fig. 4. Log-transformed response ratios ln(RR) of pH, exchangeable acidity (EA), cation exchange capacity (CEC) and base saturation (BS) in forest floor, topsoil and subsoil layers following whole tree thinning (WTT) and harvesting (WTH) as compared to conventional stem-only thinning or harvesting (SOH). Significance levels are indicated as ns (P ≥ 0.05), * (P < 0.05), ** (P < 0.01), *** (P < 0.001). For EA, BS and CEC there were only enough observations available for testing for the WTH treatment. For pH, enough observations for testing for the subsoil were only available for the WTH treatment.
units, Ca was dominated by stock units and Al had almost the same number of observations in each unit category (Table 2). Only for K was there a significantly larger negative difference in studies which reported stocks of exchangeable nutrients, while the few (6) studies that reported concentration differences did not show significant effects following WTT or WTH. For P, Mg and Na there were not enough studies to compare effect of units.

Effects of tree species
Effects of intensive harvesting, WTT, WTH and WTH + S jointly, on SOC, TN, exchangeable base cations and pH were calculated for the forest floor for stands of Norway spruce and Scots pine (Fig. 6). Most of the investigated soil nutrients showed similar significant negative logtransformed response ratios in pure stands of both Norway spruce or Scots pine. The exception was Mg, for which the response did not differ significantly from zero in Norway spruce stands, while Scots pine stands showed a significant decrease (mean of − 14%). Negative differences in ln(RR) of the remaining elements analysed ranged from − 7% to − 11% (Fig. 6). Only forest floor pH showed a significantly different (P < 0.001) response between the tree species, indicating a very small, but significant, positive difference after intensive harvesting (mean of 1%) for Scots pine stands, and a significantly negative difference in Norway spruce stands (mean of − 2%).

Effects of growing season air temperature and precipitation
Both mean air temperature over the growing season (T May-Aug ), defined as the four warmest months in a year, and growing season precipitation (P May-Aug ) had a higher number of significant regression relationships with mean ln(RR) responses in the forest floor and topsoil than did annual temperature or precipitation (data not shown). Climate variables never showed significant correlations with mean subsoil ln (RR)s (data not shown). Forest floor (Table 3) and topsoil (not shown) ln (RR)s showed similar correlations to climate variables, but the level of significance was generally greater for the forest floor (data not shown).
Meta-estimates of SOC, TN, Ca, BS and pH in the forest floor were significantly negatively correlated with growing season temperature (T May-Aug ; ranging from 10.9 • C to 14.6 • C; Table 3), suggesting stronger effects of intensified harvesting in warmer climates. The reverse was found for P and Al; their mean ln(RR) was positively correlated to T May-Aug . Al became therefore more abundant in the forest floor in warmer climates following intensified harvesting compared with SOH. Mean ln (RR) for other soil properties did not significantly correlate with T May-Aug (Table 3).
after intensified harvesting compared with SOH with higher P May-Aug . Meta-estimates of Zn and pH, on the other hand, increased significantly with increasing P May-Aug (Table 3). Slope values for P May-Aug were very low even when significant, and always lower than for T May-Aug (Table 3).

Effect of time since harvesting
Regressions were used to investigate the effect of time elapsed since harvesting on the log-transformed response ratios, ln(RR), for the different soil properties across all intensified harvest treaments. A nonsignificant slope means that ln(RR) did not correlate with time elapsed since harvesting (maximum 35-38 years) (Fig. 7, Table 4).
For SOC, ln(RR) of the forest floor increased slightly with time since harvest (Fig. 7, Table 4), and even more for the SOC in topsoil, indicating greater reductions after intensive harvesting compared with SOH in the short term that approach zero in the longer term (Fig. 7, Table 4). In the subsoil the slope was not signficantly different from zero.
As was the case for SOC, the slopes were significantly positive for P, K, Ca, Mg, Zn and BS in the forest floor, with intercepts that were negative and significantly different from zero (Table 4). Only exchangeable acidity (EA) had a significant negative slope in the forest floor, and a positive intercept significantly different from zero (Table 4). The regression models indicated that it might take 39,42,37,34,41,37 and 34 years to reach ln(RR) = 0 (meaning no difference between effects of SOH and WTH) after intensified harvesting for P, K, Ca, Mg, Zn, EA and BS, respectively. The intercepts were significantly negative for TN and significantly positive for Al, while the slopes were not significantly different from zero, which could indicate no or slower temporal change (Table 4).
In the topsoil, EA, BS and pH had a significantly positive slope and significantly negative intercept, indicating lower values for intensive harvesting in the short term, but less so in the longer term (Table 4). The regression model predicted that it would take 34, 45 and 24 years to reach ln(RR) = 0 for EA, BS and pH, respectively. Ca in the topsoil responded oppositely to Ca in the forest floor, as it had a positive intercept (increase following intensified harvesting) and a negative slope, indicating that intensified harvest caused a significant increase of Ca in the topsoil in the short term, but not in the longer term, with the regression model reaching ln(RR) = 0 at 37 years after harvesting. Ln (RR) for K did not have a significant slope in the topsoil, but still had a significantly negative intercept (Table 4).
There were no significant effects of time since harvest on the ln(RR) for the soil parameters of the subsoil (data not shown).

Effects of intensive biomass harvesting methods on SOC, TN, exchangeable elements, pH, base saturation, cation exchange capacity and exchangeable acidity
Our results generally support a greater reduction in soil nutrient concentrations and SOC in northern European coniferous forest ecosystems after intensified harvesting compared with SOH, which supported our first two hypotheses. This is in general agreement with the results of the meta-analyses of Achat et al. (2015a,b) and Wan et al. (2018), but partly contrasts with findings of other meta-analyses that have been based on world-wide (Hume et al., 2018;James and Harrison, 2016) or biome-wide datasets (Nave et al., 2010). Achat et al. (2015a,b) found that intensive harvesting (residue removal) resulted in SOC losses in all the investigated soil horizons (forest floor and mineral soil above and below 20 cm). Results for TN were similar to those for SOC (Achat et al., 2015a). Similarly, Wan et al. (2018) found in a worldwide meta-analysis that the retention of harvest Table 3 Regression slope parameters of meta-regressions of the log-transformed response ratios, ln(RR), for SOC and TN stocks, exchangeable element concentrations, exchangeable acidity (EA), cation exchange capacity (CEC), base saturation (BS) and pH in the forest floor against mean growing season temperature or mean growing season precipitation for the period May to August (T May-Aug , • C, and P May-Aug , mm, respectively), comparing intensive harvesting (whole tree thinning, WTT, whole tree harvesting, WTH, and WTH + stump harvesting, WTH + S) with conventional stem-only thinning or harvesting (SOH). Significance levels are indicated as * (P < 0.05), ** (P < 0.01), *** (P < 0.001); NS = not significant; N indicates the number of observations included in the analysis.  Fig. 7. Meta-regressions of ln(RR) with time since harvest for SOC in forest floor, topsoil and subsoil, respectively. Numbers refer to different studies (see Table 1) and circle diameters reflect the weight of the estimate in the meta-analysis. residues in SOH led to 8.2% greater soil C stocks in the 0-20 cm layer of the mineral soils compared to WTH. This contrasts with the findings of James and Harrison (2016), who found in their world-wide meta-analysis that mineral soils had 13.3% more soil C after WTH compared to SOH. This unexpected result was attributed to reduced leaching of dissolved organic carbon and nutrients after WTH reducing the priming of mineral soil SOC mineralisation. An alternative explanation was that WTH and SOH treatments were possibly on average on different soil types, as the SOH studies were ca. 3 times as many as the WTH studies (James and Harrison, 2016). Nave et al. (2010) found no additional effects of harvest intensity on soil C in their meta-analysis when comparing WTH and SOH to a no-harvest control on forest floor in the temperate zone and neither did Hume et al. (2018) for forest floor and mineral soil C; they too mostly used studies from the temperate zone.

Effects of soil depth
Generally, there were more significant effects for the forest floor than for the mineral topsoil and subsoil. For some element concentrations and stocks (e.g. Ca or Mg after WTT) this may reflect a larger number of studies with data from the upper soil horizons. For other chemical properties, or more generally, it may reflect a decreasing effect of intensive biomass harvesting with increasing soil depth.
Of all soil layers, the forest floor showed the greatest reduction in the SOC and TN stocks with increasing intensity of biomass harvesting. The large relative SOC and TN losses in the forest floor and the correspondingly stable SOC and TN stocks in the topsoil of the WTH + S treatment suggest a potential mixing of forest floor material into the upper mineral soil during the stump extraction procedure, since WTH without stump harvesting surprisingly seemed to have a larger impact on the topsoil SOC than WTH + S. Additional explanations may be that the WTH and WTH + S studies were conducted on different soil types, with different pre-treatment management regime, including intensity of site preparation, and/or climatic conditions. In a Swedish experiment, stump and root harvesting after clearcutting resulted in a significantly lower soil C pool in the FH layer and a significantly lower annual heterotrophic respiration in the whole soil profile 20-30 years after stump harvesting compared to patch scarification . Comparing effects of WTH and WTH + S in Finland 11-12 years after final felling showed that soil C and TN pools had a tendency (not significant) to be lower following stump harvesting, while in situ SOC mineralization rates did not significantly differ between the same treatments . A study of effects of WTH and WTH + S 8-13 years after final felling in Finland found no significant effect on soil C or TN stocks, while soil pH was slightly higher after stump harvesting . A comparison of effects of WTH and WTH + S on two soil types in the UK four years after harvesting found that stump harvesting led to reductions in soil C and TN concentrations (or stocks) as well as base cations in both brown podzolised soil and peaty gley soil types, although in the podzolized brown soil some of the reductions in soil C, TN and base cations in the topsoil were compensated by increases at depth, indicating a downwards redistribution of organic matter (Vanguelova et al, 2017).
In most cases, WTH/WTT led to a significant reduction in exchangeable concentrations compared with SOH in our study, with some elements, such as P, often responding most drastically. This pattern was also found in other studies, including by Achat et al. (2015a), while such a drastic reduction in soil P after WTH was not reported from the global meta-analysis of harvesting effects by Hume et al. (2018).
Surprisingly, we found that while Ca concentrations decreased in the topsoil after WTH (-13% change), they increased by 24% following WTT. This was similar to a pattern observed for SOC and TN, although more pronounced. Possibly thinnings open the canopy, leading to higher temperatures on the forest floor and increasing decomposition rates, but without the leaching loss often seen after final felling due to the existence of a tree cover with high nutrient demands. In addition, higher moisture might increase microbial activity resulting in increased decomposition of mineralized nutrients.
The forest floors of WTH plots were more acid compared to those of SOH plots (decrease in pH and BS, increase in EA). The differences in EA and BS support the above general findings of a higher level of exchangeable Al and a lower level of exchangeable nutrients after intensive harvesting. In the topsoil, both EA and BS were significantly lower for WTH compared to SOH. However, a similar reduction was not observed for CEC. In the subsoil, significant decreases after WTH compared with SOH were observed for BS and CEC, although there was no significant effect on SOC, TN, pH, EA or exchangeable nutrients. Achat et al. (2015a) also found some soil acidification after residue removal; however, many of their data were from northern boreal forests, so similarity with our results might be expected.

Effects of choice of units (stocks vs. concentrations)
For P and K in the forest floor, the log-transformed response ratio was significantly more negative when stock units rather than concentration units were used, but for other soil chemical parameters (Ca, Mg, Zn, Mn and Al) there was no significant difference. Although choice of units for Table 4 Regression parameters of a regression of the log-transformed response ratios, ln(RR), and time since harvesting for SOC, soil exchangeable elements, exchangeable acidity (EA), cation exchange capacity (CEC), base saturation (BS) and pH for intensive harvesting (whole tree thinning, WTT, or whole tree harvesting, WTH and WTH + stump harvesting, WTH + S) in forest floor and topsoil. Significance levels are indicated as * (P < 0.05), ** (P < 0.01), *** (P < 0.001). NS = not significant. N = total number of paired plots (k) used for the analysis. Their distribution to intensive harvest types can be seen in Figs. 2  reporting results may thus in some cases affect the results obtained, the possible effect of variability due to few studies in either group needs to be considered. Stock estimates were based on measurements of soil bulk density rather than pedotransfer functions, and the calculation of soil bulk density is one of the key parameters affecting SOC stock estimates at the plot as well as the landscape level (Mehler et al., 2014;Throop et al., 2012;Walter et al., 2016). However, it is unclear how much uncertainty in the determination of bulk density might have affected the stock estimates used.
According to Nave et al. (2010), the reporting units used (concentration vs. pool size) are one of the two most significant categorical factors accounting for among-study variation. However, contrary to the findings of Nave et al. (2010), we found that the overall differences between the response ratios reported based on concentrations and stocks were minor. In the current study, stock data were only included for SOC and TN in the forest floor, where the effect of units is expected to be largest, as the thickness of the layer is inherently affected by the treatment. Less difference must be expected for the mineral layers where thickness is expected to stay the same for fixed depths, assuming the effect of the harvesting on soil compaction is the same for both treatments. The dominance of data given in concentration units for Mn and Zn in the forest floor, compared to the other elements where stock units dominate, suggests that care might need to be taken when comparing the responses of these elements in the forest floor.

Effects of tree species
Meta-estimates, ln(RR), differed significantly between Norway spruce and Scots pine only for pH. As different tree species are grown under different site conditions, it is unclear whether this difference is due to tree species per se or to other site-related differences such as soil type, texture and fertility. Since both species were conifers, they may have a similar build-up of organic matter in the forest floor, which was the only soil layer analysed for species effects. This might be the reason that we did not see a stronger tree species effect.

Effects of growing season air temperature and precipitation
The effects of both growing season temperature and precipitation on the log-transformed response ratio varied greatly, as both positive and negative correlations were found. This agrees with Thiffault et al. (2011) and Wan et al. (2018), who found no clear climate-related effect of WTH on soil nutrient pools or soil C relative to the effects of SOH. Their range of climatic conditions was even wider than ours. However, our logtransformed response ratios for SOC, TN, Ca, BS and pH in the forest floor decreased significantly with increasing growing season mean temperature, suggesting greater reductions after WTH relative to SOH in a warmer climate. This agrees with Achat et al. (2015a), who found that the decrease in SOC and TN after WTH was generally higher under temperate conditions than colder conditions, although for TN the difference between climate classes was not significant.

Effects of time since harvest
Our results indicated that some soil nutrients may recover faster with time since harvest than others. Topsoil concentrations of intensively harvested and SOH plots did not approach each other in any significant way within a time span of 37 years, indicating that a longer time is needed for this soil layer to recover. The subsoil showed no significant time effects. Significant positive slopes for P, K, Ca, Mg, and Zn in the forest floor are consistent with an initial loss and gradual replenishment as time passes after harvest. As opposed to recovery of soil C and N, soil P recovery has been reported as slower after harvesting vs. no harvesting at the global level (Hume et al., 2018), indicating a decoupling of the P cycle from those of C and N after harvesting. Morris et al. (2019) also reported from Canadian sites in the North American Long-Term Soil Productivity network, that P recovery was slower and remained lower in both SOH and WTH 20 years after harvesting, as opposed to soil C, Kjeldahl N, K, extractable Ca and extractable Mg. Our study of intensified harvesting in northern Europe suggested that after WTH P approached levels comparable to SOH after ca. 40 years in the forest floor, while there was no significantly negative effect on the rate of change in P in the topsoil in the temporal analysis (Table 4). A significant positive intercept and negative slope detected in the forest floor for EA is consistent with an initial enhancement of acidity that diminished during the next 35 years after harvesting. The positive differences in the log-transformed response ratios with time suggest that the difference between effects of SOH and WTH will decrease with time, and that these differences will be more apparent in the forest floor than the mineral soil. Compared with other studies of harvesting, as such, the additional temporal impact of intensified biomass removal appears less problematic in a stand rotation perspective. For SOC in the forest floor, a significant positive slope was observed, even if small and comparable to a similar uncertain trend found by Clarke et al. (2015). Achat et al. (2015a) found that negative differences in TN stock and exchangeable K and Mg in the forest floor tended to be larger in the first ten years after removal of harvesting residues than later, which also suggests recovery with time. Our results are generally in line with the quite fast recovery of soil elements reported from North America (Morris et al., 2019). Studies of available pools of Ca, Mg and K in coarse sandy forest soils in Denmark indicated site-specific potential nutrient deficiencies related to intensive biomass removal (Vejre, 1999). Like other metaanalyses encompassing time trends, the results do not include repeated measurements at the same locations but reflect studies at different locations with different climatic and edaphic conditions which will affect both the intercept and the slope for the various elements.
The impacts on soil parameters found in this and other studies indicate a risk of reduced productivity after WTH or WTT (Vesterdal et al., 2002), which has been confirmed from a series of whole-tree harvesting experiments in Finland, Sweden and Norway (Helmisaari et al., 2011) and from another series in Norway (Tveite and Hanssen, 2013), but not for two experiments from Denmark (Nord-Larsen, 2002).

Conclusions
Our results generally support greater reductions in nutrient concentrations, SOC and TN after WTH compared with SOH in northern European temperate and boreal forest soils, consistent with the results obtained on a worldwide scale. Effects were greater in the forest floor than in the mineral soil, and greater in the topsoil than the subsoil. Spruce-and pine-dominated stands had for most elements comparable negative relative responses in the forest floor. There appeared to be greater effects of WTH relative to SOH in a warmer climate. The differences between effects of different harvest types in the forest floor and topsoil were generally reduced with time but were likely to last for several decades.

Author statement
This study was initiated by NC, IS, LF, BDS, KA, DL and SJ. The metadatabase was prepared by HMS and the metaanalysis carried out by LPK, with support from LV, BDS and NC. Interpretation of the results was done by OJK, TGB, LV, BDS and NC. All authors provided data for the metadatabase and participated in the writing and revision of the manuscript.

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.