Rising CO2 drives divergence in water use efficiency of evergreen and deciduous plants

Wuu Kuang, et al.


Intrinsic water use efficiency (iWUE), defined as the ratio of photosynthesis to stomatal conductance, is a key variable in plant physiology and ecology. Yet, how rising atmospheric CO2 concentration affects iWUE at broad species and ecosystem scales is poorly understood. In a field-based study of 244 woody angiosperm species across eight biomes over the past 25 years of increasing atmospheric CO2 (~45 ppm), we show that iWUE in evergreen species has increased more rapidly than in deciduous species. Specifically, the difference in iWUE gain between evergreen and deciduous taxa diverges along a mean annual temperature gradient from tropical to boreal forests and follows similar observed trends in leaf functional traits such as leaf mass per area. Synthesis of multiple lines of evidence supports our findings. This study provides timely insights into the impact of Anthropocene climate change on forest ecosystems and will aid the development of next-generation trait-based vegetation models.


Climate change will likely alter future carbon and hydrologic cycles (1). These cycles are closely tied to plant assimilation of atmospheric CO2 through photosynthesis by the regulation of CO2 and water vapor exchange via small pores on the leaf surface, called stomata. CO2 uptake is necessarily accompanied by water loss through stomata, and this carbon gain to water loss metric is generally referred to as water use efficiency (2). At the leaf level, variation in the photosynthesis (A)–to–stomatal conductance (gs) ratio over a leaf life span represents a time-integrative or averaged estimate of the intrinsic water use efficiency (iWUE), operating at a common evaporative demand (2). Thus, iWUE, a form of water use efficiency, is an important measure of the potential water cost of maintaining a given rate of carbon assimilation per unit leaf area.

A primary response of plants to increasing CO2 is to increase A and is often accompanied by reducing diffusive gs to minimize transpirational water loss (3). As a result, iWUE is generally known to increase with rising atmospheric CO2 (4). However, the magnitude and direction of iWUE responses to elevated CO2 at broad ecosystem and species ranges in natural ecosystems are poorly understood. Specifically, the decadal responses of two key plant functional groups, evergreen and deciduous, are not clear; this is important given that these functional groups occur across many taxonomic groups, and their relative proportions largely define global ecosystems and ecosystem functions and services (5, 6). It is hypothesized that evergreen plants are more sensitive in their iWUE response to elevated atmospheric CO2 than deciduous plants (7). However, to date, experimental CO2 enrichment studies, which were based on limited species and ecosystem type, are equivocal (7).

Here, we assessed the impact of human-driven increases in atmospheric CO2 [~45 parts per million (ppm)] over the past ~25 years on the iWUE of deciduous versus evergreen plants (244 species; table S1). We focus on iWUE responses of woody taxa from 20 field sites spanning eight biomes between two time periods: 1988–1991 and 2013–2015 (Fig. 1A and table S2). To compare the iWUE response of contemporary (2013–2015) to historical plants (1988–1991), we used a unique georeferenced herbarium collection of C3 woody flowering species, known as Climate-Leaf Analysis Multivariate Programme (CLAMP) (8), to represent historical samples. We compared these to contemporary leaves collected 25 years later by our team from the same species at the same sites (which we will refer to as species sites) and biomes (which we will refer to as species biomes). We inferred iWUE using leaf stable carbon isotopes (δ13C). To minimize variability in leaf δ13C between historical and contemporary samples—due to possible differences in phenology, seasonality, and field protocols—we operated the same field sampling protocol as CLAMP (8) and sampled during approximately the same collection season or month as the historical leaves.

Fig. 1 Major study areas and iWUE values across biomes in historical and contemporary samples.

(A) Major study areas. (B) Historical and contemporary iWUE at 355 and 400 ppm atmospheric CO2 concentration respectively arranged by increasing averaged iWUE values. Boxplots show median (center line), mean (red dot), interquartile range (IQR), 1.5 times of IQR (whiskers), and outliers (black dots). Numbers in brackets are the number of leaves. All iWUE gains are likely to be larger than zero.


A total of 2031 historical and contemporary leaves were analyzed for leaf δ13C, leaf mass per area (LMA), carbon per mass (Cmass), and nitrogen per mass (Nmass). There is no likely difference in average total LMA and Nmass between the historical and contemporary samples [ΔLMA = − 0.4 g m−2; 95% credible interval (CI95%), –1.4 to 0.6; ΔNmass = 0.06%; CI95%, –0.12 to 0.27] and the slopes of regression between the two time periods through the origin are close to 1 (LMA slope = 0.97; CI95%, 0.96 to 0.98; r2 = 0.92; Nmass slope = 0.97; CI95%, 0.94 to 1.00; r2 = 0.93) (fig. S1). Average evergreen LMA is likely higher than deciduous within each biome in both the historical and contemporary samples (table S3).

An unequivocal increase in average iWUE (ΔiWUE) was observed in all eight biomes investigated, ranging from highest in the tropical seasonal moist forest [TSF(M)] (17.2 μmol mol−1; CI95%, 14.3 to 20.0) to lowest in the tropical rainforest (TF) (5.2 μmol mol−1; CI95%, 1.6 to 8.3) (Fig. 1B and table S4). Among the seven biomes with both evergreen and deciduous groups, evergreen species generally demonstrated a greater ΔiWUE in response to ~45 ppm rise in CO2 than deciduous plants, within cooler biomes (Fig. 2A and table S5): this trend also prevailed when data were further grouped into growth habit (tree versus shrub) or high- and low-light habitat (understory subcanopy versus open canopy) (figs. S2 and S3 and tables S6 and S7). A substantial decrease in the ratio of leaf intercellular CO2 (ci) to ambient atmospheric CO2 (ca), ci/ca, in evergreens compared with deciduous taxa resulted in a higher calculated iWUE gain (fig. S4). Our results agree well with published studies that have reported either a decrease in ci/ca (9, 10) or a near constant ci/ca (11, 12) for tree species. Differences between average iWUE gain in evergreen and deciduous taxa (ΔiWUEe-d) widened, however, with decreasing mean annual temperature (MAT) from the tropical toward the boreal biomes (slope = −0.395; CI95%, –0.770 to –0.004; r2 = 0.70; Fig. 2B).

