Significance Elevated levels of atmospheric carbon dioxide affect plants directly by stimulating photosynthesis and reducing stomatal aperture. These direct effects trigger several more subtle, indirect effects via changes in soil moisture and plant structure. While such effects have been acknowledged, they have never been assessed quantitatively, partly due to the fact they are inseparable in field experiments. Here we show that the indirect effects of elevated CO 2 explain, on average, 28% of the total plant productivity response, and are almost equal to the size of direct effects on evapotranspiration. This finding has major implications for our mechanistic understanding of plant response to elevated CO 2 , forcing us to revisit the interpretation of experimental results as well as simulations of future productivity.

Abstract Increasing concentrations of atmospheric carbon dioxide are expected to affect carbon assimilation and evapotranspiration (ET), ultimately driving changes in plant growth, hydrology, and the global carbon balance. Direct leaf biochemical effects have been widely investigated, whereas indirect effects, although documented, elude explicit quantification in experiments. Here, we used a mechanistic model to investigate the relative contributions of direct (through carbon assimilation) and indirect (via soil moisture savings due to stomatal closure, and changes in leaf area index) effects of elevated CO 2 across a variety of ecosystems. We specifically determined which ecosystems and climatic conditions maximize the indirect effects of elevated CO 2 . The simulations suggest that the indirect effects of elevated CO 2 on net primary productivity are large and variable, ranging from less than 10% to more than 100% of the size of direct effects. For ET, indirect effects were, on average, 65% of the size of direct effects. Indirect effects tended to be considerably larger in water-limited ecosystems. As a consequence, the total CO 2 effect had a significant, inverse relationship with the wetness index and was directly related to vapor pressure deficit. These results have major implications for our understanding of the CO 2 response of ecosystems and for global projections of CO 2 fertilization, because, although direct effects are typically understood and easily reproducible in models, simulations of indirect effects are far more challenging and difficult to constrain. Our findings also provide an explanation for the discrepancies between experiments in the total CO 2 effect on net primary productivity.

The leaf-level response to elevated CO 2 (eCO 2 ) is well known: at current CO 2 levels, photosynthesis of C 3 plants is not saturated, whereas, for C 4 plants, it is close to saturation (e.g., refs. 1⇓⇓–4). If acclimation is limited, leaf-level carbon assimilation of C 3 plants will increase as the CO 2 concentration increases, as shown by, among others, observations in Free-Air CO 2 Enrichment (FACE) experiments (5⇓–7). Concurrently, stomatal conductance decreases consistently with eCO 2 in most species (8⇓⇓–11). Even though the leaf-level responses are well characterized and quantifiable, the ecosystem response to eCO 2 remains considerably more uncertain and difficult to predict (12⇓⇓⇓⇓–17). This discrepancy is not simply a consequence of the uncertainty in scaling up from leaf to canopy and ecosystem but derives from indirect effects and feedbacks that may lead to an amplification or dampening of the direct leaf-level response to eCO 2 .

Indirect effects may be related to (i) modifications of plant water status through changes in soil moisture within the root zone, which occur as a consequence of stomatal closure; (ii) changes in Leaf Area Index (LAI), root biomass and depth, and canopy structure; (iii) limitations due to soil nutrient scarcity or plant incapability to take up nutrients at a rate sufficient to support enhanced carbon assimilation; (iv) changes in ecosystem composition and biodiversity; and (v) higher-order interactions of the above indirect effects. Changes in soil moisture due to the reduction in stomatal conductance with eCO 2 (water saving effects) have been observed in a series of studies, most commonly in grasslands (18⇓⇓–21) but also in other ecosystems (22, 23). Water-saving effects have been hypothesized to stimulate vegetation productivity by a magnitude comparable to or larger than the direct eCO 2 effect (21, 24), but no study has quantitatively partitioned direct and indirect effects. For instance, the water-saving effects of eCO 2 can even lead to an increase in C 4 plant abundance over that of C 3 plants, despite the absence of substantial direct effects on C 4 plant growth (25, 26). Such strong indirect effects mediated by hydrology have led to the hypothesis that the response of vegetation to eCO 2 may be more an issue of water than of carbon (27). Despite evidence from individual studies for a considerable stimulation of plant productivity through eCO 2 -derived water savings, a metaanalysis of vegetation responses to eCO 2 showed no strong relationship between the effect size of CO 2 stimulation (elevated over ambient) and annual precipitation (28). Considering sites independently showed a larger CO 2 effect with decreasing precipitation, but the authors also hypothesized that the sign could be reversed approaching very dry conditions, based on observations from a desert study (29). Other indirect effects through LAI and nutrient limitations are interlinked with changes in soil moisture (e.g., refs. 30 and 31). A more favorable water status can increase LAI, but such an increase in LAI can also lead to more rapid depletion of soil moisture. Altered soil moisture may also influence soil microbial activity and soil organic matter turnover rates, ultimately modifying nutrient availability (32).

