# Global distribution and changes of leaf-level intrinsic water use efficiency and their responses to water stress

Xiang Wang, Zheng Fu, Philippe Ciais, Lixin Wang, Nina Buchmann, Trevor F Keenan, Martin de Kauwe, Josep Peñuelas, Guo Chen, Xiaoying Gong, et al.

► To cite this version:

Xiang Wang, Zheng Fu, Philippe Ciais, Lixin Wang, Nina Buchmann, et al.. Global distribution and changes of leaf-level intrinsic water use efficiency and their responses to water stress. Nature Communications, 2026, 17 (1), pp.1530. 10.1038/s41467-025-68252-9. hal-05523722

HAL Id: hal-05523722

https://hal.science/hal-05523722v1

Submitted on 23 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# Global distribution and changes of leaf-level intrinsic water use efficiency and their responses to water stress

Received: 28 January 2025

Accepted: 22 December 2025

Published online: 07 January 2026

<figure><img src="image_3.png" /></figure>

Check for updates

Xiang Wang<sup>1</sup>, Zheng Fu<sup>1</sup>✉, Philippe Ciais<sup>2</sup>, Lixin Wang<sup>3</sup>, Nina Buchmann<sup>4</sup>, Trevor F. Keenan<sup>5,6</sup>, Martin De Kauwe<sup>7</sup>, Josep Peñuelas<sup>8,9</sup>, Guo Chen<sup>10</sup>, Xiaoying Gong<sup>11</sup>, Jingfeng Xiao<sup>12</sup>, Xing Li<sup>13</sup>, Qiaoyun Xie<sup>14</sup>, Paul C. Stoy<sup>15</sup>, David Makowski<sup>16</sup>, William K. Smith<sup>17</sup>, Han Wang<sup>18</sup>, Songhan Wang<sup>19</sup>, Fangyue Zhang<sup>20</sup> & Shuli Niu<sup>1</sup>

Intrinsic water use efficiency (iWUE) at the leaf level measures water expenditures by terrestrial plants during photosynthesis, yet its global spatiotemporal dynamics and responses to water stress remain poorly understood. Using machine-learning models and carbon isotope observations in C₃ foliage, here we elucidate global patterns, trends, and water-stress responses of leaf iWUE. We find high iWUE in cold, arid regions and lower values in warm, humid areas. From 2001 to 2020, global iWUE increases at 0.2 ± 0.02 µmol mol⁻¹ year⁻¹, with strong biome specific differences. Grasslands exhibit the highest mean iWUE but the slowest increase, whereas evergreen broadleaf forests show the lowest iWUE yet the fastest increase. iWUE rises with increasing water stress, but the rate of growth diminishes as water stress intensifies. Vapor pressure deficit influence iWUE more broadly than soil moisture. The ecological optimality model reproduces the spatial patterns of leaf iWUE and identifies vapor pressure deficit as the dominant driver, but overestimates mean iWUE and its trend. Our findings suggest that increasing water stress may slow the rate of global iWUE increase as the climate continues to warm.

Climate warming and rising atmospheric CO₂ concentrations (Cₐ) change the coupling of carbon and water in terrestrial ecosystems¹,². Intrinsic water use efficiency (iWUE), the ratio of leaf net photosynthetic rate to stomatal conductance (i.e., A/g₅), is a key indicator for assessing plant responses to climate change³⁻⁵. Improved knowledge of plant iWUE is critical to a more mechanistic understanding of land-atmosphere interactions, plant function, and ecosystem vulnerability to water stress.

In C₃ plant photosynthesis, leaf ¹³C discrimination (Δ¹³C) is linearly related to leaf iWUE⁶⁻⁹. Yet, despite providing key information for iWUE, plant carbon isotope measurements are often derived from single-site studies. Such a reliance on single-location observations

introduces substantial uncertainties when isotopic data are used at larger scales, such as in global models². Single-site datasets fail to adequately represent the diversity of Earth's ecoregions and pose challenges in translating leaf or tree-ring isotope proxies into representations of entire ecosystems¹⁰. For example, using tree-ring isotope data to infer iWUE gives a long time series but does not cover non-forested areas¹¹.

Leaf-level iWUE is influenced mainly by atmospheric CO₂, vapor pressure deficit (VPD), and soil moisture (SM). Studies using free-air CO₂ enrichment experiment¹² and gas-exchange methods¹³ have shown that rising Cₐ reduces leaf stomatal conductance (g₅), thereby promoting an increase in iWUE. Recently, an analysis of global tree-ring

<figure><img src="image_4.png" /></figure>
carbon isotope data suggested that while elevated C_a has driven a substantial increase in forest iWUE over the 20th century, the rate of change in iWUE relative to C_a has decreased by half⁴. Additionally, both theory and experimental studies suggest that increased atmospheric and soil drought may enhance iWUE¹¹,¹⁵,¹⁶. Consistent with theory, increased VPD and reduced SM have been found to reduce leaf g_s⁶,¹⁷, resulting in higher C_i/C_a and higher iWUE, with these effects becoming more pronounced over time¹⁸⁻²⁰. Despite advancements in understanding how rising C_a influences iWUE, there remains a gap in spatiotemporal insights into how global leaf-level iWUE responds to atmospheric dryness and soil dryness.

To quantify the global spatiotemporal distribution of iWUE, we compile carbon isotope observations from 5964 globally-distributed C₃ plant foliage samples²¹, spanning the time period 2001–2020 while covering all continents except Antarctica, and construct global maps of carbon isotope discrimination (Δ¹³C) in C₃ plants during 2001–2020 using a random forest explainable machine learning models ensemble (Methods). The resulting annual maps of Δ¹³C are then translated into maps of iWUE for C₃ plants using Ma's²² and Yu's²³ equations to calculate leaf-level iWUE, which incorporate the effect of mesophyll resistance g_m and nonphotosynthetic fractionations (“Methods”). Based on maps of leaf iWUE, we further investigate the influence and control regions of atmospheric dryness versus soil dryness on global iWUE. Finally, we test how ecological optimality theory²⁴,²⁵ encapsulated in the P-model²⁶ simulates global leaf-level iWUE distribution and trends, and discuss differences between optimality theory and observation-based estimates.

## Results and discussion

### Spatiotemporal dynamics of global leaf-level iWUE

We quantified the spatial distribution of iWUE in C₃ plants using Δ¹³C of natural C₃ plant foliage and well-established equations (Methods). Compared to previous studies based on in situ leaf or tree-ring isotope measurements—which are geographically sparse or biased toward forest ecosystems (unrepresentative of non-forest ecosystems¹¹)—our study provides a spatially continuous, observation-constrained framework for estimating leaf-level iWUE across biomes. This approach addresses key limitations in spatial coverage and scaling, enabling cross-biome estimation of iWUE, and offering broader applicability for large-scale analyses of carbon-water coupling. From 2001–2020, we estimated a global mean iWUE of 55.4 ± 1.6 µmol mol⁻¹ (mean ± standard deviation) (Fig. 1a, b). Spatially, iWUE values were lower in humid regions (mean: 52.5 µmol mol⁻¹) and higher in arid and cold regions (mean: 60.4 µmol mol⁻¹, Fig. 1a, b). This finding aligns with ecophysiological trade-off mechanisms through which plant foliage in cold and arid regions are subject to cold or water stress and operates at reduced g_s. This leads to declines in both photosynthetic and transpiration rates¹¹,²⁷,²⁸. The high iWUE observed in high-latitude boreal regions (cold) and across the Australian continent (dry) (Fig. 1a, b) may be due to the fact that the decline was stronger for transpiration than for photosynthesis²⁹,³⁰. Conversely, in warm and humid regions, higher water availability in soils increases air humidity and SM, leading to increased leaf g_s³¹⁻³³. With elevated g_s, leaf C_i and transpiration rates rise, resulting in a decreased leaf iWUE³⁰,³⁴,³⁵, as observed in tropical rainforests regions such as Southeast Asia, central Africa, and the Amazon (Fig. 1a, b). We also compared leaf-level iWUE values across plant functional types (PFTs) (Supplementary Fig. S1). We found that iWUE in non-forests tended to be higher than that in forests. Within forest types, deciduous needleleaf forests (DNF) had the largest iWUE, followed by evergreen needleleaf forests (ENF), deciduous broadleaf forests (DBF), and evergreen broadleaf forests (EBF) (DNF > ENF > DBF > EBF, Supplementary Fig. S1). For the non-forest types, iWUE was greater in C₃ grasslands than shrublands and savannas.

Global leaf-level iWUE had significantly increased over the past two decades (Fig. 1c–i). More than half of the C₃ vegetation-dominated

