Dynamic subcanopy leaf traits drive resistance of net primary production across a disturbance severity gradient

1. Introduction

Forests store 2.4 ± 0.4 Pg carbon (C) annually in biomass and soils (Pan et al., 2011), yet the future of this large C sink is uncertain as global change continues to reshape disturbance regimes (Williams et al., 2016; Keenan and Williams, 2018; Wang et al., 2021). Across the world’s forests, disturbance agents including introduced insects, fungal pathogens, and extreme weather are becoming more frequent and spatially extensive (Lovett et al., 2006; Cohen et al., 2016; Williams et al., 2016). Unlike severe disturbances (i.e., killing a larger fraction of trees, sensu Frelich and Reich, 1998) such as stand-replacing wildfires or hurricanes, host-specific insects and fungal pathogens frequently result in spatially heterogeneous tree mortality and partial, rather than total, canopy loss (Hicke et al., 2012). The degree to which these “moderate severity” disturbances cause declines in ecosystem net primary production (NPP) is variable, and not always proportional to disturbance severity (Stuart-Haëntjens et al., 2015; Fahey et al., 2016). Indeed, some forests have shown high production resistance (i.e., capacity to resist change, sensu Hillebrand et al., 2017) to disturbance, with post-disturbance NPP remaining at or near pre-disturbance levels up to a threshold of disturbance severity, beyond which NPP declines precipitously (Flower and Gonzalez-Meler, 2015; Stuart-Haëntjens et al., 2015; Riutta et al., 2018). However, the mechanisms underpinning such high functional resistance to disturbance remain poorly characterized across extensive gradients of disturbance severity and, because of potentially non-linear relationships, are impossible to infer from studies of severe, stand-replacing disturbances such as clear-cut harvesting or fire (Amiro et al., 2010; Nave et al., 2011; Goetz et al., 2012; Hicke et al., 2012; Gough et al., 2013).

As upper canopy disturbance progresses, shifts in community-wide subcanopy leaf functional trait responses may be crucial to NPP resistance (Seidl et al., 2017; Herben et al., 2018). In this study, we define “leaf functional traits” as morphological (e.g., leaf mass per area), biochemical (e.g., leaf N content or proxies thereof), and physiological (e.g., photosynthetic rate, conductance) properties of leaves influencing whole-plant resource acquisition and production (Reich et al., 2003; Wright et al., 2004; Reich, 2014; Niklas et al., 2023). Disturbances targeting the upper canopy stratum redistribute limiting resources such as light, nitrogen, and water previously captured by upper canopy trees (Refsland and Fraterrigo, 2017; Taylor et al., 2017; Curtis and Gough, 2018), affecting an array of leaf-to-ecosystem scale C cycling processes across entire plant communities. At the leaf level, disturbance-driven within-canopy resource redistribution can precipitate shifts in multiple leaf traits correlated with ecosystem-scale production (Poorter et al., 2009; Nave et al., 2011; Prado Júnior et al., 2015; Stuart-Haëntjens et al., 2015). Community-wide, dynamic changes in average abundance-weighted (i.e., community weighted mean, CWM) leaf traits could support compensatory plant production (Conti and Díaz, 2013; Bu et al., 2019; Happonen et al., 2022) and, consequently, whole-stand NPP resistance to disturbance. Although leaf functional trait shifts following disturbance have been characterized (Mcintyre et al., 1999; Louault et al., 2005; Bernhardt-Römermann et al., 2011; Herben et al., 2018), few studies have examined to what extent traits associated with leaf-level C fixation scale to affect the resistance of ecosystem-level C cycling processes, and whether subcanopy trait dynamics can sufficiently compensate for lost forest canopy area across a wide range of disturbance severity (van der Sande et al., 2017; Haber et al., 2020; Wei et al., 2021).

The timing and magnitude of these leaf-to-ecosystem C cycling process changes may also depend on the source and severity of disturbance. For example, severe disturbances from clearcut harvesting and fire rapidly alter CWM leaf traits by immediately changing subcanopy light availability (Refsland and Fraterrigo, 2017; Huerta et al., 2021). In contrast, changes in available light span a broader within-canopy range and lag the initiation of insect disturbance because defoliators and wood-boring pests precipitate more spatially variable and gradual tree mortality (Fotis and Curtis, 2017; Atkins et al., 2020). Better constraint of leaf trait dynamics and the degree to which they support NPP resistance across diverse disturbance types and severities is critical to disturbance-related modeling, adaptive forest management, and fundamental understanding of the biotic features promoting ecosystem functional stability (Anderegg et al., 2015; Seidl et al., 2017; Domke et al., 2018).

