Significance The metabolism of North America’s oldest boreal trees (Thuja occidentalis L.) is strongly affected by rising anthropogenic C O 2 emissions. Intrinsic water use efficiency ( i W U E ) increased dramatically, although nonlinearly, since the beginning of the industrial Era. Our study shows that while T. occidentalis L. acclimated to rising [ C O 2 ], no stem growth was observed, suggesting that trees could not benefit from the increased i W U E .

Abstract Due to anthropogenic emissions and changes in land use, trees are now exposed to atmospheric levels of [ C O 2 ] that are unprecedented for 650,000 y [Lüthi et al. (2008) Nature 453:379–382] (thousands of tree generations). Trees are expected to acclimate by modulating leaf–gas exchanges and alter water use efficiency which may result in forest productivity changes. Here, we present evidence of one of the strongest, nonlinear, and unequivocal postindustrial increases in intrinsic water use efficiency ( i W U E ) ever documented (+59%). A dual-isotope tree-ring analysis ( δ 13 C and δ 18 O ) covering 715 y of growth of North America’s oldest boreal trees (Thuja occidentalis L.) revealed an unprecedented increase in i W U E that was directly linked to elevated assimilation rates of C O 2 (A). However, limited nutrient availability, changes in carbon allocation strategies, and changes in stomatal density may have offset stem growth benefits awarded by the increased i W U E . Our results demonstrate that even in scenarios where a positive C O 2 fertilization effect is observed, other mechanisms may prevent trees from assimilating and storing supplementary anthropogenic emissions as above-ground biomass. In such cases, the sink capacity of forests in response to changing atmospheric conditions might be overestimated.

Plants are an integral component of both water and carbon cycles on Earth. This is because stomata (i.e., microscopic apertures on plant leaves) allow plants to adjust their metabolism to changing environmental conditions by controlling the tradeoff between C O 2 intake for photosynthesis and water loss through transpiration. However, since the beginning of the industrial period, anthropogenic emissions and changes in land use resulted in a 40% increase in atmospheric [ C O 2 ] (1) (hereafter referred to as c a ). Plants have not been exposed to concentrations above 290 ppm for at least 650,000 y (2, 3), so the dynamic carbon–water coupling that prevailed for several millennia could be disrupted with potential, yet uncertain, impacts on hydro-ecosystem functioning (4, 5).

One of the most striking consequences of elevated c a is an increase in plant intrinsic water use efficiency ( i W U E , i.e., the rate of carbon uptake per unit of water lost). The rate of increase of i W U E in response to rising c a depends on how the assimilation rate (A) and/or stomatal conductance ( g s ) adjusts to determine the optimal internal [ C O 2 ] ( c i ) (5, 6). Three theoretical scenarios (hereafter referred to as S1–S3) have been proposed (7) as guidelines on how to interpret possible acclimation strategies and their respective effects on i W U E (Fig. 1). In S1, plants maintain a constant c i ( i W U E increases strongly); in S2, there is a constant c i / c a ratio ( i W U E increases moderately); and in S3, there is a constant difference between c a and c i ( i W U E remains constant). Short-term flux measurements (4) and free-air C O 2 enrichment (FACE) experiments (8) have shown that S1 requires a strong and active physiological response leading to the highest rates of i W U E increase, most consistent with a strong C O 2 fertilization effect. However, the active i W U E response documented over short periods of time (yearly to decadal) may not reflect the changes occurring on longer timescales characteristic of natural forests (decadal to centennial). Over these timescales, there is no clear line of evidence that points toward a single, timescale-independent, acclimation strategy. To the contrary, the emerging concept of “optimal stomatal behavior” implies that woody plants shift along a continuum of strategies (5) and adjust g s to maximize carbon gains while minimizing exposure to drought stress. This complex interplay takes place on yearly to multidecadal timescales and therefore requires the investigation of time series that are sufficiently long to be able to track and compare present-day stomatal behavior and i W U E changes to preindustrial reference levels.

Fig. 1. Evolution of i W U E , c a , RWI, and δ 18 O since 1300 CE for old-growth white cedars at lake Duparquet, Eastern Canada. Two distinct periods emerge which have been highlighted through different coloring and shading: P1 (1850–1965, red) and P2 (1965–2014, green). (A) i W U E is calculated from δ 13 C -derived physiological parameters shown in Fig. 2. The gray shading around the i W U E curve represents the intertree variability. The three theoretical i W U E scenarios presented in the Introduction (S1, constant c i ; S2, constant c i / c a ; and S3, constant c a − c i ) are plotted as dotted lines. (B) Annually resolved RWI chronology, RCS standardized with bootstrap 95% CI (gray shading). (C) δ 18 O values.

