Mid-term Impacts of Silvicultural Treatments on Soils and Understory Plant Diversity in Temperate Hardwood Forests of Quebec, Canada.


 Background: Short-term effects of silvicultural treatments on soil properties and understory vegetation in temperate hardwood forests are well documented, but few studies have examined longer term effects of treatment intensity. We hypothesized that short-term effects of silvicultural treatments on understory plant diversity do persist over the medium-term (20 years after treatment); the magnitude of these effects would be proportional to the intensity of canopy and soil disturbance.Methods: Soil properties (pH, total C and N, extractable P, exchangeable bases) and understory community diversity indices were measured in six experimental sites along a longitudinal gradient that covered different climatic and edaphic conditions in the yellow birch-sugar maple bioclimate domain. Reference condition, i.e. control forest with no anthropogenic disturbance for at least 80 years (CON) were compared to twenty years old regeneration treaments representing a gradient of canopy and soil disturbance intensity: single-tree selection cuts (SIN); group-selection cuts (GRP); and group-selection cuts with scarification (GRPS). Results: Geographic location of sites explained more variation in soil properties and community composition than did treatments. Species richness in both group-selection treatments was higher than that in CON forests. However, understory plant equitability and beta diversity among sites in GRP and GRPS were lower than in SIN and CON forests. More intense treatments (GRP and GRPS) increased the relative occurrence of vegetatively reproducing heliophilic plants, a trait syndrome associated with adaptations to disturbed environments. These treatments also contributed to the medium-term persistence of recalcitrant competitor species (e.g., Rubus idaeus, Prunus pensylvanica) whereas soil scarification appears to have negative sustained effects on species known to be sensitive to regeneration treatments (e.g. Monotropa uniflora, Dryopteris spinulosa).Conclusion: Of the treatments studied, single-tree selection cutting appears to be the most appropriate silvicultural treatment for maintaining soil functions and heterogeneous understory plant communities with compositions and structures similar to natural forests, while more intense treatments rather maintain and expand species that are better adapted to a wider range of environmental conditions, including open environments.