Our overarching goal was to determine whether, when, and to what extent changes in subcanopy CWM leaf traits support subcanopy and total production across an experimental gradient of disturbance severity. We asked three questions: (Q1) How and when do four subcanopy CWM leaf traits known to relate to ecosystem-scale C cycling processes–leaf mass per area, the red edge normalized difference vegetation index, the light-saturated rate of photosynthesis, and stomatal conductance–respond to disturbance across an extensive gradient of severity?; (Q2) To what extent are these CWM responses driven by upper canopy structural change and associated changes in light environment?; and (Q3) Which subcanopy CWM leaf traits best predict subcanopy wood NPP following disturbance? We hypothesized that CWM leaf trait responses would be proportional to disturbance severity, with the greatest increases in all four CWM traits occurring at the upper end of the disturbance severity gradient (H1; Stuart-Haëntjens et al., 2015), and that increasing canopy area loss over one to several growing seasons would largely account for subcanopy trait responses (H2; Canham, 1988). Lastly, we hypothesized that increased CWM leaf traits linked to C fixation would predict post-disturbance subcanopy NPP responses (H3; Fahey et al., 2016). To test these hypotheses, we selected the four leaf traits listed in Q1 above, each having established relationships with disturbance (Eitel et al., 2011; Prado Júnior et al., 2015; Zhang et al., 2022) and NPP (Garnier et al., 2004; Wright et al., 2004, 2019).

2. Materials and methods

2.1. Experimental design and site description

This study took place within the Forest Resilience Threshold Experiment (FoRTE; Gough et al., 2020; Grigri et al., 2020), an ecosystem-scale disturbance manipulation at the University of Michigan Biological Station in northern lower peninsula Michigan (UMBS, 45.56°N, 84.67°W). UMBS is dominated by secondary growth mixed deciduous and coniferous forest with mean annual precipitation of 817 mm and temperature of 5.5°C (Gough et al., 2013). The roughly century-old regrown forest of UMBS and the Upper Great Lakes Region is broadly undergoing successional transition (Wolter and White, 2002; Curtis and Gough, 2018), with dominant tree species shifting from early successional aspen (Populus grandidentata and P. tremuloides) and birch (Betula papyrifera) to later successional species including northern red oak (Quercus rubra), red maple (Acer rubrum), American beech (Fagus grandifolia), and white pine (Pinus strobus). FoRTE spans four distinct landform subunits, or landscape ecosystem types, characterized by differing soils, microclimate, plant community composition, and productivity according to the site-specific classification system developed at UMBS (Table 1; Pearsall et al., 1995). The pre-disturbance within-site spatial variation in CWM leaf functional traits may interact with disturbance at different severities with variable outcomes for post-disturbance NPP. Taken together, the distinct yet regionally representative FoRTE landscape ecosystem types allow for robust depiction of community functional trait shifts, and their consequences for NPP, across representative plant functional types of the Upper Great Lakes ecoregion (Gough et al., 2020). Prior to disturbance, substantial differences in productivity and forest canopy structure across landscape ecosystem types at UMBS were well documented (Hardiman et al., 2011; Scheuermann et al., 2018).

Table 1. Landscape ecosystem characteristics of the four FoRTE experimental replicates (Figure 1) prior to disturbance initiation (in 2018).

Forest Resilience Threshold Experiment was designed to independently evaluate the overlapping impacts of successional change in forest composition–specifically, the total loss of mature aspen and birch stems–and moderately severe disturbances on the temperate forest C cycle. Simulation of these coincident regional events is achieved through experimental implementation of disturbance across a replicated gradient of severity and two contrasting treatment types, each of which targets a different stem size class. Targeted gross defoliation levels of 0 (control), 45, 65, and 85% leaf area index (LAI) loss span a gradient from moderate to severe disturbance. Additionally, contrasting “top-down” and “bottom-up” treatment types were implemented to sequentially eliminate the largest and smallest (>8 cm diameter at breast height, DBH) stems, respectively, until the desired defoliation level was reached. By targeting contrasting size classes, each treatment type effectively simulates different disturbance agents: top-down attacks the oldest, early-successional species cohort (i.e., aspen and birch) similar to a host-specific fungus or insect outbreak (Hardiman et al., 2013), while bottom-up simulates a more species-agnostic agent such as deer browse (Reed et al., 2022).

We applied factorial combinations of disturbance severity and type across four experimental replicates (Figure 1A), each of which was located on a distinct landscape ecosystem type according to the classification of Pearsall et al. (1995). Replicates consisted of four 0.5 ha circular plots varying by disturbance severity, each of which was also bisected and received the contrasting treatments (“top-down” or “bottom-up”) randomly applied to the eastern or western halves (Figure 1B).