Fig. 2 iWUE gain (∆iWUE) and its relationship with mean annual temperature (MAT) and leaf mass per unit area (LMA).

Dotplots represent mean of posterior distributions (n = 6000 samples), CI95%. Red line is the fitted regression. (A) ∆iWUE of deciduous and evergreen plants in biomes arranged by increasing MAT. (B) Differences between evergreen and deciduous ΔiWUE (ΔiWUEe-d) versus MAT, ΔiWUEe-d = 11 – 0.4MAT, r2 = 0.70. (C) ΔiWUEe-d versus average difference of evergreen and deciduous LMA (LMAe-d), ΔiWUEe-d = −2.0 + 0.14 LMAe-d, r2 = 0.80. (D) Boxplots of deciduous and evergreen LMA across biomes for combined historical and contemporary samples arranged by increasing MAT. All P(LMAevergreen > LMAdeciduous) ≥ 0.95. (E) Comparison of the rate of iWUE gain per unit of CO2 concentration (ΔiWUE/ΔCO2) for total deciduous and evergreen samples [P(ΔiWUE/ΔCO2 evergreen > ΔiWUE/ΔCO2 deciduous) = 0.87]. (F) Scatter plot of LMA versus MAT of evergreen and deciduous plants for combined historical and contemporary samples, n = 2031 leaves.

In this study, atmospheric CO2 is likely a dominant factor for iWUE gain because of the likely difference in atmospheric CO2 concentration between the two time periods (Mauna Loa station; CI95%, 43.60 to 45.89 ppm). In contrast with this, other influential climatic variables, such as air temperature and vapor pressure deficit (VPD) showed only small changes with no likely difference statistically within biomes at CI95% (table S8). Furthermore, our result demonstrated that the small changes in MAT (∆MAT) and VPD (∆VPD) between historical and contemporary periods in this study were unlikely to affect ΔiWUEe-d, as the differences in ΔiWUE between evergreen and deciduous within the same biome were not highly influenced by ∆MAT or ∆VPD (fig. S5).

In relation to leaf functional traits, ΔiWUEe-d also varied increasing tightly (r2 = 0.80) with the biome average difference between LMA in evergreen and in deciduous species (LMAe-d; slope = 0.14; CI95%, 0.05 to 0.23; Fig. 2, C and D). The total average ΔiWUE value for each deciduous and evergreen group, with all biomes combined, was quantified by normalizing iWUE with VPD, temperature, precipitation, and altitude using models developed in this study (table S9). We found that average ΔiWUE was higher in evergreen than in deciduous species [P(ΔiWUEevergreen > ΔiWUEdeciduous) = 1] with gains of ~39% (17.1 μmol mol−1; CI95%, 13.8 to 20.5) and ~15% (7.8 μmol mol−1; CI95%, 5.0 to 10.4), respectively. These correspond to an iWUE gain of 0.39 μmol mol−1 ppm−1 (CI95%, 0.30 to 0.46) in evergreen and 0.18 μmol mol−1 ppm−1 (CI95%, 0.12 to 0.25) in deciduous species [P(ΔiWUE/ΔCO2evergreen > ΔiWUE/ΔCO2deciduous) = 0.99] (Fig. 2E).

The divergence of evergreen and deciduous ΔiWUE along a MAT gradient (−1.4° to 26.7°C) parallels those observed for LMA (Fig. 2F) and Nmass (fig. S6). The LMA divergence in functional groups from warmer to colder sites (27.5° to −16°C) was observed in a previous study (13) and was associated with LMA increment with leaf life span; this divergent trend is related to the requirement of leaves with longer life spans to maximize carbon gain in shorter growing seasons, i.e., in colder biomes (14). Our results demonstrated how this well-studied trend (13, 14), in LMA divergence from warmer to colder biomes, also manifests in the differential response of evergreen and deciduous taxa to anthropogenic CO2 rise. The smaller differences in LMA between the leaf habit classes in the warmer biomes compared with the colder biomes contributed to the observed trend. High LMA generally occurs in woody evergreens because of their robust leaf structure, which can incur resistance to CO2 diffusion and, hence, lower mesophyll conductance (gm) (7, 15, 16). Therefore, evergreen leaves, in general, are likely to operate at lower gm values than deciduous leaves (16, 17).

Under elevated CO2, leaves with low gm may show a higher increase in A than high gm taxa, and their A is less sensitive to reduction in gs—this, in turn, leads to strong iWUE gain (iWUE = A/gs) (7). At a given gs, A of leaves with low gm (i.e., evergreens) is more limited by lower chloroplast CO2 concentration (cc) and, thus, responds more strongly to rising CO2. The reason for this is that the higher cc gets, the less CO2 affects photosynthesis because of the saturation of the A versus cc relationship (7). We did not measure gm, but we did observe greater average LMA and iWUE responses in evergreens than in deciduous species, suggesting increased CO2 diffusion limitations in the former. LMA and gm are inversely correlated, but the relationship is confounded by mesophyll cell wall thickness and chloroplast surface area that can vary across environmental gradients and species (15, 18). Therefore, in this study, high LMA was associated with greater iWUE response to a ~45-ppm rise in atmospheric CO2 concentration in evergreen compared with deciduous leaves (Fig. 2, C and E).