Moreover, it is still unclear how these shifts in stomatal behavior impact forest productivity, i.e., whether changes in the rate of increase in i W U E are mirrored by changes in growth rates. Based on FACE studies (8), one could hypothesize that increases in i W U E (e.g., S1 or S2) would enhance stem growth. However, at least in boreal North America, tree-ring–based and bioclimatic modeling studies found little evidence of a widespread growth of forest stands that could be attributed to c a increases (9, 10). Conversely, in this zone, growth has declined over vast portions of land, which has been attributed to the negative effects of temperature-induced droughts (10, 11). However, few of these tree-ring studies were accompanied by a thorough and time-continuous investigation of changes in gas exchange in response to increasing c a . Only one study measured changes in i W U E and found that the positive response to c a could not compensate for growth declines observed at the southernmost fringe of the boreal forest (12). Nevertheless, a high level of uncertainty persists with regard to the linkages between climate variability, i W U E , and forest productivity and their changes since preindustrial times.

To remedy this uncertainty, we report on the acclimation of a population of long-lived white cedars (Thuja occidentalis L.) growing in xeric conditions that offers a unique, long-term perspective on i W U E –growth relationships. We present a long (1300–2014), annually resolved, multiproxy tree-ring chronology in eastern North America (ring widths and stable isotopes of carbon and oxygen), to track changes in i W U E and growth attributable to rising c a and climate variability for a period that extends far beyond the onset of the Industrial Era.

Conclusions North America’s oldest boreal trees (T. occidentalis L.) exhibited an unprecedented 59% increase in i W U E since the beginning of the Industrial Era, mainly in response to rising atmospheric [ C O 2 ] ( c a ). This response was not constant in time, but a clear regime shift occurred around 1965, gradually switching the leaf gas-exchange strategy from a constant c i to a constant c i / c a acclimation. Trees maximized carbon gains between 1850 and 1965 (P1) when the local climate was relatively wet and cold. Conversely, they reduced vulnerability to drought stress when the local climate became hotter and drier after 1965 (P2), possibly also combined with a closer saturation of the photosynthetic apparatus. However, the marked i W U E increase since 1850 did not translate into enhanced stem growth. This implies that other physiological processes such as nutrient availability, changes in carbon allocation strategies, or changes in stomatal density may have offset the growth benefits awarded by the increase in i W U E . Our results thus suggest that even in favorable conditions for growth, not all trees can take advantage of elevated levels of c a and i W U E . These mechanisms are not typically taken into account by ecophysiological models (45) and DGVMs (44) which may lead to them overestimating any positive effects bestowed by a higher c a on carbon assimilation and fixation. Predictions of elevated future growth and possible alleviating effects on c a might be overly optimistic if the models do not allow for the possibility of constant or even reduced growth in the context of increasing c a . It would be advisable to define improved allocation rules and determine whether those rules can change over time (and with increasing c a ). Integrating measurements acquired at shorter time intervals (daily or weekly resolution) with long times series of tree-ring measurements would further help to track and define those allocation rules in trees exposed to a varying climate and nutrient availability. SI Appendix. SI Appendix provides complementary information on methods and results: the regional curve used for regional curve standardization (RCS), quality statistics for the RCS chronology, cohorts selected for δ 13 C and δ 18 O analysis and their junction, effect of correction methods on δ 13 C chronology, and results for whole-wood and cellulose δ 13 C and δ 18 O values comparison.