Figure 1. Forest Resilience Threshold Experiment (FoRTE) site map including experimental replicate locations (A) and sampling design (B). Four replicates each containing four plots with designated whole-plot gross defoliation levels (0, 45, 65, or 85% disturbance severity) are distributed across four distinct landscape ecosystem types at UMBS according to a local classification system [map colors in panel (A)]. Each plot was halved to generate two subplots, each of which was issued a restrictively randomized assignment of “top-down” or “bottom-up” treatment type (B). All measurements of interest were located within the experimental unit of subplots.

Baseline data collection, including tree inventory, forest structure assessment, micrometeorological measurements, and above- and belowground carbon flux measurements, began in summer 2018 prior to experimental disturbance initiation within each of 32, 0.1 ha. experimental subplots (Figure 1B). Tree inventory entailed tagging, species identification, stem diameter measurement, and mapping of all stems ≥8 cm DBH using a laser-outfitted digital caliper and GPS transponder array (Haglöf Inc., Madison, MS, USA). Before leaf-out the following spring (May 2019), we stem girdled ∼3,700 trees across 8 ha of forest. Stem girdling involves cutting and removal of a ring of bark and phloem tissue from the tree stem at approximately 1 m above the ground, via chain saw scoring and pry bar removal of tissue. Intended to simulate the effects of phloem-disrupting disturbance agents, stem girdling disrupts the transport of fixed carbon through the phloem, resulting in carbon starvation to the roots and tree death in 1–3 years following girdling (Gough et al., 2013).

2.2. Community weighted mean leaf functional traits

We measured subcanopy leaf physiological, morphological, and optical characteristics of interest annually at peak leaf area index to assess temporal changes in community-level leaf functional traits across the disturbance severity gradient. To capture an extensive and representative sample of the FoRTE subcanopy, we selected leaves from ungirdled subcanopy trees in all 32 experimental subplots. Leaves without obvious signs of damage by pathogens or herbivores were chosen, each from unique woody stems ≥1 cm DBH located at approximately 1 m height above ground. We sampled three leaves within each of four established vegetation survey plots in each subplot (inset Figure 1B) for a targeted total of 12 leaves per subplot and 384 leaf samples annually across FoRTE. Where leaves at 1 m were not present, we selected the nearest stem outside the survey plot with a leaf at the targeted height. Stems were tagged and the terminal branch from which the leaf was sampled was flagged for repeated measurements in subsequent years. In a subset of cases (39 of 433 unique stems sampled over 4 years, or 9%), stems died and we replaced them with new trees the following year. Because we sampled leaves that were present at 1 m height in a systematic manner across all subplots, our leaf sampling protocol mirrors taxon-free “trait transect” methods previously shown to generate robust estimates of CWM leaf functional traits in plant communities (Gaucherand and Lavorel, 2007; Lavorel et al., 2008).

Leaf-level physiological variables, the light-saturated rate of photosynthesis (A_net) and stomatal conductance (g_sw), were measured using a LI-6400 XT portable photosynthesis system (LI-COR Biosciences Inc., Lincoln, NE, USA) with internal chamber settings of 25 °C block temperature, 2,000 μmol m^–2 s^–1 photosynthetically active radiation (PAR), 400 ppm CO₂, and a targeted vapor pressure deficit below a maximum of 2 kPa. We sampled leaves while still attached to the tree, then removed and retained them for further analysis.

Immediately following leaf removal, we measured leaf reflectance in the field for intact adaxial leaf surfaces of all broadleaf samples using a CI-710 miniature leaf spectrometer (CID Bio-Science Inc., Camas, WA, USA). Conifer needles were too narrow to be measured with this instrument, resulting in the capture of reflectance spectra for broadleaf deciduous leaves only. Using the built-in SpectraSnap software, we computed a spectral vegetation index known to closely correlate with leaf chlorophyll content (Gamon and Surfus, 1999; Evangelides and Nobajas, 2020): the red edge normalized difference vegetation index (reNDVI), a narrow-band counterpart of the well-known NDVI:

We chose to use reNDVI instead of NDVI to avoid data loss resulting from instrument malfunction in 2020, which caused a subset of reflectance spectra in that year to fully saturate in the vicinity of 800 nm and prevented calculation of NDVI for all samples.

Leaf morphology, summarized as leaf mass per area (LMA), was determined for every sampled leaf following collection in the field. We scanned leaves while still fresh using an LI-3100C area meter (LI-COR Biosciences Inc., Lincoln, NE, USA) on the appropriate resolution setting (1.0 mm² for broadleaf and 0.1 mm² for needleleaf samples). Subsequently, we dried leaves for 72 h at 60 °C before weighing them, and then computed LMA as the ratio of leaf area to dry mass (g).

