Early warning signals of changing resilience in the biogeochemistry and biology of a northern hardwood forest

Change over time in early warning signals of stream solute concentrations and fluxes as determined with Mann–Kendall tests (see supplementary information) varied among five small watersheds monitored for over 50 years at Hubbard Brook (figure 2). In watershed 6, a biogeochemical reference watershed, most early warning signal metrics for Ca concentration and flux significantly decreased over time (figure 2; table S2), implying increasing resilience (figure 1(C)). Watershed 1, which received a calcium silicate (CaSiO₃) addition in 1999, also showed decreases in early warning signals of Ca fluxes indicative of increasing resilience (figure 1(C)) that were statistically significant for skewness and kurtosis (figure 2; table S2). By contrast, experimentally harvested watersheds generally showed increasing trends in early warning signals for Ca concentrations and fluxes that suggest declining resilience (figures 1(A) and (B)). In watershed 2, which was deforested in 1965–1966 followed by herbicide application for the subsequent three summers, all four early warning signals for Ca concentration significantly increased over time, indicating decreasing resilience. Likewise, three of the four early warning signals for Ca concentration in watershed 4, which was strip-cut by thirds in 1970, 1972, and 1974, implied decreasing resilience (i.e. significant positive trends). Watershed 5, which was whole-tree harvested during the winters of 1983 and 1984, similarly exhibited significant increases in autocorrelation, skewness, and kurtosis of Ca concentration, each indicating decreasing resilience.

Zoom In Zoom Out Reset image size

For stream NO₃ concentrations and fluxes, there was an overall pattern of increasing trends in skewness and kurtosis (implying flickering or declining resilience) in all watersheds including the unmanipulated biogeochemical reference watershed 6 (figure 2; table S2). These trends in early warning signals of stream NO₃ alone may not indicate declining resilience. Long-term declines in stream N export punctuated by periodically elevated values can manifest as increased asymmetry (skewness) and longer tails (kurtosis) in the data distribution even if resilience is decreasing, unchanged, or increasing. By way of example, figure S5 illustrates how skewness and kurtosis can increase even as mean and variance decrease.

Ecosystem state variables representative of biological organisms—soil microbial biomass carbon (C) soil microbial biomass nitrogen (N), tree growth, lepidopteran larval biomass, and avian abundance—also showed trends in early warning signals that suggested changing resilience. Note that none of these variables were measured in experimentally manipulated sites, but rather were monitored in reference areas of the forest that experienced exogenous long-term disturbances and episodic stressors (figure S1). Trends for these variables were less coherent than those for stream solute concentrations and fluxes. For example, except for a significant decrease in autocorrelation, there were no trends in early warning signals for soil microbial biomass C, an indicator of overall soil nutrient cycling. By contrast, the standard deviation of soil microbial biomass N, an indicator of N cycling, significantly decreased and both skewness and kurtosis significantly increased (table S2).

Patterns of change over time for basal area increment of two of the dominant tree species at Hubbard Brook, sugar maple (Acer saccharum) and American beech, were also mixed, with basal area increment of American beech showing significant increases in three of four early warning signals (i.e. autocorrelation, skewness, and kurtosis, indicating declining resilience) and sugar maple also showing declining resilience but only for one measure, an increasing trend in standard deviation (figure 2, table S2). For lepidopteran larval biomass, the two early warning signals representative of flickering, skewness and kurtosis, significantly increased while neither standard deviation nor autocorrelation changed over time (table S2). For migratory bird abundance, kurtosis increased for American redstart (Setophaga ruticilla), black-throated blue warbler (S. caerulescens), and black-throated green warbler (S. virens) even as skewness showed no significant change. These same species showed increasing, decreasing, or no trends in standard deviation or autocorrelation, respectively. The three other birds species considered, least flycatcher (Empidonax minimus), ovenbird (Seiurus aurocapilla), and red-eyed vireo (Vireo olivaceus), also exhibited a mix of significantly increasing, decreasing, or no change over time in early warning signals (figure 2, table S2).

