Data

Data Exploration

Overstory Structure and Composition

Boxplots of live and dead basal area, tree density, and canopy cover were examined across treatments and years to assess data variability, identify potential outliers, and evaluate the range of overstory responses prior to analysis (Figs. 5-9). This exploration revealed a few data entry errors, which were corrected before proceeding with formal analysis. The Salvage treatment was plotted separately due to its distinct scale, and overstory data for this treatment were only collected in 2008 and 2025, as clearcutting in 2009 removed the original overstory and the regenerating stand was not measured until 2025 (Figs. 6, 8).

Figures 5-6. Temporal dynamics of basal area (m2/ha) across treatment plots (2008–2025).

Figures 7-8. Temporal dynamics of tree density per hectare (tph) across treatment plots (2008–2025).

Figure 9. Temporal dynamics of canopy cover percentage across treatment plots (2008–2025).

Table 2. Total number of live trees recorded per species, treatment plots, and year across all plots.

Species-level live tree counts across treatment plots and years are summarized in Table 2. Due to the low abundance of non-pine tree species across all treatment plots, boxplots could not be produced to represent their density distributions meaningfully.

Downed Woody Material

Small and large downed woody material biomass were visualized across treatment plots and years to assess data variability and identify potential outliers prior to analysis (Figs. 10-11). Small downed woody material biomass showed a declining trend across all treatment plots by 2025. Large downed woody material showed high variability in the 100-killed treatment, particularly in 2025, consistent with the expected heterogeneity of snag fall dynamics across plots.

Figure 10. Temporal dynamics of small downed woody material (DWM) biomass across treatments (2008–2025).

Figure 11. Temporal dynamics of large downed woody material (DWM) biomass across treatments (2008–2025).

Understory Plant Community Vegetation

Figure 12. NMDS ordination of all understory plant community data (Bray-Curtis dissimilarity, k = 3, stress = 0.188) showing all 648 plot-year observations.

Ordination of all understory community specie frequencies using Non-metric Multidimensional Scaling (NMDS) based on Bray-Curtis dissimilarities (k = 3, stress = 0.188) was performed to assess data quality and identify potential outliers prior to analysis (Fig. 12). Points are broadly distributed across ordination space with no plots dramatically separated from the main cloud, indicating no obvious outliers or data errors in the understory community dataset. Several plots at the periphery of the ordination space were investigated but found to be consistently distinct across multiple years, likely reflecting local microhabitat variation rather than measurement errors. No plots were removed prior to further analysis.

To minimize the influence of rare taxa in the multivariate analysis, understory vegetation species present in fewer than 5% of the plots were filtered and excluded, following the methodology applied by Steinke (2018) in a previous analysis of this dataset. This process resulted in a group of 29 common species (Table 3).

Table 3. List of common understory species used for multivariate analysis.

Using the common species retained after filtering, a three-dimensional solution (k = 3) was selected for final ordination, as a two-dimensional solution yielded a higher stress value (0.2506). Three transformation approaches were compared to assess ordination quality: no transformation, Wisconsin double standardization, and square root followed by Wisconsin double standardization (the default transformation applied by metaMDS in the vegan package). The Hellinger transformation produced the lowest stress value (0.1807) and was therefore selected for final analysis (Table 4).

Tabla 4. NMDS transformations comparison.

Soil Nutrients

To prevent statistical bias and ensure data reliability, nutrients with supply rates that frequently fell below the minimum detection limit (absence in at least 50% of plots) were excluded from analysis, resulting in the removal of NO3-N, Cu, Pb, and Cd. The remaining 11 nutrients were retained for further analysis: Al, B, Ca, Fe, K, Mg, Mn, NH4N, P, S, and Zn.

To capture key successional stages following disturbance, nutrient data from four sampling years were selected for analysis: 2008 (pre-treatment baseline), 2010 (early red stage), 2016 (early gray stage), and 2025 (late gray stage). Since soil nutrient sampling in 2025 was limited to plot corners (n = 48 sampling points), only the 48 quadrats consistently sampled across all four years were retained to ensure spatial comparability across the time series, resulting in 192 plot-year observations across four treatments.

Boxplots of raw nutrient supply rates revealed several outliers across treatments and years (Figs. 13–16). These were retained as they represent natural variability in forest soil nutrient dynamics rather than measurement errors, consistent with the high spatial heterogeneity expected in boreal forest soils.

Soil nutrient distributions were examined for skewness prior to multivariate analysis, and individual transformations were applied to reduce distributional asymmetry (Table 5). Nutrients with skewness exceeding 1.5 were log-transformed using log1p, those with skewness between 1.0 and 1.5 received a square root transformation, and nutrients with skewness below 1.0 were retained without transformation. Following transformation, all nutrients exhibited skewness values between -1 and 1.

Table 5. Summary of individual transformations applied to soil nutrient supply rates prior to multivariate analysis.

Figure 17. Principal component analysis (PCA) of soil nutrient supply rates across treatment plots and sampling years (n = 48 plots consistently sampled). Points represent treatment centroids per year.

A principal component analysis of soil nutrient supply rates was performed to explore soil nutrient trajectories across treatment plots and years (Fig. 17). Although the PCA revealed temporal shifts in nutrient composition across all years, all four treatments followed broadly similar trajectories in PCA space, making it difficult to distinguish treatment-specific patterns in soil nutrient composition. Nutrient-community relationships were therefore examined through vector fitting onto the understory community ordination in the results section.

Data Limitations

A limitation of this study is that soil nutrient sampling in 2025 was restricted to 48 of the 108 plots due to budget constraints. Therefore, only plots sampled consistently across all years were retained for nutrient analyses, which reduces the spatial replication available compared to other response variables.

Another limitation is that tree regeneration data for the Salvage treatment plots were only collected in 2025 after the application of the treatment in 2009. The recovery trajectory between those two points remains unobserved, and the pattern shown in the figures represents a 16-year gap rather than a continuous record.

Report abuse