The actual number of leaf samples retained for physiological analysis per year averaged n = 379, as some leaves failed to stabilize during photosynthetic measurements (stability here defined as < 5% variation in reported photosynthetic rate between 180 and 300 s of measurement) or yielded unrealistically high scaled A_net values due to inaccurate leaf scanning (only in the case of some needleleaf samples).

2.3. Canopy light interception

The fraction of absorbed photosynthetically active radiation (fPAR) was modeled from canopy gap fraction using a site-specific relationship between gap fraction and fPAR (Supplementary Figure 1). We obtained gap fraction data at maximum leaf expansion in each year from 2018 to 2021 (i.e., pre- to post-disturbance) by optical imaging of the canopy. Skyward facing hemispherical images were obtained at 1 m above the forest floor with a self-leveling mount and 180° fisheye lens affixed to a Sony Alpha 6000 DSLR camera (Sony Electronics, New York, NY, USA). We sampled at least 2 h before or 2 h after solar noon at five locations within each subplot, including subplot center and at 10 m from center in each of the four cardinal directions. Images were collected under variable sky conditions over 4 years of sampling, but we processed them in WinSCANOPY (Regent Instruments, Quebec, Canada) with consistent color classification and thresholding for determination of sky vs. vegetation. We developed a site-specific linear regression model (p < 0.0001, Adj. R² = 0.89; Supplementary Figure 1) for computation of fPAR from gap fraction using 2018 and 2019 gap fraction data (collected as just described) as well as fPAR data, which we collected with an AccuPAR LP-80 ceptometer (METER Group Inc., Pullman, WA, USA) using methods fully detailed in Atkins et al. (2018b). We report computed values of fPAR using this regression model in each year.

2.4. Vegetation area index

We evaluated how disturbance-related changes in canopy structure affected subcanopy leaf traits using terrestrial lidar-derived vegetation area index (VAI), a unitless measure of one-sided canopy surface area (leaves plus the woody components of the canopy) per unit ground area. VAI is correlated with canopy light absorption and productivity (Atkins et al., 2018b). VAI was estimated from annual terrestrial lidar acquisitions taken with a portable canopy lidar (PCL) system. PCL consists of a user-mounted, upward facing near-infrared (900 nm) pulsed laser operating at 2,000 Hz (Riegl USA, Inc., Orlando, FL, USA) which collects returns as the user walks transects through subplots [full methods described in Hardiman et al. (2011)]. We analyzed lidar point cloud data using the forestr R package (Atkins et al., 2018a) with subplot-level means calculated as the mean of two 40 m transects per subplot, per year. Due to accidental omission, four subplots were excluded from measurement in 2018 and two in 2019.

2.5. Aboveground wood net primary production

Aboveground wood net primary production (ANPP_w) was computed from repeated measurements of stem diameter for a subset of ungirdled subcanopy trees (1–8 cm DBH), detailed by Grigri et al. (2020) and Niedermaier et al. (2022). Stems of this size class were censused in one quarter of each subplot (0.025 ha) in 2019 to obtain species abundances and diameter distributions. Repeated DBH measurements of eight subcanopy stems per subplot were made using digital calipers during and after the growing season in each year between 2019–2021. Using species- and region-specific allometric equations (Cooper, 1981; Gough et al., 2008), wood biomass increment was inferred for each stem in each year, then multiplied by species-specific stem densities to obtain subplot-scale ANPP_w.

2.6. Statistical analysis

All statistical analysis was performed in the R environment for statistical computing (version 4.0.2, R Core Team, 2020) using the fortedata R package (Atkins et al., 2021). With both subplot CWM leaf functional traits and subplot ANPP_w, we tested relationships with disturbance severity, type, and time using split-split plot mixed effects ANOVA with experimental replicate as a blocking factor, similar in design to those used in first-year FoRTE analyses (Gough et al., 2020; Grigri et al., 2020). In these models, disturbance severity was the fully randomized whole-plot factor while treatment type (bottom-up or top-down) was the restrictively randomized split-plot factor. In addition to treatment, time (as year) could not be fully randomized within blocks and was treated as the split-split plot factor in the models. Pairwise post-hoc comparisons using Fisher’s LSD in R package agricolae (de Mendiburu and Yaseen, 2020) were tested at alpha = 0.05 where significant effects were found, with a priori expectations about the direction of change articulated in our hypotheses. While three of four CWM leaf functional traits satisfied the assumptions of normality and homogeneous variance required for ANOVA, leaf mass per area (LMA) followed a right-skewed distribution and was transformed via 1/x transformation before analysis. Additionally, subcanopy and total ANPP_w data were log transformed prior to analysis to satisfy normality assumptions.