More than half of the early warning signals calculated (58%) indicated declining resilience. Another 27% indicated increasing resilience, while 15% showed no trend (figure 3). The significance and direction of these trends varied both within and among ecosystem attributes. Of the 31 ecosystem state variables used for early warning signals analysis, 37% exhibited temporal increases in skewness and kurtosis. Our conceptual model (figure 1(B)) suggests that these positive trends are indicative of declining resilience in the northern hardwood forest at Hubbard Brook and of flickering between alternative states. When differentiated into biogeochemical and biological state variables, increases in skewness or kurtosis likewise represented the plurality of changing resilience indicators. For stream NO₃ concentrations and fluxes, 45% of early warning signals calculated across five small watersheds showed flickering or flickering ambiguous behavior. Similarly, 41% of early warning signals derived from migratory bird abundance fit with the flickering category of our conceptual model.

Zoom In Zoom Out Reset image size

Increases in standard deviation or autocorrelation that represent critical slowing down comprised a small proportion of trends, amounting to 20% of the temporal behavior of early warning signals calculated for all 31 ecosystem state variables included in the analysis. Critical slowing down was also less symptomatic of resilience loss within specific biogeochemical and biological attributes. Only 20% of early warning signals calculated for stream Ca and NO₃ concentrations and fluxes showed temporal behavior representative of critical slowing down. In the cases of lepidopteran larval biomass and soil microbial biomass C and N, critical slowing down behavior was completely absent, and the lack of change over time in both standard deviation and autocorrelation was classified as no trend (figure 3). In multiple cases early warning signals exhibited decreasing trends, which implied 'recovery' in our conceptual model (figure 1) and in the summary of early warning signals trends (figure 3). This pattern was evident for migratory bird abundance, soil microbial biomass (C and N), and some stream Ca concentrations and fluxes. Thus, while slightly over half of all early warning signals analyzed showed temporal increases that would suggest declining resilience, slightly less than half exhibiting either unchanging or increasing resilience.

Measures of climate change and atmospheric deposition (figure S2) were significantly correlated with early warning signals in the long-term data at Hubbard Brook (figure 4, tables S3–S7). For both univariate correlations (tables S3–S6) and multivariate analyses (figure 4, table S7), early warning signals were more strongly correlated with precipitation chemistry (Total N, pH) than with climatic drivers (mean annual temperature, total annual precipitation). For stream solutes and fluxes aggregated across watersheds, changing precipitation pH and atmospheric N inputs mapped onto indicators of flickering (figures 4(a) and (b)). Skewness was more associated with mean annual precipitation pH while kurtosis (e.g. see figures 4(a) and (b)) was more correlated with total precipitation inorganic N flux. For multivariate analyses of migratory bird abundance and tree basal area increment, general patterns were not observed between changing atmospheric chemistry and early warning signals of potential regime shifts (figures 4(c) and (d)).

Zoom In Zoom Out Reset image size

We used long-term records of stream solute concentrations and fluxes at Hubbard Brook to demonstrate how resilience in ecosystem nutrient retention might vary between reference and experimentally disturbed watersheds. Declining trends in early warning signals of stream Ca concentrations and fluxes in watershed 6, the biogeochemical reference, and watershed 1, which received CaSiO₃ additions, align with prior understanding of forest recovery from the long-term legacy of acid rain while also providing new insights into the linkage between Ca and resilience. The recovery of watershed 1 following Ca additions has been well-documented [e.g. 45, 50] and fits with decreasing trends in early warning signals of stream Ca concentrations and fluxes that we report here. By contrast, declining early warning signals of stream Ca in watershed 6 were somewhat surprising given changing growth rates and species composition at this site [45] that suggest resilience loss. In fact, temporal decreases in early warning signals of stream Ca concentrations and fluxes from watershed 6 suggest that the system is not approaching a regime shift and may, in fact, be slowly recovering from decades of acid deposition. Early warning signals analyzed from experimentally harvested watersheds indicated critical slowing down and flickering (table S2), or evidence of declining resilience. These increases in early warning signals of Ca concentrations and fluxes may represent inherently low rates of nutrient cation weathering in tandem with acid deposition [54] from the three experimentally cut watersheds (watersheds 2, 4, and 5) that experienced additional local perturbations in Ca loss associated experimental forest cutting [38]. This combination of pulse (e.g. forest cutting) plus the legacy of press (e.g. acid deposition) disturbances may have rendered these watersheds more vulnerable than reference locations in the forest to the long-term impacts of acidic precipitation [44]. While previous research has also shown that forest regeneration dynamics in experimentally cut areas may be attributable to Ca losses from both acid rain and harvesting [58], these dynamics have not been explicitly considered from a resilience perspective.

