Partial Retention of Legacy Trees Protect Mycorrhizal Inoculum Potential, Biodiversity, and Soil Resources While Promoting Natural Regeneration of Interior Douglas-Fir

Clearcutting reduces proximity to seed sources and mycorrhizal inoculum potential for regenerating seedlings. Partial retention of legacy trees and protection of refuge plants, as well as preservation of the forest floor, can maintain mycorrhizal networks that colonize germinants and improve nutrient supply. However, little is known of overstory retention levels that best protect mycorrhizal inoculum while also providing sufficient light and soil resources for seedling establishment. To quantify the effect of tree retention on seedling regeneration, refuge plants, and resource availability, we compared five harvesting methods with increasing retention of overstory trees (clearcutting (0% retention), seed tree (10% retention), 30% patch retention, 60% patch retention, and 100% retention in uncut controls) in an interior Douglas-fir-dominated forest in British Columbia. Regeneration increased with proximity to legacy trees in partially cut forests, with increasing densities of interior Douglas-fir, western redcedar, grand fir, and western hemlock seedlings with overstory tree retention. Clearcutting reduced cover of ectomycorrhizal refuge plants (from 80 to 5%) while promoting arbuscular mycorrhizal plants the year after harvest. Richness of shrubs, herbs, and mosses declined with increasing harvesting intensity, but tree richness remained at control levels. The presence of legacy trees in all partially cut treatments mitigated these losses. Light availability declined with increasing overstory cover and proximity to leave trees, but it still exceeded 1,000 W m−2 in the clearcut, seed tree and 30% retention treatments. Increasing harvesting intensity reduced aboveground and belowground C stocks, particularly in live trees and the forest floor, although forest floor losses were also substantial where thinning took place in the 60% retention treatment. The loss of forest floor carbon, along with understory plant richness with intense harvesting was likely associated with a loss of ectomycorrhizal inoculum potential. This study suggests that dispersed retention of overstory trees where seed trees are spaced ~10–20 m apart, and aggregated retention where openings are <60 m (2 tree-lengths) in width, will result in an optimal balance of seed source proximity, inoculum potential, and resource availability where seedling regeneration, plant biodiversity, and carbon stocks are protected.