Although it is essential to improve our understanding of ecosystem response to climate change, a detailed quantification of indirect effects of eCO 2 and the comparison with direct physiological effects is fundamentally impossible in field experiments such as FACE (3, 5, 15, 33). Field experiments can only estimate the total response to a given external treatment or combination of treatments. The different components of the total response cannot be quantitatively separated in most cases, rendering the discussion of indirect effects speculative (e.g., refs. 20, 21, 24, 29, 34, and 35). A precise separation of the various effects is possible, however, using terrestrial biosphere models. Despite numerous limitations, these tools are capable of shedding light onto complex environmental issues with multiple interacting feedbacks (e.g., refs. 36 and 37). Here, we use the state-of-the-art ecohydrological model Tethys−Chloris (T&C) (38, 39) to disentangle direct from indirect CO 2 effects on productivity and evapotranspiration (ET). The model was used to simulate the response to an eCO 2 treatment in a series of ecosystems spanning a wide range of climates and biomes (Materials and Methods). The specific questions addressed here are as follows: (i) What are the relative contributions of direct and indirect effects to the total productivity and ET response to elevated CO 2 and do these depend on ecosystem water availability? (ii) In which biomes should we expect the strongest response? (iii) Does the total response to eCO 2 correlate with wetness? We hypothesize that indirect effects will be generally significant and potentially comparable to direct effects in water-limited grassland ecosystems, especially in sites containing C 4 species. We further expect a negative correlation between indirect effects on net primary productivity and wetness index for wet and mesic sites (28).

Although quantifying direct leaf-level effects is relatively easy, untangling the indirect effects at the ecosystem scale is far more difficult and requires a dedicated approach as presented here. Quantitative knowledge of the role of indirect effects is expected to shed light on the observed differences between sites in the response to eCO 2 and to improve our understanding of the interannual variability of the eCO 2 effects. A mechanistic explanation of indirect effects also suggests avenues for improvement of global terrestrial biosphere models, which have limitations in capturing ecosystem-level responses to eCO 2 (40).

