Tree-ring δ 13C tracks flux tower ecosystem productivity estimates in a NE temperate forest

We investigated relationships between tree-ring δ13C and growth, and flux tower estimates of gross primary productivity (GPP) at Harvard Forest from 1992 to 2010. Seasonal variations of derived photosynthetic isotope discrimination (Δ13C) and leaf intercellular CO2 concentration (c i ) showed significant increasing trends for the dominant deciduous and coniferous species. Δ13C was positively correlated to growing-season GPP and is primarily controlled by precipitation and soil moisture indicating that site conditions maintained high stomatal conductance under increasing atmospheric CO2 levels. Increasing Δ13C over the 1992–2010 period is attributed to increasing annual and summer water availability identified at Harvard Forest and across the region. Higher Δ13C is coincident with an enhancement in growth and ecosystem-level net carbon uptake. This work suggests that tree-ring δ13C could serve as a measure of forest GPP and be used to improve the calibration and predictive skill of ecosystem and carbon cycle models.


Introduction
Our ability to project the future carbon cycle is limited by our lack of understanding of terrestrial carbon (C) cycle dynamics and the feedbacks that constrain C budgets. Projections of Cflux and C-sequestration from current coupled terrestrial carbon cycle models give widely divergent results (Friedlingstein et al 2006). The lack of agreement among projections is, in part, related to poorly constrained model parameters.
Constraining carbon budget projections is an important issue given the current mitigation of anthropogenic emissions by terrestrial ecosystems sequestering C (Le Quéré et al 2009). Terrestrial ecosystems sequestered approximately 30% of anthropogenic emissions from 2000 to 2006, and there is interest in the potential to manage forest ecosystems to increase the strength of the carbon sink (Hurtt et al 2002, Canadell et al 2007, Woodbury et al 2007. Longer periods of observational data on C-cycle variability, in particular, would help constrain model parameters and improve model performance. CO 2 flux measurements between ecosystems and the atmosphere from eddy covariance flux towers (Baldocchi 2003) are valuable for model evaluation (Suzuki and Ichii 2010, Richardson et al 2012, Schaefer et al 2012, Keenan Environmental Research Letters Environ. Res. Lett. 9 (2014 074011 (9pp) doi:10.1088/1748-9326/9/7/074011 Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. et al 2012a, Raczka et al 2013), but they only measure CO 2 flux conditions over short periods (maximum 20 years, most ⩽10 years) in a small footprint (1 km 2 ) and the regional tower density is low (e.g., five sites in northeastern US forests with >10 years of data). Consequently, a variety of other types of measurements are needed to constrain C-cycle model parameters and to improve projections of interannual variability in C-cycle dynamics in terrestrial ecosystems (Pan et al 2006, Sims et al 2006, Randerson et al 2009, Babst et al 2013. Annual measurements of forest carbon production can be obtained from tree rings that integrate the influence of climate and site factors on forest growth (Babst et al 2013). Attempts to associate tree-ring width with net ecosystem exchange (NEE) and gross primary productivity (GPP) measurements from flux towers have proven inconclusive (Rocha et al 2006). Integrals of annual growth over a few to several years have yielded similar net carbon fluxes as continuous flux tower measurements, but annual increment and flux tower measurements do not show strong inter-annual correlations (Barford et al 2001).
The composition of carbon isotopes (δ 13 C) of an annual tree ring records the proportions of assimilated carbon and is closely linked to photosynthetic capacity and stomatal regulation during tree growth (Farquhar et al 1989, Ogée et al 2009. Therefore, the distribution of stable carbon isotopes within a tree ring reflects the carbon assimilation and environmental conditions experienced by a tree through the growing season. The δ 13 C is widely used to study the linkage between environmental conditions and the physiological processes that control tree growth at inter-and intra-annual time resolutions (Loader et al 1995). Tree rings, particularly wood cellulose δ 13 C, contain valuable records of past climate, leaf-gas exchange, and carbon allocation within trees (e.g. Barbour et al 2002). Several studies used tree-ring δ 13 C to reconstruct the intrinsic water use efficiency (iWUE) of trees (e.g. Duquesnay et al 1998), and δ 13 C within a tree ring has been shown to track the physiological changes at the canopy level as recorded by the flux towers (Walcroft et al 1997, Michelot et al 2011. This approach has shown the potential to better understand the link between canopy-level physiology, tree-ring isotopic signature and climate drivers (Offermann et al 2011).
In this paper, we identify the relationships between interannual variability in tree-ring δ 13 C, tree-ring increment, flux tower measurements of CO 2 , and climate for a temperate forest in the northeastern United States. The period of analysis spans 18 years (1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010) and is from the longest available record of continuous flux tower measurements (Harvard Forest, Petersham, Massachusetts;Urbanski et al 2007). We analyzed the δ 13 C (of α-cellulose) from tree-ring latewood (LW) because it is less dependent on carbon stored during the previous growth year and therefore corresponds to the isotopic signal of recently assimilated C during the growing season. We conducted this analysis for one broad-leaf deciduous and one evergreen needle-leaf co-dominant species in the flux-tower footprint to contrast species response. Specifically, we: (1) identified the relationship between tree-ring δ 13 C, GPP, and tree growth, and, (2) assessed the responses of tree physiology and growth to rising atmospheric CO 2 mole fraction and climate variation. We also discuss how changes in climate conditions and atmospheric CO 2 concentration over the period of record may have modulated growth and forest productivity over the last 18 years.

Study site and tree-ring chronologies
The forest within the flux-tower footprint (1 km 2 ; Harvard Forest-EMS tower) is mainly deciduous and is dominated by Quercus. rubra (L.) (northern red oak; 52% of basal area), Acer rubrum L. (red maple; 22% of basal area) and Tsuga canadensis (L.) (eastern hemlock; 17% of basal area). Increment core sampling focused on canopy dominant trees. To identify interannual variability in tree growth and tree-ring δ 13 C, we cored Q. rubra and T. canadensis with a 5 mm increment borer at 1.37 m stem height (table 1). Two cores were taken from each tree to build ring-width chronologies. Additionally, we cored five and four trees respectively of each species at the same height with a 12 mm increment borer. The large cores were used to analyze the δ 13 C in the LW of the tree rings. A sample of four trees is sufficient to achieve good precision of the δ 13 C sample mean (McCarroll andLoader 2004, Leavitt 2008).
Ring-width chronologies were developed for each species using standard dendrochronological techniques (Speer 2010). For the site, Q. rubra mean age was 97 years (range 71-115 years) and T. canadensis mean age was 145 years (range 89-221 years). For isotope cores, Q. rubra mean age was 103 years (range 87-113 years) and T. canadensis mean age was 142 years (range 93-188 years). All ages were estimated from inner ring dates.
Basal area increment (BAI) was calculated for each tree and species using the dplR package in R (Bunn 2008). We used the 'outside in' function to convert raw ring-width measurements to BAI based on the diameter of the tree and the width of each ring moving towards the pith of the tree. The method assumes a circular growth pattern. BAI was used instead of ring width as a surrogate of radial growth and carbon gain because it represents more accurately tree annual biomass increment without the need for standardization (Biondi and Qeadan 2008). A mean BAI value was computed for each species at an annual resolution, averaged over all the trees cored for both carbon isotope analyses and ring-width chronologies (i.e., 16 trees/42 cores for Q. rubra and 22 trees/ 55 cores for T. canadensis).

Stable istope analyses
We analyzed the LW portion of each tree ring from individual trees over the 1992-2010 period (table 1). The samples were milled using an ultra-centrifugation mill (Qiagen TissueLy-serII) and α-cellulose was extracted from each wood sample following the Soxhlet method elaborated by Green (1963) and modified by Leavitt and Danzer (1993). δ 13 C ratios were measured on the CO 2 produced by α-cellulose combustion in a Costech elemental analyzer coupled with a Thermo Delta V-IRMS. The sample accuracy was determined to be ±0.08‰ (1 σ standard deviation calculated from the average difference between measured and true internal standard, n = 13). The isotopic value is expressed in the delta (δ) notation relative to the VPDB (‰ VPDB).

Calculation of Δ, c i and iWUE
We used δ 13 C to determine carbon isotope discrimination (Δ) by the plant against atmospheric δ 13 C, and variation in plant iWUE. The Δ describes the isotopic difference between the δ 13 C of air (δ 13 C air ) and the plant (δ 13 C plant ) and results from the preferential use of 12 C over 13 C during photosynthesis. Δ is calculated using: where a is the fractionation during CO 2 diffusion through the stomata (4.4‰: O'Leary 1981); b is the fractionation by RuBP carboxylase (27‰: Farquhar and Richards 1984); and c i and c a , are the leaf intercellular space and ambient CO 2 concentrations (μmol Mol −1 ), respectively. We used the data of atmospheric CO 2 concentration (c a ) measured at 29 m height on the eddy-covariance tower to calculate c i from equation (2). The fractionations due to diffusion and carboxylation are constant but additive; therefore, the δ 13 C records variations in c i as regulated by two main processes: stomatal conductance (c a /c i ) and photosynthetic assimilation rate (A). The iWUE is the ratio of the net photosynthetic assimilation rate (A) and water vapor conductance (g H2O ) and is described by Ehleringer and Cerling (1995) as: The iWUE derived from plant isotope data is used to compare photosynthetic properties independent from evaporative demand (Osmond et al 1980), and is therefore, often applied as an indicator of long-term trends in the internal regulation of carbon uptake and water loss of plants (Seibt et al 2008). We used summer (June, July and August) values of δ 13 C air and c a to calculate Δ, and to reconstruct c i and iWUE for the period corresponding to LW formation.

Meteorological and CO 2 flux data and data analyses
We used hourly gap-filled meteorological data from the eddycovariance tower measured by sensors above the canopy at 29 m height (Urbanski et al 2007, Harvard Forest Data Archive HF004) to compute monthly mean temperature, vapor pressure deficit (VPD), incident photosynthetically active radiation, relative humidity as well as monthly total precipitation from 1992 to 2010. Hourly gap-filled GPP was used to provide monthly cumulative GPP. We used the Palmer Drought Severity Index (PDSI) for the central Massachusetts climate region from the National Climatic Data Center.
The association between climate variables and tree-ring Δ chronologies were assessed using correlation analyses. The GPP time series during the growing season (from May to October) were compared to tree-ring Δ chronologies for each species. Correlations were calculated between monthly GPP and LW-Δ using the Pearson (pairwise) product-moment correlation. Trends over the 1992-2010 period were estimated using linear regression.

Results
We assessed the number of trees required to extract a common reliable climate signal (thus, the degree to which isotopic composition series vary in parallel) using the expressed population signal (EPS, Robertson et al 1997, McCarroll andPawellek 1998). The EPS values calculated using the 1992-2010 isotopic compositions are ⩾0.85 for two trees from each species. An EPS value of ⩾0.85 suggests that the sample size is adequate (Wigley et al 1984). At the high frequency investigated in this study (year to year), the δ 13 C series display a high common variance within each species (table 1). We assessed the confidence interval (CI) around annual mean values to account for the differences in the absolute isotopic values (McCarroll and Loader 2004). The sample size of 4 T. canadensis and 5 Q. rubra provided a CI 95% = 0.3‰ and 0.18‰ respectively. The mean values of 1992-2010 individual series were calculated for all tree cores to produce one main series for each species (Shi et al 2011) hereafter called δ 13 C TSCA for T. canadensis and δ 13 C QURU for Q. rubra.
δ 13 C measurements were significantly different between species (t = 2.02, P < 0.05). The mean δ 13 C TSCA display higher ratios than δ 13 C QURU , this difference is probably linked to hydraulic conductivity in conifers which is less efficient than in ring-porous deciduous species (Stuiver andBraziunas 1987, McCulloh et al 2010) or other factors such as differences in the maximum photosynthetic capacity or leaf morphology (Barbour et al 2002). The species specific time series of annual δ 13 C variations are positively correlated with each other (r = 0.61 P < 0.001). The derived Δ and c i series for each species reveal similar and significant increasing trends over time, although the increase was more pronounced for T. canadensis (figure 1). The iWUE over the period of record remained stable which translates into a constant c a − c i or an increasing c i /c a . A few years did display a higher iWUE in comparison to most years; in particular, iWUE spikes in 1999 in both species records.
The mean BAI also showed a significant increasing trend over the period of record for Q. rubra (R 2 = 0.64, P = 0.0001) but the trend was less pronounced for T. canadensis (R 2 = 0.32, P < 0.05). Overall, the trends analysis using linear regression shows that Δ, c i and BAI have increased for each species. The changes were significant (p < 0.05) and the slopes were positive (table 2).
The Δ was strongly correlated with July (0.55, P < 0.05) and October (0.53, P < 0.05) GPP for T. canadensis, and with July (0.49, P < 0.05), August (0.58, P < 0.01), and October (0.55, P < 0.05) GPP for Q. rubra. Online supplementary table S1 (available at stacks.iop.org/ERL/9/074011/mmedia) summarizes the Pearson correlation obtained between monthly climate variables and Δ QURU and Δ TSCA. Overall, the strongest positive correlations were found with PDSI from May to August for T. canadensis and for May for Q. rubra. A positive correlation was found between Δ of both species and early-growing season precipitation (April and May).

Relationship between Δ, growth, and flux tower estimates of GPP
The carbon contribution to stem growth during LW formation should primarily originate from the immediate product of The Δ-GPP correlations are consistent with the seasonality of carbon uptake from eddy flux measurements and C storage estimates at the deciduous forest-EMS tower site and the Hemlock Forest tower site located 0.5 km from the EMS tower (Hadley and Schedlbauer 2002, Hadley et al 2009). Net C uptake at the deciduous forest site occurs between the end of May and mid-October. At the hemlock site, the peak in net C uptake and storage occurs in April and May. This period is not influencing the δ 13 C in the LW, although it shows the effect of conifers on the annual pattern of the C exchange. Conifers currently represent 17% of the total basal area at the deciduous forest site. After C uptake declines in June, it increases in July, and falls in August before a second uptake peak in October, consistent with the Δ-GPP correlations for T. canadensis.
The species seasonality of carbon allocation is further reflected in the relationship between LW-Δ and growth. The correlation for Q. rubra LW-Δ and mean BAI was significant (r = 0.72, P < 0.001), and was slightly higher with the mean Q. rubra BAI of all sampled trees (r = 0.75, P < 0.001). To illustrate the relationship between aboveground carbon gain and photosynthetic discrimination we regressed BAI against Δ, and a strong positive trend was obtained for Q. rubra (R 2 = 0.39 P < 0.01). The growth increase of Q. rubra could be interpreted as the result of a greater photosynthetic rate induced by increased c i concentrations (von Caemmerer and Farquhar 1981, Schubert and Jahren 2012), and a substantial allocation of carbohydrates to LW formation and radial growth (Keel et al 2006, Palacio et al 2011. However stage of stand development may also be contributing to this increase in BAI (e.g. Foster et al 2014, figure S1). A relationship between Δ and BAI was not found for T. canadensis, probably because carbon sequestered after June contributes little to the current year radial growth or that carbon allocation differs by species (Richardson et al 2013). The trend in the BAI growth enhancement for Q. rubra is consistent with changes in the measured annual increment of aboveground biomass in the flux-tower footprint from 1993-2010 (∼20% increase, Data Archive: HF069), and provides a quantification of carbon stored in aboveground woody biomass. Interestingly, the aboveground biomass over the last decade accounted for ∼50% of the total carbon sequestered (Keenan et al 2012a) confirming the relationship we found between Q. rubra BAI and Δ, and growing season GPP (figure 2). The remaining carbon uptake could be attributed to litter or soil pools; however, these carbon proportions and stocks are not accounted for in the tree-rings δ 13 C.
The inter-annual variability and fraction of carbon uptake and sequestration explained by the tree-ring derived Δ (figures 2(c), (d)) provides a quantitative proxy that can be used to constrain models that estimate forest productivity. Process-based models fail to accurately reproduce the observed inter-annual variability of C-fluxes (Keenan et al 2012b) mainly because of inaccurate model allocation structure and lagged effects of climate variability on tree growth and physiology (Gough et al 2009). In contrast, the Δ integrates biotic factors which modulate the δ 13 C fractionation and the isotopic composition of the total pool of carbon fixed during photosynthesis at seasonal and annual cycles (Leavitt 1993).

Climate drivers and physiological implication of the Δ trend
We examined the climatic drivers of Δ using meteorological data from the eddy-covariance tower and PDSI. The climate signal within both species Δ series is dominated by monthly precipitation and PDSI during the growing season (table S1), although they have a stronger influence on T. canadensis due to its particularly drought-sensitive nature caused by a relatively shallow rooting-depth (Cook and Cole 1991).
The carbon isotope discrimination properties can be used to assess how trees are responding to increasing atmospheric concentration c a and changes in climate. Three physiological response scenarios described by Saurer et al (2004) identify possible changes in iWUE and c i following c a increase: (1) c i remains constant such that c i /c a decreases and iWUE increases; (2) c i increases proportional to c a such that c i /c a remains constant and iWUE increases; (3) c i increases at the same rate as c a, c i /c a increases and iWUE remains constant. In this study, Δ and c i increased and there was no increase in iWUE for both species ( figure 2, table 2). The c i /c a calculated for Table 2. Trend statistic, per species, for isotopic discrimination (Δ), intercellular CO 2 concentration (c i ), intrinsic water use efficiency (iWUE) and basal area increment (BAI). All units are given per year (yr). The slopes are estimated using linear regression.
both species increases indicating c i follows the increase in c a at a rate of 2 ± 0.21 ppm yr -1 for T. canadensis, and 1.7 ± 0.17 ppm yr -1 for Q. rubra. The c i /c a increase and no change in iWUE suggest a weak stomatal response to c a increase and is described as a passive response to changes in c a (McCarroll et al 2009) where neither stomatal conductance nor photosynthetic rate change. The most common response documented in trees is an active response as stomatal conductance is reduced (Bert et al 1997, Saurer et al 2004, Peñuelas et al 2011. The iWUE trends derived from the treerings δ 13 C in our study show a divergence from the instantaneous WUE (the ratio of carbon assimilation to transpiration, Farquhar and Richards 1984) trend documented at the ecosystem level using continuous measurements of CO 2 and water vapor fluxes at HF (Keenan et al 2013). The canopyintegrated water use efficiency calculated as the ratio between gross ecosystem photosynthesis (GEP) and Ecosystem evapotranspiration (E) shows a significant increase over the measurements period (3.5% change yr −1 , P = 0.01). The latter was calculated taking into account the atmospheric evaporative demand and therefore is referred to as 'inherent' water use efficiency (W ei ). The 'intrinsic' water use efficiency as derived from tree-rings δ 13 C is often used as an indicator of long-term trends in the internal regulation of carbon uptake and water loss independently from evaporative demand (Osmond et al 1980). Therefore iWUE inferred from δ 13 C may not be representative of instantaneous or inherent WUE. Although both are affected by similar processes (i.e. stomatal conductance), they can vary independently because they are influenced by additional factors like mesophyll conductance, leaf N-content and C-respiratory losses (Griffiths et al 1999, Seibt et al 2008. We note that when iWUE is calculated from the flux tower measurements (as the ratio between GEP and canopy water conductance), the magnitude of the trend is lower and the statistical strength is no longer significant compared to those calculated for W ei (1% change yr −1 , P = 0.2, Keenan et al 2013).
The Δ trend found in this study during the last 18 years is unusual but not unique. Results from literature often reported a decrease in Δ inferred from tree-ring records, and thus a strong improvement of iWUE during the last 100 years under increased atmospheric CO 2 (Arneth et al 2002). However, Marshall and Monserud (1996) showed that iWUE has remained static in sites from western US and a switch in iWUE from an active to a passive response around AD 1970 has been noted in several European sites (Gagen et al 2011). The spatial variability of stomatal regulation and iWUE response to increasing c a is explained by the strong dependence of water and carbon usage on moisture stress experienced at a specific site and species sensitivity to moisture availability (Warren et al 2001). The increased discrimination is a reflection of the availability of CO 2 . When the level of CO 2 (c i ) is high, large discrimination is observed reflecting the RuBisCO enzymatic preference for 12 CO 2 (Stewart et al 1995). Clearly, the Δ chronologies at HF indicate that despite higher c a concentrations, factors like physiology and site condition contributed to stomatal openness and conductance.
The strong correlation between Δ and growing season PDSI and precipitation for both species suggest a link between water availability and stomatal conductance (Dupouey et al 1993). The climatic data for central Massachusetts indicates that precipitation and moisture (PDSI) are high at HF and are rarely limiting tree growth (Voelker 2011). In fact, the region became wetter over the period of analysis. The mean annual precipitation from 1898 to 2010 was 107.45 cm and increased 0.28 ± 0.06 cm yr -1 (p < 0.0001). The trend in precipitation is also reflected in a trend of increasing PDSI of 0.03 ± 0.005 units yr -1 (P < 0.0001) over the same period. Since 1992, precipitation and PDSI have increased at rates of 1.30 ± 0.67 cm yr -1 (P < 0.07) and 0.15 ± 0.05 units yr -1 (P < 0.006), respectively. The increase in regional scale water availability (Wang et al 2013) is also evident at the HF site. Annual precipitation and growing season precipitation have increased by 10% and 15% over the last 18 years, respectively (Keenan et al 2013).
We conclude on the basis of the photosynthetic discrimination sensitivity to water availability that the long-term increase in moisture which controls Δ drove the increase in stomatal conductance and carbon photosynthetic discrimination. A positive correlation between Δ and mean annual and summer precipitation has been confirmed for modern leaf tissues from multiple biomes in North America (Diefendorf et al 2010). This correlation is explained by the effect of the mean annual and summer precipitation levels on both soil water status and VPD (Hartman and Danin 2010). It is noteworthy that despite the mesic conditions of the HF site, precipitation and soil moisture levels remain very strong predictors of CO 2 discrimination. Similar to our findings, drought sensitivity has been demonstrated for radial growth of T. canadensis and Q. rubra growing in mesic sites located in the New York City watershed (Pederson et al 2013).
The sensitivity to water availability is further illustrated by iWUE data for the year 1999. Both species exhibited a higher iWUE reflecting lower Δ values (figure 1). The summer of 1999 was very dry with decreased soil moisture levels resulting in the restriction of the CO 2 supply to the leaf (Brugnoli et al 1988) because of stomatal closure, and in a weaker carbon uptake observed in the summertime NEE record (Urbanski et al 2007).

Conclusion
An enhancement in photosynthetic isotopic 13 C discrimination and c i was observed over the last 18 years for two codominant species at the Harvard Forest. The combined effect of increased water availability and higher atmospheric CO 2 concentration lead to increased plant CO 2 assimilation. The higher net photosynthesis can in part explain the growth enhancement of Q. rubra and above-ground GPP (e.g. Linares and Camarero 2012) and is consistent with the observed higher CO 2 uptake and C storage documented in the flux-tower measurements over the same period. The tree-ring Δ-GPP relationship identified in this study can be used as a quantitative proxy to reconstruct and interpret past forest productivity, as driven by climate variability and in response to long-term atmospheric CO 2 increase.
The photosynthetic carbon isotopic discrimination model is embedded in most of the models of forest carbon cycling (Richardson et al 2012) that are coupled with earth-system models to project terrestrial carbon cycle and feedbacks to climate change (Sitch et al 2008). We show that the data provided by tree-ring δ 13 C, recording the seasonal and interannual patterns of plant carbon assimilation, can be used to inform forest ecosystem model parameterization (Bodin et al 2013), to improve simulations of monthly scale ecosystem-function (e.g. Medvigy et al 2013), and to constrain longer term terrestrial carbon dynamics in response to climate fluctuations and forest management (e.g. Brooks and Mitchell 2011).
Whether the documented isotopic trends in the present study are merely temporary or a site specific phenomenon remains to be seen. Further investigations of tree-ring δ 13 C and ecosystem physiology measurements from other fluxtower sites in northeastern US forests need to be explored to assess physiological responses at the regional scale.