Increasing the intensity of regeneration treatments decreased beta diversity of temperate hardwood forest understory 20 years after disturbance

Abstract

Context

In temperate hardwood forests, increased intensity of soil and canopy disturbances tends to increase species richness due to the establishment of numerous early-successional plant species. However, while competitive pioneer species from early stages of succession can become recalcitrant and alter patterns of natural regeneration, very few studies have examined longer-term effects of these treatments on plant biodiversity.

Aims

In this study, we investigated mid-term (ca. 20 years) effects of different regeneration treatments with varying soil and canopy disturbance intensities. We compared understory plant communities in temperate hardwood forests from all the South of Quebec (Canada).

Methods

Using circular experimental plots of 1962.5 m² (radius = 25 m), we measured taxonomic and functional diversity indices and soil properties using four levels of disturbance intensity in six temperate hardwood forests of Quebec distributed along a longitudinal gradient. Reference forests, i.e. control forests with no silvicultural treatment known for ≥ 80 years, were compared to 20-year-old single-tree selection cuts, group-selection cuts and group-selection cuts with soil scarification.

Results

Species richness in both group-selection treatments was higher than that in reference forests. Plant equitability and beta diversity among sites in both group-selection treatments were lower than in single-tree selection cuts and control forests. More intense treatments contributed to the mid-term persistence of recalcitrant competitor species (e.g. Rubus idaeus L., Prunus pensylvanica L.f.) whereas soil scarification appears to have negative sustained effects on species known to be sensitive to regeneration treatments (e.g. Monotropa uniflora L., Dryopteris spinulosa Kuhn).

Conclusions

In temperate hardwood forests of Southern Quebec, silvicultural treatments of higher intensities resulted in detrimental effects on soil properties, especially in the surface horizon, 20 years after disturbance. This legacy, in turn, affected the composition and diversity of understory plant communities. The more intense silvicultural treatments contributed to the persistence of pioneer species better adapted to a wider range of environmental conditions and resulted in a decrease in understory plant community heterogeneity among sites. Conversely, single-tree selection cutting appeared to be the most appropriate silvicultural treatment for maintaining soil functions and heterogeneity of understory plant communities after 20 years; composition and structure being similar to long-undisturbed forests.

1 Introduction

Within managed temperate forest ecosystems, a primary focus has long been the development of methods that enhance post-harvest tree regeneration and stand productivity (Webster and Lorimer 2005; Walmsley et al. 2009; Bilodeau-Gauthier et al. 2020) assuming that forests were an unlimited source of raw material supply for wood products (Puettmann et al. 2009). Large canopy openings resulting from high disturbance intensity can substantially decrease the structural complexity of forest stands (Chaudhary et al. 2016; Roy et al. 2021), alter wildlife micro-habitats (Work et al. 2010) and change the water balance and thermal regime of the sites (Siemion et al. 2011), modifying patterns of ecological succession (Ellum 2009; Duguid and Ashton 2013).

The assemblage of specific plant communities varies according to climate and biophysical patterns, which are generally defined as ecological units or regions (Saucier et al. 2010). Forest management practices may reduce plant diversity and result in local extinction of the more sensitive species to canopy disturbances or changes in environmental conditions such as late-successional species with large seeds, low seed production and limited dispersal (Paillet et al. 2010; Kermavnar et al. 2019). At the same time, by changing the abiotic conditions of the forest floor, such as surface temperature, repeated disturbance such as tree harvest and site preparation can provide opportunities for competitive pioneer species to establish and spread over the short term (Royo and Carson 2006). These species often characterized by prolific vegetative reproduction and a persistent seed bank can in turn interfere with the regeneration of commercial species (Shields and Webster 2007) and ultimately create shifts in plant community composition (Paillet et al. 2010; Duguid and Ashton 2013). Yet, understory plant communities are essential to maintain key ecological functions such as forest biodiversity, energy flow and nutrient cycling (Gilliam 2007).

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 treatment and disturbance intensity, 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 et al. 2004; Wyatt and Silman 2010). Also, while emulation of strong soil disturbances is 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 properties. For example, forest canopy and soil disturbances can significantly reduce C concentrations up to 25% in the first 5 years after high-intensity harvest treatments such as clearcut (Neurath et al. 2010). Strong soil disturbance can also increase soil compaction and decrease the organic horizon thickness due to lower fresh litter inputs (Siemion et al. 2011; Chaudhary et al. 2016). Finally, forest harvesting followed by soil preparation treatments can lead to greater losses in nutrient reserves (N, P, base cations) located in the upper part of the forest floor (Thiffault et al. 2011).