Background
Within temperate forest ecosystems, the diversity of understory plant communities is a key determinant of forest dynamics. These communities play important roles in the occurrence and maintenance of several key ecological processes (e.g., nutrient recycling, regeneration; Gracia et al. 2007; Moisan-De Serres et al. 2018). Yet, these communities exhibit increased environmental sensitivity that could be exacerbated by repeated harvesting pressures that disturb both the soil and forest cover (Roberts and Zhu 2002;Ellum 2009).
Since the 1990s, researchers have been interested in the relationships between silvicultural treatment intensi cation (i.e., increased disturbance intensity or frequency; Foley et al. 2005), forest productivity and biodiversity (Gilliam and Roberts 1995;Brockerhoff et al. 2008; Paquette and Messier 2010). Some have suggested that intensi cation required of forest management practices threatens forest landscapes (e.g., through decreases in structural complexity; Chaudhary et al. 2016) and biodiversity (e.g., changes in species composition, decreases in heterogeneity and functional diversity; Lindenmayer 2012, 2016; Gauthier et al. 2016;Yeboah et al. 2016;Messier et al. 2019). Others have suggested that by maintaining forest harvesting, more services in fact are produced (e.g., wood production) than could be achieved without management (Puettmann et al. 2009; O'Hara 2015).
In temperate hardwood forests, silvicultural treatments can have contrasting effects on biotic and abiotic conditions (Siemion et al. 2011;. Different levels of soil and canopy disturbance intensity, relative to the size of canopy openings and the nature of machinery that is used, can alter plant community composition and soil properties (Roberts 2004). For example, species richness and relative abundance of shade-intolerant species increases with the percentage of canopy openings (Keenan and Kimmins 1993). In contrast, small canopy openings, such as those that are created by natural disturbances, i.e., senescence and death of an individual tree or small group of individuals, tend to maintain the native ora of closed natural environments (Angers et al. 2005;Poznanovic et al. 2013).
Large canopy disturbances also provide opportunities for competitive pioneer species to establish and spread over the short-term (Webster and Lorimer 2005), which can interfere with the regeneration of commercially important species (Royo and Carson 2006;Shields and Webster 2007;Powers and Nagel 2009) and locally decrease the occurrence of late-successional species (Paillet et al. 2010; Duguid and Ashton 2013). Also, while emulation of strong soil disturbances that are caused by windthrow in the case of the natural disturbance regime, or following scari cation, are implemented to create microsites that favour desired commercial species, e.g., Betula alleghaniensis Britton (yellow birch) (Erdmann 1990), the application of these techniques can also substantially modify soil physico-chemical properties, such as reductions in C and N concentrations, soil compaction and decreases in the organic horizon thickness (Siemion et al. 2011;Chaudhary et al. 2016).
Management of North American hardwood forests has evolved towards alternative regeneration processes in order to maintain forests that are resilient to global change and to address problems that are related to biodiversity conservation. These alternative forms of management include partial cuts, which involve lower levels of disturbance that are closer to the natural disturbance regime, i.e., the ecosystem approach (D'Amato et al. 2018). Maintaining a permanent cover and uneven-aged structure allows the natural dynamics of regeneration to be maintained through the accumulation of woody biomass and the development of a complex structure with various tree ages and sizes (Rogers et al. 2018). Forest management that aims to simulate the natural disturbance regime also promotes the presence and maintenance of taxonomic groups that are sensitive to abrupt changes in abiotic conditions (e.g., largeseeded spring geophytes, with poor dispersal; Aubin  To date, the relationships between understory plant community diversity, soil properties, and the intensity of canopy and soil disturbances remain unclear. In the short-term, hardwood forests that are disturbed by logging tend to have higher species richness than unmanaged forests, but exhibit a lower degree of heterogeneity (Falk et al. 2008;Markgraf et al. 2020). This response is particularly due to the establishment of competitive pioneer species that specialize in openings within communities that are composed primarily of species that are associated with late-successional forests (Archambault et al. 1998;Naaf and Wulf 2007;Moola and Vasseur 2008). During ecological succession, the coexistence of pioneer and late-successional species is generally considered to be transient, i.e., lasting less than 20 years (Moola and Vasseur 2008). However, for a given disturbance, the functional traits of competitive pioneer species and their spatio-temporal capacities to exploit available resources may allow them to remain in the community for longer term and, ultimately, slow down successional trajectories (Roxburgh and Martin 2020). Consequently, substantial variations in understory plant diversity not only remain very di cult to measure, but also di cult to predict along a disturbance intensity gradient on the basis of species richness alone, for example, using meta-analyses (Gilliam and Roberts 1995 Nolet et al. 2018). In contrast, the functional approach provides more information on the processes of species responses to biotic and abiotic conditions (Aubin et al. 2007;Lavorel et al. 2007). The study of diversity according to these components thus permits a more detailed understanding to be gained regarding ecological processes and understory dynamics within disturbed systems (Hooper et al. 2005;Cadotte et al. 2011).
The objective of this study was to assess medium-term effects of different regeneration processes on soil properties and on the speci c and functional composition of understory plant communities in temperate hardwood forests. We hypothesized that short-term effects of silvicultural treatments on soil properties and understory plant communities persist in the medium-term (i.e., 20 years after treatment); the magnitude of these effects would be proportional to the intensity of both forest canopy and soil disturbance. We measured soil chemical properties (pH, total C and N, extractable P, exchangeable bases), together with several indices and metrics of understory plant community diversity, in six experimental sites that had been established primarily in yellow birch-sugar maple (Acer saccharum Marsh.) and yellow birch-balsam r (Abies balsamea (L.) Mill.) stands in southern Quebec, Canada. We compared a gradient of disturbance intensity treatments: 1) single-tree selection cutting; 2) group-selection cuts; and 3) group-selection cuts with soil scari cation, to reference condition, i.e. control forests with no silvicultural treatment known for ≥ 80 years.