Trends in early warning signals of stream NO₃ told a different story than Ca. Early warning signals for stream Ca concentrations differed between unharvested and harvested watersheds, but this did not occur for stream NO₃. All of the watersheds in the analysis except for watershed 1, which received CaSiO₃ additions, showed significant increases in both the skewness and kurtosis of stream NO₃ fluxes and concentrations that imply flickering within the system (figures 1 and 3). This coherence in early warning signals of stream NO₃ affirms previous research reporting that the site is transitioning into an alternative biogeochemical regime, from eutrophic, where N is abundant, toward oligotrophic, where N is relatively scarce [59]. The flickering behavior of early warnings signals of soil microbial biomass N (figure 3; table S2) suggests the ubiquity of this transition at Hubbard Brook, with implications for future forest productivity that are poorly understood [59]. Beyond confirming prior documentation of N oligotrophication at the site, the flickering behavior we report suggests loss system of ecosystem resilience over time. The earliest harvest experiment at Hubbard Brook, the clear cut and herbicide application at watershed 2, showed incredible resilience of the N cycle to harvesting and removal of plant biomass [37]. Increasing trends in skewness and kurtosis from harvested and biogeochemical reference watersheds shown in this study indicate declining resilience such that a repeat experiment at watershed 2 might result in very different outcomes with regards to N biogeochemistry.

Early warning signals of ecosystem attributes representing soil nutrient cycling, net primary productivity, primary consumers, and secondary consumers exhibited some symptoms of declining resilience such as critical slowing down and flickering. They also showed decreasing trends or lack of change over time, implying increasing or unchanging resilience, respectively. The positive trends in skewness and kurtosis of lepidopteran larval biomass (figure 3; table S2) imply declining resilience and flickering for caterpillars that match recent reports of decreasing insect populations at Hubbard Brook [60] and in many other areas [61]. This pattern is also consistent with previous studies at the site, which have described 20- to 40-fold interannual variation in lepidopteran larval biomass that may be driven by a combination of extreme winter cold and long, warm summers [62, 63]. Interannual variation in phytochemistry, particularly N, may also have influenced early warning signals of caterpillar resilience. We showed that total rates of annual NO₃ deposition were negatively correlated with both the skewness and kurtosis of lepidopteran larval biomass (tables S5 and S6). Investigations at Hubbard Brook and elsewhere have affirmed the linkages among N inputs, soil N availability, and lepidopteran abundance [64, 65]. At Hubbard Brook, observed decreases in soil N availability that arise from N oligotrophication [59, 66] may result in lower foliar N concentrations that lead to decreased caterpillar biomass [64]. Here, for the first time, we suggest declining resilience in lepidoptera in response to changing N biogeochemistry and other stressors at the site, as shown by positive trends in skewness and kurtosis.

Our results also indicate flickering in American beech basal area increment and critical slowing down for sugar maple basal area increment (figure 2; table S2). Both are foundational, long-lived, late successional species of northern hardwood forests, and sugar maple is a key component of the regional forest economy and cultural identity. Outside of the experimental watersheds, biomass accumulation for both of these species has plateaued since ∼1980 [67]. This slowing growth may be the natural result of forest maturation [68]. However, the relationships we observed between changing precipitation chemistry and early warning signals of basal area increment, either alone (tables S3–S6) or in multivariate analyses (figure 4; table S7), also fit with current understanding of how shifts in acid deposition, soil nutrient status, and soil nutrient depletion interact to influence the growth of these species [69]. For sugar maple, decades of acid deposition and subsequent soil Ca depletion may account for declining growth rates [50]. In the case of American beech, the arrival of beech bark disease at Hubbard Brook in the 1970s has led to lower growth rates and increased mortality in canopy trees [45]. When combined with current understanding of the interactive effects of sugar maple decline and beech bark disease, indicators of declining resilience suggest a possible transition in forest structure with implications for biogeochemical functions such as soil nutrient cycling [70] and possible impacts on biological diversity and the abundance of particular wildlife species [71]. Early warning signals of decreasing resilience of sugar maple and American beech also suggest that these species may not be able to recover quickly from future disturbances.