In recent decades, management of North American hardwood forests has evolved towards alternative regeneration processes in order to maintain forest resilience to global change and to better integrate biodiversity conservation objectives (Messier et al. 2019). For instance, harvest practices that promote permanent forest cover and uneven-aged structure can have negligible negative impacts on both abiotic conditions and vegetation diversity (Chaudhary et al. 2016; Raymond et al. 2018), allowing the natural dynamics of regeneration to be maintained (Rogers et al. 2018). These forest practices that aim to simulate the natural disturbance regime (i.e. small openings initiated by old individuals death or light windfall) also promote the presence and maintenance of native flora of closed natural environments (e.g. large-seeded spring geophytes, with poor dispersal; Aubin et al. 2007).

Previous studies that have assessed floristic diversity within different silvicultural treatments have mainly focused on the early response (1 to 5 years after harvest) of vegetation to different levels of disturbance intensity (Roberts and Zhu 2002; Swanson et al. 2011; Bell et al. 2014; Tinya et al. 2019). Only a few studies have investigated this phenomenon in the medium term (~ 20 years after cutting, i.e. Roy et al. 2021). While there is much evidence that understory vegetation benefits from less intense treatments such as single-tree selection cut compared with clearcutting (Lindenmayer et al. 2012), the relationships among understory vegetation diversity, soil properties and the intensity of canopy and soil disturbance remain unclear. Moreover, many studies refer to species richness only, without considering abundance (Paillet et al. 2010; Chaudhary et al. 2016). Thus, it is important to verify whether the intensification of silvicultural treatments contributes both to longer persistence of the pioneer competitive species and to the decrease in taxonomic and functional heterogeneity of understory plant communities (Gauthier et al. 2016) in comparison to mature forests that have not been managed for some decades (Paillet et al. 2010; Duguid and Ashton 2013).

In order to capture different dimensions and provide a more synthetic view of the response of understory plant communities to different intensity levels of soil and canopy disturbance, we selected a set of metrics covering taxonomic and functional diversity (Mayfield et al. 2010). By characterizing plant species with a set of biological attributes (Violle et al. 2007) that determine the organisms’ response to its abiotic and biotic environment (Wellstein et al. 2011), the functional approach provides more information on the processes of species responses to anthropogenic disturbances (Aubin et al. 2007; Lavorel et al. 2007). Moreover, the taxonomic approach can detect the presence of rare species, but has low analytical power to characterize biodiversity and the effects of disturbance on ecosystem functions (Daly et al. 2018; Willis and Martin 2020). Consequently, substantial variations in plant diversity remain very difficult to measure, but also difficult to predict along a disturbance intensity gradient on the basis of species richness alone (Paillet et al. 2010; Duguid and Ashton 2013; Nolet et al. 2018). 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 main objective of this study was to assess medium-term response of understory plant communities to different regeneration treatments applied 20 years ago representing four levels of soil and canopy disturbance intensity, in temperate hardwood forests of southern Quebec, Canada. We hypothesized that the short-term effects of the regeneration treatments on understory plant communities and soil properties persist in the medium term; the magnitude of these effects would be proportional to the intensity of both forest canopy and soil disturbance. We used the results from understory vegetation surveys and soil analysis to answer three important questions relative to biodiversity conservation in managed forests: (1) Have soil properties in forests disturbed 20 years ago recovered to the levels of the reference forests (i.e. control forests that have not been logged for at least 80 years)? (2) Is alpha diversity still greater in forests disturbed more intensively? (3) Does the composition of understory plant communities in 20-year-old managed forests differ significantly from that of reference forests?