Study sites and experimental layouts
Six experimental sites were selected within the Great Lakes-St. Lawrence Forest Region (Rowe 1972) of southern Quebec, Canada, along a longitudinal gradient (Fig. 1) that covered different climatic and edaphic conditions (Table 1). Five sites are located in the yellow birch-sugar maple bioclimatic domain (Saucier et al. 2009), where the latter two species dominated the canopy, with contributions of American beech (Fagus grandifolia Ehrh.), American basswood (Tilia americana L.), and hophornbeam or ironwood (Ostrya virginiana [Mill.] K.Koch). Site six is located in the yellow birch-r bioclimatic domain. Consequently, it is dominated by balsam r and yellow birch, with contributions from trembling aspen (Populus tremuloides Michx.), paper or white birch (Betula papyrifera Marsh.), and white spruce (Picea glauca [Moench] Voss). The creation of small canopy gaps following tree senescence and death is characteristic of the domain's natural disturbance regime, with occasional larger scale disturbances such as windthrows and freezing rain (Runkle 1985). Soils that have developed on the sites are Brunisols and Podzols (Soil Classi cation Working Group 1998; 7th Approximation: Inceptisols and Spodosols).
Legend: The experimental have a surface area of 1962.5 m 2 comprising 52 oristic inventory points and 9 soil sampling points distributed along the four inventory transects. CON: controls; SIN: single-tree selection cuts; GRP: group-selection cuts; GRPS: group-selection cuts with scari cation).
The sites were selected with the aim of comparing Control Forest that was not logged for at least 80 years (CON, n = 21) to three levels of disturbance intensity stemming from the following regeneration treatments. 1) Single-tree Selection Cuts exempli ed shelterwood cutting (SIN, n = 13), with the goal of removing about 30% of average basal area. 2) Group-Selection Cuts (GRP, n = 15) resulted in openings ranging in area from 1500 m 2 to 2500 m 2 . 3) Group-Selection Cuts with Scari cation (GRPS, n = 17) created areas that are equivalent to those of GRPs, with the addition of soil disturbance (Table 1).
Scari cation of the Lac Marcotte and Saint-Michel-des-Saints sites was carried out immediately after GRP using a harrow. At Escuminac, Kipawa and Woburn, scari cation was carried out the year following the cut using a mechanical shovel, thereby creating an average 400 pits (ca. 2 x 3 m) per hectare. We have assumed that the two scari cation treatments that were considered in this study both generated substantially greater disturbance to the soil and herbaceous layer compared to the unscari ed portions.
We therefore combined these two scari cation approaches into a single treatment for subsequent statistical analyses. Each site consisted of 3-6 randomized complete blocks in which three to four regeneration treatments were compared (Table 1).  2 The value in parentheses is the standard deviation.
3 Treatments: CON is the control, SIN corresponds to single-tree selection cut, GRP corresponds to group-selection cut, and GRPS corresponds to group-selection cut with scari cation.

