Harvest Intensity Effects on Carbon Stocks and Biodiversity Are Dependent on Regional Climate in Douglas-Fir Forests of British Columbia

Simard, Suzanne W.; Roach, W. Jean; Defrenne, Camille E.; Pickles, Brian J.; Snyder, Eva N.; Robinson, Alyssa; Lavkulich, Les M.

doi:10.3389/ffgc.2020.00088

ORIGINAL RESEARCH article

Front. For. Glob. Change, 24 July 2020

Sec. Forest Management

Volume 3 - 2020 | https://doi.org/10.3389/ffgc.2020.00088

Harvest Intensity Effects on Carbon Stocks and Biodiversity Are Dependent on Regional Climate in Douglas-Fir Forests of British Columbia

$\r\nSuzanne W. Simard*$ Suzanne W. Simard^1*

W. Jean Roach¹

Camille E. Defrenne^2†

Brian J. Pickles^3†

Eva N. Snyder¹

Alyssa Robinson¹

Les M. Lavkulich⁴

¹Department of Forest and Conservation Sciences, Faculty of Forestry, The University of British Columbia, Vancouver, BC, Canada
²Environmental Sciences Division and Climate Change Science Institute, Oak Ridge National Laboratory, Oak Ridge, TN, United States
³School of Biological Sciences, University of Reading, Reading, United Kingdom
⁴Faculty of Land and Food Systems, The University of British Columbia, Vancouver, BC, Canada

Temperate forests provide crucial ecosystems services as living sinks for atmospheric carbon (C) and repositories of biodiversity. Applying harvesting at intensities that minimize losses offers one means for mitigating global change. However, little is known of overstory retention levels that best conserve ecosystem services in different regional climates, and likewise as climate changes. To quantify the effect of harvest intensity on C stocks and biodiversity, we compared five harvesting intensities (clearcutting, seedtree retention, 30% patch retention, 60% patch retention, and uncut controls) across a climatic aridity gradient that ranged from humid to semi-arid in the Douglas-fir (Pseudotsuga menziesii) forests of British Columbia. We found that increased harvesting intensity reduced total ecosystem, aboveground, and live tree C stocks 1 year post-harvest, and the magnitude of these losses were negatively correlated with climatic aridity. In humid forests, total ecosystem C ranged from 50% loss following clearcut harvest, to 30% loss following large patch retention harvest. In arid forests this range was 60 to 8% loss, respectively. Where lower retention harvests are sought, the small patch retention treatment protected both C stocks and biodiversity in the arid forests, whereas the seedtree method performed as well or better in the humid forests. Belowground C stocks declined by an average of 29% after harvesting, with almost all of the loss from the forest floor and none from the mineral soil. Of the secondary pools, standing and coarse deadwood declined in all harvesting treatments regardless of cutting intensity or aridity, while C stocks in fine fuels and stumps increased. The understory plant C pool declined across all harvesting intensities in the humid forests, but increased in arid forests. Shannon’s diversity and richness of tree and bryoid species declined with harvesting intensity, where tree species losses were greatest in the humid forests and bryoid losses greatest in arid forests. Shrub and herb species were unaffected. This study showed that the highest retention level was best at reducing losses in C stocks and biodiversity, and clearcutting the poorest, and while partial retention of canopy trees can reduce losses in these ecosystem services, outcomes will vary with climatic aridity.

Introduction

Forests store approximately 80% of total global aboveground carbon (C) and 70% of terrestrial soil C due to photosynthetic capture by trees and plants (Sedjo, 1993; Dixon et al., 1994; Goodale et al., 2002; IPCC, 2018), and have sequestered the equivalent of 20% of global fossil fuel emissions over the past three decades (Le Quéré et al., 2018). In addition, they are home to more than three-quarters of the world’s vertebrate and invertebrate species (Lindenmayer and Franklin, 2002; Food and Agricultural Organization of the United Nations [FAO], 2018), and this biodiversity is correlated with productivity and C storage (Thompson et al., 2012; Liang et al., 2016; Lecina-Diaz et al., 2018). Because of the high C sink strength and habitat quality of forests (Goodale et al., 2002; Pugh et al., 2019), improved forest stewardship, including reduced deforestation, enhanced reforestation, and afforestation, is considered the most effective approach for mitigating global change (Smith et al., 2016; Field and Mach, 2017). Griscom et al. (2017) argues that improved stewardship has the potential to provide 37% of the mitigation needed for stabilizing global warming below 2°C by 2039. Of the stewardship approaches, reduced deforestation is the most widely effective and acceptable strategy (Keith et al., 2014; Field and Mach, 2017; Buotte et al., 2019a). However, primary forests continue to be clearcut and converted to tree plantations, agricultural crops or other land uses at a rapid rate, accounting for up to 40% of historic global anthropogenic CO₂ emissions worldwide (Houghton, 2010) and detrimental effects on biodiversity (Betts et al., 2017).

Clearcutting removes all trees and causes immediate release of aboveground C (Liu et al., 2011; Goetz et al., 2012; Wang et al., 2012) as well as up to 30% of the forest floor pool (Nave et al., 2010; James and Harrison, 2016). The magnitude of C losses from clearcutting increases with forest productivity, which is negatively correlated with climatic aridity (Buotte et al., 2019a). Partial cutting, where a portion of the trees are retained, however, has the potential to reduce aboveground C loss, increase understory C storage, reduce or eliminate losses from the forest floor, and reduce forest recovery time (Zhou et al., 2013). The retention of some trees also increases structural and species diversity (Franklin et al., 1987; Gustaffson and Perhans, 2010). Partial cutting involves removing individual trees or groups of trees from the forest in different intensities and patterns for various objectives; these can include provision of a seed source or shelter for regeneration, enhancement of growth of the remaining trees, management of species composition and stand structure, improvement in wildlife habitat, reductions in disturbance to the forest floor and erosion, or reductions in vulnerability to disturbances such as fire or pests (Weetman and Vyse, 1990; Kneeshaw et al., 2002). The appropriate harvesting intensity for conservation of carbon stocks or biodiversity in different forest ecosystems is poorly understood, however, and is expected to vary with local site and climatic conditions, and to shift along with climate change. While considerable attention is being paid to global tree restoration as a global change mitigation strategy (Bastin et al., 2019; Friedlingstein et al., 2019; Lewis et al., 2019; Luedeling et al., 2019; Veldman et al., 2019), ongoing widespread clearcutting is of serious concern, and there is an urgent need to find alternatives to help avoid dangerous positive feedbacks to global change and species extinctions.