2 Materials and methods

2.1 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). Five sites were located in the yellow birch-sugar maple bioclimatic domain (Saucier et al. 2009), with contributions of American beech (Fagus grandifolia Ehrh.), American basswood (Tilia americana L.) and ironwood (Ostrya virginiana [Mill.] K. Koch). Site six was located in the yellow birch-fir bioclimatic domain and is dominated by balsam fir and yellow birch, with contributions from trembling aspen (Populus tremuloides Michx.), 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 Classification Working Group 1998; 7^th Approximation: Inceptisols and Spodosols).

Fig. 1

Geographic location of experimental sites. Information on the bioclimatic domains and silvicultural treatments present in each experimental site. Diagram of the experimental units. The experimental plot has a surface area of 1962.5 m² comprising 52 floristic inventory points and 9 soil sampling points distributed along the four inventory transects. CON control, SIN single-tree selection cut, GRP group-selection cut, GRPS group-selection cut with scarification

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 silvicultural treatments representing a gradient of canopy and soil disturbance intensity (Jaeger et al. 2022). These are as follows: (1) single-tree selection cuts (SIN, n = 13), an uneven-aged low-intensity silvicultural treatment removing about 30% of average basal area, (2) group-selection cuts (GRP, n = 15), a moderate-intensity treatment alternative to clear-cuts and resulting in openings ranging in area from 1500 to 2500 m², (3) group-selection cuts with soil scarification (GRPS, n = 17), a moderate-intensity alternative treatment with the additional purpose of creating soil microsites for some commercial tree species (see Appendix Table 6 for more informations about site characteristics). Scarification of sites Lac Marcotte and Saint-Michel-des-Saints was carried out immediately after GRP using a harrow (Fig. 1). At sites Escuminac, Kipawa and Woburn, scarification was carried out the year following the cut using a mechanical shovel, thereby creating an average 400 pits (ca. 2 × 3 m) per hectare. We have assumed that the two scarification treatments that were considered in this study both generated substantially greater disturbance to the soil and herbaceous layer compared to the unscarified experimental plots. We therefore combined these two scarification approaches into a single treatment for subsequent statistical analyses. Each site consisted of 3 to 5 randomized complete blocks in which three to four regeneration treatments were compared.

2.2 Vegetation survey

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, using presence/absence survey for each inventory point (Fig. 1). We did not count spring ephemerals, such as Erythronium americanum Ker Gawl (Yellow trout-lily). All individuals were identified to species using the identification key for vascular plants of the Flore Laurentienne (Marie-Victorin et al. 2002). Given the difficulty of identifying taxa at the vegetative stage in both the Carex genus and the Poaceae family, they were respectively grouped under Carex spp. and Poaceae spp. Taxonomic nomenclature was validated according to the Database of Vascular Plants of Canada (VASCAN; https://data.canadensys.net/vascan).

Each experimental plot covered an area of 1962.5 m² and included 52 inventory points, each with a radius of 15 cm, and separated by 1.5 m (Fig. 1). These were systematically distributed along four 25 m transects following Aubin et al. (2007). 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. When all 52 inventory points were inventoried, we walked the experimental plot for another 10 minutes to count species present in the experimental plot but never encountered at any of the 52 inventory points. They 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 an experimental plot by dividing its occurrence value by 52.

2.3 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 in 1:2 soil to water (deionized) ratio 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) on a TruMac CNS analyser (LECO, St. Joseph, MI, USA). Phosphorus (P) concentrations were determined colorimetrically (as molybdate blue) on Mehlich III soil extracts (Zidia and Tran 2008) with a flow-injection analyser (Lachat Instruments, Milwaukee, WI, USA). Exchangeable base (Ca, Mg, K, Al) concentrations were measured on BaCl₂ (0.1 M) soil extracts by atomic absorption spectrophotometry (Varian 220 FS, Agilent Technologies, Palo Alto, CA, USA) (Ziadi and Tran 2007).

