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-fir bioclimatic domain. Consequently, it is dominated by balsam fir 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 Classification Working Group 1998; 7th Approximation: Inceptisols and Spodosols).
Legend: The experimental have a surface area of 1962.5 m2 comprising 52 floristic 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 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 levels of disturbance intensity stemming from the following regeneration treatments. 1) Single-tree Selection Cuts exemplified 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 m2 to 2500 m2. 3) Group-Selection Cuts with Scarification (GRPS, n = 17) created areas that are equivalent to those of GRPs, with the addition of soil disturbance (Table 1). Scarification of the Lac Marcotte and Saint-Michel-des-Saints sites was carried out immediately after GRP using a harrow. At 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 x 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 portions. We therefore combined these two scarification 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).
Table 1
Environmental, ecological and edaphic characteristics of the study sites that were located in southern Quebec, Canada.
Municipality | ZEC Kipawa | Lac Gagnon (Papineau-Labelle Wildlife Reserve) | Saint-Michel-des-Saints | Lac Marcotte (Mastigouche Wildlife Reserve) | Woburn | Escuminac |
Latitude | 46°52'N | 46°06'N | 47°01'N | 46°77'N | 45°21'N | 48°09'N |
Longitude | 78°42'W | 75°08'W | 74°20'W | 73°11'W | 70°48'W | 66°31'W |
Mean annual temperature (°C) | 2.7 | 4.6 | 3.1 | 4.1 | 3.9 | 4.0 |
Annual precipitation (mm) | 956 | 1090 | 933 | 1070 | 1367 | 951 |
Ecological type 1 (sub-domain) | YBSM (west) | YBSM (west) | YBSM (west) | YBSM (east) | YBSM (east) | YBBF |
Sand (%) | 53 | 57 | 53 | 38 | 53 | 21 |
Clay (%) | 5 | 33 | 37 | 9 | 38 | 57 |
pH (2:1 soil:water) | 4.5 (0.2) 2 | 4.9 (0.3) | 4.5 (0.2) | 4.7 (0.2) | 4.7 (0.2) | 4.7 (0.2) |
Treatments 3 | CON, SIN, GRPS | CON, SIN, GRP | CON, GRP, GRPS | CON, GRP, GRPS | CON, SIN, GRP, GRPS | CON, SIN, GRPS |
Number of blocks | 3 | 6 | 6 | 3 | 3 | 3 |
Year of cutting | 2001 | 2006 | 2000 | 1998 | 1997 | 1999 |
Year of scarification | 2002 | -- | 2000 | 1998 | 1998 | 2000 |
1 Ecological type: YBSM, yellow birch-sugar maple domain, YBBF, yellow birch-fir domain. |
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 scarification. |
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 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 genus Carex genus and the family Poaceae, they were respectively grouped under Carex spp. and Poaceae spp. Taxonomic nomenclature was standardized according to the Database of Vascular Plants of Canada (VASCAN; https://data.canadensys.net/vascan).
Each experimental plot (Fig. 1) covered an area of 1962.5 m2 and included 52 inventory points, each with a radius of 15 cm, and separated by 1.5 m. 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. 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).
Table 2
Individual environmental performance traits that were included in the analysis.
Trait | Description of trait |
Biological type | Qualitative variable: evergreen, deciduous, shrub, bush, sporophyte, monocot, graminaceae, asterales, other herbaceous plant. |
Reproductive mode | Qualitative variable: mainly asexual, mainly sexual, asexual and sexual. |
Shade tolerance | Qualitative variable: shade-intolerant, intermediate tolerance, shade-tolerant |
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) on a TruMac CNS analyzer (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 analyzer (Lachat Instruments, Milwaukee, WI, USA). Exchangeable base (Ca, Mg, K, Al) concentrations were measured on BaCl2 (0.1 M) soil extracts by atomic absorption spectrophotometry (Varian 220 FS, Agilent Technologies, Palo Alto, CA, USA) (Tran and Simard 2008).
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 fixed 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 (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).
Taxonomic response to treatment
To assess the effects of the different treatments on the understory plant community, we first tested the response of potentially sensitive and potentially recalcitrant species. Among the species in present in our dataset, eight were identified in the literature as being potentially sensitive to soil and canopy disturbance 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 2004; Flinn 2007). Other species that are known to be sensitive, such as spring geophytes (e.g., Trillium spp., Streptopus lanceolatus (Aiton) Reveal [= S. roseus Michx.]; Aubin et al. 2007) and orchids (e.g., Goodyera spp.; Turcotte 2008), had too low relative occurrence to be considered in the analysis. 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. We tested the differences in relative occurrence between different silvicultural treatments for potentially susceptible 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.
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 (pi) in the respective communities, which was calculated as: H' = - ∑Si pi ln pi (Shannon 1948). 3) The "effective number" of species (i.e., true diversity) was calculated from the Shannon entropy exponent formula: 1D = e(H') (Jost 2006). 4) The equitability index (E) was calculated according to the formula: E = 1D/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 first 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 (Bwithin_Treat) was calculated as the alpha diversity (1D) that was measured at the plot level, divided by the average 1D, which was calculated for each of the three plots in the treatment (Anderson et al. 2011; Royer-Tardif et al. 2018). To compare beta diversity between regeneration treatments, we considered Bbetween_Treat to be the multivariate dispersion of vegetation composition (Anderson et al. 2011) by calculating the average distance to the centroid for all plots within a treatment. To compare beta diversity between sites for each treatment, we considered Bbetween_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 significant 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 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.
We used a measure of multivariate functional dispersion
(FDis; Lalibert é and Legendre 2010) to compare trait
diversity between baseline and future conditions using data
on all traits simultaneously. FDis calculates the average
distance of individual species to their group centroid in
multivariate trait space that has been defi ned by an appropri-
ate distance measure. Tighter clustering of species in multi-
variate trait space corresponds to lower functional diversity
and thus lower FDis values. Although a range of functional
diversity indices exist we chose FDis because, unlike other
techniques, it is largely unaff ected by species richness
(sample size).
We used a measure of multivariate functional dispersion
(FDis; Lalibert é and Legendre 2010) to compare trait
diversity between baseline and future conditions using data
on all traits simultaneously. FDis calculates the average
distance of individual species to their group centroid in
multivariate trait space that has been defi ned by an appropri-
ate distance measure. Tighter clustering of species in multi-
variate trait space corresponds to lower functional diversity
and thus lower FDis values. Although a range of functional
diversity indices exist we chose FDis because, unlike other
techniques, it is largely unaff ected by species richness
(sample size
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). Beta diversity partitioning was performed with the function beta.pair() from the betapart package (Baselga 2010, 2012). All statistical analyses were performed in R (R version 3.5.2; R Core Team 2017).