Discussion This study aimed to determine the relative importance of direct and indirect effects of eCO 2 on ecosystem NPP and ET. Our results suggest that, in xeric environments, indirect effects can be comparable to or even larger than the direct, photosynthetic effect of eCO 2 on NPP. On average, indirect effects accounted for 28% of the total stimulation of NPP. The indirect/direct effect ratio ranged from less than 0.1 for tropical and moist sites to more than 1 for semiarid C 4 grasslands. The hypothesized decrease of effect size with extremely dry conditions (28) was not supported by our simulations, which represent integrated responses across multiple years. However, our results should be regarded as potential responses, in the absence of nutrient limitations. Suppression of eCO 2 effects on NPP by severe water deficit remains a possibility for explaining interannual variation in response within sites. Note that, for arid or semiarid sites characterized by herbaceous species and relatively fast biomass turnover rates, eCO 2 stimulation of NPP does not necessarily translate to an increase in standing biomass, even after several years (64, 65), but may be detected in an increase in soil organic carbon (66). This result is only partially captured in model simulations that still show a positive effect on biomass also in the most arid ecosystems, even though the effect is considerably smaller than for NPP. Limitations in nutrient uptake exacerbated by water stress can also dampen the biomass response of the most arid sites. Additionally, eCO 2 may stimulate rhizodeposition, potentially explaining the discrepancy between NPP and biomass responses. We found that changes in ET due to eCO 2 were smaller than what a pure, direct response of stomatal conductance would suggest, even at the ecosystem level (i.e., direct effect of −5 to −15%), because indirect effects tend to compensate partially or totally for the direct effect. Further, water “saved” via reduced stomatal conductance is likely to be consumed in water-limited systems, either immediately via increased LAI or by extension of the growing period if LAI is unaffected by eCO 2 (Fig. 2). Some of the effects might be due to changes in root biomass, which were included in the model; however, changes in rooting depth in response to eCO 2 were not considered, and it is therefore possible that indirect effects of eCO 2 may increase beyond those simulated here, if development of deeper roots were able to access water not otherwise available. The overall difference between the two CO 2 scenarios (375 ppm vs. 550 ppm) in terms of water fluxes (ET) was typically less than 8% and mostly constrained between −5% and 0. Changes in water use of this magnitude would rarely be observable, due to a combination of measurement uncertainty (e.g., ref. 67) and interannual variability (e.g., ref. 39). Over the large number of sites we simulated, the total change in NPP with the increase in CO 2 concentration was mostly in the order of 20 to 35%. These values are very similar or slightly larger than observations in FACE experiments when nutrient limitations do not play a role (28, 68). In fact, our results should be considered as the potential response of NPP to eCO 2 in the absence of sink limitations (e.g., ref. 58). The variation in the NPP response as a function of the wetness index is quite impressive, because these are numerical simulations from a mechanistic model rather than observations from real experiments. The large scatter in intermediate wetness conditions suggests that differences in phenology, temperature, short-term meteorological variability, biome, and soil type, all of which were accounted for in the simulations, play a significant role in the NPP response to eCO 2 . Contrary to the situation with ET, the sum of indirect effects tends to enhance the response of NPP to eCO 2 because it adds to the direct physiological response. This is especially evident in semiarid sites, which are responsive to eCO 2 even when C 4 species are predominant, as supported from observations (25). Our results demonstrate that mechanistic models of terrestrial ecosystems, despite known limitations (e.g., refs. 46, 58, 69, and 70), do provide substantial insights on ecosystem response to eCO 2 that are impossible to obtain with field experiments alone. Model limitations and structure may affect the magnitude of some of the estimates but are unlikely to change the prevailing patterns, with the important exception of nutrient limitation. Furthermore, T&C generated total responses to eCO 2 that closely matched observations. For instance, the average modeled eCO 2 effect size of NPP, ET, and EWUE is consistent for the Duke-FACE and for the first 7 y of the ORNL-FACE experiments. Regardless of inherent shortcomings of simulation models, ecosystems at the dry end of the climate spectrum, which experience repeated water stress, are expected to be the most responsive to eCO 2 in terms of productivity. When indirect LAI effects are removed, mimicking a lack of stimulation in LAI increase (Fig. 2), productivity in these sites responds more strongly to eCO 2 . Further, the significant positive relationships between VPD, a measure of atmospheric dryness, and total NPP response to eCO 2 (SI Appendix, Fig. S4B) supports the idea that the drier sites are where the most significant effects of eCO 2 on NPP should be expected. This finding agrees with modeling studies based on optimality principles (71, 72) and is supported by global patterns of positive response of semiarid ecosystems to CO 2 fertilization (73, 74), forcing the reevaluation of the role of semiarid ecosystems in the land carbon sink (75, 76). All this evidence corroborates our results and suggests that projections of eCO 2 effects at local and global scales are substantially affected by mechanisms and feedbacks contributing to indirect effects, which are inherently more challenging to model than the direct effect on carbon assimilation. Information on indirect effects derivable from conventional field experiments is necessarily limited. This issue demands both novel experiments specifically designed to target indirect effects and mechanistic solutions in models that do not strongly depend on empirical results. In this context, particular focus should be devoted to addressing the representation of water stress effects on the response of ecosystem productivity.