Understory vegetation inventories
From the end of June to mid-August 2019, herbs, ferns and woody plants up to 2 m in height (i.e., the understory vegetation) were inventoried in the 66 experimental plots of the study. We did not count spring ephemerals, such as Erythronium americanum (trout lily). All individuals were identi ed to species using the identi cation key for vascular plants of the Flore Laurentienne ( . We assigned an occurrence value of 1 for each species that was present at an inventory point, with a maximum value of 52 for that same species within an experimental plot. Species present in the plot, yet never encountered at any of the 52 inventory points were scored 0.5 to account for the total species richness of the plot. The total number of recorded occurrences provides an estimate of species abundance. We considered the relative occurrence (F, as %) of a species within a plot by dividing its occurrence value by 52.

Environmental performance traits
The environmental performance traits (sensu Violle et al 2007) that were used in the study had been obtained from the TOPIC database (Aubin et al. 2020). These were selected for the analyses because of their links with competitive abilities (biological type) and potential for colonization following disturbance (reproductive mode, shade tolerance; Table 2).

Soil sampling and laboratory analyses
We collected a composite sample of the organic (FH) and mineral (0-20 cm) soil horizons in each experimental plot. Cores (8 cm dia.) were taken at nine sampling points and bulked for each of the plot (Fig. 1). Organic horizon (FH) thickness (cm) was also measured at the nine sampling points. Samples were air-dried and sieved to pass a 2 mm mesh. Bulk pH was measured using a soil:distilled water ratio of 1:2 for mineral soil and 1:10 for organic soil (Hendershot et al. 2008). Total nitrogen (N) and carbon (C) concentrations were measured by high-temperature combustion (1450 °C) followed respectively by infrared detection (C) and by thermal conductivity (N

Data analyses
Soil properties response to treatment The effects of silvicultural treatments on each soil property were measured using random-effects generalised linear mixed models (glmm), where sites and blocks that were nested within sites were considered as random effects and silvicultural treatments were treated as xed effects.
PERMANOVA was also used to assess the variance partitioning of the set of soil properties and tests of their multivariate means among treatments were determined by permutational multivariate analysis of variance (PERMANOVA, Anderson 2017). PERMANOVA (with 999 permutations) tested the effects of silvicultural treatments, sites and their interaction, while constraining permutations within blocks to reproduce random effects.

Relationships between soil properties and plant communities
Relationships between soil properties and plant community composition and understory species distributions were analyzed using redundancy analysis (dbRDA) that employed Sorensen's index for distance matrix creation ( Several univariate measures of alpha diversity were calculated to account for the distributions and patterns of species within understory plant communities. 1) Species richness (S) was measured as the number of species that were present in each plot.
2) The Shannon index or information measure (H') combined species richness (S) and equitability (E) using the geometric mean of proportional abundances of i species (p i ) in the respective communities, which was calculated as: H' = -∑ S i p i ln p i (Shannon 1948).
3) The "effective number" of species (i.e., true diversity) was calculated from the Shannon entropy exponent formula: 1 D = e (H') (Jost 2006). 4) The equitability index (E) was calculated according to the formula: E = 1 D/S (Tuomisto 2010). Effects of silvicultural treatments on alpha diversity indices were analyzed using random-effects generalized linear mixed models (glmm) that were equivalent to those used for soil properties.

Functional response to treatment
In order to compare the functional diversity between different silvicultural treatments, we rst calculated the functional dispersion index (Fdis; Laliberté and Legendre 2010). We also calculated the functional diversity of each individual trait using the Rao index (Rao 1982), to compare the variation of species traits composition within the communities. Finally, we determined the community-weighted means (CWM) for three environmental performance traits, namely 1) biological type, 2) reproductive mode, and 3) shade tolerance, to compare functional redundancy between silvicultural treatments.

Beta diversity
Beta diversity was measured to quantify the extent of heterogeneity between plots within the same silvicultural treatment, between plots of different treatments, and between sites (Appendix A). Within a treatment, beta diversity (B within_Treat ) was calculated as the alpha diversity ( 1 D) that was measured at the To compare beta diversity between sites for each treatment, we considered B between_Site as the multivariate dispersion of vegetation composition by calculating the average distance between centroids of sites for the same treatment. Both distance calculations were based on the Euclidean distance matrix corresponding to distance to centroid (dcen) (Anderson et al. 2011).

Species composition patterns
We assessed compositional differences between regeneration treatments using PERMANOVA (based on 999 permutations; Anderson 2017) based on distance matrices that were calculated from Hellinger distances. PERMANOVA tested for differences in species assemblages among silvicultural treatments and among sites, together with their interaction (Treatments x Sites). When signi cant differences were detected for the interaction and main effects, we performed multiple comparison tests; P-values were adjusted using the Holm-Bonferroni sequential method. A signi cant result that is obtained by PERMANOVA may originate from mean differences in species assemblages between the treatments, may indicate differences in variation within treatments (i.e., heterogeneity in multivariate scatter within groups), or may be a combination of both. To provide the appropriate interpretation of signi cant results, we used the function PERMDISP to test the homogeneity of the multivariate spread. PERMDISP is a permutation-based multivariate extension of Levene's test of homogeneity of variance (Anderson 2017). When these tests detected signi cant differences, we used the function TukeyHSD() to perform pairwise means comparisons of the different regeneration treatments.
We used a measure of multivariate functional dispersion