Given that harvesting at some level to meet human needs in the future is inevitable, a critical pair of questions remains: Which harvesting practices will minimize stand-level losses of C and biodiversity, and will these responses change as climate changes? Meta-analyses indicate that on-site C stocks in aboveground trees decrease linearly with cutting intensity (Johnson and Curtis, 2001; Zhou et al., 2013; Hume et al., 2017), but variable retention harvesting can leave residual trees in a wide range of patterns that could affect losses and rate of recovery. Soil C also declines with harvesting, with more lost from forest floors than mineral soils, with magnitude of loss dependent upon forest type, climatic zone and soil order (Johnson and Curtis, 2001; Nave et al., 2010; Zhou et al., 2013; James and Harrison, 2016; Hume et al., 2017). Little is known, however, of how partial canopy retention may mitigate these losses. A synthesis of 100 studies concluded that selection harvesting positively affected understory species richness, but richness was not affected by even-aged management treatments in the first 50 years following logging (Duguid and Ashton, 2013). These meta-analyses have the advantage that datasets are large so that trends masked by variability in smaller datasets can be identified. However, they have less value for understanding responses regionally or locally because the data were derived from large geographic areas (e.g., worldwide temperate forests).

Researchers have compared the impacts of various harvesting treatments on C stocks within a climatic region using either empirical (Olsson et al., 1996; Strong, 1997; Powers et al., 2012; Strukelj et al., 2015) or modeling approaches (Nunery and Keeton, 2010; Scheller et al., 2011). Modelers have also projected the impacts of a single harvesting treatment on C stocks across a climatic gradient (Wang et al., 2014). However, few empirical studies exist where the interactive effects of harvesting intensity and regional climate have been examined, especially where more than one C pool is measured. Likewise, many studies have explored the effects of climate change on biodiversity and several major reviews have been published (Walther et al., 2002; Root et al., 2003; Parmesan, 2006; Bellard et al., 2012). Estimates have also been made of harvesting intensity impacts on biodiversity in a single climatic region (Carreño-Rocabado et al., 2012; Ding et al., 2017, 2019; Pyles et al., 2018). However, we are not aware of any field studies that have examined how climate and harvesting intensity interact to affect biodiversity, although there have been investigations of interactions between climate change and land use (Monge-González et al., 2019; Peters et al., 2019).

Here, we conducted an extensive field experiment examining how harvesting intensity interacts with regional climatic aridity to affect C stock and biodiversity components in the Douglas-fir (Pseudotsuga menziesii Mirb. (Franco)) forests of British Columbia. Douglas-fir is commercially important, has the widest natural latitudinal distribution of any tree species in North America (19–55 ^oN), and its range is predicted to expand with climate change (Wang et al., 2016). It is also extensively planted on other continents as a productive exotic timber species (Hermann and Lavender, 1999). We compared C stock and biodiversity components following five levels of harvest intensity applied in six climatic regions, from humid to semi-arid, of the Douglas-fir forests in British Columbia. We tested the following hypotheses: (1) Total ecosystem and aboveground C stocks will be reduced by harvesting proportional with harvesting intensity. (2) Total C losses from harvesting will decline with increasing climate aridity, but losses will be proportional to pre-treatment C stocks regardless of climatic region. (3) C stocks of the forest floor but not the mineral soil will decline with harvesting intensity, and losses will occur regardless of climatic region. (4) Diversity and richness of understory plants and bryoids will decline with harvesting intensity.

Materials and Methods

Study Sites Across the Climatic Gradient

The study took place in mixed Douglas-fir-conifer forests of British Columbia (BC), Canada, with five interior and one coastal site covering a 900 km-long climate gradient (Figure 1 and Table 1). The interior locations are within the current range of interior Douglas-fir (P. menziesii var. glauca) in BC (mean annual temperature (MAT) 2.3–7.7 ^oC, mean annual precipitation (MAP) 532–1059 mm); Interior Douglas-fir (IDF), Interior Cedar-Hemlock (ICH), and Sub-boreal Spruce (SBS) biogeoclimatic (BEC) zones), and are located from south of Cranbrook (49.21 ^oN, 115.37 ^o W) to Fort St James (54.65 ^oN, 124.43 ^oW). The coastal location (MAT 8.0 ^oC, MAP 2701 mm; Coastal Western Hemlock (CWH) zone) (49.32°N, 22.54°W), situated near Maple Ridge and about 60 km east of Vancouver, BC is within the range of the coastal variety of Douglas-fir (P. menziesii var. menziesii).

FIGURE 1

Figure 1. Map of the six research locations and the current distribution of interior Douglas-fir in British Columbia.

TABLE 1

Table 1. Characteristics of the six study locations In British Columbia, Canada (n = 3 for all locations, except Narrows Creek where n = 1, and Alex Fraser where n = 2).

