For this study, we incorporated a multiproxy approach, using tree ring basal area increment (BAI), ∆ 13 C, and δ 15 N to explore the drivers of recent growth and physiological changes of red spruce trees ( Picea rubens Sarg.), an economically and ecologically important conifer species in temperate forests of the eastern United States and a key indicator species of forest disturbance by acidic air pollution (Hamburg & Cogbill, 1988 ; Johnson & Siccama, 1983 ; Siccama, Bliss, & Vogelmann, 1982 ). Our objectives were to (1) identify those factors most related to proximate increases in red spruce growth in the Central Appalachian Mountains; (2) examine the relationships between tree growth, reconstructed seasonally integrated leaf physiology, including i WUE, and environmental factors; and (3) explore the use of tree ring δ 15 N as an integrator of historical ecosystem N cycling. To meet these objectives, we examined red spruce tree growth using tree rings at three forest locations spanning a 100 km transect in the Central Appalachian Mountains where we have previously examined drivers of ecosystem N availability (Smith, Mathias, McNeil, Peterjohn, & Thomas, 2016 ). We used Kendall's rank correlation analyses and generalized linear mixed models (GLMMs) to disentangle the multiple environmental effects and their sensitivities on red spruce tree growth into discreet components.

Tree ring growth and isotope chronologies have become important proxies used to examine historical changes in productivity, plant physiology, and ecosystem processes in a rapidly changing environment (Dye et al., 2016 ; McLauchlan et al., 2017 ; Thomas et al., 2013 ). While tree ring chronologies provide an opportunity to examine changes in forest productivity (Dye et al., 2016 ), these datasets are strengthened by linking observed changes in growth to underlying physiology (Belmecheri et al., 2014 ; Thomas et al., 2013 ). In particular, carbon isotope signatures of tree rings provide information regarding the diffusive and biochemical processes of photosynthesis (Farquhar, Ehleringer, & Hubick, 1989 ; Farquhar, O'Leary, & Berry, 1982 ). Nitrogen isotope signatures (δ 15 N) of tree rings record changes in ecosystem N status, although the exact mechanisms have yet to be fully reconciled (Evans, 2001 ; Gerhart & McLauchlan, 2014 ). It has been postulated that tree ring δ 15 N reflects pollutant emissions and deposition of oxidized N compounds (Elliott et al., 2007 ; Gurmesa et al., 2017 ; Hietz, Dünisch, & Wanek, 2010 ); however, there is also abundant evidence that δ 15 N reflects ecosystem N cycling dynamics (Burnham, McNeil, Adams, & Peterjohn, 2016 ; Craine et al., 2009 , 2015 ; Evans, 2001 , 2007 ; Hobbie & Högberg, 2012 ; McLauchlan et al., 2017 ; Nadelhoffer & Fry, 1994 ). These proxies become more informative when an integrative approach is taken, combining multiple indicators of ecosystem status in order to track long‐term changes in ecosystem function (Choi, Lee, Chang, & Ro, 2005 ; Kwak, Lim, Chang, Lee, & Choi, 2011 ; Leonelli et al., 2012 ; Zeng, Liu, Xu, Wang, & An, 2014 ).

In addition to the declines in acidic air pollution since the Clean Air Act, there have been substantial increases in atmospheric CO 2 concentrations (Figure S1 c), along with significant changes in climate (Figure S1 d). Increasing atmospheric CO 2 concentrations stimulates forest productivity, thereby enhancing carbon storage (Fernández‐Martínez et al., 2017 ; Norby et al., 2005 ); however, these responses are often modulated by changes in temperature or precipitation (Granda, Rossatto, Camarero, Voltas, & Valladares, 2014 ; Jennings, Guerrieri, Vadeboncoeur, & Asbjornsen, 2016 ), in addition to acid pollutants (Fernández‐Martínez et al., 2017 ). It is, therefore, critical to understand how each of these factors contribute to the recovery of trees from acid pollution and the extent, and nature, of any interactive effects between tree recovery and environmental change.

For decades, temperate forests in the eastern United States have experienced anthropogenic disturbance, including pervasive levels of acidic air pollution from fossil fuel combustion (Greaver et al., 2012 ). Before environmental legislation designed to curb pollutant emissions became effective, national levels of SO 2 and NO x emissions were as high as 29.6 and 25.6 million metric tons annually (Lefohn, Husar, & Husar, 1999 ; EPA, 2015 ) (Figure S1 a,b). Since the Clean Air Act and subsequent amendments, SO 2 and NO x emissions have fallen from their highest levels by 84% and 56%, respectively (Lefohn et al., 1999 ; EPA, 2015 ). Although significant progress has been made in understanding C cycling following many types of disturbance in forest ecosystems (Liu et al., 2011 ), mechanisms that govern recovery from acid deposition following reductions in acidic air pollution are still poorly understood. It has been projected that the recovery of terrestrial ecosystems following long‐term acid deposition takes decades (Driscoll et al., 2001 ; Likens, Driscoll, & Buso, 1996 ), and now over four decades after the Clean Air Act of 1970, there remains scant evidence of forest ecosystem recovery from the long history of disturbance by acid deposition in the eastern United States. (Engel et al., 2016 ; Lawrence et al., 2015 ; Li et al., 2010 ; Thomas et al., 2013 ).

Recent observations of enhanced tree growth in temperate forest ecosystems of the northern hemisphere have raised questions about the underlying causal mechanisms (Boisvenue & Running, 2006 ; Engel et al., 2016 ; Fang et al., 2014 ; Johnson & Abrams, 2009 ; McMahon, Parker, & Miller, 2010 ; Thomas, Spal, Smith, & Nippert, 2013 ). Is this increased productivity a consequence of forest regrowth following disturbance, or is tree growth being enhanced by environmental changes, such as CO 2 fertilization and climate change (Joos, Prentice, & House, 2002 ; Schimel et al., 2000 )? These are critical questions because forest ecosystems are a significant part of the terrestrial carbon (C) sink that removes nearly 30% of anthropogenic C emissions each year, providing a negative feedback that moderates climate change (Le Quéré et al., 2018 ). Yet, a large amount of uncertainty exists about the continued capacity of forests to sequester C emissions (Friedlingstein et al., 2006 ). Initially, forest regrowth after disturbance provides a large C sink, but this sink may not persist as forest productivity slows during the late stages of stand development, thereby limiting the capacity of the C sink (Caspersen et al., 2000 ; Fang et al., 2014 ). The persistence of the effects of increased CO 2 and climate change on forest C sequestration are harder to predict and are likely to interact with the dynamics of forest recovery, possibly resulting in either positive or negative feedbacks on the C sink (Anderson‐Teixeira et al., 2013 ). In order to partition the relative contributions by these drivers, and how their complex interactions affect forest C sinks, we need a deeper understanding of how simultaneous changes in key environmental variables affect long‐term C cycling in forest ecosystems.

The average model was used to identify the relative effect each parameter had on BAI for the 1989–2014 period (Fernández‐Martínez et al., 2017 ). To do this, we first used the average model to predict BAI for a given year. We then ran the model to predict BAI holding a given environmental parameter constant at its respective median value for the region over 1989–2014 (atmospheric CO 2 : 371.5 ppm, NO x emissions: 20.52 million metric tons, SO 2 emissions: 14.04 million metric tons, mean April temperature: 10.08°C), while letting all other parameters change as observed. The difference in slope between BAI predicted from the original model and BAI predicted when holding a particular environmental parameter constant was the contribution of that parameter to predicted BAI. The sensitivity of BAI to each parameter was then determined by dividing the contribution of given parameter to BAI by its own trend over 1989–2014. Finally, we determined the total contribution each environmental parameter had on the total change in BAI over 1989–2014 by multiplying the sensitivity of BAI to each respective environmental parameter by that parameter's total change over 1989–2014 (atmospheric CO 2 : +46.8 ppm, NO x emissions: −13.78 million metric tons, SO 2 emissions: −19.50 million metric tons, mean April temperature: +1.82°C). Standard errors were propagated throughout all calculations using the R package ‘propagate’ (Spiess, 2017 ).

To examine the influence and interactions of environmental parameters on red spruce BAI from 1940 to 2014 and from 1989 to 2014, we utilized generalized linear mixed models (GLMMs) and model averaging using the R packages ‘nlme’ and ‘MuMIn’ (Bartoń, 2017 ; Pinheiro et al., 2018 ). Models differing by less than four AICc units, with respect to the best model, were averaged to identify the nature of the relationship between predictor variables and BAI (Fernández‐Martínez et al., 2017 ). The environmental parameters included in this analysis were selected a priori based on their known importance to tree growth and physiology and informed by the results of the initial Kendall's rank correlation and bootsrapped Pearson's correlations. For this analysis, we included atmospheric CO 2 , national emissions of SO 2 , national emissions of NO x , and mean April temperatures as predictor variables, exploring all possible combinations of models, including interactions. We also used combinations of other environmental factors in the GLMM averaging, but the four factors listed above consistently provided the best model fit, capturing the greatest variation in BAI while minimizing the unidentified biotic or abiotic contributions to the temporal variation in BAI, determined as an ‘unknown’ residual effect from the consensus model. In the GLMMs, we accounted for temporal autocorrelation using an AR(1,0) structure, and included environmental parameters as fixed effects, while site was included as a random effect.

Although trends in pollutant emissions and CO 2 were clear across the 75‐year chronology, changes in climate variables were less clear due to the high interannual variability in temperature and precipitation that occurs in this region of the Central Appalachian Mountains. Thus, to carefully examine the potential influences of temperature and precipitation on the chronologies of BAI, ∆ 13 C, A , and g c , we first examined the linear trends of annual, growing season, and monthly climate variables to determine if significant changes across 1940–2014 had occurred. Second, we used analysis of covariance (ANCOVA) to determine whether trends in temperature or precipitation were different before and after 1989. Third, we used analysis of variance (ANOVA) to determine if climate variables during the period 1940–1989 were different from 1989 to 2014. Finally, ANCOVA was used to determine whether changes in the responses of BAI, A , and g c to temperature or precipitation were different before and after 1989.

Kendall's rank correlation, a nonparametric correlation analysis (Sokal & Rohlf, 2011 ), was used to examine the relationships between potential explanatory environmental parameters and the dependent variables BAI, ∆ 13 C, A , and g c over the 1940–2014 chronology, as well as the period 1989–2014. For these analyses, environmental parameters included US emissions of SO 2 and NO x (Lefohn et al., 1999 ; EPA, 2015 ), atmospheric CO 2 concentrations (Keeling et al., 2015 ), and climate variables (annual, growing season, and monthly mean temperatures, and precipitation) (WV Climate Division 4, NOAA, 2017 ). We also examined the potential effects of wet deposition of SO 4 2− and NO 3 − in precipitation (WV18; NADP, 2015 ); however, since NADP measurements only partially covered our 1940–2014 chronology, and because the trends of pollutant emissions and trends of atmospheric wet deposition were strongly related (Figure S1 a,b), we only used pollutant emissions in these analyses. For δ 15 N, we used Pearson product‐moment correlation to examine the relationships between wood δ 15 N, US NO x emissions (EPA, 2015 ), and NO 3 − deposition (NADP, 2015 ), as well as the relationships between δ 15 N and Δ 13 C, seasonally integrated A , g c , and BAI of the red spruce trees.

Breakpoint analyses were performed using the R package ‘segmented’ (Muggeo, 2008 ) with the red spruce chronologies of BAI, ∆ 13 C, and δ 15 N to determine specific critical years where significant directional shifts in the 1940–2014 chronologies of BAI, ∆ 13 C, and δ 15 N occurred. For all further analyses used to examine environmental effects on these proxies of tree growth, photosynthetic physiology and ecosystem N status, we used 1989 as the critical year where this directional shift occurred, which was determined as the average breakpoint of BAI, ∆ 13 C, and δ 15 N chronologies (Table 1 ). We used linear mixed effects (LME) models using the ‘nlme’ package in R (Pinheiro, Bates, DebRoy, Sarkar, & Team, 2018 ) with year as the single predictor variable, accounting for temporal autocorrelation with an AR(1,0) structure, and including site as a random factor, to examine temporal trends in BAI, A , g c , and δ 15 N for before (1940–1989) and after (1989–2014) the critical year. We also used breakpoint analysis to identify a critical directional shift that occurred in the i WUE chronology and used the same LME model structure to determine changes in i WUE before and after this year.

For the calculation of A and g c , we measured A ‐C i response curves from at least three red spruce trees per site during June, July, and August 2014. Between 0900 and 1600 EST, we harvested small branches containing intact sun needles from the upper canopy and immediately placed them into water picks to maintain the transpiration stream. A‐ C i curves were measured using an open‐flow gas exchange system with red/blue LED lights (LI6400XT, Li‐Cor, Inc., Lincoln, NE) by placing intact needles into the cuvette at saturating light (1,500 μmol photons m −2 s −1 ). After ca . 5 min equilibration in the cuvette, rates of photosynthesis were measured at nine ambient CO 2 concentrations between 50 and 1,200 μl CO 2 L −1 air beginning at ambient CO 2 (400 ppm) at the time of the study. After gas exchange analyses, all needles were taken to the laboratory and scanned to determine projected leaf area, and gas exchange of red spruce was expressed on a total leaf area basis.

In order to simulate changes in seasonally integrated photosynthesis (), stomatal conductance (), and intrinsic water use efficiency (WUE) across the tree ring chronology from 1940 to 2014, field measurements of the relationship betweenand leaf intercellular CO(C) (–Ccurves) were coupled with the ∆C chronology, following the methods of (Thomas et al.,). First, ∆C was converted to C/Cthe ratio of leaf intercellular COto atmospheric CO, usingwhereis the diffusional fractionation of COthrough the stomata (4.4‰) andis the biochemical fractionation of COby Rubisco (27‰) (Farquhar et al.,). For each individual year of the chronology, the respective C/Cwas multiplied by the atmospheric COconcentration (Keeling, Piper, Bollenbacher, & Walker,) to determine C. Values of isotopically derived Cwere then used to predict the seasonally integratedusing‐Crelationships derived from measurements at each of our three sites. Stomatal conductance was calculated using methods from (Farquhar & Sharkey,) usingwhereis seasonally integrated photosynthesis, Cis the atmospheric COconcentration and Cis the leaf intercellular COconcentration.

A subset of five randomly chosen tree cores was selected from each location to develop C and N isotope chronologies for the red spruce trees. As with BAI, C and N isotopes were averaged for similar years and forest site was a replicate (= 3). Samples were analyzed for their stable C or N isotope composition with a ThermoFisher Delta V+ isotope ratio mass spectrometer interfaced with a Carlo Erba NC 2500 Elemental Analyzer. For these analyses, tree cores were dissected under microscope by scalpel into individual growth rings at the boundary of earlywood and latewood for years 1940–2014. Tree ring whole wood from each individual year was then finely chopped to homogenize the sample, and 1 ± 0.2 mg was packed into tin capsules for C isotopic analysis. Carbon isotope composition was measured as δC (‰) usingwhereis the ratio ofC:C in the wood sample andis the ratio ofC:C in the standard PeeDee belemnite (PDB) from the PeeDee River Formation in Hemingway, South Carolina. The conversion of leaf carbohydrate to wood was then accounted for by subtracting 3‰ from the wood δC (Leavitt & Long,; Thomas et al.,). Leaf‐corrected δC values were converted to isotope discrimination (∆C) to account for the progressive depletion in atmospheric δCOdue to fossil fuel burning usingwhere δand δare the stable carbon isotope compositions of the atmosphere and plant, respectively (Farquhar et al.,). Atmospheric δC values for a given year up to 2004 used to calculate ∆C were interpolated using the high precision records of atmospheric δC obtained from Antarctic ice cores (Francey et al.,) by McCarroll and Loader (). After 2004, atmospheric δC was calculated as the average change per year (−0.0199‰) during the 1940–2004 period (Keeling, Piper, Bollenbacher, & Walker,). For this study, needle structure, mesophyll conductance, and mesophyll COfractionation of red spruce trees were assumed to have not changed across the 75‐year tree ring chronology (Seibt, Rajabi, Griffiths, & Berry,).

At each site, 20–25 red spruce trees, greater than 10 cm diameter at breast height (1.4 m aboveground), were cored using a 5.15 mm increment borer (Haglöf Inc., Madison, MS). Trees were randomly selected to avoid any potential biases during growth analyses, and cores were taken perpendicular to the slope to avoid compression and expansion wood. Twenty cores from MCG were harvested at the end of the 2013 growing season, with five additional cores harvested after the 2014 growing season. At SOR and CGL, all tree cores were collected at the end of the 2014 growing season. Cores were returned to lab where they were processed according to standard dendrochronological techniques, as described by Stokes and Smiley (Stokes & Smiley,). Chronologies were crossdated using WinDENDRO (Regent Instruments, Inc.), and statistically confirmed using dplR (Bunn et al.,), after which an individual calendar year was assigned to each growth ring. The average age of sampled red spruce trees from the three locations was 109 ± 27 years (mean ±). The early portion of the tree ring chronologies represented growth release of saplings after logging in the early 1900s and therefore, the focus of this study was the time period after juvenile growth, 1940–2014. The expressed population signals (EPS), a measure of the common variance in a chronology, for each red spruce site across years 1940–2014 were all greater than 0.85, indicating that our sample size adequately represents the tree population (Wigley, Briffa, & Jones,). Raw ring widths of each year were converted to basal area increment (BAI) to better characterize tree growth and to minimize the effects of tree size and age on annual growth trends usingwhereis the radius of the tree at a given year of ring formation,. At each site, BAI was averaged for each year using the 20–25 red spruce trees at each location to develop site‐level chronologies, which we treated as replicates (3) in this study.

Three red spruce forest stands were chosen along a north–south transect of 100 km along the Central Appalachian Mountains. McGowan Mountain (MCG; 38˚58′N, 79˚41′W) is positioned within the Appalachian Plateau physiographic province, while Span Oak Run (SOR; 38˚37′N, 79˚46′W), and Cranberry Glades (CGL; 38˚12′N, 80˚17′W) are located in the Allegheny Mountain Section physiographic province of West Virginia (USGS, 1946 ). Each site is ≥900 meters in elevation and has a southwest aspect with slopes ranging from 0% to 10%. The mean annual temperature for this region in West Virginia was 9.5°C over the period of 1940–2014 (WV Climate Division 4, NOAA, 2017 ). The annual total precipitation across this time period was 1229.6 mm/year, with 47% of the precipitation during the growing season from May to September (WV Climate Division 4, NOAA, 2017 ). Soils at each of the three sites consist of a thick organic layer with litter contributed by red spruce ( Picea rubens Sarg.), eastern hemlock ( Tsuga candadensis ), yellow birch ( Betula alleghaniensis ), and red maple ( Acer rubrum ) as dominant canopy species. These sites are mature closed canopy stands with no evidence of anthropogenic or natural disturbance, such as logging or fire, since logging events in the early 1900s. Overall tree density at MCG, SOR, and CGL was 534, 617, and 675 trees/ha, with red spruce contributing ca . 23%, 25%, and 48% of the trees, respectively. Soils at MCG are classified as ernest and are moderately well to poorly drained, very deep, and have parent material that is colluvium derived from acid shale, siltstone, and sandstone (USDA Soil Survey Geographic (SSURGO), 2015). SOR and CGL are snowdog series soils and are very deep, moderately well drained, and colluvium derived from acid sandstone, siltstone, and shale (USDA Soil Survey Geographic (SSURGO), 2015). The underlying geology at MCG and CGL consists of Pennsylvanian Pottsville sandstone, while SOR is underlain by Devonian Chemung Shale (WVGES, 1968 ). Thus, the soil buffering capacities are different at each site with MCG having the lowest buffering capacity and CGL having the greatest buffering capacity.

(a) Contribution of atmospheric CO, NOemissions, SOemissions, and mean April temperatures to the change in BAI each year predicted by the GLMM average model for 1989–2014, and (b) the contribution of CO, NOemissions, SOemissions, and mean April temperatures to the total change in BAI over 1989–2014. Numbers in brackets indicate the direction and magnitude of changes in environmental parameters. Numbers in brackets in panel (a) represent the trend in each respective environmental parameter over 1989–2014, while numbers in brackets in panel (b) represent the total change in each respective environmental parameter over 1989–2014. Units for COand April Tare ppm and °C, respectively, while units for NOand SOare 10metric tons. Unknown contributions in (a) were calculated as the difference between the observed change and all known contributions. Models with ∆AICc < 4 were included in the average model (see Table S3 for model structure). Asterisks (*) indicate significance at the α = 0.05 level

In order to specifically examine the recent increases in BAI from 1989 to 2014, model averaging of the best models (∆AICc < 4 from the best model; hereafter referred to as ‘average model’) determined by GLMMs indicated the importance of mean April temperatures, a negative interaction between SO 2 and NO x emissions, a positive interaction between atmospheric CO 2 and SO 2 emissions, and a positive interaction between atmospheric CO 2 and NO x emissions (Table S3 ). We found that over the 1989–2014 period the average model predicts an increase in BAI of 0.29 ± 0.02 cm 2 /year ( t = 15.60, p < .001) (Figure 6 a), in close agreement with observed BAI ( R 2 = .85, p < .001) (Figure S5 ). Of the 0.29 cm 2 /year increase in BAI, increasing atmospheric CO 2 contributed to a 0.17 ± 0.02 cm 2 /year increase in BAI ( t = 96.05, p < .01). Reductions in national emissions of SO 2 contributed to a 0.11 ± 0.02 cm 2 /year increase in BAI ( t = 17.18, p < .01), while decreasing emissions of NO x reduced BAI by 0.06 ± 0.02 cm 2 /year ( t = −7.68, p < .01). Higher April temperatures contributed to a small increase in BAI by 0.04 ± 0.03 cm 2 /year ( t = 3.65, p < .01) (Figure 6 a). Unknown contributions due to unidentified biotic or abiotic factors were calculated as the difference between the observed change and all known contributions and accounted for 0.04 ± 0.05 cm 2 /year ( t = 1.08, p = .28). Sensitivity analysis of GLMM results revealed that over 1989–2014 average BAI was the most sensitive to increases in April temperatures and reductions in SO 2 emissions, followed by reductions in NO x, and increases in atmospheric CO 2 (Table 3 ).

Kendall's rank correlation revealed that BAI of red spruce trees from 1940 to 2014 declined as atmospheric CO 2 increased (τ = −0.14, p = .002), as pollutant emissions of NO x (τ = −0.34, p < .0001) and SO 2 increased (τ = −0.12, p = .01) (Table S1 ). GLMMs identified only a negative relationship between NO x emissions and BAI over 1940–2014 ( p < .001, −0.353 ± 0.05 slope estimate ± standard error), with no other models with a ∆AICc < 4. Kendall's rank correlation showed increases in BAI from 1989 to 2014 were correlated with increases in atmospheric CO 2 (τ = 0.41, p < .0001), increases in April temperatures (τ = 0.26, p < .001) and reductions in pollutant emissions of SO 2 (τ = −0.42, p < .0001) and NO x (τ = −0.36, p < .0001). Additionally, the temperature response of BAI to April temperatures was greater during 1989–2014 than during 1940–1989 (ANCOVA; F = 6.96, p = .009) (Figure 4 b).

In this study, there was a 51% increase in seasonally integrated i WUE across the 75‐year red spruce chronology ( p < .0001; Figure 3 c), where i WUE was positively correlated with A ( r = .25, p = 0.0002) and negatively correlated with g c ( r = −.83, p < .0001). There was, however, a breakpoint in i WUE averaged across the three sites at 1995 ± 2.5 years, where i WUE declined 13% after that year ( p < .0001; Figure 3 c). Across the 75‐year chronology, BAI of red spruce trees was not correlated with i WUE ( p = 0.10), but BAI was positively correlated with A ( r = .52, p < .0001) and g c ( r = .39, p < .0001). Kendall's rank correlation revealed that i WUE during 1940–1995 was positively correlated with atmospheric CO 2 (τ = 0.71, p < .0001), NO x (τ = 0.61, p < .0001), mean March temperatures (τ = 0.14, p = .006), and mean November temperatures (τ = 0.13, p = .013). There was a negative correlation between i WUE during 1940–1995 with June precipitation (τ = −0.16, p = .002), mean May temperatures (τ = −0.13, p = .015), and mean June temperatures (τ = −0.15, p = .006), mean October temperatures (τ = −0.13, p = .015), and mean growing season temperatures (τ = −0.10, p = .05). After the 1995 breakpoint, after which i WUE declined, Kendall's rank correlation indicated that i WUE was not significantly correlated with any of the environmental factors that we examined.

(a) Relationships between seasonally integrated photosynthesis and (b) basal area increment with mean April temperatures before and after 1989. ANCOVA was used to determine if the response of seasonally integrated photosynthesis or basal area increment to April temperatures were different before and after 1989, the average critical year for the region determined via breakpoint analysis where a directional shift in BAI, ∆ 13 C, and δ 15 N occurred. Climate data are from West Virginia Division 4, which is the region that contains the three red spruce stands in this study (WV Climate Division 4, NOAA, 2017)

(a) Chronology of seasonally integrated photosynthesis modeled using ∆ 13 C for red spruce trees ( N = 5 per site) from three forest sites in the central Appalachian Mountains from 1940 to 2014, (b) seasonally integrated stomatal conductance modeled using ∆ 13 C, and (c) intrinsic water use efficiency ( i WUE) derived from the chronologies of ∆ 13 C. The vertical dashed lines in (a) and (b) are at 1989, the average critical year for the region determined via breakpoint analysis where a directional shift in BAI, ∆ 13 C, and δ 15 N occurred. The vertical dashed line in (c) is at 1995, the critical year where a directional shift in i WUE occurred. The three red spruce forest sites are McGowan Mountain (MCG), Span Oak Run (SOR), and Cranberry Glades (CGL)

Chronologies of seasonally integrated leaf physiology modeled using ∆ 13 C revealed an 11% reduction in A ( p < .0001; Figure 3 a) in red spruce and a 47% reduction in g c ( p < .0001; Figure 3 b) from 1940 to 1989. From 1989 to 2014, there was a 43% increase in A ( p < .0001; Figure 3 a) and a 67% increase in g c ( p < .0001; Figure 3 b). Using historical monthly, seasonal, and annual climate variables (WV Climate Division 4, NOAA, 2017 ), atmospheric CO 2 concentrations (Keeling et al., 2015 ), and US SO 2 and NO x emissions (Lefohn et al., 1999 ; EPA, 2015 ) to examine environmental influences on A and g c of red spruce trees, Kendall's rank correlation revealed that A across the entire 1940–2014 chronology was negatively correlated with SO 2 (τ = −0.39, p < .0001) and NO x (τ = −0.27, p < .0001) and positively correlated with atmospheric CO 2 (τ = 0.20, p < .0001), mean annual temperatures (τ = 0.22, p < .0001), and mean April temperatures (τ = 0.17, p = .0002). Across the entire chronology, g c was negatively correlated with NO x emissions (τ = −0.66, p < .0001) and atmospheric CO 2 (τ = −0.43, p < .0001) and positively correlated with May temperatures (τ = 0.13, p < .01), June temperatures (τ = 0.15, p < .01), and June precipitation (τ = 0.14, p < .01). Increases in A from 1989 to 2014 are correlated with increases in atmospheric CO 2 (τ = 0.73, p < .0001), increases in April temperatures (τ = 0.23, p < .01) and reductions in pollutant emissions of SO 2 (τ = −0.72, p < .0001) and NO x (τ = −0.69, p < .0001). Furthermore, the temperature response of A was higher for a given April temperature during 1989–2014 than during 1940–1989 (ANCOVA; F = 8.94, p = .003) (Figure 4 a). Increases in g c from 1989 to 2014 are correlated with increases in atmospheric CO 2 (τ = 0.63, p < .0001), increases in April temperatures (τ = 0.21, p < .01) and May temperatures (τ = 0.21, p < .01), as well as reductions in pollutant emissions of SO 2 (τ = −0.62, p < .0001) and NO x (τ = −0.59, p < .0001). These environmental factors were the most important influences on A and g c and complete Kendall's rank correlation analyses are shown in Table S1 for 1940–2014 and in Table 2 for 1989–2014.

Mean annual temperatures in this region have increased significantly over the past 75 years by 0.50°C (WV Climate Division 4, NOAA, 2017 ) (Figure S2 a; R 2 = .04, p = .04). In contrast, there were no significant changes in mean annual precipitation (Figure S2 b), growing season temperature (Figure S2 c), or growing season precipitation from 1940 to 2014 (Figure S2 d) (WV Climate Division 4, NOAA, 2017 ). With the exception of April mean temperatures, temperature and precipitation in this region have not changed significantly at or around 1989 in a way that could affect BAI, ∆ 13 C, and δ 15 N. Mean April temperatures increased by 0.06°C/year from 1989 to 2014, but showed no trend prior to 1989 (Figure S1 d; ANCOVA, F = 3.90, p = .05), and therefore, mean April temperatures were 0.72°C warmer during 1989–2014 than during 1940–1989 (ANOVA, F = 4.10, p = .05).

We examined whether there were changes that occurred near 1989 for one or more environmental factors that could possibly explain the observed near simultaneous changes in BAI, tree ring ∆ 13 C and δ 15 N signatures. From 1940 to 1989, SO 2 emissions in the United States averaged 23.13 million metric tons/year, but declined by 0.73 million metric tons/year from 1989 to 2014 (Figure S1 a; ANCOVA, F = 129.7, p < .001) (Lefohn et al., 1999 ; EPA, 2015 ). From 1940 to 1989, NO x emissions in the United States increased by 0.92 million metric tons/year ( R 2 = .92, p < .01), but declined by 0.56 million metric tons/year from 1989 to 2014 (Figure S1 b; R 2 = .91, p < .001) (EPA, 2015 ). On average, atmospheric CO 2 concentration increased from 311 ppm in 1940 to 397 ppm in 2014, or 1.2 ppm/year (Figure S1 c) (Keeling et al., 2015 ). Moreover, the rate at which atmospheric CO 2 concentration increased was 1.0 ppm/year greater during the 1989–2014 period than during 1940–1989 (ANCOVA, F = 136.5, p < .01) (Keeling et al., 2015 ).

Analysis of the nitrogen isotope composition in whole wood of red spruce tree rings for the 48 years prior to 1989 revealed a mean δ 15 N signature of 0.404 ‰, but no trend between 1940 and 1989 ( p = .23; Figure 2 b). After 1989, however, the δ 15 N chronology shows an abrupt, significant depletion in the heavier N isotope in all three red spruce locations as δ 15 N progressively declined, on average, by 0.083‰ each year from the mean δ 15 N signatures of 0.404 ‰ before 1989 ( p < .0001; Figure 2 b). The average N concentrations in red spruce tree rings were 32% greater after 1989 (0.059%) compared to wood N concentrations before 1989 (0.045%) ( F = 47.37, p < .0001).

Using 1989 as the overall breakpoint year, red spruce BAI showed a 0.16 cm 2 /year decline, or a total reduction of 48.9%, between 1940 and 1989 ( p < .0001), but a 0.34 cm 2 /year increase, or a total increase of 105.8%, from 1989 to 2014, ( p < .0001) (Figure 1 ). Tree ring ∆ 13 C signatures declined by 17% between 1940 and 1989 from 20.7 ‰ to 17.7 ‰ ( p < .0001), after which ∆ 13 C signatures increased steadily by 0.070 ‰ each year from 1989 to 2014 ( p < .0001) (Figure 2 a). Growth rates following the 1989 breakpoint were not different than growth rates during the juvenile period from 1908 to 1941 (ANCOVA, F = 2.92, p = .09; Figure 1 ). In addition, tree ring ∆ 13 C signatures for the year 2013 were not different from those for the year 1941 (ANOVA, F = 0.652, p = .47; Figure 2 a).

Breakpoint analysis of the 1900–2014 red spruce BAI chronology revealed two distinct breakpoints, with the first at 1941 ± 2.1 years, marking the end of the juvenile growth phase of the red spruce trees (Figure 1 ). The second breakpoint in BAI was 1988 ± 2.4 years, where a directional shift from declining BAI to increasing BAI was observed (Figure 1 , Table 1 ). Breakpoints in 1940–2014 chronologies of ∆ 13 C and δ 15 N exhibit directional changes in the late 20th century that were contemporaneous with those in BAI. Breakpoint analysis for tree ring ∆ 13 C indicates a directional shift in tree ring ∆ 13 C signatures at 1991 ± 1.5 years (Figure 2 a, Table 1 ). For tree ring δ 15 N, breakpoint analysis shows a directional shift in tree ring δ 15 N signatures at 1987 ± 3.8 years when averaged over the three sites (Figure 2 b, Table 1 ). One‐way ANOVA indicated breakpoint years were not significantly different for a given response variable (BAI, ∆ 13 C, or δ 15 N) across sites ( F = 0.74, p = .52), nor for a given site (MCG, SOR, or CGL) across response variable ( F = 0.83, p = .53). Thus, the overall average year where directional changes were observed in the 1940–2014 chronologies of BAI, ∆ 13 C, and δ 15 N at all three sites was 1989 ± 1.5 years.

4 DISCUSSION

A synchronous change in three tree ring proxy indicators, BAI, ∆13C, and δ15N, identified 1989 as the average year that three forest stands of red spruce trees, spanning 100 km in the Central Appalachian Mountains, began a 25 year period of accelerated tree growth (Table 1). Prior to 1989, there was a 49 year period of declining tree growth, characteristic of reduced growth of red spruce trees exposed to high levels of acidic air pollution and deposition in the eastern United States. (Cook, Johnson, & Blasing, 1987; Eager & Adams, 1992; Johnson, Cook, & Siccama, 1988). Our observations of recent increases in growth of red spruce trees are similar to observations of increased growth of red spruce in the northeastern United States (Engel et al., 2016) and other tree species in the eastern United States (Fang et al., 2014; Johnson & Abrams, 2009; McMahon et al., 2010; Thomas et al., 2013). In our study, GLMM model averaging used to examine the environmental influences on growth of red spruce from 1989 to 2014 indicated recovery from reduced acidic air pollution, a fertilization effect of increased atmospheric CO 2 , and a positive influence of early spring temperatures, thus highlighting the importance of multiple environmental factors on the proximate increases in red spruce tree growth in the Central Appalachian Mountains.

Since the Clean Air Act of 1970 and subsequent amendments, pollutant emissions in the United States have dropped considerably. From 1940 to 1989, SO 2 emissions averaged 23.13 million metric tons/year, but declined by 0.73 million metric tons/year from 1989 to 2014 (Figure S1a) (Lefohn et al., 1999; EPA, 2015). Kendall's rank correlation placed the decline in SO 2 emissions as the topmost environmental factor contributing to the observed growth increases in red spruce trees since 1989 (Table 2). GLMM model averaging also pointed to a strong effect of decreasing SO 2 emissions on tree growth across 1989–2014, as well as a high sensitivity to these pollutant reductions (Figure 6, Table 3). The changes in the red spruce ∆13C chronology, a proxy indicator for photosynthetic physiology, further support recovery of these trees from acid pollution. The decrease in ∆13C in the red spruce tree rings from 1940 to 1989, corresponding to reductions in A and g c , has been shown in many studies examining the effects of acidic sulfur pollution on the C isotope signatures of trees (Boettger, Haupt, Friedrich, & Waterhouse, 2014; Kwak et al., 2016; Rinne, Loader, Switsur, Treydte, & Waterhouse, 2010; Santruckova, Santrucek, Setlik, Svoboda, & Kopacek, 2007; Savard, Bégin, Parent, Smirnoff, & Marion, 2004), including species from the eastern United States (Li et al., 2010; Thomas et al., 2013). The increase in ∆13C, corresponding to increases in A and g c , that we observed after 1989 as SO 2 emissions declined, has also been found in studies as a function of decreasing pollution (Boettger et al., 2014; Li et al., 2010; Thomas et al., 2013), and is opposite from observations of reduced tree ring ∆13C with increased tree height or as trees get older (Brienen et al., 2017; McDowell, Bond, Dickman, & Ryan, 2011). Thus, our data are consistent with the hypothesis that reductions in acidic air pollution and deposition from the historically high totals in the 1970s reached a critical level around 1989 that red spruce trees were better able to tolerate, thereby beginning the recovery of tree growth and photosynthetic physiology of red spruce trees in the Central Appalachian Mountains that has persisted for 25 years as pollution levels continued to decline.

Our study suggests that increases in growth and A after 1989 were likely enhanced by increases in atmospheric CO 2 concentrations, which have increased substantially over the last 75 years and have increased by ca. 13% from 1989 to 2014 (Keeling et al., 2015). Although numerous experiments using large square‐wave increases in CO 2 , such as FACE experiments, often show increased tree growth (Norby et al., 2005), an effect of increasing CO 2 on growth is rarely observed in tree ring studies, possibly because the annual step changes in CO 2 are small and any changes in growth are masked by other environmental factors positively or negatively contributing to growth. In our study, both Kendall's rank correlation and GLMM model averaging indicate strong effects of increasing CO 2 on red spruce growth from 1989 to 2014 (Figure 6). In addition, if A is considered a diffusive process, where A = (C a –C i )*g c , then the 43% increase in A of red spruce trees after 1989 is likely produced by the combination of the 13% increase in C a and the 67% increase in g c as acidic pollution declines. Our data suggesting a CO 2 enhancement of growth and A of red spruce trees since 1989 are consistent with flux tower observations of increased forest productivity from 23 forest sites in the United States and Europe from 1995 to 2011, where increases in atmospheric CO 2 were shown to stimulate GPP and NEP (Fernández‐Martínez et al., 2017).

Both SO 4 2− deposition and elevated CO 2 reduce g c experimentally (Ainsworth & Rogers, 2007; Borer, Schaberg, & DeHayes, 2005) and, thus, the combined effects of acid deposition and increasing atmospheric CO 2 are likely responsible for the 47% decline in g c of red spruce trees from 1940 to 1989. If we assume that the increase in g c after 1989 represents a complete recovery from acid deposition, then the 15% reduction in g c from 1940 to 2014 reflects only the effect of increasing CO 2 on g c of the red spruce trees (Figure S6a) and is more in line with observed reductions in g c per ppm CO 2 increase from elevated CO 2 experiments using trees (Ainsworth & Rogers, 2007). Under this assumption, our data indicate that acidic air pollution is a 3.7× greater forcing factor on g c of red spruce trees than increasing CO 2 (Figure S6a). Numerous experiments, as well as dendroisotopic studies, have shown that iWUE, the ratio of A to stomatal conductance of water, of many plants is improved as a result of a combination of increasing CO 2 stimulating A and causing partial stomatal closure (Ainsworth & Rogers, 2007; Franks et al., 2013; Keenan et al., 2013; Knauer et al., 2016), and this has been hypothesized to be an important mechanism leading to increased forest productivity. Over the 1940–2014 tree ring chronology, iWUE of red spruce trees increased 51%, or +0.68% per year, as A increased by 27% and g c declined by 15% from 1940 to 2014. If we use the Ball–Berry model to predict the theoretical change in iWUE if impacted by changes in CO 2 alone, using a constant C i /C a ratio of 0.72 (mean at 1940), then iWUE increases by 30% over the tree ring chronology (Figure S6b), tracking the increase in atmospheric CO 2 (Knauer et al., 2016). If we model iWUE using the decline in g c attributed to increasing CO 2 for red spruce trees in this study (−15%, Figure S6a) and an increase in A that is predicted from CO 2 enrichment studies using trees (Ainsworth & Rogers, 2007), scaled over the range of CO 2 during our study (9.3% increase in A from 1940 to 1989), then iWUE increases by 41% over the tree ring chronology (Figure S6b). Both of these examples, point to the much greater sensitivity of g c to acidic sulfur pollution than to CO 2 and the subsequent role of pollution in the observed increases of iWUE. In addition, several analyses have observed that increases in iWUE do not always translate into greater tree growth (Lévesque, Siegwolf, Saurer, Eilmann, & Rigling, 2014; Peñuelas, Canadell, & Ogaya, 2011; Silva, Anand, & Leithead, 2010). In our study, BAI declined 48.9% from 1940 to 1989 as iWUE increased 67.2%, from 1989 to 1995, BAI increased 63.6% as iWUE increased 0.6%, and from 1995 to 2014, BAI increased 37.4% as iWUE decreased 9.8%. Thus, our study agrees with these analyses that iWUE is a poor predictor of tree growth, especially in mesic environments, such as our red spruce forest locations in the Central Appalachian Mountains.

From 1940 to 1989, NO x emissions increased by 0.92 million metric tons/year, but declined by 0.56 million metric tons/year from 1989 to 2014 (EPA, 2015). GLMM model averaging indicated that tree growth of red spruce across 1989–2014 was negatively impacted by reductions in NO x emissions after 1989 (Figure 6). This result was surprising given that Kendall's rank correlation identified strong negative correlations between NO x emissions and BAI, ∆13C, A, and g c after 1989 (Table 2), as well as across the entire 1940–2014 chronology (Table S1). These contradictory results highlight the uncertainty associated with the effects of N deposition on tree growth since N deposition may act as a source of soil acidity and negatively impact tree growth or may act as a fertilizer source of N and positively impact tree growth (Jennings et al., 2016; Thomas, Canham, Weathers, & Goodale, 2010).

The negative correlations between tree ring δ15N and BAI, ∆13C, A, and g c after 1989 add to the difficulty of interpreting the NO x emissions results from our GLMM model averaging, as our δ15N chronology would suggest red spruce tree growth and photosynthesis increased as ecosystem N availability declined. However, even with the dramatic declines in N deposition in this region, an examination of N cycling at these three red spruce forest sites suggests that these ecosystems are not strongly N‐limited (Smith et al., 2016). Given the close relationship between tree ring δ15N and declines in national emissions of NO x and local wet deposition of NO 3 − (Figure S3), one may hypothesize that the δ15N signatures in the red spruce tree ring chronologies reflect changes in the δ15N signatures of pollutant sources, since NO 3 − produced from fossil fuel combustion is enriched in 15N relative to biogenic sources (Felix, Elliott, & Shaw, 2012). However, the evidence that tree ring δ15N reflects δ15N from fossil fuel sources is limited (Gerhart & McLauchlan, 2014), and although the source δ15N signature of N deposition over the tree ring chronology is unknown at our locations, our data do not support this hypothesis given the constant wood δ15N prior to 1989 when NO x emissions were increasing dramatically (Figure 2b, Figure S1).

An alternative hypothesis is that δ15N in plant tissues, including tree rings, reflects ecosystem N supply relative to demand, where a more open N cycle, defined as a scenario where N supply exceeds plant demand, typically results in greater rates of net nitrate production (a process that strongly fractionates against 15N) followed by the loss of isotopically light nitrate by leaching and/or denitrification (Gerhart & McLauchlan, 2014). The loss of 15N‐depleted nitrate enriches the residual pool of plant available N that, in turn, leads to elevated values of δ15N in plant tissue (Pardo, Hemond, Montoya, Fahey, & Siccama, 2002). Following this line of reasoning, the red spruce tree ring δ15N chronology between 1940 and 1989 may reflect a more open N cycle resulting from declining tree growth rates and, therefore, a reduced plant demand for N by red spruce trees. Likewise, the decrease in wood δ15N after 1989 may have resulted from a tighter, more closed N cycle, as increased tree growth created a greater demand for N during a time of reductions in N deposition (Figure 1, Figure S1b). In addition, N concentrations in wood were 32% greater after 1989 compared to before 1989, and this, along with greater tree growth, further suggests that plant N uptake and retention increased after 1989. Our δ15N data reflect observations of declining tree ring δ15N signatures across the continental United States in recent decades (Elmore, Nelson, & Craine, 2016; McLauchlan et al., 2017) and support the patterns found in the European NITREX studies where experimental removal of N deposition from throughfall led to rapid ecosystem recovery and a more closed N cycle as losses of gaseous N through denitrification and NO 3 − leaching declined (Bredemeier et al., 1998; Corre & Lamersdorf, 2004). Thus, the δ15N chronology found in these red spruce trees has broad implications for regional forest ecosystem recovery to pollution, as well as increased water quality and fewer occurrences of eutrophication, since it is likely that less N is lost into stream runoff as N inputs from pollution decline, as has been shown in this region by Eshleman, Sabo, and Kline (2013).

Changes in historical climate were the least influential environmental variables affecting BAI, ∆13C, and seasonally integrated A and g c of red spruce trees after 1989, possibly due to the high interannual variability in temperature and precipitation in the Central Appalachian region (WV Climate Division 4, NOAA, 2017). Precipitation at our study locations did not significantly change over the study period. Precipitation was not different before and after 1989 (Figure S2b,d), and thus, was not identified by either Kendall's rank correlation or GLMMs as having an influence on BAI, ∆13C, A, or g c from 1989 to 2014. Our results are contrary to Levesque, Andreu‐Hayles, and Pederson (2017), who concluded that increasing water availability after 1984, and not increases in CO 2 or reductions in acid pollution, has been the primary factor driving increases in BAI, A, and g c of two broadleaf species, Quercus rubra and Liriodendron tulipifera, at Black Rock, NY, where annual precipitation increased 4% and growing season precipitation increased 19% during 1984–2014 compared to 1950–1983 (PRISM Climate Group, 2004).

On the other hand, mean annual temperatures in this Central Appalachian region have significantly increased by 0.50°C over the past 75 years, which was largely contributed to by a 0.72°C increase in mean April temperatures during 1989–2014 compared to 1940–1989 (WV Climate Division 4, NOAA, 2017). Our tree ring data point to changes in early spring phenology of red spruce trees due to these warmer temperatures after 1989 that are manifested as greater tree growth. First, Kendall's rank correlation identified a strong positive correlation between mean April temperatures and BAI, A, and g c of red spruce trees during 1989–2014. Second, the GLMM average model indicated a main effect of April temperatures on BAI of red spruce during 1989–2014 that contributed to tree growth increases, and showed a very high sensitivity, suggesting that even small changes in spring temperatures can affect tree growth. Finally, ANCOVA indicated that greater mean April temperatures after 1989 translated to greater seasonally integrated A and, a subtle, but significant, increase in BAI of red spruce trees compared to 1940–1989 (Figure 4). Our observations of an accelerated early spring phenology of red spruce trees due to higher spring temperatures are similar to observations of other tree species globally (Black et al., 2000; Piao, Friedlingstein, Ciais, Viovy, & Demarty, 2007), as well as in the eastern United States (Dragoni et al., 2011; Elmore et al., 2016; Richardson et al., 2009).