Study site and sampling design
We conducted our study in 2021 in a regenerating 9.6-ha forest stand at Yale-Myers Forest (41 57’ N, 41 57’W), a privately owned 3,213-ha mixed hardwood temperate forest in northeastern Connecticut, USA. Mean annual precipitation is 125 cm, and between 2000 and 2020, the average January temperature was − 3.7℃ and the average July temperature was 21.9℃ (NOAA 2020). Dominant vegetation includes dense patches of mountain laurel (Kalmia latifolia L., an ERM shrub) and regenerating black birch (Betula lenta L., an ECM tree) in the understory, alongside several mature oak (Quercus rubra L. and Q. alba L.) and maple (Acer rubrum L.) trees retained as a part of the management strategy.
The stand was harvested in 2016, five years before our sampling, using an irregular shelterwood system designed to promote oak regeneration. Shelterwood systems retain existing mature trees for shade and as a seed source for the next crop of trees while preventing plants that will dominate in full sun (e.g., Rubus spp.) from outcompeting tree regeneration. Irregular shelterwoods are distinguished from uniform shelterwoods by having an uneven distribution of multiple age classes that are retained in the long term (Raymond et al. 2009; Ashton and Kelty 2018). This retention of standing mature trees and snags (standing dead trees) can mitigate the negative effects of timber harvesting on wildlife (Duguid et al. 2016, Mossman et al. 2019, Hanle et al. 2020) and topsoil health (Carpenter et al. 2021). As part of the management prescription, mountain laurel—the dominant ERM plant species at the forest—is systematically crushed to promote tree regeneration. We therefore anticipate that observed differences in organic matter build up and nutrient cycling at our site arise through plant-soil interactions that have occurred following the harvest.
The study design consisted of 30 50-cm x 50-cm plots split evenly between 5 areas in the stand. Each area was greater than 100 m apart and included 6 independent plots, half of which were dominated by regenerating black birch and the other half by mountain laurel. The distribution of the plots was designed to capture both local and broader-scale spatial heterogeneity within the forest stand. Black birch and mountain laurel were the two most abundant woody species in the understory of the stand and therefore captured variation in dominant vegetation across our site. Each plot included at least one stem of the designated study species, and no plots had stems of both species. We situated each plot at least 6 m from overstory trees (> 10 cm diameter at breast height (DBH)) to reduce the local influence of mature trees, and all plots in an area were ~ 6 m apart.
Soil and vegetation sampling
We sampled soil and vegetation in June 2021. We measured volumetric soil moisture to a depth of 12 cm using a HS2 HydroSense Soil Moisture Probe (Campbell Scientific, Logan, Utah, USA) probe at 3 locations within each plot (to account for fine-scale spatial heterogeneity). We collected and pooled three 15-cm deep soil cores from within each plot using a 2-cm diameter soil corer. Multiple stems and leaves (at least two) were collected from study species present within each plot.
Laboratory analyses
Freshly collected soils used for measuring N and C process rates were homogenized on a plastic tray, passed through a 4-mm sieve, and stored at 4℃ in gallon-sized airtight Ziplock bags before analysis. Gravimetric moisture content was calculated for each sample through mass loss after oven-drying at 105°C to constant mass. Water holding capacity was also calculated by inundating samples of soil with DI-water and allowing the sample to drain for 2 h, then finding the mass of the soil before drying, desiccating, and reweighing. We then assessed a subsample of soil to determine rates of potential net N mineralization and potential net nitrification, comparing the initial values of extractable inorganic N to values of extractable N after a 24-d incubation (Hart et al. 1994; Robertson et al. 1999). For each extraction, we added 25 mL 2 M KCl to each soil subsample (6-g dry weight equivalent), shook vigorously, refrigerated for 12 h, and then filtered the subsamples. We colorimetrically analyzed NH4+ and NO3− concentrations using a flow analyzer (Astoria 2, Astoria-Pacific, Clackamas, Oregon, USA). Net potential nitrification was calculated as the difference in NO3− in the incubated and initial samples, and net potential N mineralization is the difference between the initial and final sum of NH4+ and NO3−. For all these assays, soils were incubated at 20°C under a moist atmosphere with a Parafilm cover to allow air diffusion with a minimum of moisture loss. Incubated samples received weekly moisture adjustments to 65% of their water holding capacity, which is regarded as the optimal moisture for microbial activity (Robertson et al. 1999).
For C cycle processes, we first measured cumulative C mineralization (labile C) over the course of a 24-day incubation period. At set intervals (days 1, 4, 10, 16, and 24), the CO2 production rate was estimated as the area under the curve (‘AUC’) function (Bradford et al.2008) in the DescTools package (Signorell, 2023). Active microbial biomass was estimated using the substrate-induced respiration (SIR) method (Anderson & Domsch 1978), modified following West and Sparling (1986) and further following Fierer et al. (2003) to use autolyzed yeast as substrate. For both C mineralization and microbial biomass, we measured CO2 concentrations of headspaces over time using an Infra-Red Gas Analyzer (Li-COR model Li-7000, Lincoln, Nebraska, USA).
We collected coarse roots and fine roots from our soil samples to determine %C and %N in these pools in addition to the stems and leaves that we collected during fieldwork. To assess the concentrations of C and N in soil pools, we obtained particulate organic matter (POM; >53 µm fraction) and soil mineral-associated organic matter (MAOM; <53 µm fraction) pools from an air-dried soil subsample using chemical dispersal and physical separation (Paul et al. 2001). Measurement of these two soil pools permitted us to text the extent to which organic matter was rebuilding under the ECM and ERM plants. We used a Costech ECS 4010 Elemental Analyzer with Conflo III interface to analyze the amount of C and N in samples from each pool.
Statistical analysis
We used hierarchical, linear mixed-effects models for all analyses to account for potential spatial autocorrelation associated with how we grouped our plots into five areas across the forest stand. Specifically, we assumed a common slope for the random effects and included ‘area’ to represent different intercepts. All models included mycorrhizal association (0 = EcM and 1 = ErM) as a binary predictor and we tested the necessity of including continuous environmental predictors such as soil temperature and moisture. Such predictors were only retained if they influenced both the response variable and the conditional effect size of mycorrhizal association. The only predictor that did so was volumetric water content (VWC), which was included as a main effect and in an interaction with mycorrhizal type. One plot had a particular high nitrification value, and we tested the sensitivity of our conclusions to inclusion versus exclusion of this datapoint; we retained it because the effects of mycorrhizal association were unaltered. Further, when response variables strongly contravened assumptions of normality, they were loge or square-root transformed. Residuals were then checked for fit and effects checked to ensure that our results remained consistent with the underlying, non-transformed data. To permit comparison of the effect sizes of the predictor terms in each model, predictor variables were standardized by centering (the binary predictor) and by dividing by two standard deviations (the continuous predictor), following Gelman (2008). All analyses were conducted using the lmer package in the statistical freeware R (version 3.6.3; Bates et al. 2017; R Core Development Team, 2018).
Note that we did not use model selection nor rely on statistical significance and multivariate analyses for our causal inference. Instead, we relied primarily on evaluating standardized coefficient estimates (i.e., conditional effect sizes) considering published evidence of ECM, ERM and disturbance effects on soil N cycling. Our approach follows calls in ecology (e.g., Addicott et al. 2022) to reduce over-reliance on indicators of statistical significance (such as p, AIC and r2 values) for causative science. The calls in ecology reflect a broader shift (e.g., Bradford et al. 2021) advocated by entities such as the American Statistical Association to improve causative statistical inference (Wasserstein et al. 2019). We only report p and r2 values in our tables, given their value in indicating the precision of effect size estimates. Further, we follow emerging best practices in data communication and interpretation by presenting the underlying data observations in all our figures (Weissgerber et al. 2015).