Effects of silvicultural treatments on soil properties
According to the multivariate variation-partitioning test (PERMANOVA), site location explained a signi cant portion of the variation in soil properties (R 2 = 0.36, P < 0.001). Treatments alone had a small but signi cant effect on variation in these properties (R 2 = 0.05, P = 0.039) ( Table 3). Multiple comparison tests showed no signi cant differences in soil properties between each pair of silvicultural treatments (Table 4).  Univariate analyses of the effect of treatments on each of the soil properties revealed signi cant differences between treatments for FH horizon thickness, C/N ratio and K content in the organic horizon and P content in the mineral horizon (Table 5). FH-horizon thickness in CON was signi cantly higher (almost two-fold) than in the GRPS. C/N ratio of the organic horizon in the control stands were higher than in the GRPS. In the mineral horizon, P content of the GRPs was three times higher than in the CON.

Relationship between soil properties and plant communities
Redundancy analysis showed that soil properties were signi cantly associated with the composition and distribution of understorey plant communities (P < 0.001; Fig. 2). The set of soil properties that we measured explained 36% of variation in the speci c assemblage of understory plant communities (R 2 adj = 0.359, P < 0.001).
Soil properties alone explained 30% of the variation in the assemblage of biological types (R 2 adj = 0.304, P < 0.001). RDA ordination and Pearson correlation test show the positive relationships between the presence of conifers and increasing C/N ratio of the organic horizon (Pearson r = 0.60, P < 0.001), and decreasing soil pH in the organic (r = 0.46, P < 0.001) and mineral (r = 0.33, P < 0.001) horizons.
Legend: Sample treatments are indicated by symbols (CON: controls; SIN: single-tree selection cut; GRP: group-selection cut; GRPS: group-selection cut with scari cation), explanatory variables by arrows.
Percentage variance explained along each axis corresponds to R 2 .
Effects of silvicultural treatments on the alpha diversity of plant communities A total of 149 species were identi ed in the oristic surveys. We counted an average 30 species in the CON plots, 31 species in SIN plots, 34 species in GRP plots, and 34 species in the GRPS plots. Species richness in GRPs and GRPSs was higher than in CON plots (P = 0.005, Fig. 3). Equitability (E) in CON and SIN was higher than in GRPS (P < 0.001). True diversity ( 1 D) did not differ among regeneration treatments. Legend: CON: controls; SIN: single-tree selection cut; GRP: group-selection cut; GRPS: group-selection cut with scari cation. Means with different letters signi cantly differ following pairwise Tukey's tests (P < 0.05). Box-and-whisker plots in each panel display 25th and 75th percentiles (the inter-quartile range from the lower and upper edges of the box), the horizontal lines within boxes indicate the 50th percentiles (medians), and bullets within boxes indicate means; whiskers below and above boxes indicate 10th and 90th percentiles, respectively, beyond which dots indicate outliers (values > 1.5 x IQR).

