Tree rings provide early warning signals of jack pine mortality across a moisture gradient in the southern boreal forest

Recent declines in productivity and tree survival have been widely observed in boreal forests. We used early warning signals (EWS) in tree ring data to anticipate premature mortality in jack pine (Pinus banksiana)—an extensive and dominant species occurring across the moisture-limited southern boreal forest in North America. We sampled tree rings from 113 living and 84 dead trees in three soil moisture regimes (subxeric, submesic, subhygric) in central Saskatchewan, Canada. We reconstructed annual increments of tree basal area to investigate (1) whether we could detect EWS related to mortality of individual trees, and (2) how water availability and tree growth history may explain the mortality warning signs. EWS were evident as punctuated changes in growth patterns prior to transition to an alternative state of reduced growth before dying. This transition was likely triggered by a combination of severe drought and insect outbreak. Higher moisture availability associated with a soil moisture gradient did not appear to reduce tree sensitivity to stress-induced mortality. Our results suggest tree rings offer considerable potential for detecting critical transitions in tree growth, which are linked to premature mortality.


Introduction
Recent reports of increasing mortality of Earth's forests (Allen et al 2010, Anderegg et al 2013) suggest increasing vulnerability of forests to climatic stress. In the Canadian boreal forest, mean annual temperatures have increased by ∼1.5°C from 1961 to 2011 and climate models predict an additional warming of 2-3°C by the end of this century (Price et al 2013). Despite the recent wet period in the boreal forest of western Canada from 2003-2011 (Régnière and Saint-Amant 2008), increased evapotranspiration soon will likely to lead to an overall decline in summer soil moisture and more frequent and intense drought events (Hogg et al 2013, Ireson et al 2015. Continued warming, patterns of fluctuating precipitation, and increased frequency and magnitude of extreme climate events are likely to increase tree mortality in the Canadian boreal forest, reducing its capacity to provide ecosystem services. Complex mechanisms integrating multiple intrinsic and extrinsic factors lead to tree mortality under climatic stress (Pedersen 1998, Anderegg et al 2015. Predisposing or chronic influences, like poor competitive status or weak pathogens (such as Armillaria rot), can affect relative tree vigor over multiple years and increase susceptibility to acute, inciting stressors like insect defoliation or drought (Mallett and Volney 1990). Moreover, trees that die as a result of drought may have had altered growth sensitivity to climate due to lower whole crown gas exchange relative to trees that live (McDowell et al 2010, Hereş et al 2012. Trees that do not fully recover from an acute stress are likely to be more vulnerable to contributing factors, and death may occur years after the immediate impact of the inciting stress event (figure 1, Manion 1991).
Tree rings provide a history of growth that allows us to examine possible causal agents of mortality. For example, pests or pathogens may cause growth patterns prior to death that are distinct from direct climate effects (e.g., Hogg et al 2002, Anderegg et al 2015. As Canada's boreal forest is forecast to begin a period of rapid climate change, these retrospective studies are a valuable tool to help us anticipate changes in tree growth and mortality over coming decades. Dynamic systems like forests are likely to have tipping points where critical transitions to a different state can happen (Scheffer et al 2012). The proximity of a system to a critical transition may be detectable from 'early warning signals' (EWS) present in time series data (Dakos et al 2012, Scheffer et al 2012. Application of EWS analysis to tree rings offers a novel approach for detecting changing growth responses and inferring causes of tree mortality (see Camarero et al 2015).
Drought stress has been inferred as an inciting factor in recent observations of declines in tree survival in boreal forests (Peng et al 2011, Ma et al 2012. We wanted to know if (1) we could detect signals of premature mortality in individual trees prior to death, and (2) if trees growing in moister sites would be less sensitive to drought stress because of higher water availability because the role of soil moisture in mediating tree response to climate is not well understood. We also wanted to determine if EWS of mortality exist in growth chronologies of individual trees, and to determine how far into the past these signals could be detected. We searched for EWS in tree ring data from jack pine (Pinus banksiana Lamb.) populations with differing soil moisture regimes in the southern boreal forest using a model-based approach and tested whether these signals differed with moisture regime. We then also tested the presence of EWS in combination with other dendroclimatic, allometric, and site variables as candidate predictors to anticipate tree mortality.

Study area
Our study focused on pure jack pine stands surrounding a long-term research site in the network of boreal ecosystem research and monitoring sites (BERMS, Griffis et al 2003). Jack pine is a shade-intolerant boreal tree species, typically found on well-drained, sandy soils in even-aged monoculture stands regenerated after stand-replacing wildfire (Burns and Honkala 1990  . Conceptual model of tree growth in response to a discrete stress event (gray shading). Line A (green) represents a healthy tree over a discrete time period, and Line B (red) represents a dying tree with reduced vigor due to predisposing stress or competition prior to disturbance. Several possible outcomes may arise following a disturbance: (1) Recovery time is of the same magnitude as disturbance time (e.g., years).
(2) Recovery time is much longer compared to disturbance time (e.g., decades). (3) The tree fails to recover from the disturbance and vigor declines gradually, possibly resulting in death. (4) Direct mortality caused by the disturbance. (5) Contributing stressors lead to mortality of the dying tree within several years of the disturbance. Tree death is indicated by •. Modified from Pedersen (1998) and Manion (1991). 17.1°C in July. Approximately 448 mm of precipitation falls on the OJP site annually, with about 27% falling as snow (Régnière andSaint-Amant 2008, Environment Canada 2015).

Field data collection
In 2011 and 2012, we sampled three 10 × 10 m plots in each of three soil moisture regimes from dry to moist (subxeric, submesic, subhygric) in the study area (figure 2). We identified potential sampling locations along a moisture gradient defined by the normalized difference water index (Gao 1996), which we derived from Landsat enhanced thematic mapper imagery for 12 August 2001 acquired from the global land cover facility at the University of Maryland (GLCF 2012). Sampling locations within the predetermined moisture gradient were finalized in the field using understory vegetation composition. We measured height and breast height (BH) diameter (1.3 m above ground) on all live and dead trees currently or previously rooted in the plots and detectable at the time of sampling. We sampled one or two BH increment cores from all living trees and BH stem sections from dead trees (snags) and downed logs Lieffers 2009, Metsaranta andKurz 2012). Stem decomposition is relatively slow at OJP, making this area well suited to reconstruct tree mortality (Metsaranta et al 2008). However, as young jack pine stands often have densities exceeding several thousand per hectare (Kenkel et al 1997), we likely missed sampling trees that died during the early part of stand development. Thus, our data are best able to make inferences about mortality patterns in the later stages of stand development (Oliver and Larson 1996).

Laboratory procedures and basal area increments (BAI)
All samples (n = 197) were dried and sanded with progressively finer sandpaper until individual xylem cells were visible. We measured ring-widths with WinDendro (Regent Instruments, Québec, Canada) and cross-dated ring-width chronologies from each tree visually and statistically using COFECHA (Holmes 1992). We averaged ring widths from each tree into one time series and converted dead-and livetree ring-width measurements to BAI (Biondi and Qeadan 2008), as BAI is more representative of overall tree growth than ring width (Pedersen 1998, Bigler andBugmann 2003). We calculated BAI using the package 'dplR, version 1.6.0' in R (Bunn 2008). To avoid noise due to age-related rapid growth during the early life stages of the tree (Fritts 2001), we focused our analyses on the 1930-2012 period.

Statistical analyses
To understand the relationship between critical transitions and tree mortality, we examined two potential EWS of tree mortality in relation to site moisture and live/dead status, and then used these in combination with site, allometric, and dendroclimatic metrics as candidate predictors to anticipate tree mortality (figure 3). To determine when and which trees exhibited EWS, we analyzed BAI time series using (1) a discrete intervention approach (Pedersen 1998), and (2) an early warning detection model (drift-diffusionjump (DDJ), Carpenter and Brock 2011). All analyses were completed in R version 3.1.1 (R Development Core Team 2014).

Intervention detection
A potential EWS of mortality is the occurrence of sudden declines, or 'interventions,' in tree growth rates, potential resulting from inciting stresses (Pedersen 1998). Detection of interventions in tree growth data requires statistical methods tolerant of the autocorrelation commonly found in tree growth data (Fritts 2001). Utilizing time-series regression (TSR) analysis (see supporting information), we sought evidence of these interventions in individual tree BAI chronologies coincident with environmental stresses (see Pedersen 1998 for a thorough description of the methods; left panels in figure 3). The TSR model is an enhanced ordinary least-squares regression (OLS) sub-model, in which an autoregressive, moving average (ARMA) sub-model of the OLS regression residuals is used to satisfy OLS regression model assumptions-principally assumptions of serial-independence (Ostrom 1990). We first fit an OLS regression model, and then examined the resulting residuals to determine their autocorrelation structure (Pedersen 1998). To develop a complete TSR model, we then re-estimated the OLS regression parameters simultaneously with the residual ARMA sub-model. Based on the model results, we developed separate lists of intervention years for each tree within the six live/ dead-moisture groups.

Early warning signals
For our second approach (right panels, figure 3), we used a non-parametric DDJ model to detect EWS and investigate critical transitions of the individual BAI chronologies (see supplementary information). This approach has historically been used in finance to capture discontinuous behavior in asset pricing (e.g., Kijima 2002), though more recently has been used to detect abrupt shifts in ecosystem dynamics (e.g., Brock and Carpenter 2012). Drift measures the instantaneous deterministic change in the time series, diffusion measures relatively small shocks that occur at each time step, and jumps are large, one-step, positive or negative shocks that are uncorrelated in time. The total variance combines the contributions of diffusion and jumps Brock 2011, Dakos et al 2012). Essentially the DDJ model evaluates both high-and low-frequency and amplitude patterns in the time series, and estimates the probability of a regime shift (see Carpenter and Brock 2011 for a thorough description of the methods). Following Dakos et al (2012), we expected that conditional variance would increase and diffusion would decrease if a time series (e.g., tree growth) were approaching a regime shift (e.g., reduced growth prior to dying).

Dendroclimatic metrics
Our hypothesis was that the main climatic variable affecting mortality patterns was moisture availability. We represented interannual variations in moisture availability using the climate moisture index (CMI, Hogg 1997), which has been used to assess drought impacts on growth and mortality in western Canadian aspen (Populus tremuloides L.) forests (Hogg et al 2005, Hogg et al 2008. We calculated CMI from temperature and precipitation data output from BioSIM software (Régnière and Saint-Amant 2008) because modeled climate data were more representative of the conditions at the OJP site than the closest long-term weather record at Prince Albert, SK (figure 4). We assessed relations between patterns of moisture availability and tree growth and mortality across the soil moisture gradient by first correlating individual tree  growth chronologies to CMI and then using those correlations as potential predictors of mortality. Prior to estimating correlations between the two time series, low frequency trends in both tree ring growth and climate data were first removed by detrending using a cubic smoothing spline, with a frequency response of 0.5 at a wavelength of 0.67 × n years in the package 'dplR, version 1.6.0' (Bunn 2008). This detrending option is known to remove low frequency variability with little or no removal of high frequency variability (Cook and Kairiukstis 1990). We created pre-whitened standardized BAI and climate chronologies by using residuals from autoregressive modeling of the detrended series (Fritts 2001), thus making the data statistically independent (Legendre and Legendre 1998). Using the package 'bootRes, version 1.2.3' (Zang and Biondi 2013), we compared tree growth to annual mean CMI from August-July. Most of the latewood in jack pine is typically formed by late-July (Heinrichs et al 2007) and this period has been found to provide a good indicator of drought stress (Hogg et al 2013)-thus, we used 31 July as a cut-off for the 'dendrochronological' year. To avoid narrow rings that are typically formed during tree decline prior to death (Fritts 2001) and maximize the common interval shared by individual tree ring series (Bunn 2008), we chose to analyze growth-climate responses for  period.
An additional factor that may contribute to growth and mortality patterns at our sites is periodic defoliation by jack pine budworm (Choristoneura pinus pinus Freeman). The coarse resolution of budworm outbreaks precluded their inclusion as a mortality predictor in our logistic model. However, we felt it would be useful to present outbreak chronologies here to put observed BAI patterns into ecological context. To estimate historical patterns of jack pine budworm outbreaks, we used historical records of Canada's former forest insect and disease survey (figure 4, Volney 1988). The highest severity outbreak peaked in 1965, concomitant with severe drought during 1961 and 1964 (Hogg et al 2005, Quiring andPapakyriakou 2005). Canada's national forestry database (http://nfdp.ccfm.org) does not record any further defoliation by jack pine budworm in Saskatchewan after 1987.

Logistic model fitting
We used logistic regression to calculate the mortality risk of individual trees based on tree-level (growth rate before and after 1964, size, density of competing neighbors, and correlation with CMI), and site-level (subhygric, submesic, or subxeric in addition to variables related to EWS (interventions and DDJ)). For the DDJ variables, we calculated ln(BAI) at the maximum variance and first zero diffusion occurrence in order to capture the growth state of the tree prior to an abrupt shift in the BAI chronologies. The ecological reasoning is that if the DDJ metrics do indeed represent EWS in tree ring series, they should be significant predictors of mortality in the logistic model. To test if competition (i.e., a dense neighborhood) predicted mortality, we calculated competition among neighbors as the tree density within a circle surrounding the tree of interest defined by a radius determined by the third nearest neighbor (k = 3). Calculations were completed using the package 'spatstat, version 1.38-1' in R (Baddeley and Turner 2005). The general expression for logistic regression is as follows: Pr(Y = 1) is the survival probability of an individual tree expressed as a function of a matrix X of independent variables; Y is the dependent variable (tree status: live = 0 or dead = 1); and β is a vector containing the regression coefficients. We developed an initial 'full' model of mortality risk using all predictor variables (table 1) and removed variables until only significant predictors remained. Variance inflation factors (VIF) of the final model were examined to confirm that the models were not compromised by multicollinearity (Neter et al 1989). We evaluated model accuracy and fit using area under the receiver operator characteristic curves (AUC) and Akaike information criterion (AIC), respectively (Crawley 2013).

Stand characteristics
A greater proportion of trees died in subhygric (44 ± 11.5%; ±SE) and subxeric plots (48 ± 2.3%) than in submesic plots (24 ± 1.9%; table 2). Tree mortality peaked between 1975 and 1990 following a period of reduced growth around 1964, coinciding with the highest severity budworm outbreak on record (figures 4 and 5). In submesic areas, 63% of dead trees died between 1965 and 1990, compared to 74% in subxeric and 86% in subhygric plots during the same  A1). There were some variations in mean growth rates prior to 1964 among moisture groups (P < 0.001) and mortality status (P < 0.001) (figure 5, A1(a) and table A2). However, the most striking contrasts were found following 1964, between live and dead trees. Both mortality groups showed a drop in BAI around 1964, but the trees that would eventually die failed to fully recover ( figure 5 and A1(b)). Initial growth rates during the first 20 years following postfire stand initiation were slightly higher for live versus dead trees in subxeric and subhygric environments  (P < 0.021), but were similar in submesic areas (P = 0.051) (figure A1(c) and table A2).

Interventions
We detected a total of 260 significant (P < 0.05) interventions in the 197 trees examined. Most (180) experienced at least one intervention (figure 6). The number of interventions per tree was similar between live and dead trees at subxeric and both submesic and subhygric plots, but differed between submesic and subhygric plots (F 1,156 = 6.614, P = 0.003). Some trees experienced as many as four interventions, but the number of interventions per tree was independent of the length of the tree's BAI record (ρ = −0.015, P = 0.846). Among dead trees, 70 out of 84 trees experienced an intervention. Trees that died were an average of six years younger than live trees when they experienced their first intervention (F 1,174 = 5.209, P = 0.024). There was no correlation among numbers of interventions and CMI (ρ < 0.271, P > 0.310).

Correlation with moisture availability
At the stand level (growth averaged across all trees within a moisture regime and live/dead class), correlations between growth and annual CMI were significant for live trees, but not for dead trees (figure 8). At an individual level, the number of significant correlations was higher for live (35% of 113 trees) versus dead (17% of 84 trees) trees, as was the mean correlation (0.149 for live trees and 0.138 for dead trees; table 3). All but one significant correlation was positive (table 3). Collectively, these results suggest a lower sensitivity to moisture availability in trees that died relative to those that lived.

Logistic mortality models
Several factors emerged as strong predictors of tree mortality. After eliminating several metrics of tree allometry (DBH, DGE, GRA; table 1) because of high VIF (>10; Neter et al 1989), significant predictors included: ln(BAI) at the first zero diffusion occurrence, radial growth rates prior to 1964, age at first intervention, ln(BAI) at the maximum total variance, and correlation with CMI (table 4). Nearest neighbor density, distance to nearest neighbor, and number of interventions were not significant predictors of tree mortality. The model predicted live-dead status with 68-91% accuracy for dead trees, and 81-96% for live trees, based on a probability threshold of 0.5 ( figure 9). Across all moisture classes, the trees that would eventually die typically showed growth decreases prior to death that persisted for up to several decades   August-July, 1930-1980 and individual residual chronologies from live (green) and dead (red) trees at each site. Both the CMI and tree ring data were detrended and pre-whitened prior to analysis. Error bars represent 95% confidence intervals (* = significant).
( figure 9). Remarkably, long periods of relatively low, but stationary growth did not necessarily indicate an increased mortality probability (e.g., see the period from 1959 to 2011 of some living trees in figure 9(a)). Rather, the probability of high mortality was associated with a strongly negative (declining) growth trend (e.g., see the dead trees in figure 9(e)).

Discussion
Contrasting growth patterns of live and dead jack pine trees suggest that periods of combined biotic and abiotic stress pushed slow-growing trees across a threshold to eventual mortality. Dead trees showed a diagnostic EWS: a diminished growth rate that was detectable before the tree actually died. Thus it seems that if a tree grew slowly, its fate was sealed early on. In particular, we found a marked growth decline during a period of intense drought and budworm outbreak in the 1960s that probably led to mortality of trees unable to recover from the growth decline. Mortality probabilities were higher for trees with lower growth rates prior to 1964, experienced an early intervention, and had lower correlations with CMI.
We found EWS of mortality that suggested a threshold of declining growth leads to tree death. We successfully predicted individual mortality of virtually all of the now-dead trees, based on a mixture of growth, climate sensitivity, and DDJ metrics. However, the DDJ model further suggested that broad-scale mortality (i.e., detected in most individual BAI chronologies) can be predicted based solely on growth patterns (e.g., Chen et al 2008, Thorpe andDaniels 2012). Essentially, DDJ models identify the 'point of no return' in a natural system, where a shift to an alternate stable state has occurred (Dakos et al 2012, Carpenter et al 2014. In our study, that alternative state was not stable per se, but rather represented a period of reduced growth that ultimately led to mortality. Conditional variance intensity increased and diffusion decreased, as growth of now-dead trees approached a local minimum centered around 1964-commensurate with a major drought and severe outbreak of pine budworm (Volney 1988, Hogg et al 2002. That a deterministic components (drift) increased and the stochastic components (diffusion) of tree growth through time decreased in now-dead trees, invokes the idea that tree mortality is determined early in their development, and that in the face of severe disturbances trees with reduced vigor will have their lifetime shortened relative to healthy trees.
Mortality and growth declines were linked to an episode of combined abiotic and biotic stress, rather than climatic factors alone. The drought in the 1960s immediately preceded a period of intense budworm defoliation that caused widespread growth decline of jack pine. Pine growth responses to the 1960s drought were greater than responses to similar drought events in the 1940s and 1980s that did not precede an insect outbreak. Nevertheless, few trees died immediately following this inciting stress in the 1960 s, and we found little evidence for the direct mortality pathway (line 4 in figure 1). Instead, many of the dead trees in our population survived the inciting stress but their growth declined, taking years or decades to die (lines 3 and 5 in figure 1). Mortality following the inciting stress preferentially affected trees that had lower radial growth rates prior to the stress (Bigler andBugmann 2003, Bigler et al 2006). Experience of a pronounced negative intervention was insufficient to predict mortality: all but three of the trees that were alive in 2012 had experienced a significant negative intervention in the 1960s.
The sensitivity of jack pine trees in wetter, subhygric areas to the 1960s drought and insect outbreak was contrary to our expectations (e.g., Hereş et al 2012). These patterns may reflect local moisture dynamics driven by environmental characteristics of our study plots. First, the subhygric plots (figure 1) have a shallow water table and are susceptible to saturation in wet years (G van der Kamp, pers. comm.), and the presence of mottling in the upper 30 cm of the soil profile suggests these soils are normally wet for part of the year. Jack pine is a flood intolerant species that has limited capacity to produce adventitious roots and can quickly die following extended wet periods (Sinclair and Lyon 2005). Growth anomalies in the 1960s suggest that drought eventually led to mortality at subhygric plots, suggesting that acclimation to site-level conditions may leave trees vulnerable to drought even when the site is relatively moist. Conversely, jack pine in subxeric areas likely have a greater isohydric sensitivity than their subhygric counterparts, and are hydraulically better suited to taking advantage of periods when water is available (Addington et al 2006). However, the tradeoff is that pine in subxeric areas have smaller margins of safety from hydraulic failure during drought, and thus show a more sensitive stomatal closure response to soil water limitation (Sperry et al 1998). Consequently trees in subxeric and subhygric areas are most likely to die prematurely compared to submesic areas, though probably due to different factors. <0.001*** a DVA = ln(BAI) for the first zero location, GRP = mean BAI before 1964, AFI = tree age at first intervention, TVA = ln(BAI) for the maximum total variance, CMI = tree correlation with August-July CMI. b Change in AIC if the predictor was removed from the model. *p < 0.05; **, p < 0.01; ***, p < 0.001. Figure 9. Ln-transformed annual basal area increments (BAI in mm 2 ) during 1930-2012 from live and dead jack pine trees. Green lines represent trees that were predicted to be living-P(mortality) <0.5. Purple lines represent trees that logistic mortality models predicted to be dead-P(mortality) >0.5. Darker colors represent higher predicted probabilities of survival (1-P(mortality)) or mortality. Percentages are the number of trees for which the mortality models correctly predicted live-dead status.
Initial size heterogeneity may be important in driving competitive interactions that may predispose some trees to mortality. Jack pine with low initial growth rates had a greater probability of dying in subsequent decades, supporting observations that hierarchies tree size established within the first decades of stand development are good predictors of long-term survivorship (Kenkel et al 1997). Once trees are mature, diseases such as Armillaria rot accumulate in the roots and stems of susceptible individuals (Mallett and Volney 1990). Metsaranta and Lieffers (2010) suggest that while competition is one cause of mortality (through shading or wind-caused collisions with larger individuals), the ultimate trigger could be a high disease load that would reduce tree growth even in the absence of competition. Therefore we hypothesize that initial size hierarchies combined with differential disease build up and subsurface neighbor interactions (Tarroux and DesRochers 2011) that are not captured by metrics quantifying the density of nearest neighbors above ground (which were not a significant predictor of mortality) reduces the resilience of trees to interventions. Interestingly, recent evidence from western Canada suggests that community-level processes like competition are more important drivers of tree mortality than climate (Zhang et al 2015). However, our results and particularly the lack of competition significance suggests tree mortality in our study may be driven more by regional-scale processes like climate and other community-level processes such as tree pest attacks (McDowell et al 2011, McDowell andAllen 2015), rather than solely competition.

Conclusions
The emerging picture from our study is that predisposing factors interact with environmental and biotic stresses to determine both which trees will die and where. Our observations of jack pine mortality are consistent with expectations that predisposing influences can affect relative tree vigor over multiple years and increase susceptibility to acute, short-term stressors (Manion 1991, Pedersen 1998. We found that if a tree had reduced growth rates prior to a major intervention and suffered its first intervention at a young age, it was more likely to perish than other trees. Individuals situated along the marginal portions of an environmental gradient also appeared to be more vulnerable to stress-induced mortality.
We identified reduced growth during the first several decades of individual tree development as an EWS of tree mortality. At the stand-scale, mortality signals were evident as punctuated changes in both the deterministic and stochastic components of tree growth time series prior to a transition to an alternative state of reduced growth before dying. This transition was likely triggered by a combination of severe drought and budworm outbreak that increased mortality in the ensuing decades. Logistic mortality models also predicted imminent jack pine mortality in marginal environments. Site-specific variations in vulnerability to mortality for a widespread boreal species like jack pine suggests that combinations of climate stress and disturbance may increase fragmentation of the southern boreal forest. If climate change predictions are fulfilled and the global increase in average and extreme temperatures is accompanied by rainfall reductions, tree mortality events like the ones reported here may become increasingly widespread.