2.4 Environmental performance traits

Three environmental performance traits (i.e. biological type, shade tolerance, reproduction mode; sensu Violle et al. 2007) were considered in our study (Table 1). They were obtained from the TOPIC database (Aubin et al. 2020) and were selected for the analyses because of their links with competitive abilities (biological type) and potential for colonization following disturbance (shade tolerance, reproduction mode). In order to compare the functional diversity between the treatments, we first calculated the functional dispersion index (Fdis; Laliberté and Legendre 2010). We 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. We determined the community-weighted means (CWM) to compare functional redundancy between the treatments.

Table 1 Individual environmental performance traits that were included in the analysis

2.5 Biodiversity measures

Several 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 (H') combined species richness and equitability 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: ¹D = e^(H') was calculated because it measures “true diversity” instead of entropy or probability (Jost 2006). (4) The equitability index (E) was calculated according to the formula: E = ¹D/S (Tuomisto 2010). Beta diversity was measured to quantify the extent of heterogeneity among sites (B_{among_Site} ; see Appendix Fig. 5). To compare beta diversity among sites for each treatment, we considered B_{among_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; Royer-Tardif et al. 2018).

To assess the conservation potential of the different treatments, we also tested the responses of both putative disturbance-sensitive and potentially recalcitrant species. Among the species present in our dataset, eight were identified in the literature as putative disturbance-sensitive to forest management in temperate forests ecosystems: Dryopteris carthusiana (Vill.) H.P. Fuchs (= D. spinulosa [O.F. Muell.] O. Kuntze), Athyrium filix-femina (L.) Roth, Lycopodium spp., Oxalis acetosella ssp. montana (Raf.) Hultén ex Löve, Coptis trifolia Salisb., Monotropa uniflora L., Circaea alpina L. and Cypripedium acaule Aiton (Haeussler et al. 2002; Moola and Vasseur 2008; Flinn 2007). Six other species have been identified in the literature as potentially recalcitrant and competitive with commercial woody and late-successional species in eastern Canada (Jobidon 1995; Bell et al. 2011): Acer spicatum L., Corylus cornuta Marshall, Populus tremuloides Michaux, Prunus pensylvanica L.f., Pteridium aquilinum (L.) Kuhn and Rubus idaeus L.

2.6 Data analyses

To test whether the management techniques changed soil properties and to answer the first question, the effects of silvicultural treatments on each soil property were measured using random effects generalized linear mixed models (glmm). Sites and blocks that were nested within sites were considered as random effects and silvicultural treatments were treated as fixed effects. Permutational multivariate analysis of variance (PERMANOVA) was also used to assess the variance partitioning of soil properties and test their multivariate means among treatments (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, understory plant community composition and species distributions were also analysed using redundancy analysis (dbRDA) that employed Sorensen’s index for distance matrix creation (Peres-Neto et al. 2006; Legendre et al. 2009). For this purpose, species that were present in less than 10% of the plots were excluded, given that rare species provides only limited information regarding habitat preferences and factors influencing co-occurrence between species (Azeria et al. 2012). To test whether the silvicultural treatments modified the mean number of species per experimental plot and to answer the second question, we compared alpha diversity measures between silvicultural treatments using random effects generalized linear mixed models (glmm) that were equivalent to those used for soil properties. Sites and blocks that were nested within sites were considered as random effects and silvicultural treatments were treated as fixed effects. To test whether the management techniques altered plant composition and species assemblage and to answer the third question, beta diversity within treatments (B_{within_Treat}) was compared using the same method that was used for alpha diversity. Beta diversity among treatments (B_{among_Treat}) and beta diversity among sites (B_{among_Sites}) were compared with the betadisper function from the vegan package (Oksanen et al. 2017). We also tested the differences in relative occurrence between different silvicultural treatments for putative disturbance-sensitive species and potentially recalcitrant species using one-way ANOVA with permutations (n = 999) (Borcard et al. 2011; Anderson 2017). Post hoc multiple comparisons with t-tests were performed to separate the treatments. When a significant difference was observed, Bonferroni correction was applied to the P-value. Finally, 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 × Sites). When significant differences were detected, we performed multiple comparison tests; P-values were adjusted using the Holm-Bonferroni sequential method. A significant 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 significant 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 significant differences, we used the function TukeyHSD() to perform pairwise means comparisons of the different regeneration treatments. The random effects mixed models were performed with the function glmmTMB() in the glmmTMB package (Brooks et al. 2017). Tukey post hoc tests for generalized linear random effects models were performed with the functions TukeyHSD() or glht() from the multcomp package (Hothorn et al. 2008). Hellinger transformations were performed with the function decostand(); multivariate dispersion analyses were performed with betadisper(); and centroid positions were tested with adonis2(). All of these functions are from the vegan package (Oksanen et al. 2017). All statistical analyses were performed in R (R version 3.5.2; R Core Team 2017).

