Long-Term Soil Fertility and Site Productivity in Stem-Only and Whole-Tree Harvested Stands in Boreal Forest of Quebec (Canada)

: Using residual biomass from forest harvesting to produce energy is viewed increasingly as a means to reduce fossil fuel consumption. However, the impact such practices on soil and future site productivity remains a major concern. We revisited 196 forest plots that were subject to either whole-tree (WTH) or stem-only (SOH) harvesting 30 years ago in the boreal forest in Quebec, Canada. Plots were stratiﬁed by four soil regions grouped by so-called ‘soil provinces’. Soil analyses indicated that after 30 years, the forest ﬂoor of WTH sites had smaller pools of N ( − 8%), exchangeable Ca ( − 6%) and exchangeable Mn ( − 21%) and a higher C/N ratio (+12%) than that of SOH sites. Mineral soil responses to the two harvesting intensities differed among soil provinces. In the two coarse-textured granitic soil provinces, organic matter, organic carbon, and nitrogen pools over the whole solum (0–60 cm soil depth) were at least 28% smaller after WTH than after SOH. Site productivity indicators followed differences between soils and were lower after WTH than after SOH in the two granitic soil provinces. The study shows that soil characteristics greatly inﬂuence a soil’s sensitivity to increased forest biomass harvesting in the long term. 0.018), but not in the deeper soil layers ( p ≥ 0.133). Finally, no practical difference in soil pH and statistical differences in exchangeable Ca and Mg pools could be detected between harvesting treatments in the soil provinces at any soil depth ( p ≥ 0.292).