Effects of silvicultural treatments on taxonomic and functional metrics
Increasing the intensity of silvicultural treatments signi cantly affected species richness (S) and relative occurrence (F) of functional categories that were related to biological type, shade tolerance and reproductive mode (Table 6). GRPS increased species richness and relative occurrence, i.e., functional redundancy, of shrubs compared to the other treatments (P < 0.05). The relative occurrence of grasses in GRPs and GRPSs was higher than in CON and SIN. Species richness of exclusively vegetatively reproducing and shade-intolerant species in GRPs and GRPSs was higher than in CONs and SINs. Richness of species with both modes of reproduction (i.e., sexual and vegetative reproduction) and the relative occurrence of shade-intolerant species in GRPs was higher than in SINs. We did not observe any differences among treatments in terms of overall functional diversity (Fdis) (Fig. 3) or functional diversity for each environmental performance trait (Rao index) ( Table 6).   (Table 7).
Of the six potentially recalcitrant species, intensi ed disturbances of the canopy and soil had a signi cant effect on the relative occurrence of pin cherry (Prunus pensylvanica) and raspberry (Rubus idaeus). These two species occurred more frequently in the GRP and GRPS treatments compared to SIN and CON (Table 7). Mean relative occurrence values (%) are presented and two-way permutation ANOVA p-values. Signi cant treatment effects are represented by a P-value in bold (P < 0.05). Within rows, means followed by the same letter do not differ signi cantly from one another at α = 0.05 according to pairwise Tukey's tests. CON: controls; SIN: single-tree selection cut; GRP: group-selection cut; GRPS: group-selection cut with scari cation.

Effects of silvicultural treatments on beta diversity of plant communities
Silvicultural treatments explained less variation in the understorey plant species assemblage and biological types than did the site location ( Table 3). The interaction between regeneration treatments and sites was also signi cant. Based on the two-way comparison of treatments, species assemblage (P = 0.048) and biological type assemblage (P = 0.024) in the control forests differed from those in the GRPS (Table 4). Silvicultural treatments did not signi cantly affect beta diversity between silvicultural treatments (P = 0.10). Nevertheless, beta diversity between silvicultural treatments in CON and SIN tended to be higher than in GRP and GRPS (Fig. 4). Silvicultural treatments had a signi cant effect on beta diversity between sites (P = 0.018). Beta diversity between sites in CON and SIN was higher than in GRPS (Fig. 4).
Legend: Beta diversity among sites as a function of treatments (B between_Sites ) are plotted as mean Euclidean distances between centroids for the same treatments and (b) beta diversity among regeneration treatments (B between_Treat ) are plotted as Euclidean distances to the centoid as a function of the regeneration treatments. CON: controls; SIN: single-tree selection cut; GRP: group-selection cut; GRPS: group-selection cut with scari cation.

Effects of silvicultural treatments on soil properties
First of all, the silvicultural treatments considered in our study had more limited longer term effects on soil properties than we had anticipated although other studies have shown medium-term recovery of soil chemical properties following forest harvesting and scari cation treatments (Hope 2007

Relationships between soil properties and understorey plant communities
Our results show that much of the variability in understory plant community composition is associated with variation in soil properties between sites. For example, we observed that sites with a more coniferous understory were associated with higher C/N ratios in the organic horizon and lower pH values in the organic and mineral horizons. Lower pH can limit activity of the soil micro ora and fauna, leading to lower rates of litter decomposition, thereby resulting in increased organic soil accumulation and C/N ratios (Sinsabaugh et al. 2008). The species equitability index in the GRPSs was lower than in the CON and SIN forests, indicating dominance of certain species in the GRPSs. Soil scari cation treatments are frequently given little consideration from a biodiversity conservation perspective

Conclusions
Our results suggest that, in temperate hardwood forests, there is a legacy of high-intensity silvicultural treatments on soil properties, particularly in the surface horizon. This legacy, in turn, can affect the composition and diversity of understory plant communities. The most intense silvicultural treatments (i.e., GRP and GRPS) not only resulted in slight increased species richness in the medium-term, but also a decrease in understory plant community heterogeneity. The most intense silvicultural treatments resulted in an increased relative occurrence of shade-intolerant, mainly vegetatively reproducing species  Multivariate analysis of variance of (a) beta diversity among sites as a function of treatments (Bbetween_Sites) and (b) beta diversity among regeneration treatments (Bbetween_Treat). Legend: Beta diversity among sites as a function of treatments (Bbetween_Sites) are plotted as mean Euclidean distances between centroids for the same treatments and (b) beta diversity among regeneration treatments (Bbetween_Treat) are plotted as Euclidean distances to the centoid as a function of the regeneration treatments. CON: controls; SIN: single-tree selection cut; GRP: group-selection cut; GRPS: group-selection cut with scari cation.