We used linear mixed effects models to assess the relationships between fPAR, canopy VAI, CWM leaf functional traits (A_net, g_sw, LMA, and reNDVI), disturbance severity, and time. We selected linear mixed effects models for this analysis because of their applicability to hierarchically structured data, as well as their capabilities in handling missing outcome values and non-normally distributed variables. All models and summary tables were generated using R packages lme4 (Bates et al., 2015) and lmerTest (Kuznetsova et al., 2017) using restricted maximum likelihood (REML) criteria and the Satterthwaite method for t-tests (Satterthwaite, 1946). The best candidate models were selected using the second-order Akaike information criterion suitable for small sample sizes (AICc), computed with R package MuMIn (Burnham and Anderson, 1998). In these models, disturbance severity, year, and their interaction were treated as fixed effects with regression coefficients of interest, while replicate was included as a random effect. Because there is no assumption of a Gaussian distribution for outcome variables in linear mixed effects models, we did not transform LMA for this analysis. Post-hoc comparisons of regression slopes were performed using R package emmeans (Lenth, 2020).

3. Results

3.1. Canopy light interception, vegetation area index, and disturbance severity

We observed a significant and temporally lagged reduction in both fPAR and VAI with increasing disturbance severity. fPAR and VAI declined most in the highest (85% defoliation) disturbance severity treatment with increasing time, indicated by significant interactions between disturbance severity and year for fPAR in 2020 [t(110) = −2.80, p = 0.006; Figure 2B; Supplementary Table 2] and for both fPAR [t(110) = −2.29, p = 0.024] and VAI [t(139) = −2.89, p = 0.005; Figure 2A; Supplementary Table 1] in 2021. VAI additionally exhibited a statistically significant response at the 65% disturbance severity level in 2021 [t(139) = −2.47, p = 0.015]. Here and throughout our subsequent analyses, treatment type (i.e., top-down vs. bottom-up tree girdling) did not emerge as a significant predictor variable and was eliminated from final models.

Figure 2. Subplot mean fraction of absorbed photosynthetically active radiation [(A): fPAR, unitless] and canopy vegetative area index [(B): VAI, unitless] across experimental disturbance severities and years. 2018 was the year prior to disturbance, while 2021 was the third year following disturbance initiation. Asterisks denote significantly different mean values using Fisher’s LSD multiple comparisons test. For this and following figures, the experimental unit was the FoRTE subplot (32 in total). Full and final models including all statistical parameters are included in Supplementary Tables 1, 2.

3.2. Community-weighted mean leaf functional traits and disturbance severity

We observed pronounced increases over time in subcanopy CWM reNDVI and more moderate increases in CWM A_net and CWM g_sw with increasing disturbance severity (Figure 3). Significant effects of year × severity [F(9, 72) = 2.06, p = 0.044; F(9, 72) = 6.70, p < 0.0001, respectively] on CWM A_net and reNDVI demonstrate that the influence of disturbance severity on subcanopy leaf traits intensified over time (Figures 3A, D; Supplementary Tables 3, 6). In particular, subcanopy CWM reNDVI and A_net values were higher in all or most disturbance severity treatments relative to control during 2020 and 2021. Stomatal conductance exhibited a marginally significant rise with increasing disturbance severity over time [year main effect: F(3, 72) = 82.5, p < 0.0001; severity × year: F(9, 72) = 1.83, p = 0.077; Figure 3B; Supplementary Table 4]. In contrast, CWM LMA was not affected by disturbance severity [F(3, 72) = 1.33, p = 0.32; Figure 3C; Supplementary Table 5].

Figure 3. The distribution of subplot CWM values of subcanopy leaf physiological [A_net, (A); and g_sw, (B)], morphological [LMA, (C)] and optical [reNDVI, (D)] functional traits across the FoRTE manipulation, 2018–2021. Asterisks (*) denote significantly different (p < 0.05) mean values using Fisher’s LSD multiple comparisons test. All statistical parameters are included in Supplementary Tables 3–6.

3.3. Subcanopy leaf functional traits, VAI, and fPAR