Introduction
The development of bioenergy and renewable energy technologies to shift away from the use of high-emission fossil fuels is key to achieving reductions in greenhouse gas emissions. The various transition pathways explored so far to limit the rise of global mean temperature abundantly rely on bioenergy due to its multiple roles in the decarbonization of energy [1]. The bioenergy industry is still in its early stages in North America. For instance, residual biomass generated by forest harvesting in public forests in Quebec amounted to 8779 green Mg in 2018-2019, of which only 8.6% was attributed to users [2]. Nonetheless, forest biomass converted into bioenergy primarily for the targeted substitution of greenhouse gas-intensive materials or energy is part of various scenarios to mitigate the rise of atmospheric CO 2 , by the Quebec forestry sector for example [3]. Among the methods used to harvest wood, whole-tree harvesting (WTH), which encompasses several methods in which all parts of the tree above the stump are harvested, and stem-only harvesting (SOH) in which only the merchantable stems are harvested, are generally used to analyze the effects of more intense harvesting, such as those for biomass production. However, as high concentrations of nutrients are found in branches, tops, and foliage, harvesting most or all parts of the tree above the stump raises concerns about the maintenance of soil fertility and site productivity over the long term [4].
Many research articles have reported the impacts of WTH on site productivity [5][6][7][8][9]. However, most results are based on short-term analyses over 5 years or less [10]. Literature reviews and meta-analyses have shown that the increase in biomass removal has limited or only short-term effects on soil chemistry or site productivity [5,[11][12][13]. Long-term experiments in Sweden showed that WTH led to smaller pools of exchangeable base cations in soils compared to SOH [14,15]. Other studies in Sweden and the United States found no evidence that WTH reduced pools of organic matter, C, or N in soils for up to 30 years after harvesting [7,16,17]. In Canada, Morris et al. [18] found no significant differences in C, N, P, and K pools in the first 20 cm of sandy soils, compared to pre-harvest levels, for SOH and WTH, 20 years after treatment; they also found no evidence of decline in tree growth after the two harvest treatments. Short-term and long-term studies may reveal different responses, since tree nutrient requirements are relatively low during the first decade after regeneration, particularly for conifer species [19].
Soil type, climate, and tree species composition appear to be critical determinants of site sensitivity to biomass harvesting [5,20,21]. A first classification of forest ecological types considered sensitive to WTH, based mainly on the concept of critical load of acidity, was used in Quebec's state forest sustainability policy framework [22]. Critical loads of acidity depend on the long-term steady-state geochemical budget of base cations, i.e., inputs (atmospheric depositions, soil chemical weathering) vs. outputs (alkaline leaching, exports through biomass harvesting). The critical load model has been successfully used to inform and guide international negotiations aimed at reducing emissions of atmospheric acidifying pollutants across Europe and North America [23]. In the context of biomass harvesting, the critical load thresholds are considered to be exceeded when the long-term nutrient exports by WTH cannot be compensated by the inflows from atmospheric depositions and soil chemical weathering [24,25]. This concept, however, has been criticized for the static nature of its calculation and because it does not address the issue of the internal nutrient cycling of biochemical and biogeochemical ecosystems [26]. So far, the focus has been mainly in the long-term management of base cations and the prevention of soil damage by compaction or erosion. Little attention has been paid to the dynamics of ecosystem element pools, such as soil organic carbon pools [21]. Long-term changes in soil carbon after WTH or SOH are not well characterized, especially in deeper mineral soil horizons. Deeper soil horizons are rarely monitored due to the effort required for sampling, and because it is assumed that mineral soil organic C is old and stable. However, evidence exists that deeper soil layers can respond dramatically to ecosystem disturbances [27,28].
In this study, we revisited 196 forest plots harvested 30 years ago by either WTH or SOH in the boreal forest of Quebec, Canada. Our objective was to study how soil element pools and tree productivity responded to the two harvest treatments. The province of Quebec contains over 27.8 Mha of public-managed forest land, with an allowable annual harvest of 32.7 Mm 3 ·yr −1 [2]. Black spruce (Picea mariana (Mill.) B.S.P.), balsam fir (Abies balsamea (L.) Mill.), and Jack pine (Pinus banksiana Lamb.) are the main commercial tree species in this boreal forest. The following hypotheses were tested: • WTH sites have smaller soil element pools than SOH sites; • Soil responses to biomass harvesting intensity differ among soil regions; • Forest stands grow less following WTH than following SOH.

Study Area
The study area includes mainly the balsam fir-white birch bioclimatic domain and the southern part of the black spruce-feather moss bioclimatic domain, from east to west in the province of Quebec, Canada, between 47 • and 50 • of latitude N ( Figure 1). The study plots were located in four of the five soil provinces recognized in southern Quebec [29]: the Appalachians (B), the Laurentians (C), the Abitibi and James Bay Lowlands (D; hereafter referred to Abitibi Lowlands for concision), and the Mistassini Highlands (E). These soil provinces can be distinguished according to main parent material (geology, geomorphology), altitude (e.g., areas invaded by post-glacial seas), topography (slopes, landforms), soil texture, vegetation, and climate (temperature, rainfall, degree-days; see Table 1). The main types of late-successional vegetation throughout the study area are the balsam fir-white birch forest and black spruce-feather moss forest. Major tree species in the study plots are balsam fir, black spruce, jack pine, white spruce (Picea glauca (Moench) Voss.), and white the Appalachians (B), the Laurentians (C), the Abitibi and James Bay Lowlands (D; hereafter referred to Abitibi Lowlands for concision), and the Mistassini Highlands (E). These soil provinces can be distinguished according to main parent material (geology, geomorphology), altitude (e.g., areas invaded by post-glacial seas), topography (slopes, landforms), soil texture, vegetation, and climate (temperature, rainfall, degree-days; see Table  1). The main types of late-successional vegetation throughout the study area are the balsam fir-white birch forest and black spruce-feather moss forest. Major tree species in the study plots are balsam fir, black spruce, jack pine, white spruce (Picea glauca (Moench) Voss.), and white birch (Betula papyrifera Marshall) (see detailed composition of the stands sampled by soil province in the Supplementary Material, Table S1).

Experimental Design
From 1982 to 1988, 562 circular plots (0.04 ha each, 11.28 m radius) were established just before clear-cut harvest throughout the boreal forest of Quebec. In 2011, we selected 196 of these plots (sites) based on their accessibility and their similarities regarding preharvest stand composition (black spruce-balsam fir stands), surface deposit types and depth, and drainage class among the four soil provinces. Trees were harvested either by WTH (feller-buncher, mechanical harvester and cable skidder, or manual felling and cable skidder) or SOH systems (manual or mechanical felling and delimbing onsite, and cable skidder). No differences in the natural regeneration composition and abundance were observed between WTH and SOH plots up to 10 years after harvesting operations [33].

Field Sampling and Measurements
Starting 20 years after treatment, plots were measured every five years by staff of Quebec's Ministère des Forêts, de la Faune et des Parcs. Height measurements were taken on the same 20 trees in each plot. These trees were selected as representative of the observed diameter distribution of each species at the moment of the first plot survey. In all, two to four height measurements per tree were available over the years.
The selected plots were also revisited for tree measurements and soil sampling from 2011 to 2013, a time that represents a median period of 30 years (interquartile range: 2 years). For this survey, five to eight dominant trees in the study plots were measured for their diameter at breast height (DBH) and total height (Supplementary Material Table S2). Two cores per tree were also taken to assess individual tree radial growth and site quality index (SQI). Overall, these trees had a DBH of 12.7 ± 2.7 cm (mean ± SD), a height of 8.7 ± 1.8 m and were 22.4 ± 5.5 years old at DBH height. Soils were also sampled at this time. The forest floor was sampled quantitatively at 12 equidistant points (subsamples every 30 gon, or approximately 5.9 m) around each plot using a bipartite 8 cm volumetric root auger. Mineral soil was also sampled quantitatively at each of the four cardinal points (subsamples every 90 gon, or approximately 16 m) at depth intervals of approximately 15 cm, down to 60 cm. The same root auger was used for all samples, and sampling depth was recorded for each sample. When it was not possible to quantitatively sample a soil layer, we used a standard Edelman auger. We also assessed soil drainage class using provincial standards [34]. In all, 16 plots were classified as having excessive drainage (class 1), 96 as well drained (class 2), 57 as moderately well drained (class 3) and 27 as imperfectly drained (classes 4 and 5).

Laboratory Analyses
All soils were kept frozen (−18 • C) until preparation for laboratory analysis. Forest floor and soil samples were air-dried for at least four days, after which the quantitative samples were weighed individually. Forest floor samples were passed through a Wiley mill, and the 12 subsamples per plot were mixed thoroughly to form a composite sample for chemical analysis. The air-dried forest floor samples and all mineral soil subsamples were passed through a 2 mm mesh sieve. The fine fraction was analyzed for remaining humidity, organic matter by loss upon ignition, pH (soil:water 1:2 weight/volume), exchangeable cations (Ca, Mg, K, Al, acidity, cation exchange capacity (CEC), base saturation (BS); extraction with NH 4 Cl 1 N, 12 h). A subsample of the fine fraction was further ground to 250 µm for the determination of total C and N by dry combustion (LECO CR-412, LECO Corporation, St. Joseph, MI, USA). The forest floor samples were further analyzed for total P, K, Ca, and Mg concentrations by digestion in concentrated H 2 SO 4 followed by measurement by inductively coupled plasma emission spectrophotometry (ICP-AES). The fraction >2 mm of the quantitative mineral soil samples was weighed in order to determine the proportion of coarser fragments (f).
Tree cores were dried, glued to wooden holders, and sanded according to the standard procedure [35]. Tree age was determined after detection of ring boundaries under binocular magnification. Ring widths were measured using WinDENDRO ™ software (version 6.1D, Regent Instruments Inc., Quebec City, QC, Canada), and validated using signature rings to assist crossdating [36]. Visual crossdating was done through the recognition of patterns of wide and narrow rings common to sites within a given soil province [37]. Each set of raw tree ring measurements was evaluated using the COFECHA computer program [38] to ensure proper crossdating. Ring width values were converted to basal area increment (BAI, cm 2 ·yr −1 ) using the bai.out function of the dplR R package [39], version 1.6, in version 3.6.1 of the R software environment [40]. We selected BAI because it better reflects tree volume increment than does diameter increment [41].

Computations
Biomass and mineralomass (N, P, K, Ca, and Mg) for initial stands, exports, and what remained on the ground (for SOH only, as we assumed no above-ground tree biomass was left with WTH) were estimated from the pre-harvest forest surveys using the stand-scale biomass and nutrients calculator for Canadian forests [42]. The inputs for this calculator are the total and individual tree species basal areas. With respect to soils, organic matter content of the forest floor was calculated using Equation (1) where OM is organic matter content (Mg·ha −1 ), M is oven-dry sample mass (g), and A is core sampling area (50.265 cm 2 ).
Element content in the forest floor was calculated using Equation (2) where Q H is element content in the forest floor (kg·ha −1 ), and [x] is element concentration in the forest floor (mg·kg −1 ).
Element and OM contents in a given mineral soil layer were calculated using Equation (3) where Q m is element content in the mineral soil (kg·ha −1 ), [x] is element concentration (mg·kg −1 ) on an oven-dry basis, D b is bulk density (g·cm −3 ), and E e is effective thickness of the layer (cm), i.e., the corrected thickness of the soil layer without fragments (f).
Effective horizon thickness was calculated using Equation (4) where E e is effective thickness of the layer (cm), E is measured thickness of the layer (cm), and f is the coarsest fraction of the volumetric sample (>2 mm) (%/100).
We estimated the bulk density of the soil layers that had not been sampled quantitatively (n = 2114) using a quantitative relationship obtained from the measured D b and OM concentrations of the quantitative mineral soil samples (n = 1368) for each great soil texture class (sand, loam, clay) following the methodology of Federer, et al. [43] (Supplementary  Material Table S3 and Figure S1). The aqp R package v. 1.17 [44] was used to recalculate all the soil variables at fixed depths of 0-15, 15-30, and 30-60 cm (measured thickness). Element pools in individual layers were summed to represent the first 60 cm of mineral soil. All soil concentrations and pools are reported on an oven-dry mass basis.
Site quality index for each tree species, expressed as the height of the corresponding 50-year-old dominant trees, was determined using the age-height equations of Pothier and Savard [45] for Picea mariana, Abies balsamea, and Pinus banksiana, and those of Prégent, et al. [46] for Picea glauca (see Supplementary Material Figure S2). Average annual height growth (units: cm·yr −1 ) per species in each plot was also calculated from the plot surveys.

Statistical Analysis
Soil chemical properties and element pools along mineral soil depths and for the whole 0-60 cm mineral soil depth range were analyzed using a generalized nested linear mixed model. Harvesting treatments (SOH and WTH), soil provinces, soil depth and their third-order interactions were considered as fixed effects. Plots within soil provinces and location of individual soil samples within plots were considered as random effects. Properties, exchangeable and total element pools, forest floor thickness, biomass and mineralomass, as well as tree height growth variables (SQI for each tree species and annual tree height growth) were analyzed with a generalized nested linear mixed model that included harvesting treatments, soil provinces, and their second-order interactions as fixed effects. Tree height growth variables were analyzed separately by tree species. Plots and their associated drainage class within soil provinces and individual tree samples within plots were considered as random effects. Assumptions of homoscedasticity and normality of sample distributions were verified by residual plot analysis. The standardized residuals of these models were then plotted against all independent variables to detect possible heterogeneity in their variances. If present, variance heterogeneity was corrected by allowing for different variances per strata. We selected the models for which the variance function structures had the lowest Akaike Information Criterion (AIC) scores. Adjusted (predicted) means were computed with emmeans R package v. 1.4.2 [47], and statistical comparisons were made using specific a priori contrasts between harvesting treatments within soil provinces. The analyses were performed with the nlme R package v. 3.1-141 [48]. Pearson's product-moment correlations were also calculated between tree species, average SQI per plot, and all soil chemical properties and element pools.
Tree BAI growth over the period after harvesting was analyzed by generalized additive modeling (GAM) for each soil province, with harvesting treatment and tree species as fixed factors, and time after harvesting and interaction between time and the fixed factors as spline functions. The bam function of the mgcv R package, v. 1.8-31 [49] was used to perform the GAMs with the cubic regression (cr) penalized spline function and the starting basis dimension (knots) equal to 30. A first order (AR1) autoregressive error model was employed to reduce autocorrelations in the residuals. BAI values had to be log-transformed to homogenize the residues. We identified the periods during which the splines differed significantly between the two harvesting treatments (α = 0.05) using the plot_diff function of the itsadug R package v. 2.3 [50]. Simultaneous confidence intervals (95% CI) were computed for this comparison.

Harvested Biomass and Mineralomass
On average, WTH generated 32% more exported biomass than SOH, and 233% (P), 136% (K), 84% (Ca), and 112% (Mg) more mineralomass (Table 2). Differences in N exports amounted to 150 ± 10 kg·ha −1 (178%). There were also some similarities among soil provinces. Exported biomass and Ca were similar in both harvesting treatments in the granitic soil provinces (Laurentians (C) and Mistassini Highlands (E), p ≥ 0.062). In addition, in these two soil provinces similar amounts of biomass were harvested, but about half mineralomass (at the exception of Ca) were exported by SOH as compared to WTH. This is due to the larger basal area-and therefore greater bole biomass-in sites subjected to SOH compared to WTH (Supplementary Material Table S1). On average, the estimated quantities of biomass and mineralomass left on the ground (foliage and branches) by the SOH treatment (Supplementary Material Table S4) were very similar to the average calculated difference between WTH and SOH sites shown in Table 2.

Differences in the Mineral Soil after 30 Years
In the mineral soil, large differences were observed among soil provinces for all studied soil chemical properties (p ≤ 0.045; Table 3, Figure 2). In addition, several differences between the two harvesting treatments were observed across all soil provinces: (i) the soil exchangeable Mn pool was 32% and 42% larger in WTH sites (35.6 ± 2.7 kg·ha −1 for the whole soil 0-60 cm layer, and 12.2 ± 0.8 kg·ha −1 on average for a given soil layer) than in SOH sites (27.0 ± 2.1 kg·ha −1 for the whole soil 0-60 cm layer, and 8.6 ± 0.7 kg·ha −1 on av- ± 95% confidence intervals. Asterisks indicate significant differences between harvesting treatments within a given soil province: * p ≤ 0.05, ** p ≤ 0.01, *** p ≤ 0.001. Some other differences in the forest floor could also be observed between SOH and WTH sites, depending on the soil province: - In the Appalachians soil province (B), the total Mn pool was 26% smaller in WTH sites (p = 0.026; Table 4). - In the Laurentians soil province (C), BS was 16% lower in WTH sites than in SOH sites (p = 0.003; Table 3). - In the Abitibi Lowlands soil province (D), OM and organic C pools were about 25% greater in WTH sites (OM: 152 ± 8 Mg·ha −1 ; C: 72 ± 4 Mg·ha −1 ) than in SOH sites (OM: 121 ± 7 Mg·ha −1 ; C: 58 ± 3 Mg·ha −1 ; p = 0.004; Figure 2). In the forest floor, total P pool was nearly 20% larger (p = 0.021; Table 4), total Mn pool was 40% smaller (p = 0.026; Table 4), and exchangeable Mg and CEC pools were respectively 41% and 28% larger in WTH sites than in SOH sites (p ≤ 0.007; Table 3). Forest floor acidity status-as expressed by BS, exchangeable Al, exchangeable Fe, and exchangeable acidity pools-was 71% to 118% greater (19% smaller in the case of BS) in WTH sites than in SOH sites (p ≤ 0.014). - In the Mistassini Highlands soil province (E), no differences in OM and organic C pools in the forest floor were detected between WTH and SOH sites (p ≥ 0.119), but total N pool was 18% smaller in WTH sites (1.64 ± 0.12 Mg·ha −1 ) than in SOH sites (2.00 ± 0.14 Mg·ha −1 , p = 0.050; Figure 2). Total K pool in the forest floor also was 16% smaller (p = 0.018; Table 4), while total Fe pool was 35% larger in WTH sites than in SOH sites (p = 0.002). This was the only instance where total K and Fe pools in the forest floor differed between harvesting treatments.

Differences in the Mineral Soil after 30 Years
In the mineral soil, large differences were observed among soil provinces for all studied soil chemical properties (p ≤ 0.045; Table 3, Figure 2). In addition, several differences between the two harvesting treatments were observed across all soil provinces: (i) the soil exchangeable Mn pool was 32% and 42% larger in WTH sites (35.6 ± 2.7 kg·ha −1 for the whole soil 0-60 cm layer, and 12.2 ± 0.8 kg·ha −1 on average for a given soil layer) than in SOH sites (27.0 ± 2.1 kg·ha −1 for the whole soil 0-60 cm layer, and 8.6 ± 0.7 kg·ha −1 on average for a given soil layer, p ≤ 0.014), and (ii) the soil exchangeable Fe pool in the whole soil 0-60 cm layer was 15% smaller in WTH sites (94 ± 7 kg·ha −1 ) than in SOH sites (111 ± 6 kg·ha −1 , p = 0.012). Soil BS also differed between harvesting treatments across all soil provinces, but also according to soil depth, with about 33% higher values for WTH sites at depths of 15-30 cm and 30-60 cm (15-30 cm: 47% ± 4%; 30-60 cm: 60% ± 4%) as compared with SOH sites (15-30 cm: 35% ± 3%; 30-60 cm: 44% ± 3%, p ≤ 0.019). Soil BS in the 0-15 cm soil layer was similar for both harvesting treatments when averaged across all soil provinces (24.7% ± 3.1%, p = 0.866). Soil exchangeable Al and exchangeable acidity pools also differed between harvesting treatments and according to soil depth, but also among soil provinces (see below).
Some differences in soil variables could be observed between SOH and WTH sites depending on soil provinces. Major differences in OM, C, and N pools in soils were observed only in the two granitic soil provinces (Laurentians (C) and Mistassini Highlands (E); Figure 2). Over the whole 0-60 cm soil depth, the OM, C, and N pools in the Laurentians were approximately 28% smaller in WTH sites (in Mg·ha −1 , OM: 272 ± 26, C: 107 ± 12, N: 4.6 ± 0.5) than in SOH sites (in Mg·ha −1 , OM: 380 ± 19, C: 148 ± 10, N: 6.2 ± 0.4, p ≤ 0.004; Figure 3). These smaller pools in WTH sites were apparent at all three mineral soil depths, with the exception of C in the two deepest soil layers in the Laurentians (Figure 2). In the Mistassini Highlands, these pools were 33% (OM), 44% (C), and 46% (N) smaller in WTH sites (in Mg·ha −1 , OM: 170 ± 24, C: 56 ± 13, N: 2.6 ± 0.4) than in SOH sites (in Mg·ha −1 , OM: 255 ± 29, C: 100 ± 15, N: 4.8 ± 0.6, p ≤ 0.001; Figure 3) over the whole 0-60 cm soil layer. For this soil province, these pools were smaller in WTH sites only at the lower soil depths (15-30 and 30-60 cm, p ≤ 0.048; Figure 2), except for the N pool, which was smaller at all soil depths (p ≤ 0.015). In this regard, the C/N ratio in the mineral soils of Mistassini Highlands was higher in WTH sites (26.8 ± 1.9 on average for a given soil layer) than in SOH sites (19.6 ± 2.1 on average for a given soil layer; p = 0.009) throughout the soil profile. Significantly lower acidity status was also observed between WTH and SOH sites at soil depths of 15-30 cm and 30-60 cm in the Mistassini Highlands. Indeed, in this soil province, BS was 74% higher and exchangeable Al and exchangeable acidity pools were 48% to 61% smaller in WTH sites than in SOH sites at these soil depths (p ≤ 0.033, Table 3). The other granitic soil province, the Laurentians, also showed differences in acidity status; more precisely, pools of exchangeable Al and exchangeable acidity in the top 0-15 cm soil layer were approximately 40% smaller in WTH sites (Al: 252 ± 46 kg·ha −1 , exchangeable acidity: 31 ± 5 keq·ha −1 ) than in SOH sites (Al: 433 ± 36 kg·ha −1 , exchangeable acidity: 51 ± 4 keq·ha −1 , p ≤ 0.001). In the Appalachians, exchangeable K pools in the 15-30 cm and 30-60 cm soil layers were respectively 30% and 35% larger in WTH sites than in SOH sites (p ≤ 0.047), but K pool was similar for both treatments in the 0-15 cm soil layer (p = 0.086). In contrast, the exchangeable K pool in the Abitibi Lowlands was 41% smaller in the 0-15 cm soil layer of WTH sites (20 ± 8 kg·ha −1 ) as compared with SOH sites (34 ± 6 kg·ha −1 , p = 0.018), but not in the deeper soil layers (p ≥ 0.133). Finally, no practical difference in soil pH and statistical differences in exchangeable Ca and Mg pools could be detected between harvesting treatments in the soil provinces at any soil depth (p ≥ 0.292). 255 ± 29, C: 100 ± 15, N: 4.8 ± 0.6, p ≤ 0.001; Figure 3) over the whole 0-60 cm soil layer. For this soil province, these pools were smaller in WTH sites only at the lower soil depths (15-30 and 30-60 cm, p ≤ 0.048; Figure 2), except for the N pool, which was smaller at all soil depths (p ≤ 0.015). In this regard, the C/N ratio in the mineral soils of Mistassini Highlands was higher in WTH sites (26.8 ± 1.9 on average for a given soil layer) than in SOH sites (19.6 ± 2.1 on average for a given soil layer; p = 0.009) throughout the soil profile. Significantly lower acidity status was also observed between WTH and SOH sites at soil depths of 15-30 cm and 30-60 cm in the Mistassini Highlands. Indeed, in this soil province, BS was 74% higher and exchangeable Al and exchangeable acidity pools were 48% to 61% smaller in WTH sites than in SOH sites at these soil depths (p ≤ 0.033, Table 3). The other granitic soil province, the Laurentians, also showed differences in acidity status; more precisely, pools of exchangeable Al and exchangeable acidity in the top 0-15 cm soil layer were approximately 40% smaller in WTH sites (Al: 252 ± 46 kg·ha −1 , exchangeable acidity: 31 ± 5 keq·ha −1 ) than in SOH sites (Al: 433 ± 36 kg·ha −1 , exchangeable acidity: 51 ± 4 keq·ha −1 , p ≤ 0.001). In the Appalachians, exchangeable K pools in the 15-30 cm and 30-60 cm soil layers were respectively 30% and 35% larger in WTH sites than in SOH sites (p ≤ 0.047), but K pool was similar for both treatments in the 0-15 cm soil layer (p = 0.086). In contrast, the exchangeable K pool in the Abitibi Lowlands was 41% smaller in the 0-15 cm soil layer of WTH sites (20 ± 8 kg·ha −1 ) as compared with SOH sites (34 ± 6 kg·ha −1 , p = 0.018), but not in the deeper soil layers (p ≥ 0.133). Finally, no practical difference in soil pH and statistical differences in exchangeable Ca and Mg pools could be detected between harvesting treatments in the soil provinces at any soil depth (p ≥ 0.292).  Data presented are model-adjusted means ± 95% confidence intervals. Asterisks indicate significant differences between harvesting treatments within a given soil province: * p ≤ 0.05, ** p ≤ 0.01, *** p ≤ 0.001.

Site Productivity
The analyses revealed that for most tree species, SQI differed between WTH and SOH sites, but that the differences varied among soil provinces (Figure 4). For Abies balsamea, SQI was 9-20% lower in WTH sites (average: 9.5 ± 0.3 m) than in SOH sites (11.2 ± 0.1 m) in the Appalachians and the Laurentians provinces (p ≤ 0.046), but not in the two others. Picea mariana also had 8% to 22% lower SQIs in WTH sites (average: 10.3 ± 0.3 m) than in SOH sites (12.4 ± 0.1 m) in three of the four soil provinces (Laurentians, Abitibi Low-lands, and Mistassini Highlands; p ≤ 0.014). Picea glauca, which was sampled only in the Appalachians, had similar SQIs in both harvesting treatments (19.9 ± 0.4 m, p = 0.362). For Pinus banksiana, SQI was 14% lower in WTH sites (13.1 ± 0.3 m) than in SOH sites (15.2 ± 0.3 m) in the Mistassini Highlands (p ≤ 0.001), but SQI values were similar for both treatments in the Laurentians soil province (15.1 ± 0.4 m, p = 0.943). Soil C/N ratio was the variable most strongly related to SQI for the three main boreal tree species (Picea mariana, Abies balsamea, and Pinus banksiana) at various soil depths (r = −0.22 to −0.54, p ≤ 0.004; Figure 5). The second soil variable most strongly related to the SQI for Picea mariana and Abies balsamea was the pool of exchangeable acidity in the forest floor (r = −0.28 to −0.32, p ≤ 0.002). None of the studied soil variables were significantly associated with the SQI for Picea glauca at any depth (r ≤ 0.25, p ≥ 0.157).
sites, but that the differences varied among soil provinces (Figure 4). For Abies balsam SQI was 9-20% lower in WTH sites (average: 9.5 ± 0.3 m) than in SOH sites (11.2 ± 0.1 in the Appalachians and the Laurentians provinces (p ≤ 0.046), but not in the two othe Picea mariana also had 8% to 22% lower SQIs in WTH sites (average: 10.3 ± 0.3 m) than SOH sites (12.4 ± 0.1 m) in three of the four soil provinces (Laurentians, Abitibi Lowlan and Mistassini Highlands; p ≤ 0.014). Picea glauca, which was sampled only in the App lachians, had similar SQIs in both harvesting treatments (19.9 ± 0.4 m, p = 0.362). For Pin banksiana, SQI was 14% lower in WTH sites (13.1 ± 0.3 m) than in SOH sites (15.2 ± 0.3 in the Mistassini Highlands (p ≤ 0.001), but SQI values were similar for both treatments the Laurentians soil province (15.1 ± 0.4 m, p = 0.943). Soil C/N ratio was the variable m strongly related to SQI for the three main boreal tree species (Picea mariana, Abies balsam and Pinus banksiana) at various soil depths (r = −0.22 to −0.54, p ≤ 0.004; Figure 5). T second soil variable most strongly related to the SQI for Picea mariana and Abies balsam was the pool of exchangeable acidity in the forest floor (r = −0.28 to −0.32, p ≤ 0.002). No of the studied soil variables were significantly associated with the SQI for Picea glauca any depth (r ≤ 0.25, p ≥ 0.157).   Mean annual height growth showed greater variability than SQI (Figure 4). For Abies balsamea, despite large variations among and within soil provinces (p = 0.032), no difference was detected between the two harvesting treatments (p = 0.224). For Picea mariana, and only in the two granitic soil provinces (Laurentians and Mistassini Highlands), mean annual height growth was 18 to 37% less in WTH sites than in SOH sites (p ≤ 0.033). For Pinus banksiana, and only in the Mistassini Highlands, mean annual height growth was 25% less in WTH sites (22.5 ± 1.2 m) than in SOH sites (30.2 ± 1.7 m, p = 0.001).
The analysis of BAI also revealed some differences between WTH and SOH sites and among soil provinces ( Figure 6). Differences in BAI, all tree species confounded, were only observed in the two granitic soil provinces. Regardless of tree species, BAI growth was less in WTH sites than in SOH sites (p ≤ 0.013), at least after 10 years in the Laurentians, and over the 30-year time span of the study in the Mistassini Highlands. Picea mariana was the main tree species where these differences occurred (see Supplementary Material Figures S4 and S6, p ≤ 0.001). The BAI of Picea mariana and Picea glauca also differed between harvesting treatments in the Appalachians soil province. Tree growth was up to 50% greater for Picea mariana in WTH sites than in SOH sites, and up to 20% less for Picea glauca at some time during the 30 years after harvesting (see Supplementary Material Figure S3, p ≤ 0.003). In the Abitibi Lowlands soil province, BAI growth for Picea mariana seemed less in WTH sites than in SOH sites, but the difference was only marginally significant (see Figure 5 and Supplementary Material Figure S5, p = 0.075). For Abies balsamea, BAI did not differ or was only slightly greater in WTH sites during short periods in all soil provinces (Supplementary Material Figure S3 to S6, p ≥ 0.014). Mean annual height growth showed greater variability than SQI (Figure 4). For Abies balsamea, despite large variations among and within soil provinces (p = 0.032), no difference was detected between the two harvesting treatments (p = 0.224). For Picea mariana, and only in the two granitic soil provinces (Laurentians and Mistassini Highlands), mean annual height growth was 18 to 37% less in WTH sites than in SOH sites (p ≤ 0.033). For Pinus banksiana, and only in the Mistassini Highlands, mean annual height growth was 25% less in WTH sites (22.5 ± 1.2 m) than in SOH sites (30.2 ± 1.7 m, p = 0.001).
The analysis of BAI also revealed some differences between WTH and SOH sites and among soil provinces ( Figure 6). Differences in BAI, all tree species confounded, were only observed in the two granitic soil provinces. Regardless of tree species, BAI growth was less in WTH sites than in SOH sites (p ≤ 0.013), at least after 10 years in the Laurentians, and over the 30-year time span of the study in the Mistassini Highlands. Picea mariana was the main tree species where these differences occurred (see Supplementary Material Figures  S4 and S6, p ≤ 0.001). The BAI of Picea mariana and Picea glauca also differed between harvesting treatments in the Appalachians soil province. Tree growth was up to 50% greater for Picea mariana in WTH sites than in SOH sites, and up to 20% less for Picea glauca at some time during the 30 years after harvesting (see Supplementary Material Figure S3, p ≤ 0.003). In the Abitibi Lowlands soil province, BAI growth for Picea mariana seemed less in WTH sites than in SOH sites, but the difference was only marginally significant (see

Soil Organic Carbon
Our results showed smaller soil OM and organic C pools in WTH sites as compared to SOH sites, but only in the two granitic soil provinces characterized by a coarser soil texture, even though exported biomass was similar. The differences in OM after 30 years between WTH and SOH sites were apparent both in the top (0-15 cm) and deeper mineral soil layers (15-30 cm and 30-60 cm). Our results are consistent with those in the literature (Achat et al., 2015a). The importance of soil regions with distinct parent materials is highlighted in our analyses as elsewhere [51]. In coarse-textured soils, OM indeed appears to lack protection against decomposition in the mineral soil layers, partly because of their low silt and clay content. Physical protection of OM through microaggregation induced by clay minerals is a major mechanism of OM stabilization in soils [52]. In acidic forest soils such as those found in the two granitic soil provinces, soil OM is associated with minerals mainly by sorption to poorly crystalline Fe and Al oxides [53]. In some WTH sites of Quebec, upper soil layers with a coarse texture were still influenced by harvest treatment after 15 to 20 years, as shown by their lower organic C concentrations compared to SOH sites [54]. Fifteen years after WTH in a northern hardwood stand with a coarse-loamy texture at the Hubbard Brook Experimental Forest in New Hampshire, mineral soil C pools had decreased by 15%, and this decline mostly occurred at depths greater than 10 cm [55]. Ten years after repeated WTH in Finland on a sandy loam soil, total C and N pools in the combined organic and mineral soil layer were smaller than in the SOH treatment [56]. In general, SOH leads to an increase of soil OM pools at mid-soil depth as compared to non-harvested sites, while WTH does the opposite [57]. Therefore, the differences in soil OM and C pools between WTH and SOH sites that we observed in the deeper soil layers of the two granitic soil provinces can be attributed, at least in part, to an increase of these pools in SOH sites and a loss in WTH sites.
In the Appalachians and Abitibi Lowlands soil provinces, soil OM and organic C pools were not related to harvesting treatments. Thus, the fine-textured of soils may have helped to promote OM and C stabilization against disturbances in the mineral layers. However, in the Abitibi Lowlands, forest floor OM and C organic pools were found to be larger in WTH sites than in SOH sites; this seems counterintuitive, given the amount of biomass left on the ground by SOH operations. Our results suggest that WTH operations favored a greater accumulation of OM on the ground by sphagnum mosses than did SOH. This process, known to occur in the cold and humid conditions that are widespread in this flat soil province, favors the accumulation of organic matter over mineral soil after the tree canopy is cleared by harvesting operations or low-severity fires [58]. In this region, forest harvesting may lead to soil moisture saturation that can inhibit decomposition of OM in coniferous forests [59]. Our results also suggest that in all studied soil provinces, the higher C/N ratio (by 12% on average) in the forest floor of WTH sites has a negative impact on the quality of OM in this layer. As well, in the Mistassini Highlands, higher C/N ratios were found in all the mineral soil layers of WTH sites, as compared with SOH. An average increase of 2% of the C/N ratio in the forest floor of WTH sites was reported in a meta-analysis [12]. However, these researchers reported no major effect of WTH on the C/N ratio at various mineral soil depths.

Soil Nitrogen
Some general trends also emerged from our analysis regarding soil N pools 30 years after harvesting. In the two granitic soil provinces total N pool was smaller by 1.6 to 2.2 Mg·ha −1 (−25 to −45%) in WTH sites as compared with SOH sites. The additional amount of N exported by WTH was estimated at 150 ± 10 kg·ha −1 . Therefore, N losses can be attributed not only to the greater N export in the biomass by WTH, but also to leaching as suggested by the smaller mineral soil N pools with WTH in the two granitic soil provinces. The disturbance caused by forest harvesting leads to an increase in soil temperature, OM decomposition, and N mineralization. Sites harvested by WTH 12 years ago showed increased nitrification in the forest floor compared to unharvested sites located on an esker in the boreal forest in Ontario [60]. Nitrate is much more prone to leaching than ammonium in soils. The mineral soil N pool also decreased significantly by 14% over 15 years in a northern hardwood forest on coarse-loamy soil [55], likely due to increased nitrification. These observations all suggest that WTH significantly reduces soil N pool in coarse-textured soils and not on more fine-textured soils.

Soil Acidification
Signs of acidification were observed 30 years after WTH in some soil provinces. Lower BS and higher levels of exchangeable Al, exchangeable Fe, and exchangeable acidity were measured after WTH in the forest floor of the Abitibi Lowlands, consistent with the levels of mineralomass exports. These results are also consistent with those of Achat, Deleuze, Landmann, Pousse, Ranger, and Augusto [12], who reported lower BS (−8.4%) in the forest floor for WTH sites compared to SOH sites. However, no sign of acidification was found in the mineral soil layers of the Abitibi Lowlands. On the contrary, in the two granitic soil provinces, BS increased and other acidity indicators decreased in the mineral soil after WTH as compared with SOH. The increased soil acidification after SOH may have been caused by accumulated OM and C that migrated as dissolved organic carbon from the decomposing harvesting residues [61]. Base cation pools did not appear to decrease in soils after WTH as compared to SOH, with the exception of the smaller exchangeable K pool in the upper mineral soil layer after WTH in the Abitibi Lowlands. These results contradict another report of lower values of CEC (−9.9%) and BS (−17.4%) in the top mineral soil layers (<20 cm depth) in WTH sites than in SOH sites [12]. Mineral soil base cation pools may have been replenished through soil chemical weathering and conserved through decreased base cation uptake by the slower-growing secondary forest succession following WTH, as suggested by the tree growth analyses.

Specificity of Stand Response to Increased Biomass Harvesting
The lower site productivity observed after WTH in the two granitic soil provinces (the Laurentians and the Mistassini Highlands) can be attributed at least in part to the smaller N pools and higher C/N ratios 30 years after harvesting. Evidence from the Scandinavian countries showed that WTH could result in growth reductions attributable to N loss from logging residue exports [8,62]. After 31 years, the overall growth of Norway spruce stands in northern Sweden was still negatively affected by WTH compared to SOH, but this growth loss resulted only from a significant but temporary reduction in site productivity on WTH plots over a five-year period (at stand age of 8-12 years) [63]. This temporary reduction was also attributed to the N losses through biomass exportation, which hindered N nutrition for second-rotation trees during the first decade. Impaired N (or P) nutrition on WTH sites has been shown to reduce tree growth for at least 20 years in some stands [5].

Management and Policy Implications
In this study, we found that soil chemistry, site productivity, and tree growth responded differently to WTH as compared to SOH, and that the responses varied among soil provinces. Although soil base cation pools may have been reduced by WTH, the main soil properties influenced by WTH were OM, C, and N pools and C/N ratio. The current forest sustainability policy framework in Quebec relies only on soil acidification criteria to classify soil sensitivity to biomass harvesting. It targets forest ecological types that correspond to some site attributes found in this study (coarse-textured, thin, or organic soils) [64]. However, it should also take into account the findings obtained in the present study by using regional soil characteristics in a manner that is more specific than soil classification. For instance, WTH did not have the same impact on soil chemistry, site productivity and tree growth in the Appalachians as in the neighboring Laurentians, although both soil provinces have mainly Podzols. Furthermore, contrary to our findings of lower SQIs after 30 years in the two granitic soil provinces with coarser soils, [18] found no evidence after 20 years that more infertile sandy sites in neighboring Ontario were more sensitive to increased biomass exports. Such examples of diverging results call for a better knowledge of forest soils in terms of fertility and resilience to biomass harvesting.
This study shows that 30 years after treatment, the amount of OM varies significantly in the mineral soil down to a depth of 60 cm in the two granitic soil provinces with coarse-textured soils (the Laurentians and the Mistassini Highlands). On average, the difference between WTH and SOH sites regarding OM in these two soil provinces (−96 ± 28 Mg OM ha −1 ) was equivalent to 4.5 times the difference in estimated exported biomass (+21 ± 4 Mg OM ha −1 ). Therefore, the amount of C used as bioenergy from exported OM cannot compensate alone for the loss of C in the mineral soil. Furthermore, the decrease in site productivity in WTH sites should further worsen the situation after 30 years. It is clear that WTH in these two soil provinces is not a winning solution for decarbonizing energy uses. It remains to be seen whether the decreases in SQI and tree growth following WTH, associated with the increase in soil C/N ratio in these soil provinces, are transient or persistent, perhaps even throughout the entire life of the forest stands.
In the two other soil provinces with more fine-textured soils (the Appalachians and Abitibi Lowlands), biomass harvesting may result in an even or positive OM and C balance. Mineral soil OM and C did not appear associated with harvesting methods, and no decrease in mineral C/N ratio was observed. However, the cause of the observed decreases in SQI following WTH for these two soil provinces need further investigation.

Conclusions
Forest biomass is becoming more and more available for energy production in eastern North America. The development of this resource requires guidelines to remain sustainable. Our hypothesis that forest stands subjected to whole-tree harvesting in the past now displayed smaller nutrient pools and less desirable soil composition than those subjected to stem-only harvesting is supported. Our study indicates that after 30 years, biomass exportation through whole-tree harvesting had negative impacts mainly on soil OM, C, and N pools and soil C/N ratio from the forest floor down to 60 cm in the mineral soil. Our hypothesis that site productivity was also negatively affected by whole-tree harvesting is supported by our observations. However, the magnitude of these changes strongly varied among soil provinces. Greater effects were observed in the two granitic soil provinces (Laurentians and Mistassini Highlands), which have a coarser soil texture: for these, SQI, height growth, and BAI were all lower in WTH sites than in SOH sites. The Appalachians was the soil province where soils of WTH and SOH sites differed the least, as shown by similar tree growth, particularly for Picea glauca and Picea mariana. In general, SQI was negatively related to the soil C/N ratio. Sensitivity to biomass exportation by WTH was greatest in the two granitic soil provinces, followed by the Abitibi Lowlands and the Appalachians. The specificity of soils must be recognized when developing forest sustainability criteria and policies on forest biomass harvesting and proper silviculture in order to protect soil fertility and to maintain soil carbon reserves.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/f12050583/s1, Figure S1. Relationship between soil bulk density (Db) and organic matter (OM) concentration according to great soil texture groups. Lines show model predicted values and 95% confidence interval predictions, Figure S2. Expected tree height as a function of tree age (at DBH height) for various site quality index (SQI) values. Data is derived from the models of Pothier and Savard [45] for Picea mariana, Abies balsamea, and Pinus banksiana, and of Prégent [46] for Picea glauca, Figure S3. Left panels: Predicted basal area increment (BAI) of individual regenerating tree species (ABA = Abies balsamea, PIM = Picea mariana, PIG = Picea glauca) over time after harvesting treatment (stem-only harvesting [SOH] and whole-tree harvesting [WTH]) in the soil province B (Appalachians). Data presented are model averages (lines) and simultaneous 95th centile confidence intervals (bands) (95% CI). Numbers of samples are displayed for each treatment. Right panels: Estimated difference in BAI (log values, simultaneous 95% CI) between WTH and SOH (WTH minus SOH) smoothed curves. Periods during which treatments differ significantly are highlighted with a red line on the X axis,