To validate our results from the two time periods, we used published tree ring δ13C datasets (1970–2013) and leaf δ13C datasets (1981–2005) (1921) containing continuous recent sampling points to track iWUE trends along a rising atmospheric CO2 gradient (ΔiWUE/ΔCO2). The meta-analysis of tree ring iWUE data showed higher average iWUE response in evergreen (0.29 μmol mol−1 ppm−1; CI95%, 0.27 to 0.33) than deciduous (0.21 μmol mol−1 ppm−1; CI95%, 0.18 to 0.24) trees (Fig. 3A, fig. S7, and table S10). Evergreen trees in the boreal-temperate region(s), which were all gymnosperms in the published datasets (seven species), showed a greater average rate of iWUE gain (0.33 μmol mol−1 ppm−1; CI95%, 0.30 to 0.36) than their angiosperm and gymnosperm deciduous counterparts (four species) (0.14 μmol mol−1 ppm−1; CI95%, 0.11 to 0.17), but in the tropics, this disparity was not observed (Fig. 3B). This result corroborated with published studies that showed the average gm of temperate evergreen gymnosperm was onefold lower than temperate deciduous angiosperms (15, 16). Furthermore, a tree ring study at 23 sites across Europe showed that evergreen gymnosperm trees (four species) increased their iWUE substantially more than deciduous angiosperm trees (two species) in the last c. 100 years at ~22 and ~14%, respectively (10). Our meta-analysis of published leaf δ13C data from woody angiosperm species showed the same trend of higher collective iWUE increase (ΔiWUEc/ΔCO2) in evergreen (0.76 μmol mol−1 ppm−1; CI95%, 0.62 to 0.91) than in deciduous (0.51 μmol mol−1 ppm−1; CI95%, 0.32 to 0.70) leaves (Fig. 3C and fig. S8). These results confirm our original observations from the two time periods: There is an overall stronger iWUE gain in evergreen compared with deciduous species (Fig. 2, A and E) in response to rising atmospheric CO2.

Fig. 3 Average slopes of iWUE responses to rising atmospheric CO2 concentration (ΔiWUE/ΔCO2) from published data.

Dotplots represent mean of posterior distributions (n = 6000 samples), CI95%. (A) ΔiWUE/ΔCO2 from published tree ring δ13C data for the various time intervals between 1970 and 2013 for evergreen (n = 29 trees) and deciduous trees (n = 23 trees). (B) Result from (A) separated into bioclimatic zones showing higher average iWUE gain in evergreen (n = 24 trees) than in deciduous trees (14 trees) in the boreal-temperate zone, but the opposite in the tropical zone (deciduous n = 9 trees; evergreen n = 5 trees) [P(ΔiWUE/ΔCO2 deciduous > ΔiWUE/ΔCO2 evergreen) = 0.95]. (C) ΔiWUEc/ΔCO2 calculated from published leaf δ13C data collected between 1981 and 2005 for deciduous (n = 470 species sites) and evergreen (n = 1053 species sites) species.

To further test this differential evergreen/deciduous response to ~45-ppm rise in CO2, we used data from a field infrared gas exchange analysis (IRGA) experiment conducted in situ on a subset of the same leaves used for this δ13C study. Leaf A and gs responses to ~355- and ~400-ppm cuvette CO2 concentration were measured, referencing values for the historical and contemporary period, respectively. The responses measured with the gas analyzer were instantaneous responses to CO2 concentration rather than long-term responses (decadal) that are most likely influenced by acclimation. This experiment showed that average gain in leaf iWUE in evergreen leaves (0.22 μmol mol−1 ppm−1; CI95%, 0.20 to 0.25) was likely higher than that in deciduous leaves (0.20 μmol mol−1 ppm−1; CI95%, 0.17 to 0.23) [P(ΔiWUE/ΔCO2evergreen > ΔiWUE/ΔCO2deciduous) = 0.92] (Fig. 4A). Results from our in situ gas exchange study showed that an increase in A can largely contribute to an increase in iWUE under a ~45-ppm CO2 rise with higher average A gain in evergreen (22.4%; CI95%, 19.1 to 25.7) than in deciduous leaves (16.7%; CI95%, 13.4 to 20.1) (Fig. 4B). However, gs instantaneous responses showed no likely change in both groups (evergreen: –0.2%; CI95%, –2.3 to 1.8; deciduous: 1.0%; CI95%, –1.1 to 3.2) (Fig. 4C). Evergreen ci/ca showed a likely decrease, but no change was observed in deciduous leaves (evergreen: −0.015 Pa; CI95%, –0.019 to –0.010; deciduous: −0.001 Pa; CI95%, –0.003 to 0.006).

Fig. 4 In situ field infrared gas exchange analyses experiments showing the change of in vivo measurements of iWUE (∆iWUE), photosynthesis (∆A), and stomatal conductance (∆gs) in response to a shift in cuvette CO2 concentration of ~355 ppm (simulated historical) to ~400 ppm (contemporary ambient).

Dotplots represent means of posterior distributions (n = 6000 samples), CI95%. Evergreen n = 135 leaf samples (33 species); deciduous n = 119 leaf samples (31 species). (A) Dotplots of ∆iWUE in evergreen and deciduous leaves. (B) Dotplots of ∆A in evergreen and deciduous leaves. (C) Dotplots showing average ∆gs in evergreen and deciduous are unlikely to be higher than zero at CI95%.

Currently, these experimental results (Fig. 4C) do not account for possible anatomical adaptions in stomatal density and/or size that could influence gs. Stomatal density in most plant species is well known to decrease with increasing atmospheric CO2 concentration that could lead to a general decrease in maximum stomatal conductance (22). Work is therefore ongoing to assess anatomical adaptations at the species and functional group level to test these conclusions further. Results from the in situ IRGA measurements, which estimate the instantaneous responses to CO2, lend support to the long-term observations from our extensive biome-level field-based δ13C study and suggest that the magnitude of iWUE change observed here is due to a substantial increase in A coupled with little or no change in gs. Together, these results suggest that notable adjustment of photosynthetic biochemistry has occurred in woody vegetation with ~45-ppm CO2 rise.