Materials and Methods Partitioning Direct and Indirect Effects of Elevated CO 2 . The contributions of the four identified effects (E1, E2, E3, and E4) were quantified by running a series of six simulations with the T&C model (38, 39, 49, 52, 77) (SI Appendix, Text S1). The first two simulations were used to compute the total eCO 2 effect (E1 + E2 + E3 + E4) and simply represent a simulation with CO 2 concentration prescribed at ambient level (375 ppm) and one with elevated CO 2 (550 ppm), where all of the identified effects cooccur as in reality. The effect magnitude was computed as (eCO 2 − aCO 2 )/aCO 2 . The other three simulations were then necessary to partition the four effects (because there were four unknowns in four equations). In these simulations, atmospheric CO 2 concentration was kept at 550 ppm, and either soil moisture or LAI or both were externally prescribed to be the same as obtained from the ambient or eCO 2 simulations, rather than being prognostic variables. We ran four additional simulations to have redundancy on the estimate and keep the simulation with the total eCO 2 effect as a counterproof. These were (i) a simulation with eCO 2 and prescribed ambient LAI and soil moisture, where the direct effect only remained (i.e., only E1); (ii) a simulation with eCO 2 but prescribed ambient LAI, where all of the indirect effects mediated by LAI were absent (i.e., E1 + E2 occurred); (iii) a simulation with eCO 2 and prescribed ambient soil moisture and eCO 2 –LAI, where all of the indirect effects related to soil moisture were eliminated (i.e., E1 + E3 occurred); (iv) a simulation with eCO 2 and prescribed eCO 2 soil moisture and ambient LAI, where only the E3 effect was eliminated (i.e., E1 + E2 + E4 occurred). Opportune combinations of the two basic simulations with the four additional simulations with prescribed soil moisture or LAI were able to provide a distinct estimate of the quantitative contributions of the four effects (E1, E2, E3, and E4) to the total (combined) effect. Note that prescribing either soil moisture or LAI externally rather than allowing its prognostic computation in the model violates, to some extent, the water and/or carbon budget in the specific simulation; however, this was the only way to separate the four effects. As a final check, the sum of the four effects estimated with the additional simulations corresponded almost perfectly to the total eCO 2 effect (E1 + E2 + E3 + E4), testifying to the correctness of the procedure (SI Appendix, Fig. S7). An example of the results obtained with the adopted methodology is illustrated by the time series of NPP simulated imposing ambient CO 2 concentration (aCO 2 ), eCO 2 concentration, and a case with eCO 2 but with LAI and soil moisture fixed to ambient values (SI Appendix, Fig. S8 and Text S2). In the article, only long-term averaged responses over the entire simulation period (SI Appendix, Table S1) are shown, which are the results of effects occurring from the hourly to the multiannual scale. Climate Forcing and Vegetation. We selected 44 locations corresponding to sites of flux towers, manipulation experiments, and experimental stations covering different climates and biomes across the globe (SI Appendix, Table S1). For each site, the six described simulations were used to partition the four effects at an hourly scale. Values averaged over the entire length of the simulation were then reported in the results. The length of meteorological time series depended on the length of the available, good-quality hourly data for each location and ranged from a minimum of 2 y to a maximum of 31.7 y with a median of 7.8 y (SI Appendix, Table S1). Eight sites were also characterized by a nonnegligible fraction of C 4 species. The broad range of climate and vegetation types allowed for robustness in the investigation of eCO 2 effect and how it is partitioned, minimizing the risk of idiosyncratic results related to parametrization of a given biome or climate in a single location. At the same time, running the model locally rather than globally allowed us to avoid generic plant functional type parameterizations and large-scale climate forcing that may lead to large biases in the ecosystem response at a given site (49, 60, 78, 79).

Acknowledgments The helpful comments of two anonymous reviewers on an earlier version of the manuscript are acknowledged. The authors thank the organizers of the INTERFACE workshop, Using Results from Global Change Experiments to Inform Land Model Development and Calibration, held in 2014 in Beijing, China, where the discussion leading to this article began. S.F. received a travel grant from the Institute for Applied Ecology New Zealand, Auckland University of Technology.