3 Results

3.1 Soil properties

Geographic site location explained a significant portion of the variation in soil properties (R² = 0.36, P < 0.001). Treatments alone had a small but significant effect on soil properties variation (R² = 0.05, P = 0.039). Treatments had significant effects on FH-horizon thickness, C/N ratio and K content in the organic horizon and P content in the mineral horizon (Table 2). FH-horizon thickness in CON was significantly higher (almost two-fold) than in the GRPS. C/N ratio of the organic horizon in the control stands was higher than in the GRPS. In the mineral horizon, P content of the GRPs was three times higher than in the CON. Redundancy analysis (dbRDA) showed that soil properties were significantly associated with the composition and distribution of understory vegetation (P < 0.001; Fig. 2).

Table 2 Soil properties of the organic and mineral horizons, as a function of regeneration treatment

Fig. 2

Tri-plot of redundancy analyses (RDA) showing existing relationships between soil properties, sites, treatments and understory plant communities. Panel a represents taxonomic groups and plots that are grouped by treatments based on soil properties contributing significantly to the model. Panel b represents the ordination of plots that are grouped by sites based on soil properties contributing significantly to the model. Sample treatments are indicated by symbols (CON: control; SIN: single-tree selection cut; GRP: group-selection cut; GRPS: group-selection cut with scarification), explanatory variables by arrows. Percentage variance explained along each axis corresponds to R²

Soil properties explained 36% of variation in the specific assemblage of understory plant communities (R²_adj = 0.359, P < 0.001). Soil properties also explained 30% of the variation in the assemblage of biological types, i.e. variation of each biological type in relation to the variation of the others (R²_adj = 0.304, P < 0.001). dbRDA 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.

3.2 Alpha diversity

A total of 149 species were identified in the floristic 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 (S) in GRPs and GRPSs was significantly 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 (¹D) did not differ among regeneration treatments.

Fig. 3

(a) Species richness (S), (b) true diversity (¹D), (c) equitability index (E), and (d) functional dispersion index (Fdis) of understory plant communities in southern Quebec, as a function of regeneration treatments. CON: controls; SIN: single-tree selection cut; GRP: group-selection cut; GRPS: group-selection cut with scarification. Means with different letters significantly 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 × IQR)

Increasing the intensity of silvicultural treatments significantly affected species richness (S) and relative occurrence (F) of environmental performance traits that were related to biological type, shade tolerance and reproductive mode (Table 3). 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 3).

Table 3 Functional diversity (Rao index), mean species richness (S) and mean relative occurrence (F) of environmental performance traits as a function of regeneration treatments

3.3 Beta diversity