regions of the world (e.g., cold and arid regions of North America, Amazon and Congo rainforests, and the Far East region of Russia) showed a significant increase in iWUE between 2001 and 2020 (78.9% of the grid cells, *p* < 0.05, Fig. 1c–e). Only a small percentage of regions showed a significant decrease (0.9% of the grid cells, *p* < 0.05). We found a rate of increase of 0.20 ± 0.02 µmol mol⁻¹ year⁻¹ per year (*p* < 0.001) (Fig. 1i). This trend was similar to the calculated trend based on independent tree-ring isotope time series from 1966 to 2000 (0.20 ± 0.004 µmol mol⁻¹ year⁻¹)³⁶. However, this trend was significantly lower than that estimated by the JULES model³⁷ (0.28 ± 0.01 µmol mol⁻¹ year⁻¹, 1979–2016), representing a difference of about 40% compared to the modeled trend. There is a widespread consensus that the sustained increase in global WUE at various scales is primarily driven by the rise in C_a, i.e., the “CO₂ fertilization effect”²⁴,²⁷,³⁸⁻⁴⁰. This is supported by plant physiological studies indicating that rising C_a promotes CO₂ uptake and enhances photosynthesis in vegetation⁴¹, while leading to a reduction in g_s, which in turn reduces transpiration³⁹. This was illustrated by the 12% increase in C_a between 2001 and 2020 (i.e., from 370 ppm to 414 ppm), corresponding to the rise in iWUE of global C₃ plants of about 8.5 ± 0.11%. However, atmospheric nitrogen (N) deposition, showing large spatiotemporal trends globally⁴², was also shown to increase WUE of temperate forests (based on stable isotope data in tree-rings)⁴³, interacting with increasing C_a and global warming.

The increasing trend of leaf-level iWUE was observed in all C₃ biomes, both in terms of the rate of increase (Fig. 1g) and the proportion of growth (Fig. 1h). We found that the rate of iWUE growth in the five forest biomes was generally higher than in non-forests types. Broadleaf forests (often classified as angiosperms) in particular had the highest rate of increase (EBF iWUE: 0.30 ± 0.04 µmol mol⁻¹ year⁻¹, DBF iWUE: 0.29 ± 0.04 µmol mol⁻¹ year⁻¹), which was 40% larger than that of conifers (ENF iWUE: 0.18 ± 0.03 µmol mol⁻¹ year⁻¹, DNF iWUE: 0.24 ± 0.02 µmol mol⁻¹ year⁻¹). In contrast, C₃ grasses had a relatively slower iWUE growth, with grasslands increasing by only 0.15 ± 0.02 µmol mol⁻¹ year⁻¹ (Fig. 1g). Previous indoor experimental studies have shown that angiosperm stomata are more sensitive to elevated C_a compared to conifers and other PFTs. As a result, angiosperm g_s tends to be more sensitive to rising C_a, leading to a more rapid increase of iWUE⁴⁴. Our findings aligned with previous studies by Saurer et al.⁴⁵ and Wang et al.⁴⁶, which observed similar iWUE trends based on Northern Hemisphere tree-ring and flux tower observations, respectively. However, they differed from the results of Frank et al.⁴ and Mathias et al.⁴⁷, who found that conifers exhibited a higher rate of iWUE growth than broadleaves during the 20th century based on tree-ring data. It is worth noting that discrepancies in absolute iWUE trends between these studies and ours may arise from differences in location (e.g., climate, N deposition), spatial scales (e.g., in situ observations versus global gridded data), study period, and methodologies (e.g., foliage, tree rings, flux towers).

### Impacts of soil dryness versus atmospheric dryness on leaf-level iWUE

Our global leaf-level iWUE maps provide a crucial understanding of how water stress impacts plant functioning. We examined how means and trends of iWUE and their sensitivity to C_a varied with SM, VPD and aridity index gradients worldwide. Our findings revealed that iWUE generally increased as VPD rose and SM decreased (Fig. 2a), overall being, consistent with established theories on plant leaf ecophysiology from field observations or experiments¹⁷,⁴⁸,⁴⁹. Previous studies suggested that water stress-induced increases of iWUE are primarily driven by changes in g_s, as plants prioritize water conservation over carbon uptake. Stomatal closure reduces transpiration more than it reduces photosynthesis¹¹, thereby raising iWUE.

Importantly, however, the rates of increase in global leaf-level iWUE from 2001 to 2020 were generally higher in regions with higher
<figure><img src="image_5.png" /><figcaption>**Fig. 1 |** Global distribution and trends of leaf-level iWUE for C₃ plants during 2001–2020. Spatial distribution of mean iWUE (a). Latitudinal patterns of mean iWUE, with solid lines indicating mean values and shaded areas representing 95% confidence intervals (b). Spatial distribution of iWUE trends (c). Areas with significant iWUE trend (Mann–Kendall test, *p* < 0.05) (d). Proportion of significant trend changes in iWUE (e). Latitudinal patterns of iWUE trends, with solid lines indicating mean values and shaded areas representing 95% confidence intervals (f). Trends in iWUE for eight plant functional types (PFTs), with statistical significance assessed using two-sided *t* tests (***, *p* < 0.001, g). Proportion of change in iWUE</figcaption></figure>

trends for eight PFTs (h). Global iWUE anomalies, with statistical significance assessed using two-sided *t* tests (i). PFTs include ENF (evergreen needleleaf forests), EBF (evergreen broadleaf forests), DNF (deciduous needleleaf forests), DBF (deciduous broadleaf forests), MF (mixed forests), Shrublands, Savannas and Grasslands. N: negative trend; N*: significant negative trend; P: positive trend; P*: significant positive trend. Non-natural vegetation areas and deserts with less than 10% vegetation coverage were masked, as were grid units dominated by C₄ vegetation.

SM and lower VPD (humid regions) compared to arid regions (Fig. 2b). Interestingly, the sensitivity of iWUE to Cₐ (diWUE/dCₐ) followed the same pattern, with diWUE/dCₐ being higher in humid regions and lower in arid ones (Fig. 2c). This pattern was consistent when mean and trend values of iWUE and diWUE/dCₐ were analyzed using the aridity index (Fig. 2d–f). In particular, both iWUE trends and diWUE/dCₐ as a whole increased with decreasing aridity (i.e., high aridity index; Fig. 2e, f). However, this was in stark contrast to the results of a regional study

that concluded that the frequency of extreme heat waves and water stress events resulted in faster iWUE growth in arid ecosystems in the southwestern U.S. than in other parts of the globe¹¹. In fact, a number of factors influence iWUE trends across climate zones, especially their interactions with water stress⁴⁷. Several tree-ring studies showed that the positive effect of rising Cₐ on iWUE was attenuated by increasing water stress¹⁴,⁴⁷,⁵⁰,⁵¹. Therefore, plant iWUE in the humid zone benefited more from the increase in Cₐ than plants in the arid zone, resulting in a
<figure><img src="image_6.png" /><figcaption>**Fig. 2 | Response of iWUE and diWUE/dC<sub>a</sub> to water stress.** The response of mean values of iWUE to vapor pressure deficit (VPD) and soil moisture (SM) (a). The response of iWUE trends to VPD and SM (b). The response of diWUE/dC<sub>a</sub> (the sensitivity of iWUE to C<sub>a</sub>) to VPD and SM (c). Mean iWUE along the aridity index (d;</figcaption></figure>

lower index indicates drier conditions). iWUE trends along the aridity index (e). diWUE/dC<sub>a</sub> along the aridity index (f). The solid lines indicate the mean values, and the shades indicate the 95% confidence interval calculated in each aridity bin through bootstrapping (n = 5000).

faster rate of increase in iWUE in this region. This was also verified by our analysis, with the sensitivity of iWUE to C<sub>a</sub> decreasing gradually with increasing water stress (Fig. 2c, f), which also implied that the rate of increase in iWUE will decrease in more regions in the context of rapid global warming and aridification.

We also identified regions where spatiotemporal changes in iWUE were primarily controlled by either SM or VPD, using partial correlation analyses ("Methods"). Globally, the regions more affected by VPD were mainly located in the southwestern United States, the Amazon and southern South America, the Congo Basin in Africa, Siberia in Russia, and eastern Australia, while the regions affected by SM were more scattered (Fig. 3a). Overall, our results highlighted that VPD-induced atmospheric dryness played a dominant role in driving spatiotemporal changes in leaf water and carbon coupling in C<sub>3</sub> vegetation across much of the globe (Fig. 3a). Specifically, VPD dominated iWUE changes across a significantly larger area (65%) compared to SM (35%) (Fig. 3a). Partial correlation analysis based on residual VPD and SM, with temperature linear and non-linear effects on VPD and SM removed, also confirmed this result, identifying VPD and SM as the dominant drivers in 59% and 41% of C<sub>3</sub> vegetation areas, respectively (Supplementary Fig. S2). Previous studies have indicated that VPD had a more substantial effect on vegetation carbon and water cycling compared to SM<sup>52-54</sup> because VPD regulates the balance between water supply and demand by controlling plant leaf water potential, thereby directly influencing water demand<sup>55</sup>. Similarly, the impact area of VPD on iWUE change remained greater than those of SM, both in all PFTs and in different climatic zones (Fig. 3c, d). This was in line with recent observations showing a sustained increase in VPD, which drove atmospheric moisture demand and contributed to rising iWUE by triggering stomatal closure<sup>15,16</sup>. This trend of VPD-induced iWUE increases has been observed in various ecosystems, including arid shrublands in the southwestern United States (via leaf carbon

isotopes) and forests in Europe (via tree-ring carbon isotopes)<sup>11,15,56</sup>. Overall, our findings suggest that changes in iWUE for global C<sub>3</sub> plants were more strongly driven by VPD than by SM (Fig. 3).