Our biome-wide field study of iWUE responses to a mere 45-ppm CO2 rise between 1988 and 2015 suggests greater average iWUE gain in evergreen than in deciduous species, particularly in the cooler climate biomes. The diverging trend in iWUE gain highlights a strong link between LMA, MAT, and plant-CO2 responses in woody evergreen and deciduous taxa: This is strongly associated with the more distinct differences in LMA and leaf phenological traits observed between evergreen and deciduous taxa in colder biomes than in warmer biomes. This knowledge has the potential to enhance development of new-generation trait-based vegetation models, of which temperature, photosynthetic water use, and LMA are important components. That the differential response of evergreen and deciduous leaf habits in natural ecosystems has been given little attention to date is unexpected given that such a profound physiological response occurring at a continental scale could incur a substantial shift in natural forest and woodland ecology (e.g., forest fraction of evergreeness and deciduousness) and alter seasonal energy, water, and carbon balance and dynamics. Our results indicate that future increases in atmospheric CO2 may confer a competitive advantage to woody angiosperm evergreens over their deciduous neighbors to a greater extent in cooler biomes than in warmer biomes. Therefore, understanding of the differential physiological response induced by climate change in evergreen and deciduous taxa will improve our ability to build more mechanistic and predictive models on vegetation response to future climate change. While our field study covered a substantial number of woody angiosperm species, and was supported by published tree ring δ13C data that included gymnosperm species (seven evergreen and two deciduous species), future research may benefit by including more gymnosperm species to confirm the differential response of leaf habits within this group to rising atmospheric CO2, particularly in the conifer-dominated boreal biome. Further profound increases in atmospheric CO2 are projected by the year 2050 under all representative concentration pathway (RCP) scenarios [RCP 2.5 = 443 ppm; RCP 4.5 = 487 ppm; RCP 6.0 = 478 ppm; RCP 8.5 = 541 ppm (23, 24)]. In this context, higher iWUE under elevated CO2 atmospheres may have contributed to evergreen expansion in past greenhouse intervals such as the Eocene (ca. 55 million years ago), particularly in seasonally dry areas of the mid latitudes (25), rather than to elevated temperatures alone, which is the current paradigm (26).


Study sites, herbarium samples, and fieldwork

Historical herbarium samples from the CLAMP collected using the same protocol and person (Wolfe) (8) in 1988–1991 were recollected in 2013–2015 by our team (W.K.S., M.M., and J.C.M.). This yielded contemporary leaf samples of the same species from the same sites/biomes. The same standard collection protocol was used for both historical and contemporary samples. This approach was used to minimize variability of leaf δ13C. To our knowledge, CLAMP, a unique georeferenced global inventory of C3 woody angiosperm leaf physiognomic data (8, 27), is the only herbarium archive that was collected by the same person (Wolfe) using the same protocol over several biomes with each including many species (average, 25 species per site). In this study, field sites in each biome were selected from the CLAMP archive. Of the original 173 sites sampled by Wolfe (8), we selected 20 to represent eight of Whittaker’s vegetation biomes (28): boreal forest (BF), temperate rainforest, temperate deciduous forest (TDF), Mediterranean (MED), subtropical desert, tropical seasonal dry forest [TSF(D)], TSF(M), and TF (table S2). We restricted selection to sites below 700 m above sea level to limit the influence of lower CO2 partial pressure and atmospheric pressure on leaf traits and carbon isotope composition (δ13C) at higher altitudes. Site selection was based on individual site accessibility within the planned data gathering schedule and acquisition of the required scientific collection permits. Where possible, we selected three sampling sites in each biome, except in the TF of Fiji (two sites) and TSF(D) in Puerto Rico (one site). As a result of using CLAMP herbarium samples, sites in the boreal and temperate biomes are restricted to Northern America, with tropical biomes in Puerto Rico [TSF(D) and TSF(M)] and Fiji (TF). Although all the tropical biome sites are situated on islands, the plants species sampled here are from areas that experience tropical climate. We are confident that our tropical sites are representative tropical biomes as there is no evidence to suggest that the physiology of tropical island vegetation differs from that on a tropical mainland, especially at the leaf level. For instance, one of the best studied tropical forests in the world is Barro Colorado Island in Panama. Only evergreen plants were sampled in Fiji, and therefore, this biome was not used to quantify ΔiWUEe-d. To obtain a representative sample of C3 woody angiosperm species within the BF, which is usually dominated by conifers, our sampling was conducted within the interior BF zone of Alaska, which has extensive areas of open and closed deciduous forests (29). Regarding our BF sites, deciduous trees make up virtually all of the native angiosperm tree population, while the gymnosperms are mostly evergreen trees. Since we are making a direct comparison of the historical CLAMP samples with contemporary samples of exactly the same species from the same locations, we were prohibited from including gymnosperms. As a result, our fieldwork study on BF only covered angiosperms of three leaf habit and growth habit groups without evergreen trees. These included deciduous trees, deciduous shrubs, and evergreen shrubs.

Contemporary leaf samples were collected in the field between 2013 and 2015 from the same species as those in the historical CLAMP herbarium collected between 1988 and 1991 from the same sites or biomes. All fieldwork was carried out in the growing season (table S2), corresponding as closely as possible to the collection month of historical samples. Tree and shrub growth habits were sampled in all biomes and were largely represented in both evergreen and deciduous plant groups. Our sampling focused on outer-canopy leaves, meaning “sun” leaves for plants growing in relatively open environments, and leaves exposed to sun flecks when sampling naturally shade-dwelling species. We sampled fully expanded leaves, the developmental stage at which many leaf traits are relatively stable. In one aspect of the statistical analyses in this study (see section on Statistical analysis), we divided our dataset into two broadly defined habitat groups based on our field observations to reflect high- and low-light habitat: open canopy and understory subcanopy. For this study, open canopy refers to plants that are located either in open areas or at the forest canopy edge and receiving direct sunlight. By contrast, understory subcanopy refers to plants occurring within the forest canopy, in shade but receiving sun flecks. In all biomes, we sampled both the open-canopy and understory-subcanopy habitats for evergreen and deciduous plants, except for the BF and TDF biomes, there were no evergreen plant samples in the open-canopy habitat, and in the subtropical desert biome, all habitats were classified as open canopy. In the historical CLAMP samples, sun-exposed twigs were collected that may be directly exposed to the sun or sun fleck subjected to a species natural habitat. On each herbarium specimen, we had carefully selected leaves that were fully expanded (i.e., visually mature) and thick to increase the chance of including mature sun-exposed leaves.

To minimize the potentially confounding influence of height on leaf δ13C and LMA, leaves from tall trees were collected at basal-exterior canopy level within arm’s reach, up to 3 m in line with CLAMP historical collection methods. This protocol standardized collection height with historical samples. Before collection, the leaves gathered for trait analysis were also used for physiological measurements (see section on “In situ field IRGA experiments”). Our sampling protocol is in accordance with the collection methods used by Wolfe (8) following the CLAMP protocol. That is, our protocol standardizes historical and contemporary sampling methods, with the aim of reducing trait variability caused by sampling method and relevant biotic and abiotic factors that may have differed between contemporary and historical sampling periods.