Silvicultural treatments had a significant effect on beta diversity among sites (P = 0.018). B_{among_Site} in CON and SIN was higher than in GRPS (Fig. 4).

Fig. 4

Multivariate analysis of variance of beta diversity among sites as a function of treatments (B_{among_Site}). Beta diversity among sites as a function of treatments (B_{among_Site}) are plotted as mean Euclidean distances between centroids for the same treatments. CON: controls; SIN: single-tree selection cut; GRP: group-selection cut; GRPS: group-selection cut with scarification

3.4 Differences in composition

The PERMANOVA analyses revealed significant multivariate differences in the species assemblage and the biological type assemblage between silvicultural treatments and geographic site location (Table 4). Geographic site location explained greater proportion of understorey plant species assemblage and biological types than did the treatments or the interaction. Based on the multiple comparison of treatments (Tukey HSD), species assemblage (P = 0.048) and biological type assemblage (P = 0.024) in the control forests differed significantly from those in the GRPS.

Table 4 Permutational Multivariate Analysis of Variance (PERMANOVA) of species assemblage and biological type assemblage across treatments, sites and their interaction followed by multiple comparisons (Tukey HSD tests) of species assemblage and biological type assemblage among treatments

Our results showed a significant decrease in the relative occurrences of the sensitive species toothed wood fern (Dryopteris spinulosa; P = 0.0154) and ghost pipe (Monotropa uniflora; P = 0.0016), together with a marginally significant decrease for northern wood sorrel (Oxalis acetosella ssp. montana; P = 0.085) in GRPS compared to the other silvicultural treatments (Table 5). Of the six potentially recalcitrant species, intensified disturbances of the canopy and soil had a significant effect on the relative occurrence of pin cherry (Prunus pensylvanica) and raspberry (Rubus idaeus). These two species occurred more frequently in both GRP and GRPS treatments compared to SIN and CON (Table 5).

Table 5 Taxon-specific responses of putative disturbance-sensitive and potentially recalcitrant species to regeneration treatments

4 Discussion

4.1 Have soil properties in forests disturbed 20 years ago recovered to the levels of the reference forests?

First of all, we found site location to explain more of the variation in soil properties than did silvicultural treatments, which has been observed in other studies (e.g. Heuvelink and Webster 2001; Nave et al. 2010). Indeed, soil properties are largely influenced by natural inter-site variability in surface deposits, and biotic and abiotic conditions that are associated with regional and local variation in climate and microclimate, together with vegetation type and the phenological and phenotypic characteristics of the most abundant species (Walmsley et al. 2009; Jang et al. 2016). 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 scarification treatments (Hope 2007). However, in our higher-intensity treatment (GRPS), substantial decrease in FH-horizon thickness, C/N ratio and exchangeable K were still observed 20 years after treatment compared to control forests. This decrease can be attributed to the extraction of forest biomass and subsequently lower inputs of fresh litter (Saint-Laurent et al. 2000; Thiffault et al. 2011; Clarke et al. 2015), higher litter decomposition rates from early successional species (i.e. with lower C/N ratios) compared to late-stage species (Kazakou et al. 2009), and the burial of surface organic matter in the soil mineral horizon (Henneb et al. 2019). In addition, an increase in light availability and soil temperature in the short term following GRPS can enhance rates of decomposition and C-mineralization in the soil organic horizon (Nieminen 2004; Bekele et al. 2007; Diochon et al. 2009). Our results indicate that higher-intensity silvicultural treatments in temperate hardwood forests can leave a legacy effect on soils that persists in the medium term, especially in the surface horizon. Indeed, other studies suggest that it may take several decades for soil properties to recover from the effects of high-intensity disturbance (Prest et al. 2014; Zhou et al. 2015).

4.2 Is alpha diversity still greater in stands disturbed more intensively?