The loss of upper canopy VAI with increasing disturbance severity and time shifted subcanopy leaves from shade- to sun-acclimated trait profiles. We found strengthening negative relationships between subcanopy CWM leaf A_net and both upper canopy VAI and fPAR in each year following disturbance, which differed significantly (p < 0.0001 for all contrasts) from the pre-disturbance weakly positive relationship in 2018 for VAI (p = 0.02) and lack of significant relationship for fPAR (Supplementary Figures 2A, E). Likewise, inverse relationships between subcanopy CWM g_sw and both canopy VAI and fPAR emerged in the third year following disturbance [VAI: t(149) = −2.31, p = 0.02, Supplementary Figure 2B; fPAR: t(118) = −2.42, p = 0.02, Supplementary Figure 2F]. LMA lacked a significant relationship with VAI in any year but exhibited a significant negative response to fPAR [t(118) = −2.35, p = 0.02, Supplementary Figures 2C, G]. In contrast, subcanopy CWM reNDVI related significantly to canopy VAI and fPAR in some, but not all, years, with varying directionality (Supplementary Figures 2D, H).

3.4. Leaf functional trait relationships with subcanopy aboveground wood net primary production

The release of subcanopy vegetation following disturbance enabled full resistance of total (i.e., subcanopy plus canopy) ANPP_w, even at the highest disturbance severity level, and was associated with an increase in subcanopy CWM reNDVI in advance of VAI change. Subcanopy ANPP_w was strongly stimulated by increasing disturbance severity in the second and third years following phloem-girdling [year: F(2, 48) = 28.9, p < 0.0001; year × disturbance: F(6, 48) = 5.81, p = 0.0001; Figure 4; Supplementary Table 7]. This increase in subcanopy ANPP_w over time at higher disturbance severities sustained total ANPP_w even at the 85% gross defoliation level. Increases in subcanopy ANPP_w occurred before disturbance-related declines in VAI (Figure 2A), suggesting that limiting resources other than light may have stimulated subcanopy growth. The subcanopy trait most strongly coupling ANPP_w with disturbance severity was reNDVI [F(4, 91) = 9.09, p < 0.0001, Supplementary Table 8] which increased and became more spatially variable after phloem-girdling. In a multivariate context, reNDVI, year, disturbance severity, landscape ecosystem type, and the interaction between reNDVI and disturbance severity explained nearly half of the variation in subcanopy production (p < 0.0001, R²_adj = 0.46, AICc = 283.6; Figure 5; Supplementary Tables 9, 10). Variance in CWM reNDVI was found to differ significantly across disturbance severities, but not years (Levene’s test, p = 0.02 and p = 0.25, respectively).

Figure 4. Aboveground net primary production of wood (ANPP_w) across the FoRTE disturbance manipulation for the subcanopy stratum (opaque fill) and the subcanopy plus canopy (i.e., total, transparent fill) for the years 2019–2021. Asterisks denote subcanopy NPP distributions significantly different from the control in each year as computed in post-hoc tests on split-split plot ANOVA analysis (p < 0.0001, Fisher’s LSD; Supplementary Table 7). No significant differences from control were found for total NPP in each year. Total NPP analysis and statistical comparisons are detailed in Niedermaier et al. (2022).

Figure 5. Subcanopy community weighted mean leaf-level red edge normalized difference vegetation index (reNDVI) and aboveground wood net primary production (ANPP_w) in FoRTE, 2019–2021. Black line illustrates a non-linear response with increasing CWM leaf reNDVI (delinearized from log-linear regression model; p < 0.001). For full statistical model, see Supplementary Table 10.

4. Discussion

In the first 3 years following disturbance, a dynamic cascade of leaf-to-subcanopy responses enabled full NPP resistance across a range of disturbance severities and an array of landscape ecosystem types. As canopy trees senesced following girdling-induced phloem disruption and gradual consumption of their stored non-structural carbohydrates, ungirdled subcanopy trees exhibited physiological stimulation that fully compensated for declining canopy NPP (Figure 4). Supporting our first and second hypotheses (H1 and H2), we found strong evidence of declining VAI and fPAR with rising disturbance severity, causing CWM subcanopy leaf functional traits to shift toward higher photosynthetic rates, stomatal conductance, and reflectance consistent with higher sun acclimation of subcanopy leaves. However, changes in canopy structure and light environment lagged disturbance initiation by over a year, consistent with the gradual senescence of girdled trees. While three of four examined CWM subcanopy leaf functional traits–A_net, g_sw, and reNDVI–increased over time with more severe disturbance, reNDVI was most strongly related to subcanopy ANPP_w stimulation following disturbance. This subcanopy trait-ecosystem functional linkage (H3) was conserved across a spectrum of disturbance severities and ecosystem types, and also suggests dynamic mechanisms underpinning post-disturbance functional stability.