At three out of six forests, Douglas-fir made up more than 50% of basal area (BA), with exceptions being the wet interior (33% BA) and coastal locations (16% BA). Each forest featured other common tree species that varied among zones: western larch (Larix occidentalis Nutt.) in the IDF zone; hybrid white spruce (Picea glauca x engelmannii) in the SBS; western larch, western red-cedar (Thuja plicata (Donn ex D Don), western hemlock, (Tsuga heterophylla (Raf.) Sarg.), and grand fir (Abies grandis (Dougl. ex D. Don) Lindl.) in the ICH; and western redcedar and western hemlock in the CWH. Several other tree species were present in minor amounts, with up to eight tree species at one location. Prior to logging, the understory was dominated by grasses at the driest site (IDF zone), shrubs and herbs at sites with moderate moisture, and was sparse at the two wettest locations.

All sites have mesic soil moisture regime (relative to their respective BEC zone), warm aspects, gentle slope gradients (< 30 percent) and mid-slope positions. Elevations vary from 540 to 1080 m. The soil texture at the interior locations tended to be coarsest in the ICH zone, followed by the IDF, and finest in the SBS. The dominant soil orders were Podzols in the ICH and CWH, and Luvisols in the IDF and SBS, and are representative of these zones (Soil Classification Working Group, 1998). Climate data for each site was obtained from ClimateNA, based on their latitude, longitude and elevation, using the 1981–2010 climate normal dataset (Wang et al., 2016).

Harvesting Treatments

The interaction between regional climate and harvesting intensity treatment was tested in a randomized block design. The regional climate factor was represented by six forest locations, each in a unique regional climate. Four of the locations included three 20-ha replicate sites while Alex Fraser contained two and Narrows contained one (total 15 sites), and each site was divided into five 4-ha treatment blocks (total 75 blocks) with harvesting intensity treatments randomly assigned to the blocks. The harvesting intensity treatments were applied in 2017/18 prior to spring green-up. They included: (i) clearcut (no tree retention); (ii) single tree retention (retention of 25 large Douglas-fir stems ha^–1 at 25 m spacing); (iii) small patch retention (retention of 30% of the stand area in 30-m wide unconnected patches, with all trees cut in the remaining 70% of the stand); (iv) large patch retention (retention of 60% of the stand area with all trees cut in the remaining 40% of the stand); and (v) uncut control (full retention). Harvesting was carried out using feller-bunchers, but trees too large for the machine were hand-felled. The high retention blocks were thinned from below by reaching into the uncut patches with the feller-buncher and removing the smaller stems. These harvesting treatments represented a graduated series of tree retention, from the lowest in the clearcut, with increasing amounts in the single tree, small patch and large patch retention treatments, roughly corresponding with regeneration systems described as seed-tree, group selection and single tree retention with thinning from below. Trees were limbed and topped but not debarked before the boles were skidded to roadside landings, leaving needles, twigs and branches distributed relatively evenly throughout the forest after harvesting.

Measurement and Sampling Methods

Field data was collected in one permanent circular 0.04 ha (r = 11.3 m) National Forest Inventory (NFI) plot per treatment block in July/August prior to logging, and these plots were re-assessed the summer following harvesting (i.e., 1 year-post logging). The NFI plot was located as close as possible to the center of the treatment plot and positioned to best represent the treatment. In the case of the 30% and 60% retention treatments, the NFI plots were positioned to represent the same proportion that had been felled and retained. In the clearcut, seedtree and control treatments, the treatments were relatively uniform and required little repositioning of the NFI plot from the center location.

Data collection followed the NFI ground sampling guidelines (Canadian Forest Inventory Committee [CFIC] and Canadian Forest Service [CFS], 2008). All live and dead standing trees ≥9 cm DBH and > 1.3 m height within the plot boundary were tagged at the base and species, diameter at breast height (1.3 m), height, crown length, and condition were recorded. At the pre-logging assessment, increment cores were collected to determine the age of at least one undamaged codominant or dominant live tree of each major tree species. Nested within each 0.04 ha plot was one circular 50 m² subplot (r = 3.99 m) in which assessments were made of density, diameter and height of live and dead trees and shrubs < 9 cm DBH and > 1.3 m tall, and height, diameter and decay class of stumps > 4 cm in diameter and ≤1.3 tall. Visual assessments were made of percent cover of all live tree and shrub species occurring within a 314 m² subplot (r = 10.0 m) and of all herbaceous, ground-level shrub and bryoid species occurring within a 100 m² (r = 5.64 m) subplot. Coarse woody debris (CWD) (> 7.5 cm diameter) was measured along two perpendicular 30-m line transects and small woody debris (SWD) (1.1–7.5 cm diameter) pieces were counted by size class along 10 m of each transect, and included all CWD and SWD remaining after logging. Ground substrate type and depth of organic and decayed wood were recorded every 2 m along each CWD transect (total 30 stations per NFI plot).

Samples were collected for laboratory analysis at one circular 1 m² (r = 0.56 m) microplot at each end of the CWD transects (four microplots per NFI plot). Post-logging microplots were established 1.5 m clockwise on the 0.04 ha boundary from the pre-logging microplot locations. All bryoids, herbs, trees/shrubs ≤1.3 m in height, and fine woody debris (FWD) (< 1 cm diameter), including needles, twigs and branches remaining after logging, occurring within the microplots were collected separately. A 20 × 20 cm area of forest floor extending from the ground surface to the underlying substrate was measured for thickness and collected. Mineral soil, including rocks ≤7.5 cm diameter and roots, was collected from all four microplots in 10-cm diameter holes at 0–15 cm depth. The excavated holes were lined with plastic and the volume of water that filled the hole was measured with a graduated cylinder for use in calculating bulk density (Walter et al., 2016; Al-Shammary et al., 2018). Cobbles (rocks > 7.5 cm across) were weighed in the field and discarded.

Site data was collected to describe the ecological properties of each NFI plot. A soil pit was dug to ≥60 cm depth, mineral soil horizons were delineated, and their depth, texture and coarse fragment content were estimated. Mineral soil and forest floor were classified to Order (Green et al., 1993; Soil Classification Working Group, 1998). The biogeoclimatic variant of each NFI plot was determined from field maps (British Columbia Ministry of Forests Lands Natural Resource Operations and Rural Development [BC MoFLNRORD], 2019a) and soil moisture regime, nutrient status and site series were estimated using vegetation, soil and site features (Lloyd et al., 1990; Braumandl and Curran, 1992; Delong et al., 1993; Green and Klinka, 1994; Steen and Coupe, 1997).

Laboratory Analysis

The plant and FWD samples were oven dried and weighed then discarded. The mineral soil samples were oven dried, sieved and separated into roots, particles < 2 mm, and gravel (2–7.5 cm) components and each fraction was weighed. The forest floor samples were oven dried and sieved into < 8 mm and ≥8 mm fractions which were weighed. Subsamples of each < 2 mm mineral and < 8 mm forest floor fractions were sent to the Ministry of Environment laboratory in Victoria, BC for determination of C and nitrogen (N) concentrations.

Data Analysis

Descriptive statistics for C storage in each major pool and the biodiversity components were calculated for each harvesting intensity treatment in each of the climatic regions by calculating an average of the replicates. Carbon storage (Mg ha^–1) in aboveground dead and live trees > 1.3 m tall, stumps (#x2264;1.3 m tall), downed coarse, small, and fine woody debris, understory plants, forest floor, and mineral soil were calculated according to the Canadian Forest Service (CFS) protocol (Canadian Forest Inventory Committee [CFIC] and Canadian Forest Service [CFS], 2008). Root biomass of live trees was calculated as aboveground tree biomass multiplied by 0.29 or 0.20 for coniferous forests with an aboveground biomass of 50–150 Mg ha^–1 or > 150 Mg ha^–1, respectively, based on the work of Mokany et al. (2006) as suggested in Table 4 of the Guidelines for national greenhouse gas inventories (IPCC, 2006). Root biomass of each dead tree was calculated as the root biomass when the tree was alive minus the fine root proportion of total root biomass (P), where P = 0.072 + (0.354)^∗(e^–0.060RB) where RB = total root biomass (Li et al., 2003). Based on visual characteristics such as lack of needles and sloughing bark, all dead trees were assumed to have died several years prior to their measurement, in which case all fine roots would be dead, based on a fine root turnover rate of 64.1% per year (Kurz et al., 2009). Root biomass of trees cut during harvesting was calculated as root biomass when the tree was alive minus the dead fine roots (0.641^∗P). Biodiversity of trees and plants were estimated as: (i) richness (the number of species occurring in a 0.0314 ha plot for trees and shrubs, and in a 0.01 ha plot for herbs and bryoids), and (ii) Shannon’s diversity index (Shannon and Weaver, 1949), calculated as:

H^{'} = 1 - Σ [p {}_{i}\ln (p {}_{i})]

where H’ is Shannon’s index and p_i is the proportion of total density of species i.

Statistical analyses were conducted in R version 3.6.1 (R Core Team, 2018) and results were considered statistically significant at P < 0.05. We investigated how Douglas-fir-dominated forest carbon stocks and biodiversity responded to harvest intensity (clearcut, single tree, small patch retention, high retention, and control (full retention)) and climatic aridity (i.e., annual heat moisture index = (mean annual temperature + 10)/(mean annual precipitation/1000)) by fitting linear mixed-effects models (LMMs). Before analyses, most of the response variables were transformed (log₁₀ or square root; see Tables 1, 2) to meet the assumptions of the models. Aridity and harvesting intensity were added as interacting fixed effects while location and replication were added as nested random effects (the plot level was not added as a nested random effect because the number of observations at this level was insufficient). We used the ‘lmer’ function for LMMs from the package lme4 (Bates et al., 2018). Models were fitted by restricted maximum likelihood. We used the function ‘step’ from the lmerTest package to eliminate non-significant effects of models based on the Akaike information criterion (AIC; Venables and Ripley, 2002). We checked model adequacy with the ‘plot_model’ function from the sjPlot package (Lüdecke, 2018) and model fit by estimating marginal R² (R²m) values (variance explained by the fixed effects) following Nakagawa and Schielzeth (2013). Model significance was tested with a likelihood ratio test and significance of fixed effects with a type II Wald χ2 test. Standardized beta coefficients and their 95% confidence intervals were extracted with the ‘plot_model’ function from the sjstats package.

TABLE 2

Table 2. Effect of harvesting intensity and climatic aridity (1981–2000) on total ecosystem and belowground carbon storage of mixed Douglas-fir forests (n = 3 for all locations, except Narrows Creek where n = 1, and Alex Fraser where n = 2).

Results

C Stocks

Before harvest, total ecosystem C stocks ranged from 478 ±47 (SE) Mg ha^–1 at the most humid extreme of our climate gradient to 229 ±27 Mg ha^–1 in the most arid forests (Figures 2, 3 and Supplementary Table S1). On average, 62% occurred in aboveground pools and 38% belowground (mineral soil sampled to 15 cm depth), with little variation in above versus belowground distribution along the climate gradient (Figure 2 and Supplementary Table S4). The distribution of C among individual above and belowground pools was also similar across the climate gradient (Supplementary Table S4). The average proportion in each pool in the humid forests was: live trees (50.0%), dead trees (8.6%), coarse woody debris (5.5%), small and fine woody debris (1.1%), understory plants (0.2%), tree roots (12.3%), forest floor (16.3%) and mineral soil to 15 cm (7.3%). In arid forests, the distribution was as follows: live trees (46.5%), dead trees (5.8%), coarse woody debris (6.8%), small and fine woody debris (1.2%), understory plants (0.2%), tree roots (12.8%), forest floor (14.6%) and mineral soil to 15 cm (12.1%) (Figure 2).

FIGURE 2

Figure 2. Distribution of carbon stocks among above and belowground pools of mixed Douglas-fir forests (n = 3 for all locations, except Narrows Creek where n = 1 and Alex Fraser where n = 2). (A) presents mean carbon pools for each location across the climatic aridity gradient (averaged across harvesting intensities). (B) presents mean carbon pools for each harvesting treatment across the harvesting intensity gradient (averaged across locations). The aridity index, AHM, was calculated as [mean annual temperature + 10]/[mean annual precipitation/1000].

FIGURE 3

Figure 3. Effects of harvesting intensity and climatic aridity on the ecosystem, aboveground and belowground carbon pools of mixed Douglas-fir forests (n = 3 for all locations, except Narrows Creek where n = 1 and Alex Fraser where n = 2). (A,C,E) present the relationship between harvesting intensity and carbon stocks. (B,D,F) present standardized beta coefficients of the linear mixed-effect models. Beta coefficients illustrate the effect of factors on the total, aboveground and belowground carbon pools in terms of their standardized effect sizes. Circles indicate average estimates and lines indicate 95% confidence intervals. The aridity index was calculated as [mean annual temperature + 10]/[mean annual precipitation/1000]. Pre-harvest carbon stocks were added as covariates in the linear mixed-effect models. For harvesting intensity, each coefficient represents the difference between the uncut control and a given treatment. ***P-value < 0.001; **P-value < 0.05.

Harvesting intensity interacted with climatic aridity in reducing total ecosystem C (R²m = 0.70), total aboveground C (live and dead trees + coarse woody debris + stumps + understory plants) (R²m = 0.70), and live tree C (R²m = 0.77; Figures 2–4 and Table 2). After harvesting, these C pools were still greater in the more humid than arid regional climates, particularly for the clearcut, seedtree and small patch retention treatments (except for the live tree C pool; Figures 3, 4). In the most humid forests, total ecosystem C stocks were reduced to 171 Mg ha^–1 with clearcutting and 359 Mg ha^–1 with small patch retention, for example, while in the most arid forests the stocks were reduced to 95 Mg ha^–1 and 131 Mg ha^–1 in clearcut and small patch retention treatments, respectively (Table 2 and Supplementary Table S1). The magnitude of loss paralleled that in remaining total C, aboveground C, and live tree C stocks, where the calculated magnitude of the losses increased with harvesting intensity and decreased with climatic aridity, with significant interactions between the two factors (R²m = 0.58–0.75, Table 2). In support of this, total C, aboveground C, and live tree C stocks were significantly correlated with pre-harvest values (Table 2). As expected, the ratio of total C and aboveground C from before compared to after harvesting increased with harvesting intensity (R²m = 0.52 and 0.57, Table 2 and Supplementary Tables S1, S2).

FIGURE 4

Figure 4. Effect of harvesting intensity and climatic aridity on live tree carbon (A,B) and understory plant carbon pools (C,D) of mixed Douglas-fir forests (n = 3 for all locations, except Narrows Creek where n = 1 and Alex Fraser where n = 2). (B,D) present standardized beta coefficients of linear mixed-effects models. Beta coefficients illustrate the effect of factors on the live tree carbon pool (B) and on the understory plant carbon pool (D) in terms of their standardized effect sizes. Circles indicate average estimates and lines indicate 95% confidence intervals. The aridity index was calculated as [mean annual temperature + 10]/[mean annual precipitation/1000]. For harvesting intensity, each coefficient represents the difference between the uncut control and a given treatment. ***P-value < 0.001; **P-value < 0.01; *P-value < 0.05.

Harvesting shifted the distribution of C amongst pools with greater shifts as harvesting intensity increased (Figures 2–5 and Supplementary Table S4). In humid forests, the average proportions following clearcutting versus large patch retention were: live trees (0.0% vs. 45.2%), dead trees (0.0% vs. 0.0%), coarse woody debris (23.7% vs. 11.6%), small and fine woody debris (8.8% vs. 1.0%), stumps (2.5% vs. 2.2%), understory plants (0.1% vs. 0.1%), tree roots (23.9% vs. 19.1%), forest floor (11.4% vs. 4.2%) and mineral soil to 15 cm (29.7% to 15.2%). In arid forest, these averages were: live trees (0.4% vs. 60.8%), dead trees (0.0% vs. 2.2%), coarse woody debris (18.2% vs. 2.3%), small and fine woody debris (13.7% vs. 3.2%), stumps (1.9% vs. 0.3%), understory plants (0.2% vs. 0.8%), tree roots (21.4% vs. 14.4%), forest floor (21.0% vs. 8.4%) and mineral soil to 15 cm (23.2 vs. 7.7%). Of the secondary aboveground pools, C stocks in total woody debris (coarse + small + fine) stayed relatively stable with harvesting intensity. We note that in the most humid forests (Malcolm Knapp), the most downed wood remained in the clearcut followed by the seedtree then the two patch retention treatments (Figure 5 and Supplementary Table S2). There was no effect of aridity on total woody debris nor was there a harvesting by aridity interaction, but pre-harvest values were correlated with post-harvest values (R²m = 0.70; Figure 5 and Table 2).

FIGURE 5

Figure 5. Effect of harvesting intensity and climatic aridity on carbon pools in total woody debris (A,B), forest floor (C,D), and tree roots (E,F) of mixed Douglas-fir forests (n = 3 for all locations, except Narrows Creek where n = 1 and Alex Fraser where n = 2). (B,D,F) present standardized beta coefficients of linear mixed-effect models. Beta coefficients illustrate the effect of factors on total woody debris (B), the forest floor (D), and tree roots (F) in terms of their standardized effect sizes. Circles indicate average estimates and lines indicate 95% confidence intervals. The aridity index was calculated as [mean annual temperature + 10]/mean annual precipitation/1000]. For harvesting intensity, each coefficient represents the difference between the uncut control and a given treatment. ***P-value < 0.001; **P-value < 0.01; *P-value < 0.05.

Harvesting intensity and aridity significantly interacted to affect total understory plant C storage (R²m = 0.31; Figure 4 and Table 2), with low harvesting intensity treatments protecting understory C stocks in more arid forests (Supplementary Table S2).

Total belowground C stocks slightly declined with harvesting intensity, although not significantly, and were unaffected by aridity (R²m = 0.53) (Figures 2, 3 and Table 3), with the greatest magnitude losses across the entire climate gradient in the clearcut followed by the seedtree, then the small and large patch retention treatments (R²m = 0.23). C stocks in the forest floor decreased with harvesting intensity, and were correlated with pre-harvest values, but they were unaffected by climatic aridity, with no significant interaction between the harvest intensity and aridity factors, and a model marginal R² of 0.30 (Figure 5 and Table 3). Any amount of harvesting reduced C stocks in the forest floor by 61% on average, with no significant difference between the three partial retention treatments (Figure 2, Table 3, and Supplementary Table S3). Mineral soil C was unaffected by harvesting intensity or aridity, but it was correlated with pre-harvest total mineral soil C stocks (R²m = 0.25; Table 3 and Supplementary Table S3).

TABLE 3

Table 3. Effect of harvesting intensity and climatic aridity (1981–2010) on belowground carbon storage of Douglas-fir dominated forests (n = 3 for all locations, except Narrows Creek where n = 1, and Alex Fraser where n = 2).

Biodiversity

Shannon’s diversity index of the whole plant community declined with harvest intensity, but it was unaffected by climatic aridity (model R²m = 0.76; Figure 6 and Table 4). Harvesting and aridity interacted to affect tree species diversity, where losses increased with harvesting intensity and decreased with aridity (R²m = 0.42; Table 4). Shannon diversity of bryoids declined with harvesting intensity, but that of shrubs or herb groups were unaffected by either harvesting or aridity (Table 4 and Supplementary Table S5). Pre-harvest diversity and richness were significant predictors of post-harvest diversity and richness, respectively, in all functional groups (Table 4).

TABLE 4

Table 4. Effect of harvesting intensity and climatic aridity (1981–2010) on Shannon’s diversity index and species richness of Douglas-fir dominated forests (n = 3 for all locations, except Narrows Creek where n = 1, and Alex Fraser where n = 2).

FIGURE 6

Figure 6. Effect of harvesting intensity and climatic aridity on whole community Shannon’s diversity index (A,B) and plant species richness (C,D) of mixed Douglas-fir forests (n = 3 for all locations, except Narrows Creek where n = 1, and Alex Fraser where n = 2). The counts for plant species richness include all vascular plants and bryoid species. (B,D) present standardized beta coefficients of linear mixed-effect models. Beta coefficients illustrate the effect of factors on Shannon’s diversity index (B) and on plant species richness (D) in terms of their standardized effect sizes. Circles indicate average estimates and lines indicate 95% confidence intervals. The aridity index was calculated as [mean annual temperature + 10]/[mean annual precipitation/1000]. For harvesting intensity, each coefficient represents the difference between the uncut control and a given treatment. ***P-value < 0.001; **P-value < 0.01; *P-value < 0.05.

Total species richness decreased with harvesting intensity and increased with climatic aridity, with no interaction between the factors (model R²m = 0.77; Figure 6 and Table 4). Total richness was greatest in the uncut control (average 25) and lowest in the clearcut (average 15), with similar intermediate levels (18–20) among the partial retention treatments (Figure 6). Harvesting intensity significantly reduced species richness in the bryoid functional group followed by the tree group, with greater tree species losses in more diverse, wetter regions; however, harvesting had no effect on richness of either the shrub or herb groups (Table 4 and Supplementary Table S6). Herb species richness increased with climatic aridity, but neither richness of shrub nor bryoid groups were affected by regional climate (Table 4).

Discussion

Our study revealed significant losses in total ecosystem C stocks and biodiversity with increasing harvesting intensity, and the magnitude of losses was greater in more humid, productive ecosystems than semi-arid interior Douglas-fir forests. The greatest losses occurred with clearcutting, generally followed by the seedtree, small patch, then large patch retention treatments, with nuanced differences observed between the intermediate retention treatments among climatic regions.

Total Ecosystem and Aboveground C Stocks

Clearcutting is known to dramatically reduce on-site C stocks but partial cutting that aims to protect the most vulnerable components and accelerate recovery could alleviate some of these losses. We found that clearcut harvesting resulted in the largest reductions in total ecosystem C stocks, equating conceptually with large pulse C losses from the forests. Reducing harvesting intensity, however, mitigated total ecosystem C declines, with the large patch retention resulting in the lowest loss, supporting our first hypothesis. These findings agree with Zhou et al. (2013), whose meta-analysis of 81 field studies worldwide showed that losses in aboveground C storage declined linearly with reduced tree cutting intensity. Expressed as a ratio, the change in pre- to post-harvest total C stocks ranged from 2.6 following clearcutting where all trees were removed to 1.6 following the large patch retention treatment where over half of the canopy trees were retained and intermediate-sized trees were thinned from below. With 1 year between our pre-harvest and post-harvest measurements, these losses can be considered primarily a “chain-saw” effect; in other words, they represent the immediate loss from cutting the live and dead trees, reducing the deadwood pools, and disturbing the forest floor with the logging equipment.

The magnitude of loss from clearcutting harvesting is a critical issue for global C stocks. In British Columbia, where this study was conducted, C losses from ecosystems as a result of clearcutting and other forestry practices in a single year (2017) were estimated at 42 million Mg C, in addition to the lost opportunity cost of forest carbon capture, estimated at an additional 26.5 million Mg C, for a total of 68.5 million Mg C lost annually (British Columbia Ministry of Forests Lands and Natural Resource Operations [BC MoFLNRORD], 2017). These losses exceeded the 65 million Mg C emitted the same year (2017) from fossil fuel use. In the last 3 years, forest C loss has more than tripled to a total of 203 million Mg C when we account for the 2.7 M hectares burned by wildfire from 2017 to 2019 (British Columbia Ministry of Forests Lands Natural Resource Operations and Rural Development [BC MoFLNRORD], 2019b). Even without fire, the 3.6 M hectares of forests clearcut across BC between 2005 and 2017 (British Columbia Ministry of Forests Lands and Natural Resource Operations [BC MoFLNRORD], 2017) would have accounted for approximately 0.8 billion tons of C based on our study average of 170 Mg C ha^–1 total ecosystem C loss with clearcutting. At this rate of harvest, and with other associated C fluxes in BC, Canada will not be able to meet its commitments to a 40% reduction in carbon emissions by 2030, nor to 80% by 2050, even with successful decarbonization of the energy sector. This also highlights the contrast between C sequestration outcomes with forestry in Canada versus among countries (Le Quéré et al., 2018).

In forests of the Pacific Northwest, it takes one to three decades for ecosystems to return to being a CO₂ sink from a source (Law et al., 2004; Turner et al., 2004). This recovery period is greatly reduced with partial tree retention (Juodvalkis et al., 2005; Zhou et al., 2013). We do not present modeling data on long-term C storage in this paper, but we calibrated the CBM-CFS3 (Carbon Budget Model of the Canadian Forest Sector) using our C stock data to project harvesting intensity effects 50 years into the future (Kurz and Apps, 2006; Kull et al., 2019). We found that the partial retention treatments are projected to shift being sources of C emissions to sinks within one to two decades of harvest, compared to two to three decades for the clearcutting treatment. Moreover, C stocks will reach only two-thirds of pre-harvest values in clearcuts in this timeframe compared to fully recovering and even exceeding pre-harvest levels in the large patch retention treatment. This acceleration of C stock recovery suggests that partial retention harvesting, rather than clearcutting, may be an important stewardship strategy to slow C emissions from forestry and mitigate global change in the coming critical decades (IPCC, 2018).

Losses in total ecosystem and aboveground C stocks with harvesting generally declined with climatic aridity, in keeping with our second hypothesis. The magnitude of C losses was proportionate to pre-treatment C stocks, as expected (Buotte et al., 2019a). The greatest absolute losses were in the most productive ecosystems and the least in the most arid. Total ecosystem C lost from the forests at the humid end of our climate gradient ranged from approximately 250 Mg ha^–1 following clearcutting to 150 Mg ha^–1 following large patch retention, representing a decline from a 50% loss to 30% loss of initial stocks (average 480 Mg ha^–1), respectively. In the most arid forests, the losses with harvesting ranged from 138 to 18 Mg ha^–1, representing a decline in C stocks from 60% with clearcutting to 8% with large patch retention compared with the pre-harvest levels of 230 Mg ha^–1. Thus, while low intensity harvesting greatly reduced total ecosystem C losses in the most arid forests, losses were still considerable in the large patch retention treatment in the humid climatic region.

The response of C pools to harvest treatments showed variability corresponding to the forest location on the climate gradient. For instance, while across the whole climate gradient clearcut harvesting left the lowest total ecosystem C stocks and large patch retention left the greatest, in arid forests the seedtree retention reduced C stocks more than the small patch retention harvest, and vice versa in the humid forests, where small patch reduced C stocks more than seed tree retention. This suggests, from a C stocks perspective, that small patch retention is more appropriate in dry forests, whereas the seedtree method can be as effective in more humid Douglas-fir forests. This difference across the climate gradient can be explained by the variation in the size and spatial distribution of trees. In the more humid forests where water is less limiting to tree growth, the dominant and co-dominant trees are more evenly spaced and larger than in the dry forests (Simard, 2009). Selecting the largest trees as seed-trees is easy for fallers, and these dominants represent a disproportionate amount of the C stocks in the humid versus arid forests. However, delineating small thrifty patches in these even-aged forests where there are few natural boundaries is not as straightforward, with selection of small patches based mainly upon achieving uniform spatial distribution of the retained clumps. In dry forests, by contrast, trees grow naturally in uneven-aged clumps where resources are most replete, with grassy gaps providing natural definition to the tree patches (Huggard et al., 2005; Vyse et al., 2006), making selection of retention patches fairly easy. As in humid forests, dominants in arid forests are easy for fallers to select as seed trees due to height stratification within the clumps.

Secondary Aboveground C Stocks

The understory C pool was minor across our climate gradient, storing an average of 3.5 Mg ha^–1, or less than 1% of total ecosystem C prior to harvest. In terms of the C budget, protecting these stocks is therefore of minor concern when designing harvesting treatments. The understory pool was lowest at our humid sites where stands were dense and understory light availability low and greatest at the arid sites where stands were more open and canopy gaps common. Despite the smaller pool size, the understory was comprised of many more species than the tree layer. Composition changed across the aridity gradient, shifting from moss-dominated in the productive, closed canopy forests to shrub- and herb-dominated in the middle of the climate gradient, and to grass-dominated at the arid extremes. We found that C stocks of understory plants were reduced to an average of 0.1 Mg ha^–1 1 year following clearcutting across the whole climate gradient, due to disturbance by logging machinery. In the humid and intermediate climate regions, losses also occurred in the partial cutting treatments with no clear signal that tree retention at any level protected the understory in these productive forests. In the most arid forests, by contrast, any level of tree retention helped protect the pinegrass (Calamagrostis rubescens Buckley) dominated plant community, which is naturally resilient due to the dense root network. Within 1 year of harvesting at the two most arid forest regions, C stocks of the understory plant community surpassed pre-harvest levels by 0.5–2.0 Mg ha^–1. Here, partial removal of the overstory stimulated understory C via the increase in light, water and nutrients, with pinegrass benefiting from greater soil water potential (Huggard et al., 2005; Vyse et al., 2006) and soil nitrate availability (Hope et al., 2003). We expect the understory at the other climatic regions will rebound within a few years as evidenced by lush plant communities in clearcuts throughout BC (Haeussler et al., 1990), and the increase in understory C storage by 393% in the meta-analysis of 81 forests by Zhou et al. (2013). Nevertheless, the C increase in the understory was such a small portion of ecosystem stocks, as is true in forests in general (Gilliam, 2007), that it could not come close to compensating for the C loss in other pools.

Deadwood was an important C pool in our forests, with coarse woody debris storing up to 7% and snags 9% of total ecosystem C prior to harvest, falling within ranges reported for temperate forests (Turner et al., 1995; Kranabetter, 2009; Matsuzaki et al., 2013). Gains and losses in the deadwood C pool varied considerably across climate regions, where the most humid forests gained 21.0 Mg ha^–1 following clearcutting, while the most arid forests lost 3.3 Mg ha^–1. Humid forests experienced the greatest losses in C stocks from snags, but at the same time gained more C in coarse and fine woody debris than drier forests. C stocks in snags and coarse woody debris generally declined with harvesting intensity (average −17.1 Mg ha^–1 and −5.9 Mg ha^–1, respectively), while C in stumps (+ 2.8 Mg ha^–1) and small + fine woody debris (+ 5.2 Mg ha^–1) increased. This agrees with other studies, in which thinning has resulted in large reductions in coarse woody debris (Powers et al., 2012; Zhou et al., 2013; Achat et al., 2015) while stumps and fine woody debris have increased (Bond-Lamberty and Gower, 2008). Overall, three times as much C in deadwood was lost as gained with harvesting in our forests. This shift would likely reduce habitat structure for cavity nesting birds and mammals (Aitken and Martin, 2004), alter rates of decomposition (Lassauce et al., 2011; Prescott et al., 2017), and decrease the richness of saproxylic species, which account for a considerable part of biodiversity in forests and play key roles in decomposition and nutrient cycling (Lassauce et al., 2011). Deadwood also plays important roles in providing nurse sites for regeneration of plants and tree seedlings, structure for reducing erosion and regulating stream flow, and cycling of nutrients and water (Franklin et al., 1987; Harmon et al., 2004).

Belowground C Stocks

Total belowground C stocks declined by 29% with harvesting, with 98% of the loss from the forest floor, and no significant changes occurring in the mineral soil pool, in accordance with our third hypothesis. Harvesting reduced forest floor C stocks from an average of 58 Mg ha^–1 to 21.9 Mg ha^–1, or by an average of 38.5 Mg ha^–1, or 61%, which is double the roughly 30% C losses found in the meta-analyses of harvesting effects on temperate forest floors around the world (Zhou et al., 2013; James and Harrison, 2016). This loss, which falls widely outside of logging impacts on forest floor elsewhere in the world, highlights an important area where skidding practices can be improved in Douglas-fir forests of BC. We found that increasing the intensity of tree removal did not increase the magnitude of C loss from the forest floor; losses were similar whether forests were clearcut or partially harvested. This lack of harvesting intensity impact agrees with previous reviews (Zhou et al., 2013; James and Harrison, 2016), but contrasts with studies examining high thinning intensities where losses were greater (Vesterdal et al., 1995; Jurgensen et al., 2012; Powers et al., 2012; Achat et al., 2015), and it is possible that our sampling intensity was too low to distinguish among harvesting intensities. Most of the losses in forest floor C would have occurred during harvesting, as well as thinning in the high retention treatment, with substantial mixing and displacement of soil organic matter (Cambi et al., 2015), and some of the material ending up in roadside slash piles that were later burned. Increases in temperature at the soil surface with disturbance and changes in litter quality with the growth of grasses and forbs could also have contributed to microbial activity and decomposition of soil organic matter (Prescott, 2010). That the mineral soil had no C loss in any harvesting treatment or climate region, agrees with other studies (Johnson and Curtis, 2001; Powers et al., 2012; Noormets et al., 2015; Strukelj et al., 2015; Holub and Hatten, 2019). Mineral soils are protected by the forest floor from the direct physical effects of harvesting, and increases in soil moisture with logging would mitigate changes in temperature and thus decomposition rates (Grand and Lavkulich, 2011).

The magnitude of loss in total belowground C stocks with harvesting was correlated with pre-harvest belowground C stocks, which declined with climatic aridity. The soils underlying the humid forests were coarse textured, carbon-rich Podzols and those in the arid forests were fine-textured Luvisols with thinner forest floors. These differences in soils across the climate gradient suggest there is potential for greater absolute C loss from belowground pools with harvesting in humid forests than arid forests (Saidy et al., 2013). The large magnitude difference in both above and belowground C stocks of the humid compared to arid forests, and the greater vulnerability of arid ecosystems to fire and insect disturbances in arid ecosystems (Buotte et al., 2019b), may help inform forest conservation decision-making. Buotte et al. (2019a) argue that humid coastal forests with high carbon sequestration capacity, low vulnerability to fire, and high emission potential following harvesting ought to be considered priorities for preservation.

Biodiversity

Species diversity and richness declined with harvesting intensity across the entire climatic gradient, supporting our fourth hypothesis. Losses in diversity and richness were lowest in the large patch retention treatment and increased with reductions in residual tree basal area, to where Shannon’s diversity index in clearcuts was approximately half of that in control forests. An average of 10 species disappeared with clearcutting and 5–7 with partial cutting. These immediate declines in species diversity and richness with harvesting agree with other studies and meta-analyses (Clark and Clovey, 2012; Schwenk et al., 2012; Thorn et al., 2018). Where conservation of forest-adapted plant diversity and richness is a high priority in managed forests, our results show that partial cutting with high canopy retention is the most successful strategy to meet this objective across Douglas-fir forests.

The declines in Shannon’s diversity and richness were driven more by reductions in trees and bryoids than understory shrubs and herbs. Overall, shrubs did not decline in diversity or richness with harvesting because they re-sprouted following mechanical damage. Richness and diversity of herb species were also unaffected by harvesting intensity, but we nevertheless observed general reductions in certain forest adapted species including Goodyera oblongifolia, Viola glabrum, Aralia nudicaulis and Streptopus amplexifolius. Most other species losses occurred among bryoids, declining from six to one species with clearcutting in the most humid forest and from 14 to three in the most arid. Declines in bryoid diversity and richness paralleled intensity of harvest, with the greatest reductions in the clearcuts and the least in the large patch retention treatment. Bryoids were damaged with scraping and mixing by the machines and skidding of logs, and by desiccation following removal of overstory trees. Mosses and lichens play crucial roles in water retention, nutrient cycling, and erosion mitigation, and losses of species performing these functions will be particularly important in the more arid forests (Richardson, 1981). Bryoids also provide sites for germination of seeds, habitat for vertebrates and invertebrates, and stimulate growth of cyanobacteria, algae, and lichens (Richardson, 1981), and preventing reductions in these functions with partial retention of the canopy will be important across the climate gradient.

Conclusion

Any level of harvesting reduced ecosystem C stocks and biodiversity across the entire aridity gradient. Overall, the magnitude of ecosystem C and tree species diversity losses were greatest in the most humid, productive forests, whereas loss of understory bryoid species was greatest in the arid forests, with all of these losses increasing with harvest intensity. Clearcuts performed the worst at conserving ecosystem services and the large patch retention treatments the best, with intermediate responses in the low retention treatments (seedtree and small patch retention). Where moderate levels of harvesting are sought, retention of small intact patches was better at conserving C stocks and biodiversity in arid regions, whereas the seed tree strategy generally performed as well as small patches in the more humid forests. These results suggest that as climate changes and aridity increases, clearcutting is not the treatment of choice if protection of ecosystem C stocks and biodiversity is a high priority. One of the important overall findings of this analysis was that C stocks in forest floors were highly vulnerable to loss with any level of harvesting (−61% on average) compared with mineral soils (no significant change), and the organic matter losses were greater than in studies elsewhere. Further research is needed to determine ways to reduce losses to the forest floor.

Where decisions in governments are made to conserve C stocks and biodiversity in the interest of mitigating global change at national or provincial levels, our results support the argument of Buotte et al. (2019a), who wrote that high productivity forests ought to be a priority for conservation given their high C stocks and biodiversity. Our findings support the premise that partial canopy retention can alleviate on-site losses of C stocks and biodiversity associated with clearcut harvesting, but the appropriate level of tree retention must consider interactions with regional climatic aridity and the ecosystem components (such as C and biodiversity) being conserved. Our study investigated the effect of different harvesting treatments on just two factors (C stocks and biodiversity), whereas operational treatments must be designed to balance a multitude of other factors as well, such as harvesting costs per volume of wood extracted and the ability of the treatment to provide a suitable microclimate for regeneration. Furthermore, we monitored ecosystem responses at the stand level only, but larger scale factors such as watershed hydrology and large mammal habitat quality may also need to be considered in landscape level plans. Further research is needed to integrate C stock and biodiversity changes with harvesting at different scales, and the trade-offs with losses and gains in other ecosystem services and societal values.

Data Availability Statement

The datasets generated for this study are available on request to the corresponding author.

Author Contributions

SS, BP, WR, and LL conceived the ideas and designed methodology. SS, WR, and ES collected the data. CD and WR analyzed the data. AR carried out the modeling. SS and WR led the writing of the manuscript. All authors contributed critically to the drafts and gave final approval for publication.

Funding

Site selection was funded by a Natural Sciences and Engineering Research Council of Canada Strategic Project Grant STPGP 478832 - 15 to SS, BP, and LL. Field data collection, lab analysis and preparation of this article was funded by the Forest Enhancement Society of BC (grant to WR). Funding for the field truck was provided by donation from the Michael and Jena King Foundation.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank the staff of Brinkman and Associates, the BC Ministry of Forests, Lands, Natural Resource Operations and Rural Development, the University of British Columbia, and each of the local partners (Alex Fraser Research Forest, Canfor, Harrop-Procter Community Forest, John Prince Research Forest, Kalesnikoff Lumber Co., and Malcolm Knapp Research Forest) for their contributions to locating the study sites, planning and harvesting, and planting. We also thank the many members of the Simard Lab who assisted in collecting the field data and processing the samples, the BC Ministry of Environment for soil carbon analysis, and the Canadian Forest Service for assistance in carbon budget calculations.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/ffgc.2020.00088/full#supplementary-material

References

Achat, D., Fortin, M., Landmann, G., Ringeval, B., and Augusto, L. (2015). Forest soil carbon is threatened by intensive biomass harvesting. Sci. Rep. 5:15991. doi: 10.1038/srep15991

PubMed Abstract | CrossRef Full Text | Google Scholar

Aitken, K. E. H., and Martin, K. (2004). Nest cavity availability and selection in aspen–conifer groves in a grassland landscape. Can. J. For. Res. 34, 2099–2109. doi: 10.1139/x04-086