The level of disturbance intensity and resource availability and heterogeneity are key factors structuring plant communities (Decocq et al. 2004). Indeed, the increase in the intensity of soil and canopy disturbance would have led to an increase in both solar radiation in the forest floor and available water and mineral nutrients, especially for early successional fast-growing nontree species (e.g. grasses, bushes and shrubs; Keenan and Kimmins 1993; Wardle et al. 2008). While this increased resource availability could have promoted competitive exclusion among species, which generally results in a species diversity decline in the short term (Tilman 1984), the more intense treatments considered in our study favoured a greater species richness compared to both the CON and the SIN forests (Deconchat and Balent 2001; Hilmers et al. 2018). Furthermore, species richness could have been even greater if Carex and Poaceae genera had been identified to the species level, as many of the species belonging to these two groups are known to be favoured in open environments with increased light levels (Aikens et al. 2007). The greater species number associated with higher forest disturbance intensity may result from the creation of more micro-habitats and longer-lasting canopy openings, leading to greater niche partitioning and higher ressource availability in the medium term. These results accord with other studies that have measured negligible effects of single-tree selection cuts on both environmental conditions and understory plant diversity (Smith et al. 2008; Raymond et al. 2018).

The species equitability index in the GRPSs was lower than in the CON and SIN forests, indicating a dominance of certain species in the GRPSs. Thus, we observed higher relative occurrence values in GRP and GRPS for early successional fast-growing nontree species that benefit from an ability to initiate efficient differential physiological responses through a wide range of functional trait values (Violle et al. 2007; Mayfield et al. 2010) compared to control forests and SINs. Nontree species often quickly invade degraded habitats. However, their ecological role in forest succession is still unclear. While the dominance of pioneer species following the intensification of forest disturbances is likely to disappear over time (Bergeron et al. 2014; Yeboah et al. 2016), it is known that soil disturbances generated by scarification or harvest activities can lead to species dominance effects by species known as recalcitrant from the literature (e.g. Rubus spp, pin cherry (Prunus pensylvanica)) (Jobidon 1995; Royo and Carson 2006) and a time lag in the recovery and recolonization processes of species adapted to natural forests (Bergstedt et al. 2008). Indeed, severe soil disturbance may constrain dispersal and recruitment processes for the regeneration of forest interior species that are known to be sensitive to habitat change (e.g. bryophytes, lichens) (Bergstedt et al. 2008; Caners et al. 2013; Venier et al. 2015), thereby favouring the longer-term persistence of competitive and vegetatively reproducing pioneer species through interspecific competition and exclusion processes (Aschehoug et al. 2016). Moreover, presence of nontree species may buffer harsh abiotic conditions and facilitate tree recruitment, thereby setting the stage for the next step in forest succession (Duncan and Chapman 2003).

4.3 Does the composition of understory plant communities in 20-year-old managed forests differ significantly from that of reference forests?

As expected, our results showed that much of the variability in understory plant community composition is associated with variation in soil properties between sites. Indeed, the combination of edaphic and microclimatic conditions within a bioclimatic region that undergoes the natural disturbance regime would generate a region-specific assemblage of species involved in the maintenance of beta diversity between geographically distant sites at the landscape scale (Vellend 2010; Saucier et al. 2010). Although control forests and SINs in our study had fewer understory species on average, they were associated with higher beta diversity between sites compared to GRPS, indicating more heterogeneous plant communities across forests and sites. Several studies have attributed similar results to increased resource heterogeneity and understory structural diversity in unevenly managed and control forests compared to the most disturbed forests, which may favour the development of biotic homogenization processes (Falk et al. 2008; Brewer et al. 2012; Markgraf et al. 2020; Roy et al. 2021). However, other studies have shown limitations to uneven-aged management that initiates only small disturbances, including lack of early successional habitats and landscape-scale homogenization of stand structure and composition (Werner and Raffa 2000; Decocq et al. 2004; Angers et al. 2005; Roy et al. 2021).

