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

Domain: biology
Published: 2026-01-07
Editorial outcome: major_revision

Abstract
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 C3 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-1 year-1, 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.

Introduction
Climate warming and rising atmospheric CO2 concentrations (Ca) change the coupling of carbon and water in terrestrial ecosystems1,2. Intrinsic water use efficiency (iWUE), the ratio of leaf net photosynthetic rate to stomatal conductance (i.e., A/gs), is a key indicator for assessing plant responses to climate change3–5. 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 C3 plant photosynthesis, leaf 13C discrimination (Δ13C) is linearly related to leaf iWUE6–9. 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 models2. 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 ecosystems10. For example, using tree-ring isotope data to infer iWUE gives a long time series but does not cover non-forested areas11.

Leaf-level iWUE is influenced mainly by atmospheric CO2, vapor pressure deficit (VPD), and soil moisture (SM). Studies using free-air CO2 enrichment experiment12 and gas-exchange methods13 have shown that rising Ca reduces leaf stomatal conductance (gs), thereby promoting an increase in iWUE. Recently, an analysis of global tree-ring carbon isotope data suggested that while elevated Ca has driven a substantial increase in forest iWUE over the 20th century, the rate of change in iWUE relative to Ca has decreased by half14. Additionally, both theory and experimental studies suggest that increased atmospheric and soil drought may enhance iWUE11,15,16. Consistent with theory, increased VPD and reduced SM have been found to reduce leaf gs6,17, resulting in higher Ci/Ca and higher iWUE, with these effects becoming more pronounced over time18–20. Despite advancements in understanding how rising Ca 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 C3 plant foliage samples21, spanning the time period 2001–2020 while covering all continents except Antarctica, and construct global maps of carbon isotope discrimination (Δ13C) in C3 plants during 2001–2020 using a random forest explainable machine learning models ensemble (Methods). The resulting annual maps of Δ13C are then translated into maps of iWUE for C3 plants using Ma’s22 and Yu’s23 equations to calculate leaf-level iWUE, which incorporate the effect of mesophyll resistance gm 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 theory24,25 encapsulated in the P-model26 simulates global leaf-level iWUE distribution and trends, and discuss differences between optimality theory and observation-based estimates.

Results
Spatiotemporal dynamics of global leaf-level iWUE
We quantified the spatial distribution of iWUE in C3 plants using Δ13C of natural C3 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 ecosystems11)—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−1 (mean±standard deviation) (Fig. 1a, b). Spatially, iWUE values were lower in humid regions (mean: 52.5 μmol mol−1) and higher in arid and cold regions (mean: 60.4 μmol mol−1, 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 gs. This leads to declines in both photosynthetic and transpiration rates11,27,28. 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 photosynthesis29,30. Conversely, in warm and humid regions, higher water availability in soils increases air humidity and SM, leading to increased leaf gs31–33. With elevated gs, leaf Ci and transpiration rates rise, resulting in a decreased leaf iWUE30,34,35, 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 functionality 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 C3 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 C3 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−1 year−1 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−1 year−1)36. However, this trend was significantly lower than that estimated by the JULES model37 (0.28±0.01 μmol mol−1 year−1, 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 Ca, i.e., the “CO2 fertilization effect”2,4,38–40. This is supported by plant physiological studies indicating that rising Ca promotes CO2 uptake and enhances photosynthesis in vegetation41, while leading to a reduction in gs, which in turn reduce transpiration39. This was illustrated by the 12% increase in Ca between 2001 and 2020 (i.e., from 370 ppm to 414 ppm), corresponding to the rise in iWUE of global C3 plants of about 8.5±0.11%. However, atmospheric nitrogen (N) deposition, showing large spatiotemporal trends globally42, was also shown to increase WUE of temperate forests (based on stable isotope data in tree-rings)43, interacting with increasing Ca and global warming.

The increasing trend of leaf-level iWUE was observed in all C3 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−1 year−1, DBF iWUE: 0.29±0.04 μmol mol−1 year−1), which was 40% larger than that of conifers (ENF iWUE: 0.18±0.03 μmol mol−1 year−1, DNF iWUE: 0.24±0.02 μmol mol−1 year−1). In contrast, C3 grasses had a relatively slower iWUE growth, with grasslands increasing by only 0.15±0.02 μmol mol−1 year-1 (Fig. 1g). Previous indoor experimental studies have shown that angiosperm stomata are more sensitive to elevated Ca compared to conifers and other PFTs. As a result, angiosperm gs tends to be more sensitive to rising Ca, leading to a more rapid increase of iWUE44. Our findings aligned with previous studies by Saurer et al.45 and Wang et al.46, 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.4 and Mathias et al.47, 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 Ca 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 experiments17,48,49. Previous studies suggested that water stress-induced increases of iWUE are primarily driven by changes in gs, as plants prioritize water conservation over carbon uptake. Stomatal closure reduces transpiration more than it reduces photosynthesis11, 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 SM and lower VPD (humid regions) compared to arid regions (Fig. 2b). Interestingly, the sensitivity of iWUE to Ca (diWUE/dCa) followed the same pattern, with diWUE/dCa 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/dCa were analyzed using the aridity index (Fig. 2d–f). In particular, both iWUE trends and diWUE/dCa 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 globe11. In fact, a number of factors influence iWUE trends across climate zones, especially their interactions with water stress47. Several tree-ring studies showed that the positive effect of rising Ca on iWUE was attenuated by increasing water stress14,47,50,51. Therefore, plant iWUE in the humid zone benefited more from the increase in Ca than plants in the arid zone, resulting in a faster rate of increase in iWUE in this region. This was also verified by our analysis, with the sensitivity of iWUE to Ca 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 C3 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 C3 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 SM52–54 because VPD regulates the balance between water supply and demand by controlling plant leaf water potential, thereby directly influencing water demand55. 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 closure15,16. 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)11,15,56. Overall, our findings suggest that changes in iWUE for global C3 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 (98% of the global C3 covered area) (Supplementary Figs. S3c, d, S4). The globally increasing rate of iWUE estimated by the P-model was 0.49±0.02 μmol mol−1 year−1, 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 to Ca changes (Supplementary Fig. S5). Despite this, the model’s response of iWUE to variations in VPD and SM gradients generally aligned with observed patterns, showing that water stress led to an increase in iWUE (Supplementary Fig. S6). Moreover, model-based analysis also indicated that VPD played a more important role than SM in influencing iWUE 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 Ci/Ca ratio in plant leaves24–26. As such, the optimal Ci/Ca 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 photosynthesis57–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 Ca (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 ξ=√β(K+Γ*)/1.6(1+gs/gm), 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 Ci/Ca)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 trees60. However, a single β value is currently used for all PFTs in the P-model developed based on the optimality theory26,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 Ci/Ca uncertainty61. It was observed that the use of β estimated from compiled data to predict Ci/Ca resulted in a mismatch with observations. Predicted Ci/Ca values were generally lower than the median observed values at most sites61. 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 changes in more than 75% of the global C3 vegetation cover, which may be caused by the fact that the P-model only considered the effect of atmospheric water demand on Ci/Ca and ignored the effect of soil water stress. Previous studies has shown that as soil dryness, indicating that the cost of plant water transport increases as SM declines60. This bias highlights the need to carefully evaluate Earth system models that adopts similar assumptions and parameterizations, as they may misrepresent future 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 CO2 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 C3 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−1 year−1. Grasslands showed the slowest iWUE growth rate (0.15±0.02 μmol mol–1 year−1), whereas evergreen broadleaf forests exhibited the fastest growth (0.3±0.04 μmol mol−1 year−1). 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 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 C3 plant foliage worldwide21 (Supplementary Figs. S7, S8a). The Δ13C was then calculated from leaf carbon isotope ratios (δ13CLeaf, ‰) and atmospheric carbon isotope ratios (δ13Ca, ‰, δ13Ca corresponds to the leaf sample sampling year) (Eq. 1):

Δ13C = (δ13Ca - δ13CLeaf) / (1 + δ13CLeaf / 1000) (1)

The dataset collected and compiled 5964 datasets of global carbon isotope composition (δ13C) of C3 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: “C3 plant”, “leaf”, “carbon isotope” and “leaf δ13C”. 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 C3 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 C3 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 δ13C 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 processes21. Random tests in earlier studies showed no significant differences between the extracted values and the true values, indicating that the software had good accuracy62. 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 Δ13C isoscape map for C3 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 timescale of the constructed isoscapes mapping of Δ13C of C3 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.

CO2 concentrations and δ13Ca
Due to the lack of a reliable and high-resolution global map of atmospheric CO2 concentrations (Ca, ppm) and carbon isotope ratios (δ13Ca, ‰), we used monthly Ca and δ13Ca 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/m2/μm/sr) and global spatially contiguous solar-induced chlorophyll fluorescence dataset (CSIF, W/m2/μm/sr), since SIF is closely related to plant photosynthesis. GOSIF was developed and made publicly available by Li and Xiao63 in 2019 with a spatial resolution of 0.05°. CSIF was developed and made publicly available by Zhang et al.64 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 GloElev_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 variables65. 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 Nix66 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, Wm−2) used in this study were obtained from the TerraClimate dataset developed by Abatzoglou et al.67, 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 Buck68. 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 C4 vegetation percentage maps with a spatial resolution of 0.5° were resampled to 0.05°, and grids with >50% C4 vegetation were labeled as C4 dominated areas and masked. Global C4 distribution data were developed by Luo et al. (2024)69 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 evapotranspiration 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 Δ13C isoscapes of C3 plants
Currently, only the leaf carbon isotope discrimination (Δ13C) of plants with the C3 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 Δ13C 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 Δ13C isoscape mapping, we used the collected satellite spectral reflectance (seven spectral bands), atmospheric CO2 concentration (Ca), 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 CO2 concentrations and carbon isotopes, annual means of SIF, and elevation information. The δ13Ca of atmospheric CO2 are used to calculate the Δ13C (Eq. 1), and the atmospheric CO2 concentration was used to predict the Δ13C.