Clearcutting reduces proximity to seed sources and mycorrhizal inoculum potential for regenerating seedlings. Partial retention of legacy trees and protection of refuge plants, as well as preservation of the forest floor, can maintain mycorrhizal networks that colonize germinants and improve nutrient supply. However, little is known of overstory retention levels that best protect mycorrhizal inoculum while also providing sufficient light and soil resources for seedling establishment. To quantify the effect of tree retention on seedling regeneration, refuge plants, and resource availability, we compared five harvesting methods with increasing retention of overstory trees (clearcutting (0% retention), seed tree (10% retention), 30% patch retention, 60% patch retention, and 100% retention in uncut controls) in an interior Douglas-fir-dominated forest in British Columbia. Regeneration increased with proximity to legacy trees in partially cut forests, with increasing densities of interior Douglas-fir, western redcedar, grand fir, and western hemlock seedlings with overstory tree retention. Clearcutting reduced cover of ectomycorrhizal refuge plants (from 80 to 5%) while promoting arbuscular mycorrhizal plants the year after harvest. Richness of shrubs, herbs, and mosses declined with increasing harvesting intensity, but tree richness remained at control levels. The presence of legacy trees in all partially cut treatments mitigated these losses. Light availability declined with increasing overstory cover and proximity to leave trees, but it still exceeded 1,000 W m −2 in the clearcut, seed tree and 30% retention treatments. Increasing harvesting intensity reduced aboveground and belowground C stocks, particularly in live trees and the forest floor, although forest floor losses were also substantial where thinning took place in the 60% retention treatment. The loss of forest floor carbon, along with understory plant richness with intense harvesting was likely associated with a loss of ectomycorrhizal inoculum potential. This study suggests that dispersed retention of INTRODUCTION Clearcut harvesting, or removal of all trees from forests, reduces mycorrhizal inoculum (Perry et al., 1987;Jones et al., 2003;Cline et al., 2007) because the symbiotic fungi are dependent on living trees and plants for photosynthate (Smith and Read, 2008). Without a supply of carbon from their hosts, mycorrhizal fungi cannot grow a mycelial network to explore the soil for nutrients and water, and in turn supply them to their partner trees and plants (Smith and Read, 2008). Low levels of mycorrhizal inoculum can be especially problematic for seedlings regenerating in clearcut forests, where early colonization with mycorrhizal fungi is crucial to survival and growth of germinants (Barker et al., 2013). In temperate forests, where most trees form obligate associations with ectomycorrhizal fungi, species composition of the fungal community also changes with clearcutting, where there is a shift toward a low diversity of "early stage" fungi that have small carbon demands and generalist strategies for nutrient acquisition. This simplification of the fungal community is related to the change in host plant species composition and age, along with losses or redistribution of the forest floor by logging and site preparation, as well as modification of soil physicochemical properties and changes in the soil foodweb (Jones et al., 2003). The loss of fungal biodiversity can lead to reduced nutrient and water uptake and resilience of regenerating forests (Tomao et al., 2020).
Along with reductions in mycorrhizal fungal inoculum, clearcutting also reduces host plant diversity as well as carbon and nutrient stocks (Halpern and Spies, 1995;Johnson and Curtis, 2001;Zhou et al., 2013;Hume et al., 2017;Buotte et al., 2019;Ding et al., 2019;Curzon et al., 2020). Clearcutting effects on C stocks result directly from removal of trees and damage to the forest floor and understory plants , but also indirectly through loss of mycorrhizal fungi, which are the primary vectors of plant carbon into the soil (Talbot et al., 2008). Extensive mechanized clearcutting is of major concern globally because forests store 80% of terrestrial carbon and are home to more than three-quarters of the world's vertebrate and invertebrate species (Lindenmayer and Franklin, 2002;FAO, 2018), and failures to regenerate the clearcuts can create persistent deficits in these ecosystem goods and services. Clearcutting and regeneration of the world's most productive forests are of particular concern because these forests have the greatest biodiversity and C storage (Thompson et al., 2012;Liang et al., 2016;Lecina-Diaz et al., 2018;Buotte et al., 2019).
Partial cutting, or variable retention harvesting, has become more commonplace in western Canada and the Pacific Northwest of the USA in response to the negative ecological outcomes of extensive clearcutting and forest simplification, as well as the intensely negative effects these have had on public perceptions of forest management practices (Beese et al., 2019;Franklin and Donato, 2020). Variable retention harvesting aims to mimic natural disturbances and processes that provide the structural heterogeneity and ecosystem functions associated with mature and old growth forests (Franklin and Donato, 2020). Legacies, including trees, plants, snags, downed wood, understory plants or forest floor that survive the disturbance, provide continuity in gene pools and sustain complexity while influencing succession and promoting biodiversity, recovery, and resilience (Michel and Winter, 2009;Lindenmayer et al., 2012;Baker et al., 2013).
Compared with clearcutting, retention of a portion of legacy trees during harvest can mitigate losses in fungal inoculum, host plant diversity and C pools, and in so doing, promote seedling regeneration and forest recovery. 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 natural regeneration, protection of plants, trees and forest floor that support mycorrhizal fungi and soil foodwebs, particularly where custom hand falling instead of feller-bunchers are used, enhancement of growth of the remaining trees, and improvement in habitat for forest dwellers (Kneeshaw et al., 2002;Lindo and Visser, 2003;Hannam et al., 2006). Natural regeneration of partially cut forests rather than planting of clearcuts is attractive not only due to lower costs, but because the recovering forests are often more species rich and the seedlings and mycorrhizal symbionts are better adapted to local conditions (Swanson et al., 2011;Kranabetter et al., 2015). Seedlings can be rapidly colonized by a diversity of fungi when they germinate in close proximity to the residual trees, whose ectomycorrhizal roots and hyphal networks contact and infect the fledgling germinant roots.
The low degree of host specificity of many ectomycorrhizal fungi (Molina et al., 1992;Bruns, 1995;Horton and Bruns, 1998) means that mycorrhizal networks of different species of residual trees and plants can colonize regenerating seedlings. Studies show that mycorrhizal fungal colonization and diversity drop precipitously with distance of host seedlings from trees along forest edges, with the best outcomes within 5 m of residual trees, and the worst in the center of clearcut openings Hagerman et al., 1999;Kranabetter et al., 1999;Teste and Simard, 2008;Bingham and Simard, 2011). The retention of legacy trees and refuge plants with harvest can not only provide fungal inoculum, but also increases host plant species diversity and structural heterogeneity of the recovering forest (Franklin et al., 1987;Gustafsson and Perhans, 2020), while reducing aboveground C loss, increasing underground C storage, and reducing regeneration time (Zhou et al., 2013). The appropriate harvesting intensity for protecting these values in different forest ecosystems is poorly understood, however, and is expected to vary with local conditions. Interior Douglas-fir is a moderately shade tolerant, early to late successional tree species and is able to regenerate under a full or partial canopy in the southern interior of British Columbia. It grows throughout the territory of the Secwepemc Nation, who know it as cq•ctp, and who use its strong wood and boughs for flooring and pit cooking, the seeds and sugars for food, the needles for tea, and the pitch for medicinal salves (Turner et al., 1992). The species is a foundational tree in pure and mixed forests, where it dominates the upper canopy for centuries, and is a highly valued timber species to the forest industry (Parish et al., 1996). It is a legacy tree because of its great size, resistance to fire, ability to regenerate prolifically from seed, and foundational role in the structure and function of forests (Simard, 2009).
Interior Douglas-fir is host to thousands of ectomycorrhizal fungal species (Trappe, 1977), with many shared in common with the other ectomycorrhizal tree and shrub species in mixed forests (Molina et al., 1992;Hagerman and Durall, 2004). Some hostspecific fungi such as Rhizopogon vinicolor and R. vesiculosus that exclusively colonize interior Douglas-fir will colonize trees of all ages (Molina and Trappe, 1994). In mapping the Rhizopogon networks of multi-aged interior Douglas-fir forests, Beiler et al. (2010) found that seedlings germinated in the network of the older trees. The oldest trees served as hubs or nuclei of the mycorrhizal network, where they formed links with most other trees and seedlings in a scale-free pattern. Douglas-fir germinants had significantly greater survival and growth where they could tap into the network than where they could not, and this was associated with greater fungal diversity and carbon, nitrogen, or water transferred from the old trees (Teste et al., 2009;Bingham and Simard, 2011). Moreover, Douglas-fir were able to recognize relatedness of the nearby regenerating seedlings and provide more carbon to relatives to facilitate greater establishment success (Pickles et al., 2017;Asay et al., 2020). The nurturing effect of the hub trees on seedling regeneration has since inspired the colloquial naming of these trees as "mother trees" (Simard, 2017).
The Mother Tree Project was established in British Columbia to compare recovery of the interior Douglas-fir forests, including seedling regeneration, plant diversity, and carbon pools, following clearcutting and partial cutting practices . Five harvesting intensity treatments were applied in a highly productive forest that retained Douglas-fir mother trees in different densities and configurations. These included: (1) clearcutting, where none of the trees were retained, (2) seed tree retention, where about 10% were retained in a dispersed pattern to provide seed and mycorrhizal fungal inoculum, (3) 30% patch retention, where small patches of trees were retained to provide seed and fungal inoculum for regeneration in the surrounding open gaps, (4) 60% patch retention with thinning from below, where seedlings could establish under the protection and within the networks of the old trees, and (5) control, where the forests were left to develop in the absence of logging disturbance.
We tested the following hypotheses 3 years following the harvest of this experiment. (1) Total ingress of natural regeneration will increase with overstory tree retention. (2) Douglas-fir regeneration will increase with proximity to legacy trees. (3) Diversity and richness of natural regeneration will be greatest at intermediate retention levels because of the environmental heterogeneity. (4) Retention of legacy trees will protect plant richness and diversity, whereas clearcutting will reduce ectomycorrhizal shrubs and promote dominance of arbuscular mycorrhizal pioneers. (5) Light and soil water availability will decline but forest floor substrate will be better protected with increasing tree retention. (6) Carbon pools, particularly live trees, understory plants, and forest floor, will decline with increasing harvesting intensity. In testing these hypotheses, we sought to determine the tree retention levels at which natural regeneration densities are balanced with resource availability and mycorrhizal fungal inoculum potential (i.e., the potential for the mycorrhizal fungal inoculum of refuge plants to colonize the root tips of regenerating seedlings).