Only broadleaf woody C3 angiosperm species were sampled for this study (gymnosperms, grasses, and crops were not included). A total of 1550 contemporary leaf samples, each from individual plants, were collected in the contemporary fieldwork. A total of 481 historical leaf samples were subsampled from the CLAMP herbarium collection. The entire dataset used in this study comprises 244 matching historical and contemporary woody angiosperm species from 64 families (table S1). All specimens were identified to species level. Taxonomic nomenclature was updated using the online Taxonomic Name Resolution Service v 4.0.

Climate data

Mean monthly precipitation, mean monthly air temperature, maximum monthly air temperature, and vapor pressure over time periods (1988–1991 and 2013–2015) for each study site were obtained from 0.5° × 0.5° resolution Climate Research Unit data (CRU TS v.4.0) (30) gridded dataset via The Royal Netherlands Meteorological Institute (KNMI) Climate Explorer. Monthly saturated vapor pressure was calculated from maximum monthly air temperature. These were then subtracted with monthly vapor pressure to obtain monthly VPD (31) and used to infer leaf-to-air VPD. MAT and mean annual precipitation (MAP) were calculated from the monthly data.

Analysis of leaf traits

Leaf samples were oven dried at 50° to 60°C for 2 days. One half of each dried leaf blade was used for LMA analysis and the other half for δ13C, carbon (C), and nitrogen (N) elemental analyses. To standardize LMA data collection from both historical and contemporary leaves, all leaves were rehydrated. Leaf area shrinkage from drying can be reversed by rehydration (32). LMA was determined by dividing the dry leaf mass by the rehydrated leaf area. For the δ13C, N, and C elemental analyses, dried leaf fragments were placed with a tungsten bead in Eppendorf tubes and finely ground in a mixer mill (Tissue Lyser, Qiagen Inc., Valencia, CA, USA). Each sample (~3 mg) was then enclosed in a tin capsule using a crimper plate. Samples were analyzed for δ13C, C, and N using a PDZ Europa ANCA-GSL elemental analyzer interfaced with a PDZ Europa 20-20 isotope ratio mass spectrometer (Sercon Ltd., Cheshire, UK) at UC Davis Stable Isotope Facility, University of California, Davis, USA. Instrumental error was ±0.18‰ (per mil) for 13C (SD). Carbon isotope composition was calculated as

δ13C ()=(RsampleRstandard)/Rstandard×1000
(Eq. 1)where Rsample and Rstandard are the 13C/12C ratio of the sample and the international standards Vienna Pee Dee Belemnite, respectively. Carbon isotopic discrimination (∆plant) is given as

(Eq. 2)

In relation to the intercellular CO2 (ci) and ambient CO2 (ca) partial pressures, ∆plant in C3 leaves is given as follows (33)

Δplant=a+(ba) (ci/ca)

(Eq. 3)where a is the fractionation due to diffusion in air (4.4‰) and b is the net fractionation caused by carboxylation (27‰). Equation 3 is widely used and assumes that the effects of boundary layer, internal conductance, photorespiration, day respiration, and allocation are negligible. Atmospheric CO2 concentration (ca) and δ13Cair information were taken from a published instrumental dataset (1980–2015) from the Mauna Loa station (3436) corresponding to the historical and contemporary collection months (table S2). The full equation of ∆plant includes several elements such as photorespiration, day respiration, and the CO2 mole fractions in the ambient air, at the leaf surface, in the intercellular air spaces, and at the chloroplast (cc) (37, 38). Photorespiration and cc are known to influence ∆plant (38), and therefore, it is desirable to include these traits. However, we did not measure photorespiration and gm; the latter is required for estimating cc. In this study, we were interested in quantifying the differences between evergreen and deciduous iWUE (ΔiWUEe-d) rather than their absolute values. On the basis of this reasoning, the use of the simplified linear model of Farquhar et al. (33) (Eq. 3) as an approximation to ∆plant is appropriate for the purpose of this study.

iWUE can be expressed as the ratio of photosynthesis (A) and leaf conductance to water vapor transfer (g) in Eq. 4 below (33) using ci/ca calculated from Eq. 3 and ca


(Eq. 4)

iWUE inferred from δ13C is an average estimate of iWUE over a leaf life span, i.e., time integrated.

Statistical analysis

All statistical analysis was undertaken using JAGS 4.1.0. (39) and R statistical software (40). Bayesian models using JAGS, through the R package rjags (41) interface, were used: Inference of each parameter was made from Markov Chain Monte Carlo (MCMC) sampling from 6000 samples of the posterior distribution from three chains, each with 10,000 iterations with a burn-in of 2000 and a thin rate of 4 (42). Normal distribution priors with mean zero and variance 100 were used for intercept and slope parameters, while a uniform (0, 10) prior was used for the SD on the variance terms. Convergence was checked by visual assessment of MCMC chains and using the Gelman-Rubin statistic (42). Mean of trait or group was calculated from posterior distributions. CI95%s of parameter estimates were calculated as the 2.5 and 97.5% quantile of posterior distributions. The 50% credible interval (CI50%) of parameter estimates were calculated as 25 and 75% quantile of posterior distributions. The CI95% represents the interval that captures 95% of the posterior distribution, e.g., when the CI95% for a statistics score is between a and b, this means that we have a 95% chance of having a score between a and b (note that credible interval is different from confidence interval). A CI50% statistics score between a and b implies a 50% chance of having a score between these two values. Therefore, the extent of CI overlapping with zero determines how likely a value is close to zero. Statistical comparisons between groups were made by examining value of CI95% and/or by probability of group differences bigger than or smaller than zero, e.g., “P(x > y) = z” denotes that the probability of variable x being bigger than variable y, given the data, is z.