## Comparison with P-model

Last, we also investigated whether the P-model developed based on ecological optimality theory could reproduce the observed spatiotemporal dynamics of global leaf-level iWUE and its response to water stress ("Methods"). We found that overall, the P-model successfully captured the global spatial distribution pattern of iWUE and its variation across different PFTs when compared to observation-based results (Fig. 4a and Supplementary Fig. S3a). However, the model overestimated global iWUE values by about 23% (Fig. 4b, c and Supplementary Fig. S3b). Trend analysis revealed an unrealistically large and significant increase in iWUE<sub>mod</sub> (98% of the global C<sub>3</sub> covered area) (Supplementary Figs. S3c, d, S4). The globally increasing rate of iWUE<sub>mod</sub> estimated by the P-model was 0.49 ± 0.02 µmol mol<sup>-1</sup> year<sup>-1</sup>, which was much higher than the rate derived from carbon isotope observations. This mismatch was partially explained by the model's overestimation of the sensitivity of iWUE<sub>mod</sub> to C<sub>a</sub> changes (Supplementary Fig. S5). Despite this, the model's response of iWUE<sub>mod</sub> to variations in VPD and SM gradients generally aligned with observed patterns, showing that water stress led to an increase in iWUE<sub>mod</sub> (Supplementary Fig. S6). Moreover, model-based analysis also indicated that VPD played a more important role than SM in influencing iWUE<sub>mod</sub> changes in most regions worldwide. However, the area dominated by VPD in the P-model was significantly larger than that observed in the isotope data. Therefore, our findings suggest that the P-model may overestimate iWUE patterns, trends, and responses to atmospheric dryness compared to what is inferred from leaf carbon isotope data.

The core of the P-model is to minimize the total costs of carboxylation and transpiration by predicting an optimal C<sub>i</sub>/C<sub>a</sub> ratio in plant
<figure><img src="image_7.png" /><figcaption>**Fig. 3 | Effects of vapor pressure deficit (VPD) and soil moisture (SM) on iWUE.**
Spatial distribution of partial correlation coefficients of VPD and SM controlling changes in iWUE (a). Proportion of the area where iWUE is predominantly controlled by VPD or SM (b). Proportion of the area controlled by VPD or SM across eight plant functional types (PFTs, c). Proportion of the area controlled by VPD or SM along an aridity gradient (d). PFTs include ENF (evergreen needleleaf forests), EBF (evergreen broadleaf forests), DNF (deciduous needleleaf forests), DBF (deciduous broadleaf forests), MF (mixed forests), Shrublands, Savannas and Grasslands. Non-natural vegetation areas and deserts with less than 10% vegetation coverage were masked, as were grid units dominated by C₄ vegetation.</figcaption></figure>

As such, the optimal $C_i/C_α$ ratio serves as both a benchmark for photosynthesis models based on the P-model and a key variable for estimating iWUE (Supplementary Fig. S15). Although the P-model performs well in predicting global vegetation photosynthesis$^{57-59}$, discrepancies arise when assessing carbon and water coupling at the leaf scale compared to our observational data. These differences, alongside with the P-model overestimation of iWUE sensitivity to rising $C_α$ (Supplementary Fig. S5), may be related to the use of fixed computational parameters. Calculated parameters with complex background information in the model (e.g., $β$ in $ξ = \sqrt{β(K + Γ*)}/1.6(1 + g_s/g_m)$, see Eq. 10) can also lead to biased results when characterized as fixed values. $β$ is the ratio of the unit cost of plant carboxylation to transpiration at 25 °C (dimensionless). $β$ represents a critical computational parameter of the P-model, which is directly related to the accurate description of changes in plant stomatal conductance and net photosynthetic rate in response to environmental gradients (i.e., the accurate quantification of $C_i/C_α$)$^{26}$.

In fact, there are large differences in $β$ values across different biomes; for instance, angiosperms have higher $β$ than gymnosperms, while deciduous plants have higher $β$ than evergreen trees$^{60}$. However,

a single $β$ value is currently used for all PFTs in the P-model developed based on the optimality theory$^{26,58}$. The use of a single value introduces prediction bias, as demonstrated by a study on plant leaf traits along an altitudinal gradient based on optimality theory, which found that a single $β$ was the main cause of the $C_i/C_α$ uncertainty$^{61}$. It was observed that the use of $β$ estimated from compiled data to predict $C_i/C_α$ resulted in a mismatch with observations. Predicted $C_i/C_α$ values were generally lower than the median observed values at most sites$^{61}$. Furthermore, our sensitivity analysis of iWUE to $β$ showed that iWUE decreased with increasing $β$ (Fig. 5). This suggests that the default $β$ value used in the P-model may be too low for most PFTs or climatic zones. Meanwhile, the model results showed that atmospheric drought dominated iWUE$_{mod}$ changes in more than 75% of the global $C_3$ vegetation cover, which may be caused by the fact that the P-model only considered the effect of atmospheric water demand on $C_i/C_α$ and ignored the effect of soil water stress. Previous studies has shown that $β$ decreases with increasing soil dryness, indicating that the cost of plant water transport increases as SM declines$^{60}$. This bias highlights the need to carefully evaluate Earth system models that adopt similar assumptions and parameterizations, as they may misrepresent future
<figure><img src="image_8.png" /></figure>

**Fig. 4 | Global distribution of intrinsic water use efficiency (iWUE<sub>mod</sub>) estimated by the P-model.** Spatial distribution of mean iWUE<sub>mod</sub> (a). Spatial distribution of mean differences between iWUE<sub>obs</sub> and iWUE<sub>mod</sub> (b). Mean differences between iWUE<sub>com</sub> and iWUE<sub>mod</sub> across eight plant functional types (PFTs) (c). Box plots display means (white dots), medians (horizontal lines), the 25th and 75th percentiles (box edges), and the 5th and 95th percentiles (whiskers). PFTs include ENF (evergreen needleleaf forests, *n* = 124731), EBF (evergreen broadleaf forests, *n* = 404627), DNF (deciduous needleleaf forests, *n* = 13949), DBF (deciduous broadleaf forests, *n* = 107884), MF (mixed forests, *n* = 258565), SHR (shrublands, *n* = 689230), SAV (savannas, *n* = 1102465) and GRA (grasslands, *n* = 927489). *n* denotes the number of image pixels for each PFT. Non-natural vegetation areas and deserts with less than 10% vegetation coverage were masked, as were grid units dominated by C<sub>4</sub> vegetation.

carbon water feedbacks under climate change. Despite the bias in the P-model's estimation of iWUE at large scales with physiological significance, we believe that the model's performance can be improved by incorporating the observed isotope data. This would enable better estimation of the vegetation-atmosphere energy, water and carbon cycles.

Vegetation plays a critical role in removing CO₂ from the atmosphere through photosynthesis, a process inherently linked to water vapor loss. The iWUE, which measures the rate of carbon uptake per unit of water lost, serves as a key characteristic of ecosystem function and is central to the global cycles of water, energy, and carbon. However, the inconsistent application of iWUE across different scales and the lack of a comprehensive global measure of leaf-level iWUE has limited our capacity to characterize and understand the

spatiotemporal dynamics of iWUE globally. In this study, we constructed carbon isoscapes for C₃ plant leaves and utilized mechanistic equations to produce an assessment of the spatiotemporal dynamics in global leaf-level iWUE. We found that mean leaf iWUE is higher in cold and arid regions, and lower in warm and humid areas. Grasslands exhibited the highest iWUE, whereas evergreen broadleaf forests showed the lowest. Over the past two decades, global iWUE has increased significantly at a rate of 0.2 ± 0.02 µmol mol⁻¹ year⁻¹. Grasslands showed the slowest iWUE growth rate (0.15 ± 0.02 µmol mol⁻¹ year⁻¹), whereas evergreen broadleaf forests exhibited the fastest growth (0.3 ± 0.04 µmol mol⁻¹ year⁻¹). We confirmed that global leaf iWUE increased with rising water stress and further demonstrated that atmospheric dryness drove iWUE changes across a larger area than soil dryness. The rate of increase in global leaf-level iWUE was faster in
<figure><img src="image_9.png" /><figcaption>Fig. 5 | Sensitivity analysis of iWUE<sub>mod</sub> to the parameter β in the P-model. Scenario I: β varied continuously from 0 to 500, with all other parameters held constant (a). Scenario II: β fixed at six discrete values, while other parameters varied continuously (b).</figcaption></figure>