Study Area
The study took place at an interior Douglas-fir-dominated forest situated near the town of Nelson in British Columbia (B.C.), Canada (49.63 • N,117.03 • W). The site is located within the West Kootenay Dry Warm Interior Cedar-Hemlock (ICHdw1) biogeoclimatic variant, with mean annual temperature (MAT) of 6.8 • C and mean annual precipitation (MAP) of 868 mm. The climate is continental with warm, dry summers and cool, snowy winters.
The site had a medium soil moisture regime, southeast aspect, 10-30% slope gradient and occupied a mid-slope position. Elevation averaged 850 m. Soils were predominantly Humo-ferric Podzols with a sandy loam texture, and humus forms were Moders (Green et al., 1993;Soil Classification Working Group, 1998). The biogeoclimatic site series were ICHdw1-101 and -104.

Harvesting Treatments
Treatments were laid out in a completely randomized design with five harvesting intensity treatments replicated 4 times (total 20 treatment units). Harvesting treatments were randomly assigned to 5 ha treatment units and logging was conducted prior to spring green-up in 2017. The goal of the harvesting treatments was to create a gradient of legacy tree densities where large Douglas-fir were targeted for retention, but alternate species were also selected where Douglas-fir was absent. The treatments included: (i) clearcut (no tree retention); (ii) single tree retention (retention of 25 large Douglas-fir stems ha −1 at 25 m spacing, or ∼10% dispersed retention); (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) high retention (retention of 60% of the stand area with all trees cut in the remaining 40% of the stand); and (v) uncut control (100% retention). Harvesting was carried out using fellerbunchers run along random skid trails, but trees too large for the machine were hand-felled. The stands in the high retention treatment were thinned from below by reaching into the uncut patches with the feller-buncher and removing the smaller stems. Trees were limbed and topped before the boles were skidded to roadside landings.