Changes in species richness and mean relative occurrence of functional groups helped us identify the factors shaping communities. When considering all species or functional groups at once in our PERMANOVA, we observed only slightly significant differences between the different treatments. In fact, by measuring plant diversity 20 years after different types of regeneration treatments, communities had already time to naturally converge. This process occurring with increasing stand age is well-known (Paillet et al. 2010). Compared to the CON and the SIN forests, shade-tolerant species were less frequent in both the GRP and GRPS forests where they were replaced by shade-intolerant species. This result indicates that the vertical structure of the vegetation was clearly different between the more intensed treatment (i.e. GRP and GRPS) and the close-to-nature single-tree selection cuttings. This may explain why forest interior species such as Dryopteris spinulosa and Monotropa uniflora were more representative of SIN forests, compared to pioneer and light-demanding species associated with group-selection cuttings.

Another factor possibly contributing to the floristic composition differences that were observed between CON, SIN, GRP and GRPS floristic communities is the site preparation and the level of soil disturbance relative to the machinery used during the treatments. In temperate forests of southern Quebec, soil scarification significantly contributed to increase the number of species with vegetative reproduction mode and also those with vegetative and sexual reproduction mode. These results are consistent with Newmaster et al. (2007), who found shade-intolerant, pioneer species with prolific vegetative reproduction and persistent seed bank such as Prunus pensylvanica and Rubus idaeus to be more numerous under more intensive forest floor disturbances. Some of these species are strong competitors that may form a recalcitrant understory layer that limits tree recruitment, growth and survival. In this situation, and in the context of a further silvicultural treatment within recently disturbed plots, more intensive control of the vegetation may be required to ensure tree establishment and late-successional species conservation.

Heavy canopy and soil disturbances, although less frequent in the natural disturbance regime within temperate hardwood forests (Runkle 1985), may allow specific regeneration needs to be met and favour the presence of pioneer plant species important for biodiversity (Swanson et al. 2011). Finally, plant communities are sensitive to both intense silvicultural practices and ecosystem-based close-to-nature management that is characterized by a low-level intensity of canopy and soil disturbance (Schall et al. 2018; Roy et al. 2021). In both cases, biotic homogenization can occur at the landscape scale and steer ecological successional processes towards either predominantly juvenile or old growth communities (Angers et al. 2005; Kern et al. 2017; Schall et al. 2018). In this sense, silvicultural management needs to find a balance at the landscape scale that allows for both light and moderate disturbances while drawing on natural dynamics.

5 Conclusions

Our results suggest that, in temperate hardwood forests, there is a legacy of higher-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 a slight increase of species richness in the medium term, but also a decrease in understory plant community heterogeneity among sites. The most intense silvicultural treatments resulted in an increased relative occurrence of shade-intolerant species, mainly vegetatively reproducing species adaptated to disturbed environments. These treatments contributed to the persistence of pioneer species, which could be competing with later-successional species and commercial species, especially if cutting intervals are shorter than recovery time, and eventually initiate considerable changes in plant community composition. Among the studied treatments, single-tree selection cutting appeared to be the most appropriate silvicultural treatment for maintaining soil functions and understory plant communities of long-term unmanaged forests. On the other side, more intense treatments favoured species better adapted to a wider range of environmental conditions, including open environments. Therefore, at the landscape level, the guarantee of maintaining the highest biodiversity may lie in the maintenance of heterogeneity in plant communities and soil properties. This can be achieved in the search for adequate spatial representativeness of the different silvicultural treatments by considering the functions of the species that they favour over the long term.

Appendix

Table 6 Environmental, ecological and edaphic characteristics of the study sites that were located in southern Quebec, Canada

Fig. 5

Different spatial scales involved in the beta diversity calculations. Beta diversity was calculated within each site for a given treatment (B_{within_Treat}), within each site among treatments (B_{among_Treat}) and among sites for each of the treatments (B_{among_Site}), adapted from Royer-Tardif et al. (2018). CON: control, SIN: single-tree selection cut, GRP: group-selection cut, GRPS: group-selection cut with soil scarification