Materials and Methods Study Site and Sampling Strategy. Lake Duparquet is located at the southern limit of the boreal forest (48, 4 7 ° N, 79, 2 7 ° W), in northeastern Canada (SI Appendix, Fig. S1). Climate is continental, cold, and humid, with a mean annual temperature of 1 ° C and a mean annual total precipitation of 985.2 mm (climate normals for 1981–2010 from Station Mont-Brun, located 12 km south and 42 km east from Duparquet) (46). In 1987, 39 trees were sampled by ref. 40. To extend their ring-width chronology (1186–1987) until 2014, 35 trees, spread across five islands and six peninsulas of Lake Duparquet (SI Appendix, Fig. S1), were resampled in 2012 and 2014 using a 5-mm core increment borer at chest height. Sixteen trees were also part of the first chronology by ref. 40. Trees were measured, cross-dated, and detrended using the RCS approach (47) (SI Appendix). Stable Isotope Chronologies. Eleven trees were selected for carbon and oxygen stable isotope ratio analysis to cover the maximum time period with the smallest number of possible cohort changes (SI Appendix, Fig. S4). This avoids loss of a reliable mean as we go back in time, considering the systematic offsets typically encountered in the stable isotope analyses of different trees (25, 48, 49). Isotope chronologies have been shown to require few replicates (as little as four) to obtain a strong common signal (expressed population signal ≥0.85) (25). Two cohorts of six (cohort 1: 1620–2014) and five trees (cohort 2: 1295–1645) were selected, with an overlap of 25 y (SI Appendix, Fig. S4). The trees selected for stable isotope analysis had different ages that were well spread through time. This avoids incorporating noise due to nonclimatic age effects (SI Appendix, Fig. S4). The first 50 y of each sample were excluded to avoid any nonclimatic age-related trends (49, 50) (SI Appendix). δ 13 C -Derived Physiological Parameters. δ 13 C variability in tree rings depends on processes by which trees fractionate carbon at the leaf gas-exchange site (discrimination against 13C vs. 12C) before fixing it in the stem during the growing season (25, 51). Carbon fractionation (Δ, ‰) is expressed as the difference between δ 13 C values of atmospheric C O 2 ( δ 13 C a i r ) and tree-ring δ 13 C values ( δ 13 C t r e e ) as shown by Eq. 1 (7): Δ ( ‰ ) = δ 13 C a i r − δ 13 C t r e e 1 + δ 13 C t r e e / 1,000 ≈ δ 13 C a i r − δ 13 C t r e e . [1]Carbon fractionation thus occurs in two main steps (52): C O 2 diffusion through stomata (step i ≈ −4.4‰; ref. 53) and carboxylation (step ii ≈ −27‰; refs. 25 and 54). However, the overall fractionation is dependent on the ratio between leaf intercellular ( c i ) and ambient ( c a ) C O 2 concentrations: Δ ≈ δ 13 C a i r − δ 13 C t r e e = a + ( b − a ) c i c a . [2] Records of δ 13 C a i r and c a obtained from Antarctic ice cores by ref. 55 were linearly interpolated until 1295. From 1850 to 2003, values from ref. 25 were used (following refs. 2, 4, and 55) and linearly interpolated until 2014. Based on plant physiology theory and models, physiological parameters can be derived from tree-ring δ 13 C values. c i is calculated from Eq. 2, provided that all other variables are known. i W U E , a measure of the amount of carbon assimilated per unit leaf area per unit time per unit cost of water, is calculated as (18) i W U E = A g s = c a − c i 1.6 . [3]According to Eq. 3, the difference between c a and c i is proportional to i W U E and reflects the balance between the rate of C O 2 assimilation (A) and stomatal conductance ( g s ). For each of the three theoretical scenarios of leaf–gas strategy (7), the δ 13 C response to c a is shown in Fig. 2 (main text). This was achieved by fixing c i accordingly, using Eq. 3. Scenarios start in 1850, at the beginning of the industrial period, and exhibit a linear response to c a increases and assume a negligible c i response to climate variability. The c i preindustrial level was 150 ppm (the 1295–1850 mean). Dendroclimatic Analyses. Daily weather data (1901–2014 maximum and minimum temperatures) and precipitation (1901–2014 sum) and VPD (1953–2014) were retrieved from BioSIM v.10, which interpolates site-specific estimates from historical weather observations (56) as described in ref. 57. The quantity of available soil moisture was estimated for each month using the quadratic + linear (QL) formulation procedure described in ref. 58, which accounts for water loss through evapotranspiration (simplified Penman–Monteith potential evapotranspiration) and gain from precipitation. Parameter values for maximum and critical available soil water were set at 300 mm and 400 mm, respectively; the number of weather stations for interpolation was set to eight. We used values of annual mean c a recorded at Mauna Loa observatory since 1958 (59); the data were extended to 1901 using estimates from ice cores (60). Commonality analysis (61) was used to explore the climate effects of c a and climate predictors on i W U E , δ 18 O , and RWI. With a commonality analysis, one seeks to decompose linear regression R 2 into its unique and common effects. Unique effects indicate how much variance is uniquely accounted for by a single variable, taken as a single predictor. Common effects indicate how much variance is common to a set of predictors. Total effects refer to the combination of unique and common effects and describe the explanatory power of a predictor, summing up all commonalities. Commonality analysis is particularly useful to estimate effects in the presence of multicollinearity in a predictor dataset. For each analysis, we show regression beta weights, unique effects, common effects, and total effects for each variable, as well as R 2 values for the whole regression model. For simplicity, we chose to use a limited set of predictor variables: May–August (growing season) temperatures (T), SMIs, and VPDs. We performed the commonality analysis iteratively; i.e., the data were passed several times through a high-pass filter which led to incremental degradation of the original i W U E variable, as each pass removed a bigger part of the low frequencies. Hence, at the beginning of the loop i W U E contained all frequencies from 0 ≤ f ≤ 0.5 y−1 (cycle lengths between ∞ and 2 y) while at the end of the iteration i W U E was degraded to a perfectly detrended seesaw signal with f = 0.5 y−1 (cycle length = 2 y). The iterative degradation of the original i W U E signal coupled to the commonality analysis allowed us to pinpoint those frequencies above or below which c a and climate had a dominant influence on i W U E , δ 18 O , and RWI.

Acknowledgments We acknowledge the contribution of many people who helped with field and laboratory work: A. Barbe, G. Proulx, P. Leclerc, M. Gratton, and A. Adamowicz-Walczak. We thank the Natural Sciences and Engineering Research Council (NSERC) and Her Majesty The Queen in due right of Canada. C.G.-C. acknowledges financial support from both NSERC and Fonds de Recherche Québécois Nature et Technologies for her master’s thesis.

Footnotes Author contributions: C.G.-C., É.B., Y.B., I.D., and M.G. designed research; C.G.-C., É.B., and M.P.G. performed research; C.G.-C., É.B., M.P.G., and J.-F.H. contributed new reagents/analytic tools; C.G.-C., É.B., Y.B., M.P.G., I.D., L.C.R.S., and J.-F.H. analyzed data; C.G.-C., É.B., Y.B., M.P.G., I.D., L.C.R.S., and J.-F.H. wrote the paper; and É.B. provided funding.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1816686116/-/DCSupplemental.