Measurement and Sampling Methods
A one-hectare measurement plot was established in the center of each treatment unit, in which all measurements and sampling were conducted. Inside the measurement plot, a nine X nine grid of 81.10 m 2 (r = 1.78 m) subplots were established in which natural regeneration and environmental conditions were assessed 3 years after treatment. All seedlings that had regenerated since harvesting (i.e., aged 3 years or less) were counted by species, the average height by regenerating species measured, distance to the nearest Douglas-fir seed tree determined, and percent cover of the overstory trees and each regeneration substrate [mineral, organic, highly decayed wood (decay classes 4 and 5; British Columbia Ministry of Forests Range British Columbia Ministry of Environment, 2010), or recently decayed wood (decay classes 1, 2, and 3), hereafter referred to as "fresh wood"] visually estimated. The distance to the nearest Douglas-fir seed tree was not necessarily related to percent cover of the overstory trees, because the overstory contained a mix of species and ages, with Douglas-fir seed trees not always present where the overstory was dense. In the center of 10 randomly selected subplots per treatment unit, photosynthetically active radiation (PAR) and volumetric water content were also evaluated. PAR was measured in the four cardinal directions using a Sunfleck PAR Ceptometer (Model SF-80, Decagon Devices, Pullman, WA) and averaged. Volumetric water content was measured in mineral soil to 15 cm depth using a Dynamax ML2 Theta Probe (Delta-T Devices Ltd.), which measures changes to the dielectric constant of the soil as a proxy for volumetric soil water content (m 3 /m 3 ) (Delta-T Devices Ltd, 1999). Volumetric water content and PAR were measured over a 2-day period in late August, following 2 weeks of warm, dry weather.
Carbon and plant diversity data were collected in one central permanent circular 0.04 ha (r = 11.28 m) National Forest Inventory (NFI) plot per treatment unit in July/August prior to logging, and these plots were re-assessed both one and three summers following harvesting. Data collection followed the NFI ground sampling guidelines (Canadian Forest Inventory Committee Canadian Forest Service, 2008) as follows. 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. Nested within each 0.04 ha plot was one circular 50 m 2 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 2 subplot (r = 10.0 m) and of all herbaceous, ground-level shrub and bryoid species occurring within a 100 m 2 (r = 5.64 m) subplot. Each plant species was identified by its known dominant mycorrhizal functional group (ectomycorrhizal, or arbuscular, arbutoid, orchid, and ericoid mycorrhizal). 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. 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 2 (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 previous microplot locations. All bryoids, herbs, trees/shrubs <1.3 m in height, and fine woody debris (FWD) (<1 cm diameter) 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 10cm 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 then 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 at least 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, 2019) and soil moisture regime, nutrient status and site series were estimated using vegetation, soil, and site features (Braumandl and Curran, 1992).

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, roots, gravel, and charcoal removed and weighed, then sieved into <8 mm and >8 fractions, which were also weighed. All oven drying was done at 70 • C for 72 h. Subsamples of each <2 mm mineral and <8 mm forest floor fraction were sent to the Ministry of Environment laboratory in Victoria, B.C. 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 using the four replicate treatment units per treatment. Carbon storage (Mg ha −1 ) in aboveground dead and live trees >1.3 m tall, stumps (<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. 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, as suggested in Table 4.4 (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. For the regeneration and resource availability measures, the subplot measurements were averaged for each treatment unit, and these means were used to compare the parameters among the harvesting intensity treatments.
Statistical analyses were conducted in Excel and results were considered statistically significant at P < 0.05. We investigated how natural regeneration, plant diversity, mycorrhizal functional groups, water and light availability, and carbon stocks responded to legacy tree retention [0% retention (clearcut), 10% retention (seed tree), 30% retention (small patch), 60% retention (high retention), and 100% retention (control)] using one-way ANOVA or, where responses were measured more than once, ANCOVA (n = 4). The general form of the ANCOVA model was: where Yi is the response variable; µ is the general mean; τ i is the fixed effect parameter for the ith harvesting intensity treatment; βi is the estimated coefficient, Xi is the covariate, and εi is the residual (Steel and Torrie, 1980). Means were compared using the Tukey-Kramer post-hoc test. Linear relationships between variables were measured using the Pearson product-moment correlation coefficient.  Table 1.

RESULTS
Three years after harvesting, total natural regeneration density was lowest in the clearcuts and generally increased with overstory tree retention to where it was greatest in the uncut controls, with the exception of the 30% patch retention treatment where the harvested gaps behaved similarly to the clearcuts (Figure 1, Table 1). The trends were species-specific, however, where regeneration densities of coniferous species (interior Douglasfir, western hemlock, western redcedar, and grand fir) tended to be greatest in the 60 and 100% (control) retention treatments, whereas those of broadleaf tree species (paper birch, black cottonwood [Populus balsamifera ssp. trichocarpa (Torr. Et A. Gray) Brayshaw] and trembling aspen (Populus tremuloides Michx.) tended to be greatest in the 0% (clearcut) and 10% (seed tree) retention treatments ( Table 1). Height of interior Douglasfir seedlings significantly declined with increasing overstory tree retention, whereas western redcedar, western hemlock and grand fir tended to be tallest in the 30 or 60% retention treatments ( Table 1). Although the density of broadleaf tree regeneration was less than that of conifer seedlings, the broadleaves were substantially taller. Species richness of seedling regeneration did not significantly differ among the tree retention treatments ( Table 1).
Douglas-fir regeneration density declined with distance from the nearest Douglas-fir seed tree, with the greatest densities within 10 m and declining by more than half distally (Figure 2). Regeneration density continued to decline from 20 to 60 m, beyond which germinants were rare. Germinant density of Douglas-fir was unrelated to cover of the substrates except that it was negatively affected by fresh wood (r = −0.16), which was most abundant in the low tree retention treatments. Density of grand fir generally increased with organic substrate cover (r = 0.15) and decreased with fresh wood (−0.15), while western hemlock tended to regenerate more densely on highly decayed 1 | Effect of tree retention (mean ± one standard error) on overstory cover, photosynthetically active radiation (PAR), volumetric water content, density, height, and species richness of natural regeneration, and cover of seedbed substrates in the ICHdw1 variant at Redfish Creek (n = 4). wood (r = 0.15). Density of western redcedar germinants was equal among substrates. Species richness of bryoid, herb and shrub functional groups 3 years after harvesting was lowest in clearcuts, and generally increased with overstory tree retention ( Figure 3A, Table 2). Diversity of bryoids and herbs also tended to be lower in clearcuts, but differences between treatments were not significant ( Figure 3B, Table 2). Tree species richness and diversity varied little by harvesting treatment after 3 years when all ages of trees, including germinants, were included in the functional group ( Table 2).

Response
Prior to harvesting, the cover of plant species that were ectomycorrhizal averaged 80.6%, followed by arbuscular (45.2%), arbutoid (3.7%), orchid (0.9%), and ericoid (0.1%) mycorrhizal (Figure 4). One year post-harvest, cover of plants declined with harvesting in all tree retention treatments, but the most pro-found impact occurred in the clearcuts, where cover of ectomycorrhizal plants declined to 6.1%, or less than half the cover of arbuscular mycorrhizal plants (13.4%). Arbuscular mycorrhizal plants were also more prolific (42.0%) than ectomycorrhizal plants (31.2%) in the 30% patch retention treatment 1 year-post harvest, whereas ectomycorrhizal plants remained dominant in the seed tree, 60% retention and control treatments. By the third year post-treatment, ectomycorrhizal plants had recovered to 24.1% cover in the clearcuts and 28.6-40.7% in the higher retention treatments, and surpassed cover of arbuscular mycorrhizal plants in most cases. Cover of arbutoid, orchid and ericoid plants declined in all of the harvesting treatments, but ericoid plants recovered to pretreatment levels in the 60% retention treatment within 3 years. Photosynthetically active radiation significantly increased from 166 W m −2 in the uncut control to over 1,000 W m −2 in the clearcut, seed tree, and 30% patch retention treatments 3 years after harvesting ( Figure 5A, Table 1). Light availability increased exponentially with declining overstory cover, with levels exceeding 1,000 W m −2 below 30% cover ( Figure 5B). Volumetric water content tended to be greatest in clearcuts, but differences among treatments were not significant and values were low in all treatments ( Table 1).
Total ecosystem C stocks 1 year after harvesting averaged 242.6 ± 43.3 Mg ha −1 in the uncut controls, with 168.5 ± 37.5 Mg ha −1 aboveground and 74.1 ± 5.9 Mg ha −1 belowground (69 and 31%, respectively) (Figure 6, Table 3). The average portion in each pool was live trees (49.3%), dead trees (11.2%), coarse woody debris (7.0%), small and fine woody debris (1.7%), stumps (<0.1%), understory plants (<0.1%), tree roots (14.5%), forest floor (8.4%), and mineral soil to 15 cm (7.6%). Prior to treatment there were no significant differences in C stocks among the five treatments, except C in coarse and total woody debris was highest in the 30 and 60% retention treatments (p < 0.05; Table 3). Post-harvest, aboveground C stocks were lowest in the most intensive harvesting treatment (clearcutting) and highest  Table 2. for the least intensive treatment (60% retention) and controls. Post-harvest total ecosystem C stocks in the clearcut were onethird the pre-harvest value, while three-quarters of the total ecosystem C stocks was retained following the 60% retention treatment. Forest floor C declined in the clearcut, seed tree and 60% retention treatments, but not in the 30% patch retention treatments where more of the coarse woody debris was mixed into the forest floor by the logging equipment. Mineral soil C to 15 cm depth was unaffected by harvesting (Figure 6, Table 3). The total C stocks in coarse woody debris was unaffected by harvesting, but there was significantly more small woody debris in the low retention treatments, and greater fine woody debris in the harvesting treatments than the control, representing slash left from the logging operations. The understory plant community represented a minor C pool in all treatments and had recovered to only 14-20% of control levels after 3 years in the intermediate retention treatments (Figure 5, Table 3).

DISCUSSION
Our study revealed a large reduction in ectomycorrhizal fungal inoculum potential with clearcut harvesting, and this was  associated with considerably lower regeneration of interior Douglas-fir in spite of greater light availability than in treatments where various densities of overstory trees were retained. The reduction in ectomycorrhizal fungal inoculum was the direct result of removal of legacy trees and refuge plants as well as loss of forest floor C stocks. Retention of legacy trees in dispersed or small aggregate patterns mitigated most of these effects and promoted an abundance of conifer regeneration.

Regeneration
Total ingress of natural regeneration increased with tree retention in the range of harvesting treatments, in keeping with our first hypothesis. This agrees with many studies showing increased recruitment and regeneration density with retention harvesting (Coates, 2002;Newsome et al., 2010;Gottesman and Keeton, 2017) or thinning (Kuehne and Puettmann, 2008;Roberts and Harrington, 2008;Shatford et al., 2009;Dodson et al., 2012;Olson et al., 2014;Gauthier and Tremblay, 2019;  Williams and Powers, 2019). The regeneration response we observed with partial overstory retention was associated with changing resource levels, providing better growing conditions for leave trees and natural regeneration. For example, partial removal of neighboring trees increased light and water availability and this would have promoted greater seed production of leave trees (Reukema, 1982;LePage et al., 2000;Claveau et al., 2006;Kuehne and Puettmann, 2008;Roberts and Harrington, 2008;Devine and Harrington, 2009;Shatford et al., 2009;Reich et al., 2012;Urgenson et al., 2013;Zhou et al., 2016;de Montigny et al., 2018;Willis et al., 2018;Williams and Powers, 2019). Additionally, soil disturbance resulting from harvesting itself produced a range of suitable substrates for seedling establishment (LePage et al., 2000;Roberts and Harrington, 2008;Shatford et al., 2009;Dodson et al., 2012).
Interior Douglas-fir seedlings in particular established in the greatest abundance within 10 meters of retained legacy trees, beyond which density generally declined by more the half, then to very low levels at 60 m and further. This result is concordant with our second hypothesis that density of Douglas-fir regeneration would increase with proximity to seed trees because this species disperses its seed mainly by wind. It also agrees with many studies of wind dispersed species showing a rapid decline in seed fall from the seed source, with the majority of seed falling within 25 m (e.g., Ferguson and Carlson, 1990;LePage et al., 2000;Whyte and Halpern, 2019). In agreement with previous studies, Douglas-fir ingress still occurred up to 100 meters distant of leave trees, but in very low abundance (Whyte and Halpern, 2019). In all treatments, including the clearcuts, density of conifer regeneration exceeded target stocking standards 3 years postharvest (Breshears et al., 1997). If we assume 50% mortality in the first decade, it is likely these standards would still be met when at the early free-to-grow assessment date of 12 years in the seed tree and patch retention treatments, but not in the clearcut (British Columbia Ministry of Forests Range British Columbia Ministry of Environment, 2000). In the clearcuts, stocking of Douglasfir was generally below targets where distance from the nearest seed tree exceeded 60 m. Whereas, distance to seed source was a strong factor in Douglas-fir ingress density, there was little effect of substrate quality, with the exception that fresh slash inhibited germination of seed. While the harvesting produced suitable substrates for seedling establishment (Roberts and Harrington, 2008;Shatford et al., 2009;Dodson et al., 2012), Douglas-fir does regenerate on a wide range of seedbeds (Ferguson and Carlson, 1990).
Interior Douglas-fir is moderately shade tolerant, with germinants surviving best where they are protected from the sun and drought (Ferguson and Carlson, 1990). While ingress  Table 1. of this species was expected to be high in the seed tree and 30% patch treatments where light levels exceeded 1,000 W m −2 , it's high densities in the 60% retention treatment and uncut control where light intensity was below 400 W m −2 was surprising given that the species performs well in open conditions in favorable climates. Most of the recruitment in these high density treatments was concentrated in gaps that had been created by harvesting or natural disturbances such as root disease, windthrow, and bark beetles. Height growth patterns were congruent with our expectations, however, where interior Douglas-fir seedlings in the high retention treatments were half the size of those in the more open clearcuts, seed tree, and 30% patch retention treatments.
Ingress of late seral, shade tolerant species (western redcedar, western hemlock, and grand fir) also increased with overstory tree retention, in keeping with other studies (Zald et al., 2008), but the differences between treatments was only significant for grand fir. All of these species establish readily in the shade and generally on forest floor or decayed wood, although western redcedar also establishes well in full light on exposed mineral soil (Klinka et al., 2000). Ingress and height growth of western redcedar was by far greater in the 60% retention treatment, where the small openings provided the best balance of shade and substrate. Multiple studies have shown that ingress density of shade tolerant species increases with any amount of thinning (Olson et al., 2014;James, 2016;Puettmann et al., 2016), but with shifts toward later successional species at higher retention densities (Nabel et al., 2013;Crotteau et al., 2018;Gauthier and Tremblay, 2019).
The shade intolerant pioneer broadleaf species, including paper birch, black cottonwood and trembling aspen, generally germinated, sprouted, or suckered most abundantly following the clearcut and seed tree treatment, and were taller following all types of logging than in the control, as expected. Low densities of paper birch and trembling aspen also occurred in the larger gaps in the 60% retention and control treatments where seed sources were nearby. Ingress density of early-seral and gap-dependent species have been responsive to greater light conditions in larger gaps created by harvesting or natural disturbances in many studies (Coates, 2002;Kuehne and Puettmann, 2008;Roberts and Harrington, 2008;Newsome et al., 2010;Nabel et al., 2013;Urgenson et al., 2013;Olson et al., 2014;Zhu et al., 2014;James, 2016;Poirier et al., 2018). Dispersed retention with low to moderate tree density has consistently yielded abundant ingress of early-seral species such as the broadleaf species in our study (Aubry et al., 2009;Urgenson et al., 2013;Crotteau et al., 2018;de Montigny et al., 2018), with increasing seedling density with decreasing overstory density, or where gaps in the crown are present in dispersed (Ares et al., 2010;de Montigny et al., 2018) or aggregate retention treatments (Roberts and Harrington, 2008;Crotteau et al., 2018). Heavy thinning maximizes ingress density of shade intolerant trees (Zald et al., 2008;Ares et al., 2010;Otto et al., 2012), due to proficient seed production and dispersal (Kuehne and Puettmann, 2008;Urgenson et al., 2013;de Montigny et al., 2018). As crowns expand in thinned treatments, later regeneration of shade intolerant species has been inhibited (Palik et al., 2003;Maguire et al., 2007).
Species richness of seedling regeneration ingress was unaffected by harvesting intensity, but rather was reflective overall of the surrounding forest composition. We had expected to find the open and shaded conditions present in the 30% patch retention treatment to favor shade intolerant and shade tolerant species, respectively, as found in studies of aggregate retention elsewhere (Fedrowitz et al., 2014;Gottesman and Keeton, 2017;Martínez Pastur et al., 2019), but this was not reflected in our analysis. This led us to reject our third hypothesis that diversity and richness of natural regeneration would be greatest at intermediate retention levels because of environmental heterogeneity created by the harvesting treatments. Instead, sufficiently variable conditions were created across all of the harvesting treatments to support high species richness. More importantly, the species rich tree canopy of the entire forest seeded an equally rich regeneration layer as found in previous research (Clark and Covey, 2012). Many studies have found that tree species diversity and composition tend to either increase or remain unchanged with variable retention harvesting, relative to clear-cuts and mature stands (Fedrowitz et al., 2014;de FIGURE 6 | Distribution of carbon stocks among above and belowground pools in the ICHdwl variant at Redfish Creek 1 year after treatment (n = 4). Mean carbon pools are presented for each harvesting treatment across the legacy tree retention gradient. Belowground pools are represented with brown shading and aboveground pools with green shading. ANCOVA test statistics are shown in Table 3. Montigny et al., 2018). Ares et al. (2010) found that species richness increased with all thinning treatments (particularly so at lower rates of retention), with this finding corroborated in other studies (Chan et al., 2006;Otto et al., 2012;Zhou et al., 2016). Even so, we did observe an overall preponderance of early seral, shade intolerant broadleaf species in the low retention treatments, and greater abundance of late-seral shade tolerant species in the higher retention treatments, as found by Urgenson et al. (2013). It is possible that the abundant and rapid regeneration of the shade-tolerant grand fir, western hemlock and western redcedar, as well as the lower light conditions in the high retention treatments could have impeded the growth of the more shade-intolerant species (Shatford et al., 2009;Dodson et al., 2012Dodson et al., , 2014Nabel et al., 2013;Curtis et al., 2017;Crotteau et al., 2018).
The high richness of regenerating tree species in our harvesting treatments ought to promote ecosystem productivity due to niche partitioning and synergies among species (Tilman, 1999;Ares et al., 2010;Thompson et al., 2011;Forrester and Bauhus, 2016;Liang et al., 2016). The greater broadleaf densities in the lower retention treatments in particular will help protect the conifers against root pathogens according to earlier studies in these forests (Baleshta et al., 2005;Simard et al., 2005), enhance nutrient availability (Prescott et al., 2004;Zhang et al., 2012), and promote bird diversity (Aitken and Martin, 2007). The enhancement of shade tolerant species in small gaps of the higher retention treatments will also promote a shift toward later successional stand structures, which would provide habitat for more old-growth dependent plants, birds, and mammals .
Legacy interior Douglas-fir trees support mycorrhizal networks that connect with regenerating seedlings, and their presence will not only result in faster fungal colonization of the growing root tips (Barker et al., 2013), but with a greater diversity of fungi that can access heterogeneous nutrient pools (Cline et al., 2005). In the dispersed and patch retention treatments, the abundant germinants in the seed shadow of retained trees would have benefited from the mycorrhizal networks of the old trees not only by enhanced colonization, but by access to the network's extensive uptake capacity and transmission of organic nutrients and water from the old trees (Simard et al., 2012). Studies in partially cut forests show that seedlings of Douglas-fir and other ectomycorrhizal tree species have the greatest ectomycorrhizal diversity and survivorship within 10 m of residual trees (Kranabetter and Wylie, 1998;Hagerman et al., 1999;Cline et al., 2005;Teste and Simard, 2008;Bingham and Simard, 2011), corresponding closely in our study with the distribution of regenerating seedlings in proximity to legacy trees. Whether legacy trees are retained in dispersed or aggregate patterns at the low retention levels is likely important for ectomycorrhizal species with short dispersal distance, whereas these taxa may be lost in large gaps (Luoma et al., 2004).
Thinning to the degree we conducted in our 60% patch retention treatment, on the other hand, ought to have little effect on ectomycorrhizal diversity based on findings in other thinning studies (Peter et al., 2003;Buée et al., 2005;Lazaruk et al., 2005;Mosca et al., 2007;Choi et al., 2014;Castaño et al., 2018). Both Lazaruk et al. (2005) and Luoma et al. (2004) recommend a combination of dispersed and aggregate retention for preserving ectomycorrhizal fungal diversity.

Plant Communities and Refuge Plants
The reduction in richness of shrub, herb and bryoid species in the clearcut and seed tree treatments compared with the control and high retention treatments supports our fourth hypothesis and reflects the extensive ground disturbance and altered microclimatic conditions associated with the intense harvests (Bartels et al., 2018). Ectomycorrhizal shrub and broadleaf species declined from 81 to 6% cover with clearcutting, including loss of kinnickinick (Arctostaphylos uva-ursi), Bebb's willow (Salix bebbiana), alder (Alnus viridis subsp. sinuata (Regel) Ä. Löve & and D. Löve), paper birch (Betula papyrifera Marsh.), and trembling aspen (Populus tremuloides Michx.), which are known to share ectomycorrhizal fungal species in common with interior Douglas-fir (Hagerman and Durall, 2004). Loss of cover of these shrubs and trees in the clearcuts would have been accompanied by a dramatic decline in ectomycorrhizal fruiting bodies, particularly of late-stage fungi, based on studies in similar forests (Bradbury et al., 1998;Durall et al., 1999Durall et al., , 2006. Arbuscular mycorrhizal plant cover also declined from 42% preharvest to 13% in the clearcuts predominantly fireweed (Epilobium angustifolium L.), pinegrass (Calamagrostis rubscens), white hawkweed (Hierarcium albiflorum) and Hooker's fairybells (Disporum hookeri)], but they remained relatively dominant or co-dominant with the ectomycorrhizal shrubs for 4 years. Clearcutting has commonly been found to encourage an influx of arbuscular mycorrhizal pioneer species, including herbs and grasses (Bradbury, 2004;Aikens et al., 2007;Craig and MacDonald, 2009). These shade-intolerant colonizers have competitive advantages in dispersal and reproduction allowing them to dominate in areas with limited crown cover and high levels of forest floor disturbance (Beese and Bryant, 1999;Newmaster and Bell, 2002;Halpern et al., 2005;Lencinas et al., 2011;Arsenault et al., 2012;Baker et al., 2013;Haughian, 2018).
By contrast, retention of legacy trees in the 30 and 60% patch retention treatments better protected a variety ectomycorrhizal shrubs and trees as well as forest floor that are sources of ectomycorrhizal fungal inoculum for regenerating interior Douglas-fir, in support of our fourth hypothesis. By 3 years post-treatment, shrub species richness had recovered to control levels in these treatments, and cover values were reaching 40%. Ectomycorrhizal shrubs compete well in more intact ecological conditions because their fungal symbionts can access organic nutrients from the forest floor. This is supported by studies finding late-stage ectomycorrhizal fungi fruiting within 7 meters of residual trees in partially cut forests , and seedlings colonized by a variety of fungal species when they are planted within a few meters of the forest edge (Hagerman et al., 1999). They can also share organic compounds with conifer seedlings linked into their mycorrhizal networks (Simard et al., 2012). Later successional understory arbuscular species like Pacific yew (Taxus brevifolia) and Douglas maple (Acer glabrum var. douglasii) were also protected and their reproductive success was possibly enhanced by the aggregates, and the arbuscular networks of these species likely played a role in the establishment and colonization of the abundant western redcedar germinants.
Richness and abundance of bryoids were maintained in the forest patches of the 30% and 60% retention treatments, where forest floor and coarse woody debris substrates were best protected and the microclimate was wetter, cooler and shadier; by contrast, mosses were almost eliminated from the clearcut and seed tree treatments, as commonly found in other studies (Arsenault et al., 2012;Baker et al., 2013;Bartels et al., 2018;Curzon et al., 2020). Thus, our aggregated retention treatments generated sufficient resource and structural heterogeneity to support rich communities not only of tree, shrub, and herb species, but mosses as well (Clark and Covey, 2012;Haeussler et al., 2017). It is also likely that nearby source populations in adjacent intact forests and retained patches accelerated recolonization and re-establishment of these species in the small harvested gaps, as found in other studies (Halpern et al., 2005;Lencinas et al., 2011;Baker et al., 2016;Hu et al., 2018).

Resource Availability
Light availability declined with increasing tree retention, in keeping with our fifth hypothesis, where it exceeded 1,000 W m −2 in the clearcut, seed tree and 30% retention treatments, but was under 400 W m −2 , in the 60% retention treatment and uncut control. That conifer regeneration densities increased with overstory cover indicates that light was not the most important limiting factor in seedling establishment. However, growth of interior Douglas-fir and the broadleaf trees benefited from the higher light conditions in the clearcut and seed tree treatments, whereas the shady condition of the 60% retention treatment was optimal for cedar height growth.
Soil water availability did not vary significantly among the retention treatments, leading us to reject that portion of the sixth hypothesis, and conclude that soil moisture also did not limit seedling establishment. Values were low in all treatments because measurements were taken during the dry season in late August. Soil water levels did tend to rise in the clearcut where there was complete loss of overstory cover, lower rainfall interception and reduced evapotranspiration, as has occurred in previous studies (Bosch and Hewlett, 1982;Brown et al., 2005), but the increases were so variable among subplots that the response did not differ from the retention treatments. We may have been unable to discern changes in soil water content because our sample size was too small (30 samples per retention treatment) to capture the variability in below-canopy microclimates (Chen et al., 1993;Breshears et al., 1997). It is also possible that reduced evapotranspiration with the lower leaf area in the clearcuts not only increased available soil water but at the same time the greater incident light resulted in increased surface evaporation (Heithecker and Halpern, 2006). It is also possible that the legacy trees in the high retention treatments were accessing deep soil layers and redistributing the water to shallow layers through hydraulic lift (Brooks et al., 2006;Schoonmaker et al., 2007;Bingham and Simard, 2011). In a similar forest to ours, Walters et al. (2006) found that soil water increased in openings wider than one tree-length, or larger than 600 m 2 , but that these differences disappeared within 10 years due to the recovery of the plant community.

Carbon Pools
Forest floor C stocks (Mg ha −1 ) declined by 80-85% compared with pre-harvest levels one year after the clearcut, seed tree and 60% retention treatments but did not decline in the 30% retention treatment or controls, generally in support of our sixth hypothesis. Meta-analyses of harvesting effects on temperate forest floors around the world report an average 30% loss of forest floor C, regardless of whether clearcut or partial harvesting is employed (Zhou et al., 2013;James and Harrison, 2016), indicating feller-buncher and skidding practices at our site were particularly damaging. The loss of forest floor carbon represents not just a loss of C and nutrients from the site but also a substantial loss of mycorrhizal fungal inoculum. That the proportion of forest floor C lost in the 60% retention treatment was similar to the clearcut likely reflects the ground disturbance incurred during understory thinning. The retention treatment that best conserved belowground C stocks compared with pretreatment levels was the 30% patch retention. Shifting falling methods from feller-buncher to custom hand-falling would likely mitigate these losses.
As expected, the amount of C stored in aboveground live trees declined with increasing harvesting intensity, as did C stored in understory plants. Carbon stored in downed woody debris decreased in the three intermediate treatments relative to preharvest but the clearcut had a similar amount of C in woody debris as the control. It is likely that similar amounts of decayed wood were lost from the clearcut as the intermediate retention treatments, but C stocks of fine and small woody debris were elevated by logging slash. Mineral soil C stocks (to 15 cm depth) were unaffected by the treatments, agreeing with other studies (Johnson and Curtis, 2001;Powers et al., 2012;Noormets et al., 2015;Strukelj et al., 2015;Holub and Hatten, 2019).

Conclusions
Overstory tree retention increased natural regeneration of four ectomycorrhizal conifer species compared with clearcuts. Interior Douglas-fir regeneration was most dense within 10 meters of retained legacy trees, but was still abundant within 20 m, and occurred at lower yet adequate stocking within 60 m. Beyond 60 m from seed trees, there were few interior Douglasfir germinants. This distribution corresponds with patterns of abundance and diversity of ectomycorrhizal fruiting bodies and seedling root tip colonization found in previous studies in similar forests. The greatest abundance and richness of shrubs that serve as refuges of ectomycorrhizal inoculum occurred in the 30% patch retention and at higher retention densities, corresponding with conifer regeneration success. Richness of bryoids, herbs and shrubs were reduced by low retention treatments but diversity of these plant groups was unaffected. Overall, diversity and richness of regenerating tree species was not affected by retention treatment, but ingress of later seral species was favored by the treatments that retained more overstory cover, and pioneer broadleaf species were favored by treatments that created open conditions. Light availability was greatest in the low retention density treatments, remaining >1,000 W m −2 in the clearcut, seed tree and 30% patch retention treatments. Soil moisture did not significantly differ among treatments but tended to be highest in the clearcuts. Total ecosystem carbon declined with increasing harvesting intensity, and while intensive harvesting treatments were associated with large losses in forest floor carbon, the 60% retention treatment, where ground disturbance resulted from thinning, also experienced high forest floor losses. These losses would likely be mitigated by switching from feller-buncher to hand falling. We conclude that dispersed retention of overstory trees where seed trees are spaced ∼10-20 m apart, or aggregated retention where openings are <60 m, or two tree-lengths, in width, will result in an optimal balance of seed source proximity, ectomycorrhizal refuge plants, light availability, and forest floor C stocks for prolific and species-rich regeneration.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
SS and WR conceived the ideas, designed the methodology, and contributed equally to the writing of this paper. While sheltered at home during COVID, JBe, JBu, DC, DL, TS, AM-S, and AZ reviewed the literature contributing to the discussion. All authors conducted field work, contributed critically to the drafts and gave final approval for publication.