To evaluate the robustness of our sampling method in minimizing the variability between the historical and contemporary samples, we first statistically test the difference in the mean of LMA and Nmass in the two time points. Second, we plotted historical and contemporary samples through the origin each for LMA and Nmass. A regression slope that is close to 1 would indicate a general level of uniformity between the historical and contemporary samples. LMA and Nmass are well known to vary with plant height, sun and shade leaf morphotypes, and age (43, 44).

We aggregate across biomes the iWUE at each time point (historical versus contemporary) to calculate the total gain in iWUE (∆iWUE) for the deciduous and evergreen species groups, using statistical models incorporating environmental variables (environment-normalized model) (Fig. 2E). However, samples from the TF biome (Fiji) were excluded because of the absence of deciduous plant samples. The environment-normalized model standardizes the aggregated iWUE values when calculating the total gain in iWUE: Leaf δ13C or its derived variables (e.g., iWUE and ci/ca) are widely known to be confounded by latitude (20), altitude (19, 20), and site climatic variables such as VPD (45), temperature (1921, 45), and precipitation (1921). Using our own dataset, we examined the relationship between iWUE and environmental variables such as altitude, latitude, and bioclimatic variables (precipitation, temperature, and VPD). Our aim was to generate an equation that could be used to normalize iWUE values against environmental variables when aggregating data across biomes (see Fig. 2E).

For evergreen species, we averaged site monthly precipitation, temperature, VPD, and atmospheric CO2 concentration by 12 months up to and including the collection month to match the average period of photosynthetic opportunities. One meta-analysis study showed that mean annual climate parameters were more likely to match evergreen photosynthetic windows for carbon isotope discrimination of C3 plants (21). Although photosynthesis of evergreens is reduced during winter time with small winter carbon gain (46, 47), this may still influence the average carbon isotope discrimination in a leaf life span. The leaf life span of evergreen angiosperms in the boreal-temperate and tropical biomes each showed a skewed distribution with central tendencies (median) of approximately 18 and 15 months (48), respectively (fig. S9). Therefore, our approach of averaging site climatic data by a period of 12 months up to and including the collection month was a reasonable approximation of evergreen leaf life span collected at the time. This approximation took into consideration the fact that we sampled only fully expanded leaves that were neither young nor too old (i.e., visibly unhealthy). For deciduous species, we averaged these climate variables from the start of growing months up to and including the collection month.

The correlation matrix between iWUE and the foregoing environmental variables are presented in table S11. VPD shows the strongest correlation with iWUE (r2 = 0.26) followed by precipitation (r2 = 0.24), altitude (r2 = 0.20), and absolute latitude (r2 = 0.10). Temperature shows the weakest correlation with iWUE (r2 = −0.05) but is instead strongly correlated with absolute latitude (r2 = −0.93), precipitation (r2 = 0.65), and VPD (r2 = 0.53), and weakly correlated with altitude (r2 = 0.10). Therefore, temperature was not included in our model because of the extreme collinearity between covariates, which could lead to high correlation in some of the posterior parameter estimates. Last, our statistical model consists of iWUE as the dependent variable, while time (factor), altitude, averaged site VPD, and precipitation are the independent variables (Model 1). Latitude was excluded from the model because its coefficient was subsequently shown to likely contain zero at CI95% when included in the regression. To calculate the rate of iWUE change in relation to atmospheric CO2 concentration, the same model was used with time factor replaced by CO2 concentration (Model 2). In the following models, each i represents one leaf. See table S9 for coefficient values.

(Model 1)where iWUEi is the iWUE of individual i; Timei is the categorical time variable (historic and contemporary) corresponding to individual i; VPDi is the VPD corresponding to individual i; PREPi is the precipitation corresponding to individual i; ALTi is the altitude corresponding to individual i; αj(i) is the intercept of the iWUE-time relationship in categorical leaf habit j (deciduous and evergreen); βj(i) is the slope of the iWUE-time relationship in categorical leaf habit j (deciduous and evergreen), this is ΔiWUE; λ1 is the slope of the iWUE-VPD relationship; λ2 is the slope of the iWUE-PREP relationship; λ3 is the slope of the iWUE-ALT relationship; and εi is the residual of individual i.

(Model 2)where, iWUEi is the iWUE of individual i; (CO2)i is the atmospheric carbon dioxide concentration corresponding to individual i; VPDi is the atmospheric VPD corresponding to individual i; PREPi is the precipitation corresponding to individual i; ALTi is the altitude corresponding to individual i; αj(i) is the intercept of the iWUE-CO2 relationship in categorical leaf habit j (deciduous and evergreen); βj(i) is the slope of the iWUE-CO2 relationship in categorical leaf habit j (deciduous and evergreen), this is ΔiWUE/ΔCO2; λ1 is the slope of the iWUE-VPD relationship; λ2 is the slope of the iWUE-PREP relationship; λ3 is the slope of the iWUE-ALT relationship; and εi is the residual of individual i.

For βj(i), the slope of the iWUE-CO2 relationship, the actual full unit of ΔWUEi/ΔCO2 is μmol CO2 mol−1 H2O/μmol CO2 mol−1 air: For simplicity and readability, we prefer to use μmol mol−1 ppm−1. We further investigate ∆iWUE in evergreen and deciduous plants in each biome by dividing the dataset into growth habit (shrub versus tree) or habitat (understory-subcanopy versus open-canopy) categories. In each category, the probability of evergreen ∆iWUE higher than deciduous ∆iWUE was calculated.

In situ field IRGA experiments