Before constructing isoscapes of leaf carbon isotopes in C3 plants, we needed to filter the selected predictor variables. Because not all predictor variables are relevant for the construction of leaf Δ13C isoscapes, variable importance need 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 accuracy70. We used the vifstep function (multicollinearity test, Eq. 2) of the “usdm” R package71 and the recursive feature elimination (RFE, Eq. 3) of the “Caret” package72 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 1073–75, 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, Ca, Lat, Long, ELE) (2)

RFE = rfeControl(B1, B2, B3, B4, B5, B6, B7, SIF, Ca, Lat, Long, ELE) (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; Ca is the atmospheric CO2 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, Ca, 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 variables76. Here, we used the Random Forest algorithm to construct a global Δ13C isoscapes of C3 plant leaves based on Wang et al.21 (Eq. 4).

Δ13C = RandomForest(B3, B5, B7, SIF, Ca, Lat, Long, ELE) (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 R2 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 C3 plant leaves for Δ13C isoscapes. Furthermore, to remove the interference of the C4 photosynthetic pathway, we defined grid cells with pixel values greater than 50% on the C4 vegetation cover map as dominated by C4 plants and then masked these areas from the predicted Δ13C isoscapes of C3 plant leaves globally21.

Estimation of iWUE based on Δ13C of C3 plant leaves
Intrinsic water use efficiency (iWUE) was defined as the ratio of net photosynthetic rate (A) and stomatal conductance (gs), and the calculation process depended on the atmospheric CO2 concentration (Ca) and the leaf intercellular CO2 concentration (Ci) (Eq. 5):

iWUE = A / gs = Ca (1 - Ci / Ca) / 1.6 (5)

Ci/Ca can be estimated from leaf Δ13C (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.34. 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.22 and Yu et al.23 proposed an improved version of the formula (Eq. 7), by taking into account the parameters of mesophyll conductance and nonphotosynthetic fractionations:

Ci / Ca = (Δ13CLeaf - a) / (b - a) (6)

iWUE = Ca * (b - (Δ13CLeaf - dLeaf)) / (1.6 * (b - a) + gs * (b - a) / gm) (7)

where 1.6 is the diffusion ratio of water vapor to CO2 as it passes through leaf stomata; a and b are the fractionation coefficients for CO2 as it diffuses through leaf stomata (4.4‰) and Rubisco carboxylation process (29 ‰), respectively; f denotes the fractionations due to photorespiratory fractionation (11 ‰); Γ* is the carbon dioxide compensation point in the absence of mitochondrial respiration (Eq. 8)77; and gs/gm is the ratio of the stomatal conductance to the mesophyll conductance, Ma et al.22 found gs/gm 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); dLeaf (2.5‰) is the isotope correction of bulk leaf biomass to account for nonphotosynthetic effects23; am is the fractionation of CO2 during solubilization and diffusion in chloroplasts (1.8‰).

Γ* = 42.7 + 1.68(T - 25) + 0.012(T - 25)2 (8)

where T (°C) is leaf temperature, which were presented 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 Δ13C for C3 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 CO2 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 capacity24–26. Based on this theory, the iWUE equation parameter Ci/Ca (considering the effect of gm, substituting Ci/Ca for Ci/Ca) can be estimated by the P-model developed by optimality theory26, i.e., Eq. 9 and Eq. 10:

χ = Ci / Ca = Γ* / Ca + (1 - Γ* / Ca) * ξ / (ξ + VPD) (9)

ξ = sqrt(β * (K + Γ*) / 1.6) (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 Ci/Ca 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 gs/gm) based on Wang et al.26; 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.26 used unpublished data from I.J. Wright to obtain a median gs/gm of 1.4. The χ estimated based on the P-model (Supplementary Fig. S15) was consistent with the results of Wang et al.26. The calculated χ (Ci/Ca) was finally substituted into the iWUE calculation formula, i.e., Eq. 11 or Eq. 12:

iWUEmod = Ca * (1 - χ) / 1.6 (11)

iWUEmod = (Ca - Γ*) / (1.6 * (1 + sqrt(ξ) / VPD)) (12)

where the Ca unit entered in iWUE is μmol/mol, so the iWUE unit obtained is μmol/mol.

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

iWUE = (diWUE / dCa) * Ca + β1 * VPD + β2 * MAT + β3 * SM + β4 * radiation + ε (13)

where diWUE/dCa is the sensitivity of iWUE to Ca, β1, β2, β3 and β4 are the regression coefficients, and ε is the residual error term; Ca, 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 202079. 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 non-parametric trend statistics and tests that have been widely used in the analysis of spatial time series data80. 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 variable81. 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 cell82. In addition, given that the climatic effect of VPD inherently includes a temperature component, and that plants assimilation exhibit a nonlinear optimum response to temperature83,84, 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.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The study presents a large-scale analysis of leaf-level intrinsic water-use efficiency (iWUE) using a combination of machine learning, carbon isotope data, and ecological models. It provides important insights into the spatial and temporal dynamics of iWUE, as well as its response to water stress at the global scale.\nThe topic is highly relevant to global change biology, plant ecophysiology, and climate modelling. However, several major issues need to be addressed before the manuscript can be considered for publication in Nature Communications. I am suggesting major revision of the manuscript.\nComments:\n1. My main critique is that the VPD climatic effect incorporates the temperature effect. The temperature response of iWUE follows the optimum temperature curve of assimilation (e.g. 10.1007/s00704-021-03636-7; 10.32615/bp.2023.017), therefore the low iWUE of Boreal region could be due to low temperature rather than low VPD. Could you separate the temperature effect from VPD?\n2. The study claims to establish a \"novel framework for large-scale carbon-water coupling\", but does not sufficiently explain how it differs from previous approaches (e.g., tree-ring-based studies or flux-tower estimates). Strengthen the discussion by explicitly comparing the advantages and limitations of this approach over previous studies (e.g., why does this outperform tree-ring iWUE reconstructions?).\n3. The P-model overestimates iWUE trends compared to observations, but no robust explanation is provided for this discrepancy. Provide a more detailed analysis of the model-observation mismatch rather than just noting that the P-model overestimates trends.\n4. The machine learning random forest model lacks an independent validation dataset. While cross-validation is useful, an independent test dataset would provide stronger validation.\n5. The sampling distribution is not evenly spread globally—some regions (e.g., boreal forests, Africa) are underrepresented. Discuss the implications of uneven sampling distribution and how this may bias the results.\n6. The manuscript has frequent grammatical errors (e.g., \"water sress\" instead of \"water stress\"; \"trand\" instead of \"trend\") or weird sentence composition (e.g., Line 51 “Atmospheric dryness (vapor pressure deficit, VPD) increases the iWUE across a larger area than soil dryness“). Some redundant statements repeat findings across different sections. Improve language clarity and eliminate typographical errors.\n7. Line 61. I would say that net assimilation rate or net photosynthetic rate is more precise than photosynthetic rate.\n8. Line 65. The leaf 13C discrimination is also linearly related to leaf iWUE, maybe easier for reader to comprehend (https://doi.org/10.1111/ppl.14619). Moreover, I don’t understand how cited reference by Adams et al. links carbon discrimination to iWUE, it should be replaced with something more relevant like https://doi.org/10.2135/cropsci2005.0119.\n9. Line 79: Only theory? There are multiple experimental studies that showed that.\n10. Line 280: „We found that the global distribution of leaf iWUE aligned with physiological responses of plants: high iWUE values in cold and arid regions, and lower values in warm and humid areas.“ I firstly don’t understand what the authors mean by physiological responses here in connection to climate zones. Secondly, how does this align with the found positive trends? You also mentioned beforehand that the boreal forests had the lowest iWUE, now cold regions have high iWUE?", "reviewer_2": "This manuscript provides a valuable global assessment of leaf-level intrinsic water use efficiency (iWUE) across C3 vegetation, leveraging a large leaf δ¹³C dataset, machine-learning upscaling, and established physiological equations. It offers a novel product at 0.05° resolution from 2001 to 2020, revealing strong spatiotemporal patterns of iWUE and clear evidence that VPD, rather than soil moisture, dominates iWUE changes in most regions. These results will be beneficial to both carbon–water cycle modeling and ecological research.\nMajor Strengths\n1. The new global iWUE dataset covers diverse ecosystems with a coherent methodology.\n2. Incorporation of mesophyll conductance and nonphotosynthetic fractionations strengthens physiological realism.\n3. Comparisons with the P-model highlight uncertainties, underscoring the need to refine certain model parameters.\nAreas for Improvement\n1. Data Coverage and Uncertainty: Provide clearer uncertainty estimates (e.g., standard-error maps), especially in under-sampled regions (e.g., parts of Africa, high latitudes).\n2. Atmospheric δ¹³C Assumption: Justify using Mauna Loa δ¹³C globally or discuss potential biases in regions far from Hawaii or with unique CO₂ sources/sinks.\n3. Additional Confounding Factors: Briefly note how nitrogen deposition, forest age, or ontogeny might affect iWUE trends in certain regions.\n4. Methodological Clarification: Include a short sensitivity analysis (e.g., varying the gₛ/gₘ ratio) and indicate how partial correlations (VPD vs. SM) cope with variable collinearity.\nOverall, the methods are solid, and the findings have broad significance. Addressing these points in a revision will further strengthen the manuscript and clarify its conclusions."}, "author_response": "Response to Reviewers\nWe greatly appreciate the constructive comments from the reviewers invitation to submit a\nrevised version. We have thoroughly revised the manuscript following the reviewers’ suggestions\n. Please see below our point-to-point responses in blue text followin g the reviewers’comments.\nThe line nu mbers referred to are for the clean version of th e revised manuscript (non-track-change version).\n---------------------------------------------------------------------------------------------------------------------\nReviewer #1 (Remarks to the Author):\n1.1 The study presents a large-scale analysis of leaf-level intrinsic water-use efficiency (iWUE)\nusing a combination of machine learning, carbon isotope data, and ecological models. It provides\nimportant insights into the spatial and temporal dynamics of iWUE, as well as its response to water\nstress at the global scale.\nThe topic is highly relevant to global change biology, plant ecophysiology, and climate modelling.\nHowever, several major issues need to be addressed before the manuscript can be considered for\npublication in Nature Communications. I am suggesting major revision of the manuscript.\nWe greatly appreciate the reviewer’s positive comments and support of the manuscript. We have\nthoroughly revised the manuscript following your suggestions. As a consequence, our manuscript\nhas been considerably improved. Please see below our detailed responses.\nOther specific comments:\n1.2 My main critique is that the VPD climatic effect incorporates the temperature effect. The\ntemperature response of iWUE follows the optimum temperature curve of assimilation (e.g.,\n10.1007/s00704-021-03636-7; 10.32615/bp.2023.017), therefore the low iWUE of Boreal region\ncould be due to low temperature rather than low VPD. Could you separate the temperature effect\nfrom VPD?\nThank you very much for this valuable suggestion. We fully agree that the climatic effect of VPD\ninherently includes a temperature component, and that iWUE exhibits a nonlinear response to\ntemperature. We have clarified this point in the revised manuscript. Specifically, we note that\nestimated iWUE values are higher in boreal regions. We have described the spatial distribution of\nmean iWUE in Lines 104–115 and illustrated this in Fig. 1a. We also discussed the potential\nmechanisms by which high iWUE values in cold regions (e.g., boreal ecosystems) may result from\nlow-temperature stress. The revised text reads as follows:\nLines 104–115:\n“From 2001-2020, we estimated a global mean iWUE of 55.4 ± 1.6 μmol mol-1 (mean ± standard\ndeviation) (Fig. 1a, b). Spatially, iWUE values were lower in humid regions (mean: 52.5 μmol\nmol-1) and higher in arid and cold regions (mean: 60.4 μmol mol-1, Fig. 1a, b). This finding aligns\nwith ecophysiological trade-off mechanisms through which plant foliage in cold and arid regions\nare subject to cold or water stress and operate at reduced g . This leads to declines in both\ns\nphotosynthetic and transpiration rates1, 2, 3. The high iWUE observed in high-latitude boreal regions\n(cold) and across the Australian continent (dry) (Fig. 1a, b) may be due to the fact that the decline\nwas stronger for transpiration than for photosynthesis 4, 5. Conversely, in warm and humid regions,\nhigher water availability in soils increases air humidity and SM, leading to increased leaf g 6, 7, 8.\ns\nWith elevated g , leaf C and transpiration rates rise, resulting in a decreased leaf iWUE5, 9, 10, as\ns i\nobserved in tropical rainforests regions such as Southeast Asia, central Africa, and the Amazon\n(Fig. 1a, b).”\nFig. R1. Global distribution of leaf-level iWUE for C plants during 2001-2020. a, Spatial\n3\ndistribution of the mean iWUE. b, Latitudinal patterns of mean iWUE. Non-natural vegetation\nareas and deserts with less than 10% vegetation coverage were masked, as were grid units\ndominated by C vegetation.\nand soil moisture (SM) was not to explain the absolute magnitude of iWUE, but rather to evaluate\nthe relative influence and spatial dominance of atmospheric dryness (VPD) versus soil dryness\n(SM) on global iWUE. Specifically, we used 20 years of iWUE, VPD, and SM data to assess which\nvariable had stronger effect on iWUE at each grid cell. In this context, this analysis is not intended\nto precisely determine the iWUE values in given regions.\nIn addition, following your suggestion, we conducted an additional analysis using the residuals of\nfits to VPD and SM after removing the influence of air temperature, to perform a refined partial\ncorrelation analysis with iWUE. This approach tentatively isolates the confounding effect of\ntemperature and enables a clearer evaluation of the independent contributions of atmospheric and\nsoil drought to iWUE. Specifically, we applied linear regression to remove a temperature-driven\ncomponents from VPD and SM, treating the residuals as temperature-independent indicators of\nwater stress. The results from this residual-based analysis confirms that VPD is the dominant driver\nof iWUE variability in most regions, accounting for 60.4% of C vegetation-covered areas, while\nfindings reinforce our conclusion that VPD dominated iWUE changes across a larger area\ncompared to SM at the global scale.\nWe have added this supplementary analysis and Supplementary Fig. 2 into the revised manuscript\n(Lines 546–551, Lines 217–220).\nLines 546–551:\n“In addition, given that the climatic effect of VPD inherently includes a temperature component,\nwe further computed the correlations between iWUE and residuals of VPD and SM (obtained after\nremoving a linear influence of air temperature). Specifically, we first applied linear regression to\nremove temperature-driven components from VPD and SM, and then treated the residuals as\ntemperature-independent indicators of water stress, although we acknowledge the fact that non-\nlinear effects of temperature may still be included in the residuals.”\nLines 217–220:\n“Partial correlation analysis based on residual VPD and SM, with temperature linear effects on\nVPD and SM removed, also confirmed this result, identifying VPD and SM as the dominant drivers\nin 60.4% and 39.6% of C vegetation areas, respectively (Fig. S2).”\n3\nfrom temperature influence on iWUE. The spatial distribution of partial correlation coefficients of\nVPD and SM controlling changes in iWUE (a), the proportion of the area controlled by VPD or\nSM (b). Non-natural vegetation areas and deserts with less than 10% vegetation coverage were\nmasked, as were grid units dominated by C vegetation.\n4\n1.3 The study claims to establish a \"novel framework for large-scale carbon-water coupling\", but\ndoes not sufficiently explain how it differs from previous approaches (e.g., tree-ring-based studies\nor flux-tower estimates). Strengthen the discussion by explicitly comparing the advantages and\nlimitations of this approach over previous studies (e.g., why does this outperform tree-ring iWUE\nreconstructions?).\nThank you for this insightful comment. In the revised manuscript, we have strengthened the\ndiscussion by explicitly comparing the advantages and limitations of our approach relative to\nprevious studies, as suggested (Lines 99–104, 655–670).\nLines 99–104:\n“Compared to previous studies based on in situ leaf or tree-ring isotope measurements—which are\ngeographically sparse or biased toward forest ecosystems (unrepresentative of non-forest\necosystems1)—our study provides a spatially continuous, observation-constrained framework for\nestimating leaf-level iWUE across biomes. This approach addresses key limitations in spatial\ncoverage and scaling, enabling cross-biome estimation of iWUE, and offering broader\napplicability for large-scale analyses of carbon–water coupling.”\nLines 655–670:\n4\n“However, the uneven spatial distribution of C leaf isotope samples used for training our\n3\nupscaling models may affect the stability of our upscaled estimates. Therefore, the uncertainty in\nthe model construction process was assessed by performing iterative model training and different\nvalidation methods (including training–testing splits, ten-fold cross-validation, and random\nresampling validation), as well as spatial pattern comparisons (Figs. S10, S12). To better present\nthe spatial uncertainty, we also plotted the coefficient of variation (CV = (standard\ndeviation/average) × 100) of the output values predicted by the algorithms. The pattern of global\nCVs was uniform and low overall, with higher CVs in regions such as the southern edge of the\nAfrican Sahel, the Central Asian highlands, and the Tibetan Plateau, which may be due to the\nsparse distribution of sample points in these regions (Fig. S12). Importantly, this uncertainty was\nnot expected to affect the relative changes in the temporal scales of the predicted results, as\nevidenced by the range of the predicted results and the robustness of the size patterns and their\nannual variations across PFTs, as well as the consistency with the independent validation dataset\n(Figs. S9-S14). Overall, despite these limitations, the isoscape of C plant Δ13C predicted by our\n3\nmethod matches well with observed leaf isotope data and with the results from both a linear\nstatistical model11 and an Earth system model12, suggesting that our isoscapes better capture the\nglobal spatial variability of C plant Δ13C.”\n3\n1.4 The P-model overestimates iWUE trends compared to observations, but no robust explanation\nis provided for this discrepancy. Provide a more detailed analysis of the model-observation\nmismatch rather than just noting that the P-model overestimates trends.\nThank you for this helpful suggestion. While the core goal of our study is to map the global\ndistribution and temporal changes of leaf-level intrinsic water use efficiency, we agree that a more\ndetailed explanation of the model–observation mismatch is useful. Following your suggestion, we\nhave added the explanation and provided a model-observation discrepancy analysis in the revised\nmanuscript (Lines 241–242, 262–279, 279–284).\nFirstly, it is well established that the increase in atmospheric CO concentration (C ) is the primary\n2 a\ndriver of the global increase in iWUE, as elevated C reduces leaf stomatal conductance, thereby\na\nenhancing iWUE. However, through separate analyses of observational data and model\nsimulations, we found that the model-estimated sensitivity of iWUE to C was significantly higher\na\nthan that of observations. We have described this overestimation in Lines 241–242:\n“This mismatch was partially explained by the model’s overestimation of the sensitivity of\niWUE to C changes (Fig. S5).”\nmod a\na a\nareas with significance sensitivity (p < 0.05) (b) and proportion of positive and negative\nsensitivities (c) for iWUE response to C rate of change (diWUE /dC ) estimated from isotope\na obs a\nobservations. Spatial distribution (d), areas with significance sensitivity (p < 0.05) (e) and\nproportion of change (f) for iWUE response to C rate of change (diWUE /dC ) estimated by\na mod a\nthe P-model. Spatial differences of observation-based minus model-based estimates\n(diWUE /dC – diWUE /dC )(g). N: negative trend; N*: significant negative trend; P: positive\nobs a mod a\ntrend; P*: significant positive trend. Non-natural vegetation areas and deserts with less than 10%\nvegetation coverage were masked, as were grid units dominated by C vegetation.\nrepresents the ratio of unit cost of carboxylation to transpiration. In the current implementation, a\nfixed β value is applied across all plant functional types (PFTs), despite evidence that β varies\nsignificantly among biomes and species groups. For example, angiosperms tend to have higher β\nthan gymnosperms, and deciduous species higher than evergreens. This simplification introduces\nstructural bias, particularly in the prediction of C/C ratios. We have discussed the uncertainty of\ni a\nthe parameter β in the P-model and the resulting overestimation bias in the model in Lines 262–\n279:\nLines 262–279:\n“These differences, alongside with the P-model overestimation of iWUE sensitivity to rising C\na\n(Fig. S4), may be related to the use of fixed computational parameters. Calculated parameters with\ncomplex background information in the model (e.g., β in ξ = √β(K+Γ*)/1.6(1+ g /g ), see Eq. 10)\ns m\ncan also lead to biased results when characterized as fixed values. β is the ratio of the unit cost of\nplant carboxylation to transpiration at 25°C (dimensionless). β represents a critical computational\nparameter of the P-model, which is directly related to the accurate description of changes in plant\nstomatal conductance and net photosynthetic rate in response to environmental gradients (i.e., the\naccurate quantification of C/C )13.\ni a\nIn fact, there are large differences in β values across different biomes; for instance, angiosperms\nhave higher β than gymnosperms, while deciduous plants have higher β than evergreen trees14.\nHowever, a single β value is currently used for all PFTs in the P-model developed based on the\noptimality theory13, 15. The use of a single value introduces prediction bias, as demonstrated by a\nstudy on plant leaf traits along an altitudinal gradient based on optimality theory, which found that\na single β was the main cause of the C/C uncertainty16. It was observed that the use of β estimated\ni a\nfrom compiled data to predict C/C resulted in a mismatch with observations. Predicted C/C\ni a i a\nvalues were generally lower than the median observed values at most sites16. Furthermore, our\nsensitivity analysis of iWUE to β showed that iWUE decreased with increasing β (Fig. 5). This\nsuggests that the default β value used in the P-model may be too low for most PFTs or climatic\nzones.”\nmod\ncontinuously from 0 to 500, with all other parameters held constant (a), Scenario II: β fixed at six\ndiscrete values, while other parameters varied continuously (b).\nFinally, we addressed the model’s overestimation of VPD dominance in driving iWUE variability.\nAs noted in Lines 279–284, the P-model attributes iWUE changes to atmospheric dryness (VPD)\nmod\nin over 75% of global C vegetation areas, a much larger fraction than observed in our isotope-\n3\nbased empirical data. This discrepancy likely arises because the standard version of P-model\naccounts only for atmospheric water demand (via VPD), but omits the influence of soil water stress.\nWe now include:\nLines 279–284:\n“Meanwhile, the model results showed that atmospheric drought dominated iWUE changes in\nmod\nmore than 75% of the global C vegetation cover, which may be caused by the fact that the\n3\noptimality theory model only considered the effect of atmospheric water demand on C/C and\ni a\nignored the effect of soil water stress. Previous studies have shown that β decreases with increasing\nsoil dryness, indicating that the cost of plant water transport increases as SM declines14.”\n1.5 The machine learning random forest model lacks an independent validation dataset. While\ncross-validation is useful, an independent test dataset would provide stronger validation.\nThank you for this constructive suggestion. To assess the robustness of our random forest model\nin predicting the leaf Δ13C isoscapes of C plants, we have implemented three validation\n3\napproaches: (1) split-sample validation using a training set (70%) and an independent test set\n(30%), (2) 10-fold cross-validation, and (3) 100 random resampling. The independent test set\nconsidered in the first approach was not used during model training and served to evaluate model\nperformance. Results were consistent across all approaches, with the training/test split yielding an\nR of 0.77 and RMSE of 1.27 ‰, comparable to other validation methods (Lines 468-472, 600-605;\nLines 468-472:\n“We used three evaluation schemes, i) training set (70%)-test set (30%), ii) ten-fold cross-\nvalidation, iii) 100 random resampling, and four statistical criteria, namely, correlation coefficient\n(R), Nash Sutcliffe efficiency, root-mean-square error (RMSE), and mean absolute error (MAE),\nto evaluate the performance of the constructed Random Forest model for predicting the\nperformance of the C plant leaves for Δ13C isoscapes.”\n“The three validation methods (including training set-test set, tenfold cross-validation, and random\nresampling validation) and four validation indexes (including correlation coefficient (R), Nash\ncoefficient (Nash Sutcliffe), root-mean-square error (RMSE), and mean absolute error (MAE))\nproved the high predictive accuracy of our prediction model, in which the R was between 0.75-\n0.77, the Nash were both 0.59, RMSE between 1.22-1.28‰, and MAE between 0.92-0.97‰ (Fig.\nS10).”\nSupplementary Fig. 10. Analysis of predictive modeling (Random Forest) of Δ13C isoscapes of C\n3\nplant leaves based on three validation methods. Parameters encompass correlation coefficient (R)\n(a), Nash-Sutcliffe efficiency coefficient (Nash-Sutcliffe) (b), root mean square error (RMSE) (c),\nand mean absolute error (MAE) (d).\nIn addition, following your suggestion, we further validated our predicted Δ13C isoscapes using\nindependent datasets provided by Cornwell et al. (2018)11 and Lavergne et al. (2022)12. The dataset\nfrom Cornwell et al. (2018)11 is continuously updated (https://github.com/traitecoevo/leaf13C). To\nensure independence, we first removed any duplicated records shared between the two sources.\nexcluded these records and retained only pre-2001 observations (1951–2000). Third, we averaged\nΔ¹³C values for duplicate site–year combinations to match the format of our multi-year mean\npredictions. The resulting independent validation dataset spanned a wide geographic range (45.3°S\nto 75.3°N; 165.6°W to 179.4°E). Comparison with our predicted mean Δ13C for 2001–2020\nshowed a significant correlation (R = 0.69, p < 0.05), supporting the robustness of our results. In\nthe revised manuscript, we have now included this independent validation analysis as suggested.\nLines 624–628:\n“Additionally, we used leaf Δ13C data from Cornwell et al. (2018)11 and Lavergne et al. (2022)12\nthat were not included in our model development (i.e., pre-2001) as an independent validation\ndataset to compare with our model-predicted multi-year mean Δ13C (2001–2020). We removed\nduplicates and averaged values for repeated site–year records. The comparison yielded a\nsignificant correlation (R = 0.69, p < 0.05) (Fig. S13), further supporting the robustness of the\nmodel.”\nSupplementary Fig. 13. Comparison between predicted Δ13C from the optimal random forest\nmodel and observed Δ13C from the independent dataset compiled from Cornwell et al. (2018)11\nand Lavergne et al. (2022)12. Observations used in model development (i.e., from 2001–2020)\nwere excluded. Duplicate records were removed, and Δ¹³C values from the same site and year were\naveraged.\n(Cornwell et al.11 and Lavergne et al.12) in the Supplementary Text (Lines 619-624). Specifically,\nwe noted that the spatial distribution patterns were consistent, showing higher Δ13C in tropical\nrainforests and lower values in cold, arid, and high-altitude regions (Fig. S12a, b). Our predicted\nglobal average Δ13C (20.6 ± 0.3‰) also closely matches the estimate from Lavergne et al.¹² (20.9\n± 1.5‰), which was derived from the JULES land surface model that incorporates optimal\nstomatal regulation and mesophyll conductance.\nLines 619-624:\n“We found that the mapped Δ13C isoscapes displayed spatially similar patterns to the previous\nproducts11, 12, with higher Δ13C values in tropical rainforests and lower values in cold and arid\nmountainous regions and high latitudes (Fig. S12a, b). Our predicted global average Δ13C (20.6 ±\n0.3‰) closely matched the global average (20.9 ± 1.5‰) estimated by Lavergne et al.12 based on\nthe JULES model, which incorporates the Prentice stomatal model17 and accounts for the effect of\nmesophyll conductance.”\n1.6 The sampling distribution is not evenly spread globally—some regions (e.g., boreal forests,\nAfrica) are underrepresented. Discuss the implications of uneven sampling distribution and how\nthis may bias the results.\nThank you for this suggestion. In the revised manuscript, we have added a detailed discussion on\nthe implications of uneven sampling distribution and its potential impact on the model outputs per\nyour suggestion.\nLines 655–670:\n“However, the uneven spatial distribution of C leaf isotope samples used for training our\n3\nupscaling models may affect the stability of our upscaled estimates. Therefore, the uncertainty in\nthe model construction process was assessed by performing iterative model training and different\nvalidation methods (including training–testing splits, ten-fold cross-validation, and random\nresampling validation), as well as spatial pattern comparisons (Figs. S10, S12). To better present\nthe spatial uncertainty, we also plotted the coefficient of variation (CV = (standard\ndeviation/average) × 100) of the output values predicted by the algorithms. The pattern of global\nCVs was uniform and low overall, with higher CVs in regions such as the southern edge of the\nAfrican Sahel, the Central Asian highlands, and the Tibetan Plateau, which may be due to the\nsparse distribution of sample points in these regions (Fig. S12). Importantly, this uncertainty was\nnot expected to affect the relative changes in the temporal scales of the predicted results, as\nevidenced by the range of the predicted results and the robustness of the size patterns and their\nannual variations across PFTs, as well as the consistency with the independent validation dataset\n(Figs. S9-S14). Overall, despite these limitations, the isoscapes of C plant Δ13C predicted by our\n3\n11\nmethod matches well with observed leaf isotope data and with the results from both a linear\nstatistical model11 and an Earth system model12, suggesting that our isoscapes better captures the\nglobal spatial variability of C plants Δ13C.”\n3\n1.7 The manuscript has frequent grammatical errors (e.g., \"water sress\" instead of \"water stress\";\n\"trand\" instead of \"trend\") or weird sentence composition (e.g., Line 51 “Atmospheric dryness\n(vapor pressure deficit, VPD) increases the iWUE across a larger area than soil dryness“). Some\nredundant statements repeat findings across different sections. Improve language clarity and\neliminate typographical errors.\nWe appreciate your careful reading and helpful suggestions. As you suggested, we have thoroughly\nproofread the entire manuscript to correct typographical errors, improve sentence clarity, and\neliminate redundancy. For example, the sentence on Line 52 has been revised to: “Vapor pressure\ndeficit (VPD) was found to influence iWUE across a broader area than soil moisture.” We also\ndouble-checked this the manuscript. Thank you again for the careful read.\n1.8 Line 61. I would say that net assimilation rate or net photosynthetic rate is more precise than\nphotosynthetic rate.\nChanged as the reviewer suggested.\n1.9 Line 65. The leaf 13C discrimination is also linearly related to leaf iWUE, maybe easier for\nreader to comprehend (https://doi.org/10.1111/ppl.14619). Moreover, I don’t understand how cited\nreference by Adams et al. links carbon discrimination to iWUE, it should be replaced with\nsomething more relevant like https://doi.org/10.2135/cropsci2005.0119.\nThank you for pointing this out. We have modified this sentence to: “In C plant photosynthesis,\n3\nleaf 13C discrimination (Δ13C) is linearly related to leaf iWUE18, 19, 20, 21. (Line 67)”. We have\nalso replaced the references with more relevant ones, as you suggested.\n1.10 Line 79: Only theory? There are multiple experimental studies that showed that.\nThank you for this suggestion. We have modified this sentence to: “Additionally, both theory and\nexperimental studies suggest that increased atmospheric and soil drought may enhance iWUE1, 22,\n23.” (Lines 79–80).\n1.11 Line 280: „We found that the global distribution of leaf iWUE aligned with physiological\nresponses of plants: high iWUE values in cold and arid regions, and lower values in warm and\nhumid areas. “I firstly don’t understand what the authors mean by physiological responses here in\n12\nconnection to climate zones. Secondly, how does this align with the found positive trends? You\nalso mentioned beforehand that the boreal forests had the lowest iWUE, now cold regions have\nhigh iWUE?\nThank you for pointing this out. To improve clarity and avoid confusion, we have removed the\nphrase “physiological responses” and revised the sentence as follows: “We found that the global\ndistribution of mean leaf iWUE shows higher values in cold and arid regions, and lower values in\nwarm and humid areas.”\nThis revision simplifies the statement and better distinguishes between spatial patterns of mean\niWUE and its temporal trends. As noted in the manuscript, the spatial distribution (mean iWUE)\nand its rate of change over time (i.e., the trend) are distinct metrics. While cold and arid regions\ngenerally exhibit higher mean iWUE, the strongest increasing trends in iWUE do not necessarily\noccur in those regions. Instead, trends reflect the combined effects of regional climate trajectories\nand vegetation responses over time.\nAdditionally, we would like to clarify that the manuscript does not state that boreal forests have\nthe lowest iWUE values. As described in Lines 115–120, the lowest mean iWUE values are found\nin evergreen broadleaf forests located in warm and humid tropical regions.\n2.1 This manuscript provides a valuable global assessment of leaf-level intrinsic water use\nefficiency (iWUE) across C3 vegetation, leveraging a large leaf δ¹³C dataset, machine-learning\nupscaling, and established physiological equations. It offers a novel product at 0.05° resolution\nfrom 2001 to 2020, revealing strong spatiotemporal patterns of iWUE and clear evidence that VPD,\nrather than soil moisture, dominates iWUE changes in most regions. These results will be\nbeneficial to both carbon–water cycle modeling and ecological research.\nMajor Strengths\n2.2 The new global iWUE dataset covers diverse ecosystems with a coherent methodology.\n2.3 Incorporation of mesophyll conductance and nonphotosynthetic fractionations strengthens\nphysiological realism.\n2.4 Comparisons with the P-model highlight uncertainties, underscoring the need to refine certain\nmodel parameters.\nWe greatly appreciate the reviewer’s positive comments and support of the manuscript.\nAreas for Improvement\n2.5 Data Coverage and Uncertainty: Provide clearer uncertainty estimates (e.g., standard-error\nmaps), especially in under-sampled regions (e.g., parts of Africa, high latitudes).\nThank you very much for this suggestion. In the revised manuscript, we have provided clearer\nuncertainty and data coverage metrics (Figs. R2, R3) and a more detailed discussion following\nyour suggestion.\nLines 655–670:\n“However, the uneven spatial distribution of C leaf isotope samples used for training our\n3\nupscaling models may affect the stability of our upscaled estimates. Therefore, the uncertainty in\nthe model construction process was assessed by performing iterative model training and different\nvalidation methods (including training–testing splits, ten-fold cross-validation, and random\nresampling validation), as well as spatial pattern comparisons (Figs. S10, S12). To better present\nthe spatial uncertainty, we also plotted the coefficient of variation (CV = (standard\ndeviation/average) × 100) of the output values predicted by the algorithms. The pattern of global\nCVs was uniform and low overall, with higher CVs in regions such as the southern edge of the\nAfrican Sahel, the Central Asian highlands, and the Tibetan Plateau, which may be due to the\nsparse distribution of sample points in these regions (Fig. S12). Importantly, this uncertainty was\nnot expected to affect the relative changes in the temporal scales of the predicted results, as\nevidenced by the range of the predicted results and the robustness of the size patterns and their\n14\nannual variations across PFTs, as well as the consistency with the independent validation dataset\n(Figs. S9-S14). Overall, despite these limitations, the isoscapes of C plant Δ13C predicted by our\n3\nmethod matches well with observed leaf isotope data and with the results from both a linear\nstatistical model11 and an Earth system model12, suggesting that our isoscapes better capture the\nglobal spatial variability of C plant Δ13C.”\nwere grid units dominated by C vegetation.\ndistribution (b) of Δ13C for C plants over 2001-2020. Non-natural vegetation areas and deserts\n3\nwith less than 10% vegetation coverage were masked, as were grid units dominated by C\n4\nvegetation.\n2.6 Atmospheric δ¹³C Assumption: Justify using Mauna Loa δ¹³C globally or discuss potential\n15\nbiases in regions far from Hawaii or with unique CO₂ sources/sinks.\nThank you for raising this important point. We agree that using atmospheric δ13C data from the\natm\nMauna Loa station as a global proxy may introduce regional uncertainties. Due to the current\nabsence of reliable, high-resolution, spatially explicit global δ13C datasets, most global-scale\natm\nand tree-ring–based studies have similarly relied on either ice-core reconstruction (for one global\nvalue of the -industrial periods) or long-term atmospheric observations from Mauna Loa (available\nsince the late 1950s). Since our study focuses on the period 2001–2020, we adopted the Mauna\nLoa δ13C record as a practical and widely accepted surrogate, following the precedent set by\natm\nearlier large-scale studies. We have clarified this rationale in the revised Methods section.\nLines 357–362:\n“Due to the lack of a reliable and high-resolution global map of atmospheric CO concentrations\n2\n(C , ppm) and carbon isotope ratios (δ13C , ‰), we used monthly C and δ13C from the Mauna\na a a a\nLoa station (19.54°N, 155.57°W, 3397 m a.s.l., United States), data provided by the National\nOceanic and Atmospheric Administration (NOAA), for which annual averages were calculated.\nSubsequently, the annual averages were aggregated into raster data with a global 0.05° spatial\nresolution for the period 2001–2020 using the “raster” package in R.”\nAdditionally, we also discussed the potential bias of using only Mauna Loa atmospheric δ13C in\nthe revised manuscript.\nLines 714–724:\n“Finally, the use of δ13C values from the Mauna Loa station may affect the accuracy of leaf Δ13C\na\nestimates in regions with different atmospheric CO sources or isotopic compositions24. For\n2\nexample, δ13C at Mauna Loa tends to be higher than in regions like the EBF canopy of Amazonia,\na\npotentially leading to overestimated Δ13C values in these areas. A canopy-level study in French\nGuiana reported Δ13C values of ~26.4‰ in mid-canopy and 22–23‰ in the upper canopy25,\nwhereas our modeled maximum Δ13C in the Amazon reached 28.11‰—a difference of\napproximately 3–6‰, which carries over to iWUE estimation. While this introduces potential\nregional biases, it does not compromise the stability of our machine-learning upscaling framework\nor the global spatial patterns and relative trends in iWUE. We emphasize the urgent need for 4D\nglobal datasets of atmospheric CO concentration and isotopic composition (lon, lat, month,\n2\naltitude) to better constrain the discrimination maps based on plant material and δ13C .”\na\n2.7 Additional Confounding Factors: Briefly note how nitrogen deposition, forest age, or ontogeny\nmight affect iWUE trends in certain regions.\nage, ontogeny, and nitrogen deposition on observed iWUE trends.\nLines 671–683:\n“Second, we did not explicitly account for other potential factors that may influence temporal\ndynamics in leaf Δ13C and iWUE, such as vegetation ontogeny (e.g., changes in tree height and\nlight environment), forest age, and nitrogen (N) deposition. Ontogenetic changes in tree\nphysiology—including increasing hydraulic resistance and evolving light availability—can alter\nstomatal conductance and photosynthesis, thereby affecting Δ13C and iWUE. Some studies have\nfound that Δ13C values of tree rings changes significantly with tree height and light conditions,\nrather than an effect of tree age itself, specifically that iWUE increases with tree growth26.\nTherefore, many tree-ring isotope studies exclude the first 20–50 years of growth to minimize\nontogenetic effects, assuming that post-juvenile rings are less biased27, 28. However, other research\nindicates that ontogenetic effects may persist beyond early developmental stages and continue to\ninfluence isotope signatures throughout the tree's life cycle26. It is worth noting that this\nontogenetic effect does not significantly affect absolute changes in iWUE over time29. Meanwhile,\nincreased N deposition over the past decades could also affect photosynthesis and thus iWUE\ntrend30.”\n2.8 Methodological Clarification: Include a short sensitivity analysis (e.g., varying the gₛ/gₘ ratio)\nand indicate how partial correlations (VPD vs. SM) cope with variable collinearity.\nThank you very much for this suggestion. In the revised manuscript, we have added a sensitivity\nanalysis and discussed how our partial correlation analysis copes with variable collinearity\nfollowing your suggestion.\nSensitivity analysis:\nFollowing Ma et al. (2021)31, we conducted a sensitivity analysis to evaluate how variations in the\ng /g ratio affect the estimation of iWUE. We fixed all other parameters using their average values\ns m\nfor the year 2020 and varied the g /g ratio across four representative levels with a very generous\ns m\nrange (0.4, 0.79, 1.2, and 1.8). In addition, based on the sensitivity test framework proposed by\nMa et al. (2021)31, we further explored how iWUE responds to changes in g /g under different\ns m\nvalues of leaf Δ13C (ranging from 18‰ to 22‰). Results showed that iWUE consistently\nLeaf\ndecreases with increasing g /g , with the strongest sensitivity observed at low values (Fig. S16).\ns m\nThis trend aligns with theoretical expectations and confirms findings from Ma et al. (2021)31.\nIn our study, we adopted the representative value of 0.79 for g /g , as proposed by Ma et al.\ns m\n(2021)31. This value was derived from an analysis of 198 paired observations of g and g across\ns m\n32 species and 15 experimental studies under different drought stress conditions. Their analysis\n17\nshowed that the g /g ratio did not vary significantly across plant functional types or water\ns m\ntreatments, with a global mean of 0.79 ± 0.07 (95% CI). Even though both g and g declined by\ns m\n~50% under drought, their ratio remained relatively constant. The suitability of this value was\nfurther supported by an independent experimental dataset involving seven species31. Therefore,\nusing a fixed g /g ratio of 0.79 provides a practical and reasonably accurate way to account for\ns m\ng effects when directly measured g data are unavailable, and avoids the upward bias in iWUE\nm m\nestimation associated with infinite-g assumptions.\nm\nIt is worth noting, however, that as pointed out in Ma et al. (2021)31, most reported g /g values\ns m\nare based on instantaneous measurements under steady-state photosynthesis and may not fully\ncapture short-term responses to acute stress (e.g., extreme drought or high VPD). Thus, while\ncaution is needed in short-term applications, this value is appropriate and robust for our long-term,\nlarge-scale iWUE assessment.\nSupplementary Fig. 16. Sensitivity of iWUE to the g /g ratio. Modeled iWUE values under\ns m\ndifferent g /g ratios (0.4, 0.79, 1.2, and 1.6), with all other parameters fixed at 2020 means (a).\ns m\nVariation of iWUE across a range of leaf carbon isotope discrimination (Δ13C = 18‰ – 22‰)\nLeaf\nunder the same four g /g ratios (b); panel b follows the framework of Ma et al. (2021)31.\ns m\nFurthermore, we replaced the g /g ratio used in the original P-model (1.4; Wang et al., 2017)13\ns m\nwith the optimal value of 0.79 derived by Ma et al. (2021)31 to assess the P-model’s sensitivity to\nthis parameter. The results show that, irrespective of the g /g value used, the P-model consistently\ns m\noverestimated mean iWUE compared to leaf δ13C-based observations (Fig. S17), suggesting\nlimited sensitivity of model performance to this parameter.\nunder different g /g values.\ns m\nWe have added this discussion in the manuscript (Lines 692–713) and included results in\nSupplementary Figures S16 and S17.\nLines 692–713:\n“In addition, sensitivity analysis showed that iWUE decreases with increasing g /g (Fig. S16),\ns m\nconsistent with the theoretical predictions by Ma et al.31, who also demonstrated that the elasticity\nof iWUE to g /g is highest at low values of g /g and tends to level off at higher values. In this\ns m s m\nstudy, we adopted the representative g /g value of 0.79 proposed by Ma et al.31 as a model\ns m\nparameter. This value was derived from 198 paired measurements of stomatal and mesophyll\nconductance across 32 species under varying drought stress conditions, and was found to be\nconsistent across plant functional types and watering treatments, with a global mean of 0.79 ± 0.07\n(95% confidence interval). Although both g and g declined by ~50% under drought, their ratio\ns m\nremained relatively stable. The robustness of this value was further supported by an independent\nexperimental dataset comprising seven species31. Therefore, using a fixed g /g ratio of 0.79\ns m\nprovides a practical and reasonably accurate means of accounting for the influence of gm when\ndirect measurements are unavailable, and avoids the overestimation of iWUE associated with the\nassumption of infinite g . However, as noted by Ma et al.31, most reported g /g values are based\nm s m\non instantaneous measurements under steady-state photosynthesis and may not fully capture short-\nterm responses to acute stress (e.g., extreme drought or high VPD). Caution is therefore warranted\nwhen applying a fixed g /g ratio to assess short-term iWUE dynamics. Nonetheless, for the\ns m\n19\npurposes of this study—focused on long-term, large-scale iWUE estimates—the use of a\nrepresentative value of 0.79 remains appropriate and robust. To further assess model sensitivity,\nwe replaced the default g /g value of 1.4 used in the original P-model (Wang et al.13) with the\ns m\noptimal value of 0.79 from Ma et al.31. The results showed that regardless of the g /g value used,\ns m\nthe P-model consistently overestimated mean iWUE compared to leaf δ13C-based observations\n(Fig. S17), suggesting limited model sensitivity to variations in g /g .”\ns m\nVariable collinearity:\nTo investigate the influence and control regions of atmospheric dryness versus soil dryness on\nglobal iWUE, we employed a partial correlation analysis, which quantifies the strength of the\nrelationship between iWUE and each water stress variable (VPD and SM), while controlling for\nthe (linear) influence of the other.\nAlthough VPD and SM can be collinear due to shared climatic influences, partial correlation\neffectively isolates the unique contribution of each variable by removing their mutual linear effects.\nThis approach offers a robust alternative to simple correlation or multiple regression in identifying\nindependent drivers of iWUE variability. We have clarified this approach in the revised Methods\nsection.\nLines 542-546:\n“We used partial correlation analysis to assess the effects of atmospheric dryness (VPD) and soil\ndryness (SM) on iWUE means. This method controlled for the influence of one variable while\nevaluating the partial relationship of the other with the response variable 32. To identify the\ndominant driver of iWUE, we ranked VPD and SM based on the absolute values of their partial\ncorrelation coefficients in each grid cell33.”\n2.9 Overall, the methods are solid, and the findings have broad significance. Addressing these\npoints in a revision will further strengthen the manuscript and clarify its conclusions.\nWe greatly appreciate the reviewer’s positive comments and support of the manuscript. We have\ncarefully followed the reviewer’s suggestions and revised the manuscript as above. As a\nconsequence, our manuscript has been considerably improved. Please see our detailed responses\nabove.\n1. Kannenberg SA, Driscoll AW, Szejner P, Anderegg WRL, Ehleringer JR. Rapid increases\nin shrubland and forest intrinsic water-use efficiency during an ongoing megadrought. Proc\nNatl Acad Sci U S A 118, (2021).\n2. Körner C, Diemer M. In situ Photosynthetic Responses to Light, Temperature and carbon\ndioxide in herbaceous plants from low and high altitude. Funct Ecol 1, 179-194 (1987).\n3. Wang H, et al. Photosynthetic responses to altitude: an explanation based on optimality\nprinciples. New Phytol 213, 976-982 (2017).\n4. Santos CMD, Endres L, Ferreira VM, Silva JV, Rolim EV, Wanderley HCL.\nPhotosynthetic capacity and water use efficiency in Ricinus communis (L.) under drought\nstress in semi-humid and semi-arid areas. An Acad Bras Cienc 89, (2017).\n5. Yang YJ, et al. Evolution of stomatal closure to optimize water-use efficiency in response\nto dehydration in ferns and seed plants. New Phytol 230, 2001-2010 (2021).\n6. Adams MA, Buckley TN, Turnbull TL. Rainfall drives variation in rates of change in\nintrinsic water use efficiency of tropical forests. Nat Commun 10, 3661 (2019).\n7. Belmecheri S, et al. Precipitation alters the CO effect on water-use efficiency of temperate\n2\nforests. Glob Change Biol 27, 1560-1571 (2021).\n8. Diefendorf AF, Mueller KE, Wing SL, Koch PL, Freeman KH. Global patterns in leaf 13C\ndiscrimination and implications for studies of past and future climate. Proc Natl Acad Sci\nU S A 107, 5738-5743 (2010).\n9. Farquhar GD, O'Leary MH, Berry JA. On the relationship between carbon isotope\ndiscrimination and the intercellular carbon dioxide concentration in leaves. Funct Plant\nBiol 9, 121-137 (1982).\n10. Stewart G, Turnbull M, Schmidt S, Erskine P. 13C natural abundance in plant communities\nalong a rainfall gradient: a biological integrator of water availability. Funct Plant Biol 22,\n51-55 (1995).\n11. Cornwell WK, et al. Climate and soils together regulate photosynthetic carbon isotope\ndiscrimination within C3 plants worldwide. Glob Ecol Biogeogr 27, 1056-1067 (2018).\n12. Lavergne A, Hemming D, Prentice IC, Guerrieri R, Oliver RJ, Graven H. Global decadal\nvariability of plant carbon isotope discrimination and its link to gross primary production.\nGlob Change Biol 28, 524-541 (2022).\n13. Wang H, et al. Towards a universal model for carbon dioxide uptake by plants. Nat Plants\n3, 734-741 (2017).\n21\n14. Lavergne A, Sandoval D, Hare VJ, Graven H, Prentice IC. Impacts of soil water stress on\nthe acclimated stomatal limitation of photosynthesis: Insights from stable carbon isotope\ndata. Glob Change Biol 26, 7158-7172 (2020).\n15. Stocker BD, et al. P-model v1.0: an optimality-based light use efficiency model for\nsimulating ecosystem gross primary production. Geosci Model Dev 13, 1545-1581 (2020).\n16. Xu H, Wang H, Prentice IC, Harrison SP, Wang G, Sun X. Predictability of leaf traits with\nclimate and elevation: a case study in Gongga Mountain, China. Tree Physiol 41, 1336-\n1352 (2021).\n17. Prentice IC, Dong N, Gleason SM, Maire V, Wright IJ. Balancing the costs of carbon gain\nand water transport: testing a new theoretical framework for plant functional ecology. Ecol\nLett 17, 82-91 (2014).\n18. Farquhar GD, Ehleringer JR, Hubick KT. Carbon Isotope Discrimination and\nPhotosynthesis. Annu Rev Plant Biol 40, 503-537 (1989).\n19. Wang X, et al. Differences in the patterns and mechanisms of leaf and ecosystem-scale\nwater use efficiencies on the Qinghai-Tibet Plateau. Catena 222, (2023).\n20. Impa SM, Nadaradjan S, Boominathan P, Shashidhar G, Bindumadhava H, Sheshshayee\nMS. Carbon isotope discrimination accurately reflects variability in WUE measured at a\nwhole plant level in rice. Crop Sci 45, 2517-2522 (2005).\n21. Petrík P, et al. Linking stomatal size and density to water use efficiency and leaf carbon\nisotope ratio in juvenile and mature trees. Physiol Plant 176, e14619 (2024).\n22. Driscoll AW, Bitter NQ, Sandquist DR, Ehleringer JR. Multidecadal records of intrinsic\nwater-use efficiency in the desert shrub Encelia farinosa reveal strong responses to climate\nchange. Proc Natl Acad Sci U S A 117, 18161-18168 (2020).\n23. Grossiord C, et al. Plant responses to rising vapor pressure deficit. New Phytol 226, 1550-\n1566 (2020).\n24. Buchmann N, Brooks JR, Ehleringer JR. Predicting daytime carbon isotope ratios of\natmospheric CO within forest canopies. Funct Ecol 16, 49-57 (2002).\n2\n25. Buchmann N, Guehl JM, Barigah TS, Ehleringer JR. Interseasonal comparison of CO\n2\nconcentrations, isotopic composition, and carbon dynamics in an Amazonian rainforest\n(French Guiana). Oecologia 110, 120-131 (1997).\n26. Brienen RJW, et al. Tree height strongly affects estimates of water-use efficiency\nresponses to climate and CO using isotopes. Nat Commun 8, 288 (2017).\n2\n22\n27. Adams MA, Buckley TN, Turnbull TL. Diminishing CO -driven gains in water-use\n2\nefficiency of global forests. Nat Clim Chang 10, 466-471 (2020).\n28. Saurer M, et al. Spatial variability and temporal trends in water-use efficiency of European\nforests. Glob Change Biol 20, 3700-3712 (2014).\n29. Gong XY, et al. Overestimated gains in water-use efficiency by global forests. Glob\nChange Biol 28, 4923-4934 (2022).\n30. Walker AP, et al. Integrating the evidence for a terrestrial carbon sink caused by increasing\natmospheric CO . New Phytol 229, 2413-2445 (2021).\n2\n31. Ma WT, et al. Accounting for mesophyll conductance substantially improves 13C-based\nestimates of intrinsic water-use efficiency. New Phytol 229, 1326-1338 (2021).\n32. Piao S, Wang J, Li X, Xu H, Zhang Y. Spatio-temporal changes in the speed of canopy\ndevelopment and senescence in temperate China. Glob Change Biol 28, 7366-7375 (2022).\n33. Jiao W, Wang L, Smith WK, Chang Q, Wang H, D'Odorico P. Observed increasing water\nconstraint on vegetation growth over the last three decades. Nat Commun 12, 3777 (2021).\n23"}, {"round": 2, "reviewer_comments": {"reviewer_1": "Authors addressed all my comments. I am suggesting acceptance of the manuscript for publication.", "reviewer_2": "Reviewer Report\nThe authors present a global analysis of leaf-level intrinsic water use efficiency (iWUE) based on upscaling carbon isotope data with machine learning and comparing these estimates with an ecological optimality model (P-model). The topic is timely and highly relevant for understanding carbon–water coupling under climate change. The revisions have substantially improved the manuscript: the VPD–temperature confounding has been addressed, independent validation datasets have been included, and the discussion of uncertainties is stronger. Overall, the manuscript now makes a valuable contribution to global change biology.\nThat said, I still see several areas that require attention before publication:\nVPD vs. temperature effects. The residual analysis is useful, but it only accounts for linear temperature effects. Nonlinear temperature responses of assimilation (optimum curve) are not considered, which could bias results in boreal regions. This limitation should be discussed explicitly, and ideally a sensitivity test with nonlinear terms provided.\nNovelty of the framework. The claim of a “novel framework” is somewhat overstated. The main innovation is the global upscaling of foliage isotope data via machine learning, which complements but does not replace tree-ring and flux-tower approaches. The language in the abstract and conclusions should be toned down accordingly.\nP-model bias and implications. The added discussion of β parameterization and Ca sensitivity is helpful, but the broader implications remain underdeveloped. If the P-model systematically overestimates iWUE trends, Earth system models adopting similar assumptions could misrepresent future carbon–water feedbacks. This should be stated explicitly.\nSampling distribution. The new coefficient of variation maps are useful, but the uneven sampling distribution still affects confidence in biome-level conclusions (e.g., slow grassland vs. fast evergreen broadleaf responses). Please qualify which biome-level results are robust and which are more tentative."}, "author_response": "Response to Reviewers\nWe greatly appreciate the constructive comments from the reviewers invitation to submit a\nrevised version. We have thoroughly revised the manuscript following the reviewers’ suggestions\n. Please see below our point-to-point responses in blue text followin g the reviewers’comments.\nThe line nu mbers referred to are for the clean version of th e revised manuscript (non-track-\nchange version).\n---------------------------------------------------------------------------------------------------------------------\nReviewer #1 (Remarks to the Author):\n1.1 Authors addressed all my comments. I am suggesting acceptance of the manuscript for\npublication.\nWe thank the reviewer for supporting this paper and for recognizing our effort to improve it.\nReviewer #2 (Remarks to the Author):\nReviewer Report\n2.1 The authors present a global analysis of leaf-level intrinsic water use efficiency (iWUE) based\non upscaling carbon isotope data with machine learning and comparing these estimates with an\necological optimality model (P-model). The topic is timely and highly relevant for understanding\ncarbon–water coupling under climate change. The revisions have substantially improved the\nmanuscript: the VPD–temperature confounding has been addressed, independent validation\ndatasets have been included, and the discussion of uncertainties is stronger. Overall, the manuscript\nnow makes a valuable contribution to global change biology.\nWe greatly appreciate the reviewer’s positive comments and support of the manuscript. We have\nthoroughly revised the manuscript following your suggestions. Please see below our detailed\nresponses.\n2.2 That said, I still see several areas that require attention before publication:\nVPD vs. temperature effects. The residual analysis is useful, but it only accounts for linear\n1\ntemperature effects. Nonlinear temperature responses of assimilation (optimum curve) are not\nconsidered, which could bias results in boreal regions. This limitation should be discussed\nexplicitly, and ideally a sensitivity test with nonlinear terms provided.\nThank you very much for your suggestions. Following your suggestions, we have incorporated the\nnon-linear effects of temperature when calculating the VPD residuals. We regressed annual mean\ntemperature (MAT) and its squared term (MAT2) against VPD and extracted the residuals, which\nrepresent water stress signals independent of both linear and non-linear temperature effects. These\nresiduals were then used in partial correlation analyses with iWUE. This approach more\nthoroughly isolates the confounding effects of temperature and allows a clearer assessment of the\nindependent contributions of VPD and SM to iWUE variability, particularly in boreal regions\nwhere non-linear temperature responses may be pronounced. Based on this new analysis, we found\nthat VPD remained the dominant driver of iWUE variability in 59% of C vegetation areas (Fig.\neffects were removed. These results reinforce our conclusion that VPD dominated iWUE changes\nacross a larger area compared to SM at the global scale.\nSupplementary Fig. 2. Effects of vapor pressure deficit (VPD) and soil moisture (SM) removed\nfrom temperature influence on iWUE. The spatial distribution of partial correlation coefficients of\nVPD and SM controlling changes in iWUE (a), the proportion of the area controlled by VPD or\nSM (b). Non-natural vegetation areas and deserts with less than 10% vegetation coverage were\nmasked, as were grid units dominated by C vegetation.\ncompared partial correlation coefficients between iWUE and VPD derived from two models: one\ncontrolling for both linear and quadratic terms of mean annual temperature (MAT + MAT2), and\n2\nthe other controlling for linear MAT only. The difference between these two estimates\n(Difference ) provides a quantitative measure of the sensitivity of iWUE–VPD\nVPD–MAT\nrelationships to temperature nonlinearity. Results show that the difference distribution is\nsymmetrical and centered around zero, indicating that incorporating non-linear temperature\nresponses has only minor effects on most regions and thus supports the robustness of our main\nfindings (Fig. R1).\nFig. R1. Density plot of the difference in partial correlation coefficients (Difference )\nVPD–MAT\nbetween nonlinear (MAT + MAT2) and linear (MAT) temperature control models.\nWe have added discussion and updated the analysis (Supplementary Figure 2) into the revised\nmanuscript (Lines 546–552, Lines 216–219).\nLines 546–552:\n“In addition, given that the climatic effect of VPD inherently includes a temperature component,\nand that plants assimilation exhibit a nonlinear optimum response to temperature1, 2, this may bias\nthe attribution of iWUE in boreal regions. To address this, we further computed the correlations\nbetween iWUE and residuals of VPD and SM (obtained after removing both linear and non-linear\ninfluence of air temperature). Specifically, we first applied regression models including both linear\nand quadratic terms to remove temperature-driven components from VPD and SM, and then\ntreated the residuals as temperature-independent indicators of water stress.”\n“Partial correlation analysis based on residual VPD and SM, with temperature linear and non-\nlinear effects on VPD and SM removed, also confirmed this result, identifying VPD and SM as\nthe dominant drivers in 59% and 41% of C vegetation areas, respectively (Fig. S2).”\n3\n2.2 Novelty of the framework. The claim of a “novel framework” is somewhat overstated. The\nmain innovation is the global upscaling of foliage isotope data via machine learning, which\ncomplements but does not replace tree-ring and flux-tower approaches. The language in the\nabstract and conclusions should be toned down accordingly.\nFollowing your suggestion, we have revised the wording by removing “novel” from the abstract\nand “new” from the conclusion.\n2.3 P-model bias and implications. The added discussion of β parameterization and Ca sensitivity\nis helpful, but the broader implications remain underdeveloped. If the P-model systematically\noverestimates iWUE trends, Earth system models adopting similar assumptions could\nmisrepresent future carbon–water feedbacks. This should be stated explicitly.\nThank you for this suggestion. As per your suggestion, we have emphasized the insights of the P-\nmodel bias for Earth system models in the revised manuscript.\nLines 283–285:\n“This bias highlights the need to carefully evaluate Earth system models that adopt similar\nassumptions and parameterizations, as they may misrepresent future carbon water feedbacks under\nclimate change.”\n2.4 Sampling distribution. The new coefficient of variation maps are useful, but the uneven\nsampling distribution still affects confidence in biome-level conclusions (e.g., slow grassland vs.\nfast evergreen broadleaf responses). Please qualify which biome-level results are robust and which\nare more tentative.\nThank you very much for this suggestion. Following your suggestion, we performed an additional\nassessment of sampling distribution and uncertainty across different plant functional types to\nevaluate the robustness of biome-level results. Specifically, we first quantified the sample size and\nspatial coverage (latitude and longitude) of each biome to reflect sampling representativeness (the\nvegetation samples in this study were categorized into forests, grasslands, and shrublands during\ndata compilation). We then applied a bootstrap procedure (1000 resamples) to estimate the 95%\nconfidence intervals of Δ13C , thereby quantifying the uncertainty of the estimates. Finally, we\nLeaf\ncalculated the coefficient of variation (CV) as an indicator of within-biome variability.\nThe results show that forests had the largest sample size (n = 4038) and the widest spatial coverage,\nwith the narrowest confidence interval (20.6–20.7‰) and a CV of ~0.093, indicating the most\nrobust conclusions. Grasslands ranked second (n = 1107), also with broad spatial coverage and a\nsimilarly narrow confidence interval (19.2–19.4‰), yielding reliable results. By contrast,\nshrublands had the fewest samples (n = 819); although their CV (~0.092) was comparable to\ngrasslands, the more limited coverage renders these results more tentative (Supplementary Table\n3).\nSupplementary Table 3. Sampling distribution and assessment across plant functional types.\nLat. range Lon. range Mean Δ13C\nleaf\nPlant type n (samples) 95% CI (‰) CV\n(°) (°) (‰)\nForests 4038 -41.8~70.6 -140~153 20.6 20.6–20.7 0.093\nGrasslands 1107 -35.0~70.6 -140~148 19.3 19.2–19.4 0.091\nShrublands 819 -51.8~70.6 -156~151 19.7 19.5–19.8 0.092\nWe have added the Supplementary Table 3 and discussion into the revised manuscript (Lines 127–\n135 of supplementary files).\nLines 127–135 (supplementary files)\n“Additionally, differences in sample size and spatial distribution among vegetation types\ninfluenced the robustness of biocenosis-level inferences. Forests were the best represented (n =\n4038), with the broadest spatial coverage and the narrowest confidence interval (CI = 20.6–20.7‰,\nCV = 0.093), providing the strongest support for robust conclusions. Grasslands were less\nrepresented (n = 1107) but still yielded relatively reliable estimates (CI = 19.2–19.4‰, CV =\n0.091). Shrublands had the fewest samples (n = 819); although the coefficient of variation was\nsimilar to grasslands (0.092), limited spatial coverage reduced the robustness of shrubland-level\nresults compared with forests and grasslands.”\n1. Petek-Petrik A, et al. The combined effect of branch position, temperature, and VPD on\ngas exchange and water-use efficiency of Norway spruce. Biol Plant 67, 136-141 (2023).\n2. Xue L, et al. Analysis on water use efficiency of Populus euphratica forest ecosystem in\narid area. Theor Appl Climatol 145, 717-730 (2021).\n6"}, {"round": 3, "reviewer_comments": {"reviewer_2": "Thanks for your great efforts.\nMy questions have been addressed."}, "author_response": ""}]

Decision letter