Our findings demonstrate that community-wide changes in surviving (i.e., ungirdled) subcanopy trees’ leaf traits are coupled with NPP resistance in multiple forest ecosystem types, and this subcanopy response sufficiently compensated for even the highest level of disturbance severity. Even at the highest 85% targeted gross defoliation level, total ANPP_w remained stable as surviving-tree leaf traits rapidly responded to disturbance. As girdled trees in the canopy layer senesced gradually over several years following disturbance, putatively fueled by stores of non-structural carbohydrates (Nave et al., 2011), the ungirdled subcanopy stratum saw marked shifts in leaf functional traits as well as NPP. Prior work across forest biomes demonstrates that following upper canopy disturbance, subcanopy functional traits can respond rapidly to higher light availability and reduced competition for water and nutrients, thereby increasing photosynthetic capacity per leaf area (Ellsworth and Reich, 1992; Niinemets et al., 2015; Oguchi et al., 2017). Such community-wide shifts in functional traits may stabilize ecosystem C cycling processes (Refsland and Fraterrigo, 2017; Huerta et al., 2021) by accelerating subcanopy growth and compensating for declining upper canopy production (Stuart-Haëntjens et al., 2015; Castorani et al., 2021). Although a number of studies have examined community trait responses to severe disturbances such as fire (Pausas et al., 2004; Huerta et al., 2021; Pellegrini et al., 2021), the timing and magnitude of plant functional trait change across more extensive gradients of severity, as well as how these shifts may relate to production, remains unclear (Diaz et al., 2007; Biswas and Mallik, 2010; Herben et al., 2018). Recent work from FoRTE has shown that growth of surviving subcanopy and canopy trees was able to fully compensate for reductions in canopy area up to 3 years following disturbance initiation (Grigri et al., 2020; Niedermaier et al., 2022). This stimulated subcanopy growth (i.e., increased carbon uptake of this stratum) coupled with declines in soil respiration (i.e., carbon loss from the system) that were proportional to disturbance severity and sustained through time (Mathes et al., 2023) enabled the forest to maintain remarkable NPP stability at even the highest levels of disturbance.

We found that changes in subcanopy leaf traits and production preceded disturbance-driven changes in canopy structure by at least a year, suggesting that resources other than light may have stimulated initial responses to disturbance. Significant increases in CWM leaf A_net and reNDVI were measured in the second year following disturbance, prior to a substantial loss of canopy vegetation area in the third year. The timing of these responses implies that disturbance-mediated subsidies (i.e., nutrients, water, or both) likely enhanced leaf physiological capacity before extensive changes in the subcanopy light environment occurred. In forests undergoing similarly gradual mortality from insect outbreaks, high nutrient retention by understory and surviving canopy trees, and concomitant increases in foliar N, have been found in some (Griffin et al., 2011; Rhoades et al., 2017) but not all (Lovett et al., 2002; Davis et al., 2019) cases. Similarly, increases in soil water content and water availability for surviving trees have been noted in diverse forest types following both natural (Reed et al., 2014; Goeking and Tarboton, 2020) and experimental (He et al., 2013) phloem-disrupting disturbances.

Earlier experimental work from our own site found that moderate severity disturbance via stem girdling of early successional trees allowed for nutrient investment in expanded late successional canopy leaf area (Nave et al., 2011, 2014; Stuart-Haëntjens et al., 2015), thought to be a key mechanism underpinning post-disturbance C cycling stability. Our study substantially extends such prior work in temperate forests, illustrating a lagged change in subcanopy light environment. The significant decline in canopy VAI and fPAR in the second and third year following disturbance signals a transition to higher within-canopy light availability and a potential mechanistic switch from nutrient- to light-driven physiological change in subcanopy CWM functional traits. We additionally found that while disturbance influenced both canopy VAI and subcanopy ANPP_w, the two did not directly relate to one another, suggesting that rapid, sustained subcanopy ANPP_w stimulation following disturbance is likely driven by a temporally dynamic series of drivers rather than immediate changes in canopy structure.

Our results indicate that leaf traits derived from spectral data may aid in predicting forest production responses to disturbance. Though several related leaf-level functional traits examined in this study responded to rising disturbance severity and declining canopy closure, leaf-level reflectance index (reNDVI) emerged as the most useful predictor of subcanopy production response. Across the broad range of disturbance severity levels in FoRTE, reNDVI significantly responded by the second year after disturbance. Notably, the relationship between disturbance severity, reNDVI trait response, and ANPP_w was also consistent across our array of four landscape ecosystem types within FoRTE, signaling an ecologically-conserved mechanism within a compositionally diverse spectrum of regionally-representative forest ecosystem types. Like its parent index NDVI, reNDVI is strongly correlated with chlorophyll content (Gamon and Surfus, 1999) and thus can indicate changes in foliar N allocation in advance of trait shifts associated with higher light availability (Zhang et al., 2006; Yang et al., 2017). Numerous recent studies have highlighted the utility of remote sensing in unraveling disturbance-driven ecological change and forecasting forest function under changing disturbance regimes (Senf et al., 2017; Stovall et al., 2019; McDowell et al., 2020). If used in combination with appropriate mechanistic models and scaling techniques (Wieczynski et al., 2022), remotely sensed changes in spectral indices like reNDVI may provide functionally meaningful information as disturbance unfolds across a range of severities and landscape types. Particularly given the potential to rapidly and cheaply detect leaf-level optical properties using field to spaceborne instrumentation (Pettorelli et al., 2016; Shiklomanov et al., 2019), such remote sensing-forest functional links hold promise for improving ecological forecasting capabilities under global change.