Photosynthesis and photosynthetic water use were measured on 254 leaf samples from 64 of our δ13C study species. Measurements were made with a CIRAS-2 gas analyzer (PP Systems, Amesbury, MA, USA) attached to a PLC6 (U) cuvette fitted with a 1.7-cm2 measurement window and a red/white-light light-emitting diode unit. Measurements were carried out between June and August 2014 at two BF sites (16 species, Bird Creek and Kenai, Alaska, USA), one TDF site (11 species, Smithsonian Environmental Research Center, Maryland, USA), two TSF(M) sites (15 species, Cambalache and Guajataca, Puerto Rico), and one TSF(D) site (9 species, Borinquen, Puerto Rico), all from a subset of the contemporary samples. Photosynthesis (A) and stomatal conductance (gs) were assessed on an average of four individual plants per species between 9:00 am and 13:00 pm. A sun-exposed branch was sampled from each plant using a pruner and was immediately recut under water (49). Following this, a fully expanded leaf from each branch was enclosed in the cuvette of the gas analyzer, which was running at a subambient 1988–1991 averaged reference CO2 concentration of 355 ppm. Stomatal conductance at subambient CO2 concentration was recorded upon stabilization of its value, which typically took less than 15 min. Subsequently, reference CO2 was established at 400 ppm (year 2016 values), and the leaf was left to equilibrate for at least 15 min before gs at contemporary ambient atmospheric CO2 was recorded. Randomization of the sequence of the two treatments was ensured; overall, about 65% of the measurements started at 400 ppm and were reduced to 355 ppm, while the rest of measurements (35%) started at 355 ppm and were increased to 400 ppm. On several occasions, the reversibility of the CO2 effects on A and gs was tested. This was done by measuring gs at a starting CO2 concentration of 400 ppm, after which CO2 was reduced to 355 ppm for several minutes before it was returned to the initial concentration of 400 ppm. The final A and gs values at 400 ppm were the same as those initially recorded.

Meta-analysis: iWUE to atmospheric CO2 trend from tree ring δ13C studies

iWUE data calculated from tree ring δ13C were used to quantify the iWUE-CO2 response of individual deciduous and evergreen trees along a decadal time series of various time intervals between 1970 and 2013. Data were compiled from 17 published studies (5066) consisting of 52 trees from 22 species, of which 23 trees were deciduous (12 species) and 29 evergreen (10 species). Atmospheric CO2 concentration data were acquired from the Mauna Loa station data (3436). Annual δ13Cair information was obtained from published ice-core data. iWUE values were calculated from δ13C by using Eq. 3. For each study site, we obtained mean monthly precipitation, mean monthly air temperature, maximum monthly air temperature, and vapor pressure from 0.5° × 0.5° resolution CRU TS v.4.0 (30) gridded dataset for the period of δ13C for each individual tree. VPD values were calculated as per the method described in the section “Climate data.” Regression slopes (ΔiWUE/ΔCO2) for individual trees were determined by fitting a simple linear model (using the Bayesian linear regression approach, see section on “Statistical analysis”) with iWUE as the dependent variable, and atmospheric CO2 concentration, VPD, and MAP as the independent variables. In the following model, each i represents a value from a growth ring as determined in a study, from a tree, j


(Model 3)where, iWUEi is the iWUE of individual i; (CO2)i is the atmospheric carbon dioxide concentration corresponding to individual i; VPDi is the atmospheric VPD corresponding to individual i; MAPi is the MAP corresponding to individual i; αj(i) is the intercept of the iWUE-CO2 relationship in categorical individual tree j; βj(i) is the slope of the iWUE-CO2 relationship in categorical individual tree j; λ1 is the slope of the iWUE-VPD relationship; λ2 is the slope of the iWUE-PREP relationship; and εi is the residual of individual i.

By including VPD and MAP in the regression, we normalized the response slope of each tree with climatic variables, VPD and MAP. MAT is excluded from the model because of the strong collinearity with VPD (r2 = 0.72). The values for λ1 and λ2 are 5.47 (CI95%, 4.01 to 6.97) and −0.08 (CI95%, −0.09 to −0.06), respectively. On a centennial scale, a long-term iWUE fluctuation along the atmospheric CO2 gradient generally follows an exponential increase. However, we can reasonably approximate the iWUE trend with a linear model at a shorter decadal time scale. This shorter decadal time scale varies between 10 and 40 years from 1970 to 2013 depending on studies. Last, ΔiWUE/ΔCO2 values from posterior distributions of trees (6000 samples for each tree) were aggregated into deciduous and evergreen plant groups by averaging ΔiWUE/ΔCO2 values from posterior distributions. This approach therefore takes account of the uncertainty of ΔiWUE/ΔCO2 values of each tree. Further, we also aggregated deciduous and evergreen plant groups for two climatic zones: boreal-temperate and tropical.

Meta-analysis: iWUE to atmospheric CO2 trend from leaf δ13C studies

Published (1921) and unpublished angiosperm leaf δ13C data collected between 1981 and 2005 were used for meta-analysis. Year of data collection was added to the collated dataset based on original publications. Any data source without collection dates was assumed to be 2 years before the date of paper submission (~5% of datasets). Atmospheric CO2 concentration and δ13Cair information corresponding to collection year were obtained from a published instrumental dataset (1980–2015) at the Mauna Loa station (3436). For δ13C values without environmental data, we obtained MAT and MAP data from 0.5° × 0.5° resolution CRU TS v. 4.0 (30) gridded dataset. The final dataset includes 1523 species site points from 76 studies of 1000 species across eight biomes. To quantify the response of deciduous and evergreen leaves to elevated CO2, we used a linear model with iWUE as the dependent variable and atmospheric CO2 with interaction between deciduous and evergreen groups. The iWUE trend along rising atmospheric CO2 gradient across collective leaf samples from different studies in various localities may be influenced by environmental conditions of the location. To investigate the likely influential environment factor that may have contributed to the observed iWUE trend, we quantified the amount of variation contributed by atmospheric CO2 concentration, MAT, MAP, altitude, and latitude across time. We first regressed collection year against all the foregoing environmental variables and then used R package relaimpo (67) to quantify the amount of variation contributed by each environmental factor. The proportion of variance explained by the model was 99.3%, of which 98% was contributed by CO2 followed by MAT at ~1%. Therefore, we can be confident that CO2 was influential in driving iWUE trends across collection time compared with other environmental variables. [Emphasis added] We designated the iWUE gain across collective leaf samples of different species and environmental conditions/locations as ΔiWUEc to differentiate it from ΔiWUE. The latter is derived from iWUE gain of the same species composition and locality.


Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/5/12/eaax7906/DC1

Fig. S1. Historical and contemporary leaf functional trait plots through the origin.

Fig. S2. iWUE gain (∆iWUE) of deciduous and evergreen plants in biomes for growth habit, arranged by increasing MAT.