Although we observed flickering in the statistical behavior of variables as disparate as soil microbial biomass N, basal area increment of American beech, and lepidopteran larval biomass that may have implications for secondary consumers, changes in annual bird abundance at Hubbard Brook did not always show clearly interpretable patterns. Two of the six species we examined, American redstart and least flycatcher, were abundant at Hubbard Brook in the 1970s, but subsequently declined and have mostly disappeared from the area. In the analyses here, they showed significant, negative trends in standard deviation (figure 2, table S2), which, according to theory, could indicate a 'critical speeding up' that can precede local species extirpation [72], though this concept is not well-established. Declines in skewness for both species, in fact, point to a transition that has already occurred [32] due to a change in forest structure as the forest aged from mid-successional to older, more mature forest stands [73]. One bird species whose population has started declining in recent years, the black-throated blue warbler, showed positive trends in autocorrelation, skewness, and kurtosis (figure 2; table S2), which may be indicative of critical slowing down or flickering (figure 1). However, the black-throated green warbler showed similar trends in early warning signals as the black-throated blue warbler (table S2), even as its abundance has been increasing over time at Hubbard Brook. Interannual variation in avian populations is a function of complex interactions among fecundity in the previous season, adult density, nest depredation, habitat quality, insect food resources, weather conditions in overwintering grounds, and survival during migration [62, 74, 75]. These complex interactions may obscure the direct impacts of acid deposition and climate change on bird abundance [18]. Contributing factors outside of Hubbard Brook may also render birds and other migratory wildlife poor candidates for analysis of changing resilience within the site. These findings highlight the importance of detailed, ecosystem-specific knowledge in choosing variables and interpreting trends in early warning signals to avoid drawing inappropriate conclusions [21, 29].

This study used early warning signals to examine whether and how resilience has changed at Hubbard Brook, and what this change might portend for the biogeochemistry and biology of the future northern hardwood forest. Temporal trends in critical slowing down and flickering were neither coherent across ecosystem state variables, nor were they unambiguous. Early warning signals exhibited temporal behavior that implied critical slowing down or flickering (increasing trends), recovery from perturbation (decreasing trends), unchanging forest resilience (no pattern across early warning signals). This lack of synchrony among signals, even within specific state variables, illustrates well-known difficulties of using them to measure and interpret changing resilience in situ [24, 26].

In additional to their lack of coherence, another challenge to the application of early warning signals is that they do not always perform as theory predicts [22, 23]. Analysis of early warning signals in ecological data may result in false positives, suggesting an approach to a critical threshold that does not exist [23, 24]. Early warning signal assessments may also generate false negatives, failing to detect a regime shift that has already occurred [24–26]. The periodicity, length, and variability of time series data all influence the extent to which early warning signals can detect declining resilience [25, 26]. Decisions that guide the analysis, such as whether and how to detrend the data and the length of the sliding window over which to calculate early warning signals, also affect the power to detect an approaching regime shift and the likelihood of false positives [32]. Several model-based approaches have been proposed to overcome these shortcomings [76]. For example, time-varying autoregressive state-space models can detect critical slowing down [77]. Yet these and other dynamic linear models still require subjective decisions related to data transformation and detrending; they are also prone to generating both false positives and false negatives [77]. More complex models tend to be system-specific, limiting their broad applicability [27].

Even with known limitations, empirically-derived early warning signals remain a compelling means of detecting changing resilience. Our study demonstrates their utility in identifying clear resilience shifts that are both positive and negative, particularly for stream Ca and NO₃ concentrations and fluxes. Although acid deposition caused ecosystem Ca losses at Hubbard Brook, this historic Ca depletion only impacted resilience when combined with forest harvest. In contrast, the coherence of early warning signals of stream NO₃ concentrations across watersheds implies that loss of resilience—which may be due to N oligotrophication—is a widespread phenomenon regardless of management history. Our analysis also illustrates the potential of applying early warning signals analyses to long-term, in situ observations of forests. The alignment we observed between temporal trends in early warning signals calculated from the monitoring record at Hubbard Brook and our current understanding of the system suggests that they might be used in other forest ecosystems to gain insights into changing resilience.