Our study is subject to several important limitations. First, our experimental disturbance, although designed to simulate the effects of phloem-boring insects, does not exactly replicate such insect outbreaks either in its simultaneous initiation of phloem flux disruption to all stems or in its species-agnostic targeting of trees for girdling. Second, while reNDVI has been shown to tightly correlate with leaf N content in multiple forest types (de Oliveira et al., 2017; Yang et al., 2017), we were not able to directly measure foliar N for the large number of samples included in our study and thus cannot provide our own direct analysis of changes in N content mediated by disturbance. A third consideration is that leaf-level measurements may not always scale directly to ecosystem processes for a variety of ecological (Funk et al., 2017; Messier et al., 2017) and modeling (Gara et al., 2019; Berzaghi et al., 2020) reasons, underscoring the need for continued ecosystem process model development in tandem with more available data capturing pre- to post-disturbance trait responses in forests. Finally, we note that forests are highly dynamic systems, and annual measurements of processes occurring at variable timesteps from seconds to growing seasons necessarily fail to capture the nuance and potentially high temporal variability of system processes over time. Nonetheless, our work leverages a controlled and replicated manipulative field experiment to advance understanding of how disturbance-mediated trait shifts underpin forest functional stability.

5. Conclusion

We demonstrated that disturbance-driven shifts in community-wide subcanopy traits sufficiently stimulated subcanopy NPP to support complete resistance of ecosystem NPP even up to the highest levels of experimental disturbance. Lagged but significant shifts in CWM subcanopy leaf traits coincided with increases in ANPP_w across compositionally diverse, regionally representative landform types varying nearly twofold in pre-disturbance site productivity. This traits-production relationship held across a wide range of disturbance severity levels from 0 to 85% gross defoliation, suggesting a conserved ecological mechanism in forests experiencing partial to near total stem mortality. Our findings reinforce the critical role of material legacies, or the biotic and abiotic resources remaining after disturbance (Niedermaier et al., 2022). Although prior work attributed the stimulation of subcanopy production following disturbance to greater light availability (Muscolo et al., 2014; Stuart-Haëntjens et al., 2015; Kennard et al., 2020), our 4-year analysis demonstrates that multiple leaf functional traits responded to disturbance before deterioration of the upper canopy. Remotely sensed leaf reflectance (as the index reNDVI) emerged as both the earliest sentinel of subcanopy trait changes following disturbance as well as the best predictor of ecosystem-scale production response across the gradient of severity. The potential of remotely sensed traits such as reNDVI to be used as early indicators of forest functional change following disturbance is promising as spectral and hyperspectral data collection platforms become more widely available. Future work should prioritize more complete characterization of such temporally dynamic mechanisms underpinning temperate forest production resistance, particularly across realistic gradients of disturbance.

Data availability statement

All data and code supporting this project are freely available. Data collected in this and related FoRTE analyses can be accessed through the fortedata R package: https://essd.copernicus.org/articles/13/943/2021/ (Atkins et al., 2021). R scripts needed to reproduce all statistical models are available in the following GitHub repository: https://github.com/lisahaber/forte_subcanopy_leaf_functional_traits.

Author contributions

LH, CG, and BB-L conceived the ideas and designed the methodology. LH and JA collected and processed the data. LH analyzed the data. LH and CG led the writing of the manuscript. All authors contributed critically to the drafts and gave final approval for publication.

Funding

We acknowledge the National Science Foundation (DEB award 1655095) for funding this project, as well as the University of Michigan Biological Station for providing housing, equipment, and technical support throughout the 4 years of data collection involved.

Acknowledgments

We thank research technicians and former graduate students Alexandra Barry, Cameron Clay, Maxim Grigri, Laura Hickey, Kerstin Niedermaier, and Autym Shafer for their time and assistance with data collection in various field seasons from 2018–2021. We also thank Jason Tallant for his help with interpreting the site-specific landscape ecosystem classification of UMBS.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/ffgc.2023.1150209/full#supplementary-material

References