humid areas than in arid areas. Although the spatial pattern of iWUE predicted by the P-model, based on optimality theory, was similar to our observations, the model estimated a higher mean and stronger trends of iWUE changes, highlighting uncertainties in its ability to predict coupled water-carbon responses to climate change. Collectively, these results present a framework for large-scale estimation of carbon-water coupling in terrestrial ecosystems and offer a benchmark for models predicting terrestrial carbon-water cycles. Our findings also underscore the potential for intensifying water stress to further reduce the rate of increase for global iWUE as the Earth system continues to warm.

## Methods

### Global leaf carbon isotope data

We used the published carbon isotope dataset of C₃ plant foliage worldwide²¹ (Supplementary Figs. S7, S8a). The Δ¹³C was then calculated from leaf carbon isotope ratios (δ¹³C<sub>Leaf</sub>, ‰) and atmospheric carbon isotope ratios (δ¹³C<sub>a</sub>, ‰, δ¹³C<sub>a</sub> corresponds to the leaf sample sampling year) (Eq. 1):

$$ \Delta^{13}C_{Leaf} = \frac{\delta^{13}C_a - \delta^{13}C_{Leaf}}{1 + \delta^{13}C_{Leaf}/1000} \quad (1) $$

The dataset collected and compiled 5964 datasets of global carbon isotope composition (δ¹³C) of C₃ plant foliage (Supplementary Fig. S7). Specifically, we searched the literature by screening the search engines of Web of Science, Google Scholar, and China Knowledge Network using the following keywords: “C₃ plant”, “leaf”, “carbon isotope” and “leaf δ¹³C”. The literature was selected based on the following criteria: (1) leaf samples were from natural plant plots, far from anthropogenic interventions; (2) photosynthesis type was C₃ plants; (3) leaf sampling information was well documented, including sampling latitude, longitude and plant type; and (4) leaf samples were obtained between 2001 and 2020. In addition, we screened citations in this literature for other carbon isotope data from C₃ plant leaves that may not have been identified by the search engine. A total of 116 papers were selected based on these criteria. When leaf δ¹³C data were presented in body tables or supplementary tables, we manually incorporated them into our database. When presented graphically, we manually digitized each individual leaf isotope using the digitization software GetData Graph Digitizer (version 2.25). The GetData Graph Digitizer has been widely applied in digitization processes²¹. Random tests in earlier studies showed no significant differences between the extracted values and the true values, indicating that the software had good accuracy⁶². In a few cases, when data could not be obtained from the literature, we contacted the original authors to obtain the data.

In this study, the predicted Δ¹³C isoscape map for C₃ plants was generated on an annual scale with a spatial resolution of 0.05°. This choice was guided by two considerations: first, most of the leaf samples collected from the literature lack precise sampling months, often providing only broad time ranges, making the annual scale a better match for the temporal accuracy of the data. Second, a spatial resolution of 0.05° provided sufficient detail on geographical and ecosystem variations at a global scale while maintaining a manageable data size, thereby reducing the demands on processing and storage and making large-scale data analysis more feasible.

### MODIS surface reflectance

We acquired the global 0.05° (approximately equal to 5600 m) surface reflectance (MOD09CMG) product of the Moderate Resolution Imaging Spectroradiometer (MODIS) observations for the years 2001–2020 from the Terra satellite platform, encompassing seven spectral bands ranging from 459 nm to 2155 nm, with a temporal resolution of days. All bands were atmospherically corrected (data version 6.1). We excluded poor-quality image elements based on quality assurance bands (QA), minimizing residual contamination from aerosols, snow, shadows, and clouds. Since the time scale of the constructed isoscapes mapping of Δ¹³C of C₃ plant leaves was interannual, we subjected the daily MOD09CMG data to annual average synthesis. The annual average synthesis helps to reduce the noise in the spectral bands on the one hand, and on the other hand, it can better capture the overall information of each spectral band. In addition, annual averaging can improve the processing efficiency of large-volume spectral data.

### CO₂ concentrations and δ¹³C

Due to the lack of a reliable and high-resolution global map of atmospheric CO₂ concentrations (Cₐ, ppm) and carbon isotope ratios (δ¹³Cₐ, ‰), we used monthly Cₐ and δ¹³Cₐ from the Mauna Loa station (19.54°N, 155.57°W, 3397 m a.s.l., United States), data provided by the National Oceanic and Atmospheric Administration (NOAA), for which annual averages were calculated. Subsequently, the annual averages were aggregated into raster data with a global 0.05° spatial resolution for the period 2001–2020 using the “raster” package in R.

### SIF

We also use the widely used globally gridded solar-induced chlorophyll fluorescence dataset (GOSIF, W/m²/µm/sr) and global spatially contiguous solar-induced chlorophyll fluorescence dataset (CSIF, W/m²/µm/sr), since SIF is closely related to plant photosynthesis. GOSIF was developed and made publicly available by Li and Xiao⁶³ in 2019 with a spatial resolution of 0.05°. CSIF was developed and made publicly available by Zhang et al.⁶⁴ with a spatial and temporal resolutions are
0.05° and 4-day, respectively. We selected GOSIF data with interannual resolution, where annual values of grid cells were obtained by averaging eight days of data within each year. The GOSIF data spanned from 2001 to 2020. The CSIF dataset is developed using discrete OCO-2 SIF retrievals, MODIS remote sensing data combined with neural network algorithms. The CSIF data spanned from 2001 to 2022, and data between 2001 and 2020 were selected for this study. The annual CSIF values were obtained by averaging four days of data within each year. Compared to the discrete OCO-2 SIF, GOSIF and CSIF offer finer spatial resolution, continuous global coverage, and a longer temporal record. Finally, we averaged GOSIF and CSIF over the corresponding times to obtain the SIF variable for each year.

## Geographic information

We obtained the geographic location of sample sites from the literature, including latitude, longitude, and elevation data. For altitude information of the foliage sample, original records were used if they were available in the literature. If original records of elevation information were not available, we instead extracted data from the GloE-lev_30as (ELE, m) product of the Food and Agriculture Organization of the United Nations Historical Soil Map and Database. For isoscapes construction, we converted the ELE of 30 arcsec to 0.05°.

## Climate data

The mean annual temperature (MAT, °C) and annual precipitation (PRE, mm) we used are bioclimatic variables, which were derived from monthly precipitation and temperature values to generate more biologically meaningful climate variables⁶⁵. We used monthly minimum temperature, monthly maximum temperature, and monthly precipitation data (temporal resolution of months, spatial resolution of 2.5 minutes, -0.042°) from 2001–2020 provided by WorldClim. Bioclimatic variables, including temperature, were generated using the “biovars” function of the R package “dismo”. The method for generating bioclimatic variables followed the same process as BIOCLIM in the ANUCLIM framework. BIOCLIM, originally designed by Nix⁶⁶ in 1986, derives bioclimatic parameters from the provided climate surfaces based on the known habitat locations of a specific plant or animal species. MAT and PRE were used to calculate the iWUE equation based on carbon isotope estimates. Water vapor pressure difference (VPD, kPa), soil moisture (SM, mm) and downward solar radiation flux at the surface (radiation, W m⁻²) used in this study were obtained from the TerraClimate dataset developed by Abatzoglou et al.⁶⁷, with a spatial resolution of 1/24° (-0.042°). The TerraClimate data take into account the influence of topographic effects. Annual VPD, SM and radiation were composited from monthly data by taking the mean. We resampled VPD, SM and radiation data to 0.05°.

To driving the P-model, required growing season temperatures were obtained from the ERA5 (WFDE5) dataset with spatial and temporal resolutions of 0.5° and hourly scales, respectively. Also, we selected the specific humidity of WFDE5 to calculate the growing season VPD according to the method of Buck⁶⁸. The growing season for this study was defined with the average number of all days with a mean daily temperature exceeding 0 °C, where daily data were obtained by averaging over a 24 h period, and 0.5° was resampled to 0.05°. Since the WFDE5 data were only updated to 2019, we used a linear interpolation method combining data from 2001 to 2019 to obtain growing season temperatures and VPDs for 2020.

## Land cover and aridity index data