Fig. S3. iWUE gain (∆iWUE) of deciduous and evergreen plants in biomes for habitat group, arranged by increasing MAT.

Fig. S4. The changes in the ratio of leaf intercellular (ci) to ambient CO2 (ca), Δci/ca, in evergreens and deciduous species in biomes, arranged by increasing MAT.

Fig. S5. iWUE change (∆iWUE) of deciduous and evergreen plants versus MAT change (∆MAT) and VPD change (∆VPD) in biome growth habit and habitat group.

Fig. S6. Scatter plot of Nmass versus MAT for combined historical and contemporary samples of evergreen and deciduous plants.

Fig. S7. Trend of iWUE from tree ring data along increasing atmospheric CO2 concentration between the years 1970 and 2013.

Fig. S8. Evergreen and deciduous iWUE plotted against atmospheric CO2 concentration showing slope of response.

Fig. S9. Kernel density plots of leaf life span (month) of deciduous and evergreen plants in the boreal-temperate and tropical biomes.

Table S1. List of species studied, their leaf habit (evergreen, deciduous), habitat (understory subcanopy and open canopy), and growth habit (shrub and tree).

Table S2. Summary of historical and contemporary site location, vegetation type, and collection date in alphabetical order by biome and site name.

Table S3. Historical and contemporary samples showing average LMA in evergreen and deciduous group within biome and probability of evergreen LMA larger than deciduous LMA, P* = P(LMAevergreen > LMAdeciduous).

Table S4. Average iWUE change (ΔiWUE) in biome between two time points 1988–1991 and 2013–2015 with CI95% from posterior distributions in Bayesian analysis.

Table S5. Average iWUE gain (ΔiWUE) in evergreen and deciduous plants within biome with CI95% from posterior distributions in Bayesian analysis.

Table S6. Shrub and tree, average iWUE gain (ΔiWUE) in evergreen and deciduous plants within biome, with CI95% from posterior distributions in Bayesian analysis.

Table S7. Understory-subcanopy and open-canopy habitat, average iWUE gain (ΔiWUE) in evergreen and deciduous plants within biome, with CI95% from posterior distributions in Bayesian analysis.

Table S8. Average annual air temperature change and average annual VPD change of biomes between two time periods 1988–1991 and 2013–2015 with CI95% from posterior distributions in Bayesian analysis.

Table S9. Average of coefficients of Model 1 and Model 2 with CI95% from posterior distributions in Bayesian analysis.

Table S10. Slope of iWUE response to atmospheric CO2 concentration (ΔiWUE/ΔCO2) for individual trees arranged by leaf habit, species, and references.

Table S11. Pearson correlation matrix (lower half panel in gray) and significance (upper half panel) between iWUE, VPD, precipitation, temperature, altitude, and latitude.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


  1. 1.
  2. 2.
  3. 3.
  4. 4.
  5. 5.
  6. 6.
  7. 7.
  8. 8.
  9. 9.
  10. 10.
  11. 11.
  12. 12.
  13. 13.
  14. 14.
  15. 15.
  16. 16.
  17. 17.
  18. 18.
  19. 19.
  20. 20.
  21. 21.
  22. 22.
  23. 23.
  24. 24.
  25. 25.
  26. 26.
  27. 27.
  28. 28.
  29. 29.
  30. 30.
  31. 31.
  32. 32.
  33. 33.
  34. 34.
  35. 35.
  36. 36.
  37. 37.
  38. 38.
  39. 39.
  40. 40.
  41. 41.
  42. 42.
  43. 43.
  44. 44.
  45. 45.
  46. 46.
  47. 47.
  48. 48.
  49. 49.
  50. 50.
  51. 51.
  52. 52.
  53. 53.
  54. 54.
  55. 55.
  56. 56.
  57. 57.
  58. 58.
  59. 59.
  60. 60.
  61. 61.
  62. 62.
  63. 63.
  64. 64.
  65. 65.
  66. 66.
  67. 67.
Acknowledgments: We are grateful to S. Wing and staff at Smithsonian NMNH for the hospitality and for access to herbarium specimens and the loan of leaves from the CLAMP collection. We are also grateful to the following people and institutions for permissions and field assistance: Smithsonian Environmental Research Center, Maryland, USA (P. Megonigal, S. McMahon, and J. Shue), Jasper Ridge Biological Preserve, California, USA (N. Chiariello and T. Corelli); The University of the South Pacific, Fiji (M. Tuiwawa, A. Naikatini, and S. Pene), Tonto National Forest (E. Hoskins and C. Denton), California State Parks (T. Hyland and J. Kerbavaz), Alaska State Parks (P. Russell and L. Ess), and Oregon State Parks (N. Bacheller). Many thanks to S. Culhane, E. Doyle, and C. Egan for field assistance. Funding: We gratefully acknowledge funding from a Science Foundation Ireland (SFI) Principal Investigator Award (PI) 11/PI/1103. A.P. was supported by SFI Career Development Award grant 17/CDA/4695 and SFI center grant SFI/12/RC/2289_P2. R.A.S. was supported by a Natural Environment Research Council grant (no. NE/P013805/1) and an XTBG International Fellowship for Visiting Scientists. Author contributions: W.K.S. led the writing, with input from J.C.M., C.Y., and M.M. J.C.M., C.Y., M.M., I.J.W., A.P., R.A.S., T.L., and R.C. discussed and commented on the manuscript. W.K.S., M.M., C.Y., and J.C.M. designed the study and organized and conducted fieldworks. W.K.S. and M.M. sampled CLAMP historical herbarium samples and curated all leaf samples. W.K.S. contributed to the LMA, Nmass, and δ13C data. C.Y. and W.K.S. contributed to the IRGA experiment data. C.Y. processed the IRGA experiment data. W.K.S. and A.P. performed the statistical analysis. W.K.S. conducted meta-analysis for published tree ring and leaf δ13C data. I.J.W. contributed leaf δ13C data for meta-analysis. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.
This article appeared on the Science Advances website at https://advances.sciencemag.org/content/5/12/eaax7906.full]]>

Subscribe to Our Informative Weekly Newsletter Here:

  • This field is for validation purposes and should be left unchanged.