To quantify differences in terrestrial carbon and water coupling among different PFTs, we used MODIS annual land cover data (MCD12C1) as a reference land cover map (2001–2020, 0.05° × 0.05°). We summarized the International Geosphere-Biosphere Program classification (IGBP) types into eight biomes (evergreen needleleaf forest (ENF), evergreen broadleaf forest (EBF), deciduous needleleaf forest (DNF), deciduous

broadleaf forest (DBF), mixed forest (MF), shrublands, savannas, and grasslands. Shrublands included closed shrublands, open shrublands, while other land cover types were masked (Supplementary Fig. S8a). Meanwhile, the C₄ vegetation percentage maps with a spatial resolution of 0.5° were resampled to 0.05°, and grids with > 50% C₄ vegetation were labeled as C₄ dominated areas and masked. Global C₄ distribution data were developed by Luo et al. (2024)⁶⁹ and span the period 2001-2019. For the missing year 2020, data from 2019 were used instead. Finally, the ratio of mean annual precipitation to mean annual potential evaporation was used to define the aridity index (AI), thereby identifying arid (AI < 0.2), semi-arid (0.2 ≤ AI < 0.5), sub-humid (0.5 ≤ AI < 0.65), and humid (AI ≥ 0.65) regions (Supplementary Fig. S8b).

## Construction of a predictive model for Δ¹³C isoscapes of C₃ plants

Currently, only the leaf carbon isotope discrimination (Δ¹³C) of plants with the C₃ photosynthetic metabolism pathway provides a physiologically meaningful and reliable estimate of leaf-scale iWUE, which is widely used in terrestrial ecosystem carbon and water coupling and cycling studies (see the relationship between Δ¹³C and iWUE presented in the next section). The leaf carbon isotopes characterize the trade-off between leaf photosynthesis and transpiration at the in situ level and need to be extrapolated to the regional or global scale (Supplementary Text 1 and Supplementary Text 2). To enable Δ¹³C isoscape mapping, we used the collected satellite spectral reflectance (seven spectral bands), atmospheric CO₂ concentration (Cₐ), vegetation physiology (solar induced chlorophyll fluorescence, SIF), and geographic information (longitude, latitude, elevation). We used information from the location of leaf samples and the year of sampling to extract annual means in seven spectral bands, annual means of atmospheric CO₂ concentrations and carbon isotopes, annual means of SIF, and elevation information. The δ¹³Cₐ of atmospheric CO₂ are used to calculate the Δ¹³C (Eq. 1), and the atmospheric CO₂ concentration was used to predict the Δ¹³C.

Before constructing isoscapes of leaf carbon isotopes in C₃ plants, we needed to filter the selected predictor variables. Because not all predictor variables are relevant for the construction of leaf Δ¹³C isoscapes, variable importance needs to be assessed. Redundancy among variables were reduced by selecting the minimum number of variables able to provide accurate predictions. This could reduce model computation time and potentially improve prediction accuracy⁷⁰. We used the vifstep function (multicollinearity test, Eq. 2) of the “usdm” R package⁷¹ and the recursive feature elimination (RFE, Eq. 3) of the “Caret” package⁷² to filter the optimal number of predictor variables. When the multicollinearity test showed that the variance inflation factor (VIF) of a predictor variable was greater than 10⁷³⁻⁷⁵, the variable was discarded because of its high correlation with other variables. RFE constructs a model based on the recursive removal of predictor variables and identifies which variables should be kept in the model to obtain a high level of prediction accuracy (Cross-validation via ten-fold). In the process of determining the predictor variables, we first performed a multicollinearity test on the predictor variables and then performed feature selection on the predictor variables that passed the multicollinearity test.

$$VIF = vifstep(B1, B2, B3, B4, B5, B6, B7, SIF, C_a, Lat, Long, ELE) \quad (2)$$

$$RFE = rfeControl(B1, B2, B3, B4, B5, B6, B7, SIF, C_a, Lat, Long, ELE) \quad (3)$$

where VIF is the predictor variable that passed the multicollinearity test; RFE is the predictor variable that passed the feature selection; B1-B7 are the reflectance variables in the seven spectral bands observed by MODIS; SIF is the solar-induced chlorophyll fluorescence variable; Cₐ is the atmospheric CO₂ concentration; Lat, Long, and ELE are latitude,
longitude, and elevation, respectively. These predictor variables are closely related to plant structure and function and can characterize the physiological and ecological information of plants. B3, B5, B7, GOSIF, Cₐ, Lat, Long and ELE were retained based on multicollinearity test (Supplementary Table 2) and RFE filtering (Supplementary Fig. S9).

Machine learning methods have been used extensively to upscale in situ observations to regional or global. Machine learning algorithms can consider nonlinear relationships between explanatory and response variables and are capable of integrating continuous and discrete variables²⁶. Here, we used the Random Forest algorithm to construct a global Δ¹³C isoscapes of C₃ plant leaves based on Wang et al.²¹ (Eq. 4).

$$ \Delta^{13}C_a = \text{Random Forest}(\text{B3}, \text{B5}, \text{B7}, \text{SIF}, C_a, \text{Lat}, \text{Long}, \text{ELE}) \quad (4) $$

Meanwhile, in order to achieve the Random Forest model to reach the optimal state, the hyperparameters of the model needed to be optimized iteratively. We optimized the two most critical hyperparameters of random forest, namely ntree (100 - 1000, divided into 10 levels with an interval of 100) and mtry (1 - 10, divided into 10 levels with an interval of 1), and set the random forest algorithm to go through 100 iterations (10 × 10). We determined the optimal hyperparameters based on the highest R² and the lowest RMSE value, and the iterative optimization showed that the best ntree and mtry values were 1000 and 2, respectively.

We used three evaluation schemes, (i) training set (70%)-test set (30%), (ii) ten-fold cross-validation, (iii) 100 random resampling, and four statistical criteria, namely, correlation coefficient (R), Nash Sutcliffe efficiency, root-mean-square error (RMSE), and mean absolute error (MAE), to evaluate the performance of the constructed Random Forest model for predicting the performance of the C₃ plant leaves for Δ¹³C isoscapes. Furthermore, to remove the interference of the C₄ photosynthetic pathway, we defined grid cells with pixel values greater than 50% on the C₄ vegetation cover map as dominated by C₄ plants and then masked these areas from the predicted Δ¹³C isoscapes of C₃ plant leaves globally²¹.

## Estimation of iWUE based on Δ¹³C of C₃ plant leaves

Intrinsic water use efficiency (iWUE) was defined as the ratio of net photosynthetic rate (A) and stomatal conductance (g₅), and the calculation process depended on the atmospheric CO₂ concentration (Cₐ) and the leaf intercellular CO₂ concentration (Cᵢ) (Eq. 5):

$$ iWUE = \frac{A}{g_s} = \frac{C_a - C_i}{1.6} = \frac{C_a(1 - \frac{C_i}{C_a})}{1.6} \quad (5) $$

Cᵢ/Cₐ can be estimated from leaf Δ¹³C (Eq. 6), and estimation methods based on leaf carbon isotope composition can directly obtain water use efficiency (intrinsic water use efficiency, iWUE), which was first developed by Farquhar et al.³⁴. It is worth noting that the Farquhar version of the discrimination equation ignored the effects of mesophyll conductance and non-photosynthetic fractionations (Supplementary Text 2). In order to improve the accuracy of the estimation, Ma et al.²² and Yu et al.²³ proposed an improved version of the formula (Eq. 7), by taking into account the parameters of mesophyll conductance and nonphotosynthetic fractionations:

$$ \frac{C_i}{C_a} = \frac{\Delta^{13}C_{Leaf} - a}{b - a} \quad (6) $$

$$ iWUE = \frac{C_a \cdot b - (\Delta^{13}C_{Leaf} - d_{Leaf}) - f \cdot \frac{\Gamma^*}{C_a}}{1.6 \cdot b - a + \frac{g_s}{g_m} \cdot (b - a_m)} \quad (7) $$

where 1.6 is the diffusion ratio of water vapor to CO₂ as it passes through leaf stomata; a and b are the fractionation coefficients for CO₂ as it

diffuses through leaf stomata (4.4 ‰) and Rubisco carboxylation process (29 ‰), respectively; f denotes the fractionations due to photosynthetic fractionation (11‰); Γ* is the carbon dioxide compensation point in the absence of mitochondrial respiration (Eq. 8)²⁷; and g₅/gₘ is the ratio of the stomatal conductance to the mesophyll conductance, Ma et al.²² found g₅/gₘ to be 0.79 ± 0.07 for species of different PFTs under wet versus water stress conditions based on a meta-analysis (Supplementary Text 2); dₐₗₐf (2.5‰) is the isotope correction of bulk leaf biomass to account for nonphotosynthetic effects²³; aₘ is the fractionation of CO₂ during solubilization and diffusion in chloroplasts (1.8‰).

$$ \Gamma^* = 42.7 + 1.68(T - 25) + 0.012(T - 25)^2 \quad (8) $$

where T (°C) is leaf temperature, which we represented as a calculated biometeorological variable of air temperature due to the lack of measured data.

To achieve iWUE of leaves with physiological significance globally, we borrowed the idea of biogeochemical process modeling and raster pixelized each computational parameter in iWUE according to the size of the study area (global scale) and to the set spatial and temporal resolution (spatial resolution of 0.05°, and synthesized values for each year spanning from 2001 to 2020). Constant parameter raster pixels contain unique values, while variable parameter raster pixels represent continuous values. It should be emphasized that an important step in realizing this work was the construction of isoscapes of Δ¹³C for C₃ plant leaves around the globe. The calculated parameters of the rasterization were then substituted into the equation calculations to obtain the iWUE of a large scale with physiological significance (Supplementary Text 2).

## Estimation of iWUE based on the P-model

Optimality theory assumes that plants influence leaf intercellular CO₂ concentration by regulating the degree of leaf stomatal opening and closing during long-term evolutionary and short-term growth stages to minimize the total cost of maintaining carboxylation and transpiration capacity²⁴⁻²⁶. Based on this theory, the iWUE equation parameter Cᵢ/Cₐ (considering the effect of g₅, substituting Cₐ/Cₐ for Cᵢ/Cₐ) can be estimated by the P-model developed by optimality theory²⁶, i.e., Eq. 9 and Eq. 10:

$$ \chi = \frac{C_i}{C_a} = \frac{C_c}{C_a} = \frac{\Gamma^*}{C_a} + (1 - \frac{\Gamma^*}{C_a}) \frac{\xi}{(\xi + \sqrt{VPD})} \quad (9) $$

$$ \xi = \sqrt{\frac{\beta(K + \Gamma^*)}{1.6(1 + \frac{g_s}{g_m})}} \quad (10) $$

where the units of all parameters in Eqs. 9 and 10 based on the P-model need to be standardized to Pa, and the final Cᵢ/Cₐ obtained is a dimensionless value. VPD is the growing season saturated water vapor pressure deficit; β is a dimensionless constant, estimated to be 343 (considering the effect of g₅/gₘ) based on Wang et al.²⁶; K is the effective Michaelis-Menten coefficient for Rubisco enzyme at a given atmospheric pressure and temperature; Γ* is the carbon dioxide compensation point; Wang et al.²⁶ used unpublished data from I.J. Wright to obtain a median g₅/gₘ of 1.4. The χ estimated based on the P-model (Supplementary Fig. S15) was consistent with the results of Wang et al.²⁶. The calculated χ (Cᵢ/Cₐ) was finally substituted into the iWUE calculation formula, i.e., Eq. 11 or Eq. 12:

$$ iWUE_{\text{mod}} = \frac{C_a(1 - \chi)}{1.6} \quad (11) $$

$$ iWUE_{\text{mod}} = \frac{C_a - \Gamma^*}{1.6(1 + \frac{\xi}{\sqrt{VPD}})} \quad (12) $$
where the Cₐ unit entered in iWUEₘ₉₅ is µmol/mol, so the iWUE unit obtained is µmol/mol.

## Estimating the sensitivity of iWUE to Cₐ

We used multiple linear regression⁷⁸ (Eq. 13) to calculate the sensitivity of iWUE to Cₐ, the rate of change of iWUE with respect to Cₐ, for each pixel, taking into account the effects of VPD, MAT, SM, and radiation:

$$iWUE = \frac{diWUE}{dC_a} \times Ca + \beta_1 \times VPD + \beta_2 \times MAT + \beta_3 \times SM + \beta_4 \times radiation + \varepsilon \quad (13)$$

where diWUE/dCₐ is the sensitivity of iWUE to Cₐ, β₁, β₂, β₃ and β₄ are the regression coefficients, and ε is the residual error term; Cₐ, VPD, MAT, SM, and radiation are annual averages.

## Statistical analyses

We used the Theil-Sen Median method at the raster pixel scale to quantify the spatial patterns of temporal trends for iWUE at the global scale from 2001 to 2020⁷⁹. The Mann-Kendall method was then used to test the significance of the trend, which indicates a significant trend change when p < 0.05. Both methods are nonparametric trend statistics and tests that have been widely used in the analysis of spatial time series data⁸⁰. We also quantified the latitudinal patterns of global annual mean iWUE. Trends in iWUE were calculated globally and over different PFTs regions using ordinary least squares.

We conducted a statistical analysis of the differences of mean iWUE for the eight biomes based on the 2020 land cover raster data. Because we analyzed data for global raster pixels, the parameter analysis data volume was large. We adopted the Kruskal-Wallis H-test and Dunn's test, which allowed us to determine whether there was a difference in median values between two or more biomes. These two tests are non-parametric and did not require the assumption that the data conformed to normal distribution or chi-square, and the Dunn test was realized by the R package "dunn.test". In addition, we calculated the proportion of changes in the eight biomes over the past two decades based on the 2020 land cover data.

We used partial correlation analysis to assess the effects of atmospheric dryness (VPD) and soil dryness (SM) on iWUE means. This method controlled for the influence of one variable while evaluating the partial relationship of the other with the response variable⁸¹. To identify the dominant driver of iWUE, we ranked VPD and SM based on the absolute values of their partial correlation coefficients in each grid cell⁸². In addition, given that the climatic effect of VPD inherently includes a temperature component, and that plants assimilation exhibit a nonlinear optimum response to temperature⁸³,⁸⁴, this may bias the attribution of iWUE in boreal regions. To address this, we further computed the correlations between iWUE and residuals of VPD and SM (obtained after removing both linear and non-linear influence of air temperature). Specifically, we first applied regression models including both linear and quadratic terms to remove temperature-driven components from VPD and SM, and then treated the residuals as temperature-independent indicators of water stress.

## Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

## Data availability

The data used in this study are openly available in the following databases. The MOD09CMG are downloaded from the LP DAAC (Land Processes Distributed Active Archive Center) of the USGS (United States Geological Survey) (https://lpdaac.usgs.gov/products/mod09cmgv061/). GOSIF and CSIF are obtained from https://doi.org/10.11888/Ecolo.tpdc.271751 and https://globalecology.unhed/

data.html, respectively. Atmospheric CO₂ concentration (Cₐ) and carbon isotope ratios (δ¹³Cₐ) are available on the World Data Center for Greenhouse Gases (WDCGG) (https://gaw.kishou.go.jp/). GloElev_30as are downloaded from https://www.fao.org/soils-portal/soil-survey/. MOD12C1 are from https://lpdaac.usgs.gov/products/mcd12c1v006/. C₄ Vegetation Percentage data are downloaded from https://zenodo.org/records/10516423. The aridity index is available from https://cgiarcsi community/data/global-aridity-and-pet-database/. WorldClim data are from https://www.worldclim.org/. VPD, SM, and radiation are obtained from https://www.climatologylab.org. ERA5 (WFDE5) data are from https://doi.org/10.24381/cds.20d54e34. The download links for parameters are provided in Supplementary Table 1. The data underlying the results are available at: https://doi.org/10.6084/m9.fishshare.27261399.

## Code availability

The primary code used to generate the results is publicly available at https://doi.org/10.6084/m9.fishshare.27261399.

## References

1. Ballantyne, A. P., Alden, C. B., Miller, J. B., Tans, P. P. & White, J. W. Increase in observed net carbon dioxide uptake by land and oceans during the past 50 years. *Nature* **488**, 70–72 (2012).

2. Adams, M. A., Buckley, T. N., Binkley, D., Neumann, M. & Turnbull, T. L. CO₂, nitrogen deposition and a discontinuous climate response drive water use efficiency in global forests. *Nat. Commun.* **12**, 5194 (2021).

3. Niu, S. et al. Water-use efficiency in response to climate change: from leaf to ecosystem in a temperate steppe. *Glob. Change Biol.* **17**, 1073–1082 (2011).

4. Frank, D. C. et al. Water-use efficiency and transpiration across European forests during the Anthropocene. *Nat. Clim. Chang.* **5**, 579–583 (2015).

5. Peñuelas, J., Canadell, J. G. & Ogaya, R. Increased water-use efficiency during the 20th century did not translate into enhanced tree growth. *Glob. Ecol. Biogeogr.* **20**, 597–608 (2011).

6. Farquhar, G. D., Ehleringer, J. R. & Hubick, K. T. Carbon Isotope Discrimination and Photosynthesis. *Annu. Rev. Plant Biol.* **40**, 503–537 (1989).

7. Wang, X. et al. Differences in the patterns and mechanisms of leaf and ecosystem-scale water use efficiencies on the Qinghai-Tibet Plateau. *Catena* **222**, 106874 (2023).

8. Impa, S. M. et al. Carbon isotope discrimination accurately reflects variability in WUE measured at a whole plant level in rice. *Crop Sci.* **45**, 2517–2522 (2005).

9. Petrík, P. et al. Linking stomatal size and density to water use efficiency and leaf carbon isotope ratio in juvenile and mature trees. *Physiol. Plant* **176**, e14619 (2024).

10. Seibt, U., Rajabi, A., Griffiths, H. & Berry, J. A. Carbon isotopes and water use efficiency: sense and sensitivity. *Oecologia* **155**, 441–454 (2008).

11. Kannenberg, S. A., Driscoll, A. W., Szejner, P., Anderegg, W. R. L. & Ehleringer, J. R. Rapid increases in shrubland and forest intrinsic water-use efficiency during an ongoing megadrought. *Proc. Natl. Acad. Sci. USA* **118**, e2118052118 (2021).

12. De Kauwe, M. G. et al. Forest water use and water use efficiency at elevated CO₂: a model-data intercomparison at two contrasting temperate forest FACE sites. *Glob. Change Biol.* **19**, 1759–1779 (2013).

13. Li, Q. et al. Differential influence of elevated CO₂ on gas exchange and water use efficiency of four indigenous shrub species distributed in different sandy environments in central Inner Mongolia. *Ecol. Res.* **33**, 863–871 (2018).

14. Adams, M. A., Buckley, T. N. & Turnbull, T. L. Diminishing CO₂-driven gains in water-use efficiency of global forests. *Nat. Clim. Chang.* **10**, 466–471 (2020).
15. Driscoll, A. W., Bitter, N. Q., Sandquist, D. R. & Ehleringer, J. R. Multidecadal records of intrinsic water-use efficiency in the desert shrub *Encelia farinosa* reveal strong responses to climate change. *Proc. Natl. Acad. Sci. USA* **117**, 18161–18168 (2020).

16. Grossiord, C. et al. Plant responses to rising vapor pressure deficit. *New Phytol.* **226**, 1550–1566 (2020).

17. Olano, J. M. et al. Drought-induced increase in water-use efficiency reduces secondary tree growth and tracheid wall thickness in a Mediterranean conifer. *Oecologia* **176**, 273–283 (2014).

18. Fu, Z. et al. Global critical soil moisture thresholds of plant water stress. *Nat. Commun.* **15**, 4826 (2024).

19. Yuan, W. et al. Increased atmospheric vapor pressure deficit reduces global vegetation growth. *Sci. Adv.* **5**, eaax1396 (2019).

20. Liu, L. et al. Soil moisture dominates dryness stress on ecosystem production globally. *Nat. Commun.* **11**, 4892 (2020).

21. Wang, X. et al. Establishing the global isoscape of leaf carbon in C3 plants through the integrations of remote sensing, carbon, geographic, and physiological information. *Remote Sens. Environ* **302**, 113987 (2024).

22. Ma, W. T. et al. Accounting for mesophyll conductance substantially improves ¹³C-based estimates of intrinsic water-use efficiency. *New Phytol.* **229**, 1326–1338 (2021).

23. Yu, Y. Z. et al. Reconciling water-use efficiency estimates from carbon isotope discrimination of leaf biomass and tree rings: non-photosynthetic fractionation matters. *New Phytol.* **244**, 2225–2238 (2024).

24. Wright, I. J., Reich, P. B. & Westoby, M. Least-cost input mixtures of water and nitrogen for photosynthesis. *Am. Nat.* **161**, 98–111 (2003).

25. Prentice, I. C., Dong, N., Gleason, S. M., Maire, V. & Wright, I. J. Balancing the costs of carbon gain and water transport: testing a new theoretical framework for plant functional ecology. *Ecol. Lett.* **17**, 82–91 (2014).

26. Wang, H. et al. Towards a universal model for carbon dioxide uptake by plants. *Nat. Plants* **3**, 734–741 (2017).

27. Körner, C. & Diemer, M. In situ Photosynthetic Responses to Light, Temperature and carbon dioxide in herbaceous plants from low and high altitude. *Funct. Ecol.* **1**, 179–194 (1987).

28. Wang, H. et al. Photosynthetic responses to altitude: an explanation based on optimality principles. *New Phytol.* **213**, 976–982 (2017).

29. Santos, C. M. D. et al. Photosynthetic capacity and water use efficiency in *Ricinus communis* (L.) under drought stress in semi-humid and semi-arid areas. *An. Acad. Bras. Cienc* **89**, 3015–3029 (2017).

30. Yang, Y. J. et al. Evolution of stomatal closure to optimize water-use efficiency in response to dehydration in ferns and seed plants. *New Phytol.* **230**, 2001–2010 (2021).

31. Adams, M. A., Buckley, T. N. & Turnbull, T. L. Rainfall drives variation in rates of change in intrinsic water use efficiency of tropical forests. *Nat. Commun.* **10**, 3661 (2019).

32. Belmecheri, S. et al. Precipitation alters the CO₂ effect on water-use efficiency of temperate forests. *Glob. Change Biol.* **27**, 1560–1571 (2021).

33. Diefendorf, A. F., Mueller, K. E., Wing, S. L., Koch, P. L. & Freeman, K. H. Global patterns in leaf ¹³C discrimination and implications for studies of past and future climate. *Proc. Natl. Acad. Sci. USA* **107**, 5738–5743 (2010).

34. Farquhar, G. D., O’Leary, M. H. & Berry, J. A. On the relationship between carbon isotope discrimination and the intercellular carbon dioxide concentration in leaves. *Funct. Plant Biol.* **9**, 121–137 (1982).

35. Stewart, G., Turnbull, M., Schmidt, S. & Erskine, P. ¹³C natural abundance in plant communities along a rainfall gradient: a biological integrator of water availability. *Funct. Plant Biol.* **22**, 51–55 (1995).

36. Gong, X. Y. et al. Overestimated gains in water-use efficiency by global forests. *Glob. Change Biol.* **28**, 4923–4934 (2022).

37. Lavergne, A. et al. Global decadal variability of plant carbon isotope discrimination and its link to gross primary production. *Glob. Change Biol.* **28**, 524–541 (2022).

38. Barbosa, I. C. R., KÖHler, I. H., Auerswald, K., LÜPs, P. & Schnyder, H. Last-century changes of alpine grassland water-use efficiency: a reconstruction through carbon isotope analysis of a time-series of *Capra ibex* horns. *Glob. Change Biol.* **16**, 1171–1180. (2010).

39. Battipaglia, G. et al. Elevated CO₂ increases tree-level intrinsic water use efficiency: insights from carbon and oxygen isotope analyses in tree rings across three forest FACE sites. *New Phytol.* **197**, 544–554 (2013).

40. Lu, X., Wang, L. & McCabe, M. F. Elevated CO₂ as a driver of global dryland greening. *Sci. Rep.* **6**, 20716 (2016).

41. Pastore, M. A., Lee, T. D., Hobbie, S. E. & Reich, P. B. Strong photosynthetic acclimation and enhanced water-use efficiency in grassland functional groups persist over 21 years of CO₂ enrichment, independent of nitrogen supply. *Glob. Change Biol.* **25**, 3031–3044 (2019).

42. Zhu, J. et al. Changing patterns of global nitrogen deposition driven by socio-economic development. *Nat. Commun.* **16**, 46 (2025).

43. Gharun, M. et al. Effect of nitrogen deposition on centennial forest water-use efficiency. *Environ. Res. Lett* **16**, 114036 (2021).

44. Brodribb, T. J., McAdam, S. A. M., Jordan, G. J. & Feild, T. S. Evolution of stomatal responsiveness to CO₂ and optimization of water-use efficiency among land plants. *New Phytol.* **183**, 839–847 (2009).

45. Saurer, M. et al. Spatial variability and temporal trends in water-use efficiency of European forests. *Glob. Change Biol.* **20**, 3700–3712 (2014).

46. Wang, M., Chen, Y., Wu, X. & Bai, Y. Forest-type-dependent water use efficiency trends across the Northern Hemisphere. *Geophys. Res. Lett.* **45**, 8283–8293 (2018).

47. Mathias, J. M., Thomas, R. B. Global tree intrinsic water use efficiency is enhanced by increased atmospheric CO₂ and modulated by climate and plant functional types. *Proc. Natl. Acad. Sci. USA* **118**, e2014286118 (2021).

48. Rumman, R., Atkin, O. K., Bloomfield, K. J. & Eamus, D. Variation in bulk-leaf ¹³C discrimination, leaf traits and water-use efficiency-trait relationships along a continental-scale climate gradient in Australia. *Glob. Change Biol.* **24**, 1186–1200 (2018).

49. Zhao, N., Meng, P., He, Y. & Yu, X. Interaction of CO₂ concentrations and water stress in semiarid plants causes diverging response in instantaneous water use efficiency and carbon isotope composition. *Biogeosciences* **14**, 3431–3444 (2017).

50. Linares, J. C. & Camarero, J. J. From pattern to process: linking intrinsic water-use efficiency to drought-induced forest decline. *Glob. Change Biol.* **18**, 1000–1015 (2012).

51. Correa-Díaz, A., Gómez-Guerrero, A., Castruita-Esparza, L. U., Silva, L. C. R. & Horwath, W. R. Divergent responses of fir and pine trees to increasing CO₂ levels in the face of climate change. *J. Geophys. Res. Biogeosci.* **129**, e2023JG007754 (2024).

52. Lu, H. et al. Large influence of atmospheric vapor pressure deficit on ecosystem production efficiency. *Nat. Commun.* **13**, 1653 (2022).

53. Novick, K. A. et al. The increasing importance of atmospheric demand for ecosystem water and carbon fluxes. *Nat. Clim. Chang.* **6**, 1023–1027 (2016).

54. Sulman, B. N. et al. High atmospheric demand for water can limit forest carbon uptake and transpiration as severely as dry soil. *Geophys. Res. Lett.* **43**, 9686–9695 (2016).

55. Tyree, M. T. & Sperry, J. S. Vulnerability of xylem to cavitation and embolism. *Annu. Rev. Plant Biol.* **40**, 19–36 (1989).

56. Huang, M. et al. Change in terrestrial ecosystem water-use efficiency over the last three decades. *Glob. Change Biol.* **21**, 2366–2378 (2015).
57. Keenan, T. F. et al. A constraint on historic growth in global photosynthesis due to rising CO₂. *Nat. Clim. Chang.* **13**, 1376–1381 (2023).

58. Stocker, B. D. et al. P-model v1.0: an optimality-based light use efficiency model for simulating ecosystem gross primary production. *Geosci. Model Dev.* **13**, 1545–1581 (2020).

59. Zhu, Z. et al. Optimality principles explaining divergent responses of alpine vegetation to environmental change. *Glob. Change Biol.* **29**, 126–142 (2023).

60. Lavergne, A., Sandoval, D., Hare, V. J., Graven, H. & Prentice, I. C. Impacts of soil water stress on the acclimated stomatal limitation of photosynthesis: Insights from stable carbon isotope data. *Glob. Change Biol.* **26**, 7158–7172 (2020).

61. Xu, H. et al. Predictability of leaf traits with climate and elevation: a case study in Gongga Mountain, China. *Tree Physiol.* **41**, 1336–1352 (2021).

62. Li, X. et al. Differences in responses of tree-ring δ¹³C in angiosperms and gymnosperms to climate change on a global scale. *For. Ecol. Manage.* **492**, 119247 (2021).

63. Li X., Xiao J. A global, 0.05-degree product of solar-induced chlorophyll fluorescence derived from OCO-2, MODIS, and reanalysis data. *Remote Sens.* **11**, 517 (2019).

64. Zhang, Y., Joiner, J., Alemohammad, S. H., Zhou, S. & Gentine, P. A global spatially contiguous solar-induced fluorescence (CSIF) dataset using neural networks. *Biogeosciences* **15**, 5779–5800 (2018).

65. Fick, S. E. & Hijmans, R. J. WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. *Int. J. Climatol.* **37**, 4302–4315 (2017).

66. Nix, H. A. A Biogeographic Analysis of Australian Elapid Snakes, *7*, (Australian Government Publishing Service: Canberra, ACT, 1986).

67. Abatzoglou, J. T., Dobrowski, S. Z., Parks, S. A. & Hegewisch, K. C. TerraClimate, a high-resolution global dataset of monthly climate and climatic water balance from 1958-2015. *Sci. Data* **5**, 170191 (2018).

68. Buck, A. L. New Equations for Computing Vapor Pressure and Enhancement Factor. *J. Appl. Meteorol. Climatol.* **20**, 1527–1532 (1981).

69. Luo, X. et al. Mapping the global distribution of C4 vegetation using observations and optimality theory. *Nat. Commun.* **15**, 1219 (2024).

70. Gregorutti, B., Michel, B. & Saint-Pierre, P. Correlation and variable importance in random forests. *Stat. Comput.* **27**, 659–678 (2016).

71. Naimi, B., Hamm, N. A. S., Groen, T. A., Skidmore, A. K. & Toxopeus, A. G. Where is positional uncertainty a problem for species distribution modelling? *Ecography* **37**, 191–203 (2014).

72. Kuhn, M. Building predictive models in R Using the caret Package. *J. Stat. Softw.* **28**, 1–26 (2008).

73. O'Brien, R. M. A caution regarding rules of thumb for variance inflation factors. *Qual. Quant.* **41**, 673–690 (2007).

74. Stine, R. A. Graphical interpretation of variance inflation factors. *Am. Stat.* **49**, 53–56 (1995).

75. Wang, X., Ciais, P., Wang, Y. & Zhu, D. Divergent response of seasonally dry tropical vegetation to climatic variations in dry and wet seasons. *Glob. Chang. Biol.* **24**, 4709–4717 (2018).

76. Hengl, T. et al. SoilGrids250m: Global gridded soil information based on machine learning. *PLoS ONE* **12**, e0169748 (2017).

77. Brooks, A. & Farquhar, G. D. Effect of temperature on the CO₂/O₂ specificity of ribulose-1,5-bisphosphate carboxylase/oxygenase and the rate of respiration in the light. *Planta* **165**, 397–406 (1985).

78. Wang, S. et al. Recent global decline of CO₂ fertilization effects on vegetation photosynthesis. *Science* **370**, 1295–1300 (2020).

79. Sen, K. umarP. Estimates of the regression coefficient based on Kendall's Tau. *J. Am. Stat. Assoc.* **63**, 1379–1389 (1968).

80. Wang, X. et al. Spatio-temporal patterns and drivers of carbon-water coupling in frozen soil zones across the gradients of freezing over the Qinghai-Tibet Plateau. *J. Hydrol.* **621**, 129674 (2023).

81. Piao, S., Wang, J., Li, X., Xu, H. & Zhang, Y. Spatio-temporal changes in the speed of canopy development and senescence in temperate China. *Glob Change Biol.* **28**, 7366–7375 (2022).

82. Jiao, W. et al. Observed increasing water constraint on vegetation growth over the last three decades. *Nat. Commun.* **12**, 3777 (2021).

83. Petek-Petrik, A. et al. The combined effect of branch position, temperature, and VPD on gas exchange and water-use efficiency of Norway spruce. *Biol. Plant* **67**, 136–141 (2023).

84. Xue, L. et al. Analysis on water use efficiency of Populus euphratica forest ecosystem in arid area. *Theor. Appl. Climatol.* **145**, 717–730 (2021).

## Acknowledgements

This work was supported by the National Key R&D Program of China (No. 2024YFF1309000, Z.F.), the National Natural Science Foundation of China (Nos. 42471122 and 32588202, Z.F.; No. 32371683, G.C.), the Chinese Academy of Sciences (No. 2024000275, Z.F.), and the NSFC Excellent Young Scientists Fund (Overseas). We acknowledges Xiaowei Li, Edmund February, Gustavo Saiz, and Ralf Kiese for the leaf isotope data. T.F.K acknowledge support by the LEMONTREE (Land Ecosystem Models based On New Theory, Observation and Experiments) project, funded through the generosity of Eric and Wendy Schmidt by recommendation of the Schmidt Futures program. T.F.K. acknowledges additional support from a NASA Carbon Cycle Science Award 80NSSC21K1705. J.P. acknowledges support by the grant TED2021-132627 B-100 funded by the Spanish MCIN, AEI/10.13039/501100011033 European Union Next Generation EU/PRTR, and the European project CONCERTO (HORIZON-CL5-2024-D1-01). N.B. acknowledges support by the SNF project CERES (IZCOZO_220291).

## Author contributions

X.W. and Z.F. designed the study. X.W. performed the analysis. X.W. and Z.F. wrote the paper with the inputs from all co-authors. P.C., L.X.W., N.B., T.F.K., M.D.K., J.P., G.C., X.Y.G., J.F.X., X.L., Q.Y.X., P.C.S., D.M., W.K.S., H.W., S.H.W., F.Y.Z. and S.L.N. provided methodological suggestions and contributed to the interpretation of the results.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-025-68252-9.

**Correspondence** and requests for materials should be addressed to Zheng Fu.

**Peer review information** *Nature Communications* thanks Peter Petrik and the other anonymous reviewer(s) for their contribution to the peer review of this work. A peer review file is available.

**Reprints and permissions information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
**Open Access** This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

© The Author(s) 2026

<footnote n="1">¹Key Laboratory of Ecosystem Network Observation and Modeling, Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences, Beijing, China. ²Laboratoire des Sciences du Climat et de l'Environnement, LSCE/IPSL, CEA-CNRS-UVSQ, Université Paris-Saclay, Gif-sur-Yvette, France. ³Department of Earth and Environmental Sciences, Indiana University Indianapolis, Indianapolis, IN, USA. ⁴Institute of Agricultural Sciences, ETH Zurich, Zurich, Switzerland. ⁵Climate and Ecosystem Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA, USA. ⁶Department of Environmental Science, Policy and Management, UC Berkeley, Berkeley, CA, USA. ⁷School of Biological Sciences, University of Bristol, Bristol, UK. ⁸CSIC, Global Ecology Unit, CREAF-CSIC-UAB, Barcelona, Catalonia, Spain. ⁹CREAF, Cerdanyola del Vallès, Barcelona, Catalonia, Spain. ¹⁰College of Ecology and Environment, Chengdu University of Technology, Chengdu, Sichuan, China. ¹¹Key Laboratory for Subtropical Mountain Ecology, College of Geographical Sciences, Fujian Normal University, Fuzhou, China. ¹²Earth Systems Research Center, Institute for the Study of Earth, Oceans, and Space, University of New Hampshire, Durham, NH, USA. ¹³School of Geography and Planning, Sun Yat-sen University, Guangzhou, China. ¹⁴School of Engineering, The University of Western Australia, Perth, WA, Australia. ¹⁵Department of Biological Systems Engineering, University of Wisconsin-Madison, Madison, WI, USA. ¹⁶Unit Applied Mathematics and Computer Science (UMR MIA-PS) INRAE AgroParisTech Université Paris-Saclay, Palaiseau, France. ¹⁷School of Natural Resources and the Environment, University of Arizona, Tucson, AZ, USA. ¹⁸Ministry of Education Key Laboratory for Earth System Modeling, Department of Earth System Science, Tsinghua University, Beijing, China. ¹⁹College of Agriculture, Nanjing Agricultural University, Nanjing, China. ²⁰College of Water Sciences, Beijing Normal University, Beijing, China. e-mail: fuzheng@ignnrr.ac.cn</footnote>