A meta-analysis of 1,119 manipulative experiments on terrestrial carbon-cycling responses to global change

Direct quantification of terrestrial biosphere responses to global change is crucial for projections of future climate change in Earth system models. Here, we synthesized ecosystem carbon-cycling data from 1,119 experiments performed over the past four decades concerning changes in temperature, precipitation, CO2 and nitrogen across major terrestrial vegetation types of the world. Most experiments manipulated single rather than multiple global change drivers in temperate ecosystems of the USA, Europe and China. The magnitudes of warming and elevated CO2 treatments were consistent with the ranges of future projections, whereas those of precipitation changes and nitrogen inputs often exceeded the projected ranges. Increases in global change drivers consistently accelerated, but decreased precipitation slowed down carbon-cycle processes. Nonlinear (including synergistic and antagonistic) effects among global change drivers were rare. Belowground carbon allocation responded negatively to increased precipitation and nitrogen addition and positively to decreased precipitation and elevated CO2. The sensitivities of carbon variables to multiple global change drivers depended on the background climate and ecosystem condition, suggesting that Earth system models should be evaluated using site-specific conditions for best uses of this large dataset. Together, this synthesis underscores an urgent need to explore the interactions among multiple global change drivers in underrepresented regions such as semi-arid ecosystems, forests in the tropics and subtropics, and Arctic tundra when forecasting future terrestrial carbon-climate feedback. A synthesis of global change experiments that manipulated temperature, precipitation, carbon dioxide or nitrogen identifies a need to consider site-specific factors and interactions in Earth system models.

O ver the past four decades, numerous experiments have explored how terrestrial ecosystems respond to climate warming, changing precipitation regimes, elevated CO 2 (eCO 2 ) and nitrogen (N) enrichment [1][2][3][4][5] . These extensive experimental results are scattered across many publications, and thus it is impractical to compare the outcome of experiments among different ecosystems 6,7 and to assess possible nonlinear effects between global change drivers 8,9 . Although great efforts have been devoted to comparing model outputs against experimental observations with an attempt to reduce uncertainties on future projections of coupled carbon-climate system (for example, refs. [10][11][12], the global models used in the upcoming IPCC (Intergovernmental Panel on Climate Change) AR6 (the Sixth Assessment Report) have not been tested against a large set of experiments.
Here, we synthesized 2,230 peer-reviewed studies with 1,119 experiments published during 1973-2016 to compare the magnitudes of treatments achieved by experimental manipulations with the changes in global change drivers projected by the end of this century, and to quantify the responses and sensitivities of key terrestrial carbon-cycle variables, including net primary productivity and biomass as well as their above-and belowground components, root-to-shoot ratios, litter mass, gross and net ecosystem productivity, and ecosystem and soil respiration, to warming, increased or decreased precipitation, eCO 2 , N deposition and their combinations. Finally, data of multifactor experiments which factorially performed both multiple single-factor treatments and their combinations were pooled to examine whether multiple global change drivers cause multiplicative, synergistic or antagonistic responses.  Table 1 A full list of affiliations appears at the end of the paper. and Supplementary Text 1 for the distribution of experiments across countries and publication list, respectively). The majority of singlefactor (73.2%) and multifactor (76.6%) experiments were conducted in temperate ecosystems of the Northern Hemisphere (~30-60° N), in the USA, Europe and China. The climate space sampled by all the experiments ranges from 110 mm to 2,301 mm in mean annual precipitation and from −7.6 °C to 22.3 °C in mean annual temperature ( Supplementary Fig. 1). This climate space domain covers 96% of the northern temperate lands, but only 50% of the Earth's land area. Regions with scarce experiments are semi-arid ecosystems, forests in the tropics and subtropics, and Arctic tundra, despite the sensitivity of these ecosystems to climate variability and climate change being highlighted by recent studies of modelling and analyses of Earth observation data [13][14][15] .

Overview of global change manipulative experiments
The magnitudes of temperature increase experimentally applied at more than 60% of the mid-latitude warming experiments are within the projected warming range in the same latitude band by the end of the twenty-first century (Fig. 1b). However, the experimental warming magnitudes at the high latitudes were generally lower than those projected in all the IPCC Representative Concentration Pathway scenarios (RCPs) 16 . More importantly, the majority of the experiments imposed greater changes in precipitation amount and N input than the global mean rates of changes projected under different RCPs (Fig. 1c,e) 17 . By contrast, eCO 2 levels were within the projected range of RCP2.6 (475 ppm) and RCP6.0 (800 ppm), but the highest levels from RCP8.5 have not been covered by any experiment (1,313 ppm; Fig. 1d).

Local climate mediating warming effects on ecosystem carbon cycling
Across all the terrestrial ecosystems, experimental warming enhanced total, aboveground and root biomass, as well as soil respiration (Fig. 2a). Although belowground net primary productivity slightly increased under warming (P = 0.066), net primary productivity and its aboveground component remained unaltered. The root-to-shoot ratios showed a neutral warming response (Fig. 2f). Hence, carbon allocation patterns appeared to be largely unaffected by warming. In addition, gross and net fluxes of CO 2 exchange between terrestrial ecosystems and the atmosphere (gross and net ecosystem productivity, and ecosystem respiration) did not show significant changes under warming, partly due to the use of multiple sites and the considerable response variations among the experiments. The observations suggest that the slight enhancement of belowground carbon input could be offset by increased soil carbon release (enhanced decomposition of soil organic matter, accelerated root turnover or stimulated root respiration) 18,19 .
Soil moisture plays a predominant role in regulating plant photosynthesis and carbon allocation 20 as well as ecosystem respiration 21 under warming. In water-limited ecosystems, warming usually decreases carbon-cycle variables by exacerbating water limitation 22 . To test such indirect warming effects, we conducted metaregressions between carbon-cycle variable changes under warming and wetness indices, which are proxies of the water available for plants 23 . Results demonstrated that the warming responses of net primary productivity ( Fig. 3a; R 2 = 0.20, P = 0.027) and its belowground component ( Fig. 3b; R 2 = 0.44, P = 0.001), aboveground biomass ( Fig. 3c; R 2 = 0.13, P = 0.004) and ecosystem ( Fig. 3d; R 2 = 0.24, P = 0.009) and soil respiration ( Fig. 3e; R 2 = 0.59, P < 0.001) shifted from negative to positive with increasing wetness indices, suggesting negative and positive warming effects under dry and wet climate conditions, respectively. In contrast, the warming effects on plant root-to-shoot ratios changed from positive to negative with increasing wetness (Fig. 3f; R 2 = 0.23, P = 0.003), indicating that when exposed to warming, plants growing in dry areas tend to prioritize carbon allocation to roots at the expense of allocation to shoots. Alternatively, decreased root turnover in order to adjust to warming-induced drought may also contribute to this response. Together, these findings indicate that water limitations emerge as a likely indirect response to higher evaporative demand, counterbalancing the direct warming effects on plant growth and productivity 21 .
Finally, we estimated temperature sensitivities (percentage change of carbon variable per 1 °C temperature increase) for all observed variables with the exception of net ecosystem productivity. Net ecosystem productivity is the transient response of altered gross ecosystem productivity followed by lagged and warming-dependent ecosystem respiration. Theoretically, if warming duration exceeds carbon residence time in an ecosystem, gross ecosystem productivity and ecosystem respiration should balance each other, defining a new equilibrium state where net ecosystem productivity approaches zero. Across all the terrestrial ecosystems, root biomass showed the highest temperature sensitivity with an enhancement of 19.8% °C −1 ( Table 1). Despite the huge variability among the sites, our analyses suggest that the temperature sensitivity of aboveground biomass becomes negative at sites with higher mean annual temperature ( Fig. 3g; R 2 = 0.10, P = 0.007). Nevertheless, the temperature sensitivities of the other nine carbon-cycle variables were independent of mean annual temperature ( Supplementary Fig. 2). This is partly supported by an ecosystem-level study which has revealed that temperature sensitivity of ecosystem respiration is converging globally but is not mediated by mean annual temperature 24 . These observations are inconsistent with the findings of a recent study on global drylands where temperature sensitivity of soil heterotrophic respiration declines with increasing mean annual temperature 25 .

Non-uniform influences of precipitation anomalies
Increased precipitation stimulated all ecosystem carbon-cycle variables except litter mass, whereas decreased precipitation generally reduced carbon-cycle variables (Fig. 2b,c) 5 . At the ecosystem scale, we tested the first hypothesis (H 1 ) that semi-arid and temperate grasslands are more sensitive to increased precipitation than forests and deserts 26 . The responses to increased precipitation were normalized to 33% above the mean annual precipitation of each experiment (the median absolute value of precipitation treatments; see Methods for details). In contrast to the H 1 , the responses of net primary productivity and its belowground component as well as plant biomass (total, aboveground and root biomass) under a 33% precipitation increase were comparable between forests and grasslands (Fig. 2g), although the enhancements of these variables were statistically significant in grasslands only. In addition, the significant stimulations of aboveground net primary productivity under a 33% precipitation increase did not differ among forests, grasslands and deserts (between-group heterogeneity (Q B ) = 2.30, P = 0.31). However, increased precipitation enhanced soil respiration more in grasslands than in forests (Q B = 7.90, P = 0.02), but comparably between grasslands and deserts, which partly supported our H 1 .
We also proposed a second hypothesis (H 2 ), that arid ecosystems would be more sensitive to increased precipitation and less sensitive to decreased precipitation than humid ecosystems 23 . To test H 2 , we conducted meta-regressions of carbon-cycle variable responses to normalized increased and decreased precipitation (33% above and below the mean annual precipitation) to examine whether the changes of carbon-cycle variables are related to wetness indices. We found that the positive responses of aboveground net primary productivity ( Fig. 3h; R 2 = 0.04, P = 0.013) and soil respiration ( Fig. 3i; R 2 = 0.30, P < 0.001) under a 33% precipitation increase declined with increasing wetness and the negative impacts of a 33% precipitation decrease on ecosystem respiration were exacerbated with increasing wetness (Fig. 3j; R 2 = 0.53, P = 0.010). The results suggest that the effects of increased precipitation on aboveground net primary productivity and soil respiration were lower, but the impacts of decreased precipitation on ecosystem respiration were higher in humid areas-a result supporting our H 2 . No relationship of precipitation-induced changes in any other carbon-cycle variables with wetness was detected ( Supplementary Figs. 3 and 4). The non-uniform patterns in the dependence of changing precipitation effects upon wetness indices, which are proxies for the current moisture conditions at given sites, among carbon-cycle variables pose a great challenge for robust projections of carbon-cycle responses to changing precipitation regimes in global models 27,28 .
We calculated the relative precipitation sensitivity as the percentage change of carbon-cycle variables per 10 mm precipitation increase or decrease. All carbon-cycle variables except litter mass   were more sensitive to increased than to decreased precipitation across all ecosystems ( Table 1). The results generalize previous observations on ANPP responses to inter-annual rainfall variation at the Long Term Ecological Research Network sites 26 to all the carbon-cycle processes analysed here and the larger range of precipitation changes applied in manipulative experiments. All these findings indicate that most terrestrial ecosystems have resistance to drought, which might arise from intrinsic plant physiological adaptation, enhanced water-use efficiency and/or shifts in plant community composition 5,26,[29][30][31] . Increased mortality and decreased vegetation cover in response to water stress under decreased precipitation could help to explain the greater sensitivity of litter mass to decreased rather than increased precipitation. These patterns indicate a temporary mismatch of litter mass and aboveground net primary productivity in response to reduced precipitation, highlighting the non-uniform influences of precipitation anomalies on different carbon-cycle variables and the large variability in the analysed studies 27 , precluding specific interpretations. The potential disconnection between production and decomposition processes also suggests avenues for future comparative research, which would reveal novel mechanisms for ecosystem carbon cycling. Furthermore, the sensitivities of gross ecosystem productivity and ecosystem respiration to changing precipitation were generally greater than those of plant productivity, biomass and litter mass.
The greater sensitivity of aboveground than belowground net primary productivity in response to increased precipitation suggests that plants invest less carbon in water acquisition structures when relieved from water stress (Table 1), leading to reduced root-to-shoot ratios (Fig. 2f). Our observations are supported by the finding from a previous review revealing greater carbon allocation to aboveground wood production than to belowground carbon fluxes in forest ecosystems in response to increasing resource (for example, water and nutrient) supply 32 . In contrast, the enhanced root-to-shoot ratios under decreased precipitation indicate that plants allocate more carbon to water-and nutrientacquiring structures to resist drought, and thus there is a higher sensitivity of above-(aboveground net primary productivity and biomass, and litter mass) than belowground carbon-cycle variables a-e, Results on the percentage changes (mean ± 95% confidence intervals (CIs)) in ecosystem carbon-cycle variables, including net primary productivity (NPP) and its above-and belowground components (ANPP and BNPP), total, aboveground and root biomass (TB, AGB and RB), litter mass (LM), gross and net ecosystem productivity (GEP and NEP), and ecosystem and soil respiration (ER and SR) induced by W (a), increased precipitation (IP; b), decreased precipitation (DP; c), eCO 2 (d) and eN (e). f-h, The changes in rootto-shoot ratios (ΔR/S) under individual global change drivers, and the changes in carbon-cycle variables under normalized IP (NIP, a 33% precipitation increase) and eCO 2 for different ecosystems, respectively. i, The changes in AGB, RB, soil pH and SR under different eN rates. Significant (P < 0.05) Q B (between-group heterogeneity) indicates that the response ratios differ among ecosystems (g,h) or eN rates (i). The sample sizes are given along the y axis. The vertical dashed lines are zero lines. See Fig. 1 for abbreviations.   (belowground net primary productivity, root biomass and soil respiration) under decreased precipitation (Table 1). However, a previous study on carbon allocation of 14 mature forest stands has found a decreased root-to-leaf biomass ratio under drought, largely due to accelerated root turnover in a drier climate 33 . Therefore, it should be noted that changing precipitation regimes may have differential effects on carbon cycling, depending on other local climate conditions (for example, existing moisture conditions and precipitation timing and intensity). For example, precipitation increase at a dry or mesic site could stimulate biological activity, as opposed to a moist site where responses may be neutral and even negative 6,27,34,35 .
Our observations that the sensitivities of net primary productivity ( Fig. 3k; R 2 = 0.19, P = 0.042) and its aboveground component ( Fig. 3l; R 2 = 0.29, P < 0.001), and soil respiration ( Fig. 3m; R 2 = 0.13, P = 0.005) to increased precipitation shifted from positive to negative with increasing wetness indices provided direct evidence supporting this conclusion. Nevertheless, caution should be taken when extrapolating our findings to generate regional and global effects of changing precipitation amount, because rainfall frequency and intensity may be potentially important confounding factors in this meta-analysis.

More carbon allocated to plant belowground tissues under eCO 2
This synthesis demonstrated that, despite large variations among ecosystems, eCO 2 significantly stimulated plant productivity, biomass, litter mass and ecosystem carbon fluxes, with the exception of net ecosystem productivity which appeared highly uncertain (Fig. 2d). The observations could have been primarily attributed to increased leaf photosynthesis and water-use efficiency associated with CO 2 fertilization effects and water-saving mechanisms under eCO 2 (for example, decreased stomatal conductance and consequent evapotranspiration) 36 . We expected (H 3 ) that positive eCO 2 effects would be greater in drier than humid ecosystems 37,38 because enhanced water-use efficiency under eCO 2 would benefit plants especially under water limitation. However, results of meta-regressions to explore the correlations of the eCO 2 effects with mean annual precipitation and wetness indices showed no clear trends of the 12 carbon-cycle variables along mean annual precipitation. Higher wetness did not affect most of the carbon-cycle variables either, but was associated with a larger eCO 2 effect on litter mass ( Fig. 3n; R 2 = 0.53, P = 0.001), contradicting our hypothesis H 3 . Nevertheless, ecosystem respiration responses to eCO 2 significantly declined with increasing wetness indices ( Fig. 3o; R 2 = 0.40, P = 0.048), partly supporting our hypothesis H 3 . These observations suggest that CO 2 fertilization effects may be mediated by not only wetness expressed with a simple average index, but also other factors such as precipitation seasonal distributions 39,40 and extreme climate events 41 .
Further analyses indicate that ecosystem type also explains a fraction of variations in the eCO 2 effects among sites. For example, grasslands with low net primary productivity and aboveground biomass responded more weakly to eCO 2 than forests (Fig. 2h). The stronger responses of forests may be attributable to the fact that the 'forests' label in our analyses included tree seedlings with exponential growth over multiple growing seasons, which can greatly exaggerate the responses of older trees in a linear growth phase 42 . In addition, different mycorrhizal types associated with temperate forests, grasslands and tundra may help to explain variations among ecosystems in the response to eCO 2 43 . Unfortunately, mycorrhizal information was not reported in most studies.
Based on statistically significant changes in carbon-cycle variables under eCO 2 (Fig. 2d), we calculated the eCO 2 sensitivity as percentage change per 100 ppm CO 2 increase 44 . Across all eCO 2 experiments, gross ecosystem productivity and net primary productivity showed strongly positive CO 2 sensitivity (Table 1). Plants under eCO 2 allocated more carbon to roots for uptake of water and N for sustaining growth, at least in young forest/seedling experiments where an increase in plant root-to-shoot ratios was observed across ecosystems (Fig. 2f)-an observation supported by findings of a previous meta-analysis 45 . Nevertheless, another experimental study has demonstrated that enhanced root carbon allocation in a high CO 2 atmosphere generally occurs in a field (free-air CO 2 enrichment) but not in a controlled (growth chamber) environment 46 . In addition, we found that the positive effects of eCO 2 on root-to-shoot ratios declined with increasing depths of root sampling ( Fig. 3p; R 2 = 0.21, P = 0.035)-a result more consistent with the proliferation of shallow fine roots for harvesting nutrients rather than of deep roots for water uptake 47 . In contrast, a previous synthesis on forest ecosystems has revealed that eCO 2 stimulates root proliferation in deep soil, likely because the limited resource of topsoil cannot meet the increased demand for production enhancements under CO 2 enrichment 48 . Given that the depths of root sampling are usually greater in forests than grasslands, these observations suggest that root-to-shoot ratios may respond to eCO 2 more strongly in grasslands than forests, as shown by our meta-analysis results ( Fig. 2h; Q B = 5.70, P = 0.02). Similar to the observed increases in gross ecosystem productivity, ecosystem and soil respiration were also enhanced under eCO 2 , indicating that eCO 2 accelerates ecosystem carbon turnover through belowground processes (for example, microbial decomposition and soil CO 2 efflux) [49][50][51] . Root exudation and possibly root litter and organic matter decomposition might be stimulated under eCO 2 , thus a considerable fraction of increased net primary productivity is in fast-turnover pools that do not lead to carbon sequestration but are released back to the atmosphere. The coupling between net primary productivity and heterotrophic respiration under eCO 2 likely depends on the experimental duration. Under long-term exposure to eCO 2 , the two processes are expected to be closely coupled and net ecosystem productivity should approach zero after a new equilibrium is reached 52 , although none of the experiments were long enough to verify this assumption.

More carbon allocated to plant aboveground than belowground tissues under N addition
Nitrogen addition significantly stimulated plant productivity, biomass, litter mass and ecosystem carbon fluxes, but not belowground net primary productivity or net ecosystem productivity (Fig. 2e). In contrast, N addition reduced soil respiration across all the experiments assessed. Presumably, N availability limits net primary productivity in many ecosystems 53 . In contrast to the effects of eCO 2 , more carbon appeared to be allocated to above-than belowground tissues under N addition. Decreased root-to-shoot ratios in response to N addition across ecosystems suggest that the sensitivities of aboveground carbon-cycle variables to N addition are higher than those of belowground carbon-cycle variables ( Fig. 2f and Table 1). The results are consistent with the general observations that plants invest less carbon in constructing nutrient-acquiring structures (that is, roots) when more nutrients are available 32,54-56 . Our findings of substantial variability of plant carbon allocation in responses to different global change drivers are critical in evaluating Earth system models 57,58 .
The large number of N addition experiments compiled in this study (429) provides a unique opportunity to further analyse the effects of N addition on aboveground biomass, root biomass and soil respiration. Aboveground biomass increased with increasing levels of N addition rates (≤5, 5-10 and >10 g N m −2 yr −1 ; Fig. 2i; Q B = 7.20, P = 0.03). The stimulation of aboveground biomass was greater at N addition rates >10 g N m −2 yr −1 (32.1%) than ≤5 g N m −2 yr −1 (13.6%). However, there was no difference in root biomass responses among the three N addition levels, likely due to the lower sensitivity of below-than aboveground biomass and productivity in response to N addition (Table 1). By contrast, soil respiration did not respond to the low or medium levels of N addition but decreased by 10.2% at the high N addition rates of >10 g N m −2 yr −1 , which is partly inconsistent with the findings of a previous meta-analysis that both low (<5) and high (>5) N m −2 yr −1 additions suppressed soil CO 2 efflux in forest ecosystems 54 . The inclusion of grasslands and tundra in this study makes the conclusions more convincing. High N addition rates may reduce soil heterotrophic respiration due to acidificationinduced reduction in microbial activity, and thus suppress soil respiration 54,59 . Indeed, the decrease in soil pH was exacerbated with higher rates of N addition (Fig. 2i). Clearly, N addition rate is critical in estimating ecosystem carbon-cycle responses to atmospheric N deposition. Yet, given that manipulative experiments generally imposed greater amounts of N than the current global mean rates of atmospheric N deposition (Fig. 1e), the responses of ecosystem carbon cycling in these experiments are likely stronger than those caused by the forecasted levels of atmospheric N deposition.

Nonlinear effects of two-driver pairs were rare
Under current and future conditions, ecosystems will be exposed to concurrently multiple global change drivers (Fig. 1b-e). Nonlinear effects of co-occurring drivers pose challenges for understanding and projecting terrestrial carbon-climate feedbacks 8 . Two previous syntheses demonstrated that multiplicative influences were rare and ecosystem responses declined with increasing numbers of manipulated factors 60,61 , suggesting that antagonistic mechanisms are widespread when global change drivers are altered simultaneously. In contrast, other meta-analyses found that multiplicative, rather than synergistic or antagonistic, effects on woody plant biomass, ecosystem carbon storage and soil CO 2 efflux were more common [62][63][64] . These contrasting conclusions leave model predictions of such interactions in question 60 , especially given that three-and four-factor studies are too scarce to draw general conclusions. We analysed data from 80 twofactor experiments to examine whether the influence of one driver on carbon-cycle variables was modified by a second driver.
Our results showed that nonlinear effects (that is, synergistic or antagonistic) of two-driver pairs-including warming and increased precipitation, warming and decreased precipitation, warming and eCO 2 , warming and N addition, increased precipitation and N addition, decreased precipitation and eCO 2 , as well as eCO 2 and N addition-on all carbon-cycle variables were rare (Fig. 4a-g). These observations differ from the findings of Leuzinger et al. 60 and Dieleman et al. 61 -a discrepancy that probably resulted from the different analytical methods. For example, Dieleman et al. 61 assessed interactions between warming and eCO 2 on plant biomass and soil process variables (for example, soil respiration, soil carbon content and net N mineralization) by comparing the slopes of the linear correlations of the combined warming and eCO 2 influences with the sum of individual effects with 1:1 line. Leuzinger et al. 60 pooled plant biomass and soil process variables together to explore the differences in the response magnitudes under one, two and three global change drivers. When compared with the two previous reviews using Hedges' d method that concluded predominantly multiplicative effects on ecosystem carbon storage and soil CO 2 emissions 63,64 , our results, obtained with different analytical methods and database, confirm that the nonlinear effects of multiple global change drivers on the carbon-cycle processes analysed are relatively infrequent. Another study using the same methodological approach as ours found that warming and eCO 2 multiplicatively affected aboveground biomass of woody plants 62 , further supporting our conclusion. Moreover, a most recent study reporting a fourfactor (temperature, precipitation, CO 2 and N) experiment revealed predominantly multiplicative influences among the four global change drivers on ecosystem carbon cycling in a temperate grassland 51 , providing other lines of supporting evidence. All the above findings indicate that changes in carbon-cycle variables under multiple global change drivers can be under-or overestimated if more antagonistic 60,61,65 or synergistic 44,66 effects are assumed.
Nevertheless, two experimental studies on N cycling and net primary productivity by manipulating three 65 and four 67 factors demonstrated that the majority of interactions among multiple factors are antagonistic rather than multiplicative or synergistic, thereby dampening the net effects. From a global perspective, more direct evidence on whether multiplicative effects among global change drivers hold for experimental studies manipulating three or four factors is crucial for convincing projections of future global change terrestrial carbon-cycling feedbacks.

The future of global change manipulative experiments
During the past two decades, meta-analyses have reported ecosystem carbon-cycle responses to warming 2,4,5,68,69 , changing precipitation amount 5,23 , eCO 2 (refs. 1,45,[70][71][72][73] ) and N addition 3,54,73,74 , but relatively few have examined the combined effects (but see refs. [61][62][63][64]. The data synthesis presented in this study, encompassing all global change manipulative experiments published to date, revealed global-scale general patterns of the responses and sensitivity of all the terrestrial carbon-cycle processes to multiple global change drivers and their interactions. Further discussion on the discrepancy between our findings and those from previous meta-analyses is provided in Supplementary Text 2. In addition, sensitivity analyses were developed to inform model validation, thus caution should be taken when using our findings to predict the impact magnitudes of global change drivers on carbon-cycle variables at given experiments ( Supplementary Fig. 5). While global change research has made tremendous progress in developing manipulative experiments to better simulate climate change scenarios 7,75 , our synthesis highlights three key issues that need to be addressed in future experimental studies to better facilitate and assess model projection on terrestrial carbon cycling and its feedback to climate change.
First, although the existing manipulative experiments are critical to understand and disentangle mechanisms (for example, ref. 76 ), identify thresholds and tipping points of ecosystem responses to global change drivers 77 , and provide necessary information (input data and process understanding) for the modelling community 11 , the majority of the experiments have done a poor job of simulating the more complex climate change scenarios. For example, the stochastic nature of precipitation events causes substantial uncertainty in predicting future scenarios and great challenges in designing experiments to capture the full breadth of possibilities inherent in these scenarios. This is made even more difficult if experiments also attempt to capture the impacts of extreme events and exceeded thresholds 6,7,[78][79][80] .
Second, the majority of existing studies focused on single-factor effects and thus precluded the analysis of nonlinear influences on carbon-cycle feedbacks among multiple global change drivers 8 . Although the combined treatment effects of two-driver pairs on carbon cycling were mostly multiplicative (Fig. 4), three-and fourfactor experiments are still urgently needed to generate multidimensional (that is, multiple factors) response surfaces for carbon-cycle variables aiming to accurately quantify terrestrial carbon-global change feedbacks 81 . More specifically, changing plant community composition over a long-term scale may also lead to the occurrence of nonlinear effects among multiple factors 82,83 , indicating the importance of plant competition in regulating ecosystem responses to concurrent global change drivers 30 .
Finally, the influences of global change on carbon cycling are ecosystem-and climate-specific, suggesting that more experimental data from several hotspot regions are necessary to improve the global perspective. Semi-arid regions, forest ecosystems in the subtropics and tropics, and Arctic tundra are underrepresented in experiments, although they store large quantities of organic carbon and thus have a disproportionate impact on global carboncycling and climate change-carbon feedbacks (Fig. 1a) 6,13-15 . In sum, forecasting terrestrial carbon-cycle responses to ongoing global change requires future manipulative experiments aiming at knowledge useful to inform mechanistic models (for example, gradient treatments) 84 , to include multiple driving factors, and to focus on underrepresented arid/semi-arid, subtropical/tropical and tundra ecosystems. First, we identified publications from all 310,177 records through reading each title. In this step, we excluded 202,257 records whose titles showed that these studies were reviews/meta-analyses or conducted in non-terrestrial ecosystems (for example, oceans; Supplementary Fig. 6). The remaining 107,920 publications that might be relevant to our topic were included. Second, we further screened publications that may have reported global change effects on terrestrial ecosystems by reading the abstracts of all the 107,920 publications collected in the first step. Similar to the first step, we excluded 89,179 publications whose abstracts indicated that these studies were reviews/meta-analyses or conducted in non-terrestrial ecosystems. The PDF versions of the remaining 18,741 publications that might be relevant to the synthesis were downloaded. Third, we selected publications by reading the methods of the 18,741 studies to identify which of them met the following inclusion criteria:

Methods
• Studies reported results of manipulative experiments conducted in the outdoor environment and having both control and treatment groups. Each group had at least three plots as replicates. In addition, the area of each plot was larger than 1 m 2 . The treatment groups were manipulated artificially to simulate global change drivers. For example, grassland plots of treatment groups were warmed using infrared radiators 51 . • Experiments were conducted in terrestrial ecosystems. However, cropland and lab incubation studies were not included in this synthesis. • Studies examined effects of simulated global change drivers on ecological processes, including carbon (for example, net primary productivity and its above-and belowground components, total, aboveground and root biomass, root-to-shoot ratios, litter mass, gross and net ecosystem productivity, and ecosystem and soil respiration, and so on), nitrogen (for example, soil N content and N mineralization rates) and water-cycle variables (for example, evapotranspiration and water-use efficiency) as well as plant (for example, plant cover, species richness and community temporal stability) and microbial parameters (for example, microbial biomass and microbial community composition and structure).
In the third step, we excluded 17,451 publications which did not meet the above inclusion criteria. After screening this huge number of publications, 1,290 studies met the defined criteria. This list of 1,290 publications was subsequently cross-checked with references cited in the previous review/meta-analysis articles that studied global change influences on ecological processes, as well as publications used in this list to ensure that our publication survey was as inclusive as possible, and 756 additional publications were found. Finally, by searching the websites of ecology laboratories and experiment networks-for example, http:// www2.nau.edu/luo-lab/, https://mnspruce.ornl.gov/, http://www.nutnet.org/, https://www.ars.usda.gov/plains-area/fort-collins-co/center-for-agriculturalresources-research/rangeland-resources-systems-research/docs/range/prairieheating-and-co2-enrichment-phace-experiment/research/, http://www.biocon. umn.edu/, http://www.utas.edu.au/profiles/staff/plant-science/Mark-Hovenden and http://web.ics.purdue.edu/~jsdukes/Dukes.html-and checking the references of the papers downloaded from these websites, another 184 studies were collected. In total, 2,230 publications were included in the synthesis.
Information including publication year, simulated global change drivers, experimental duration, site locations and altitudes, site climate conditions (mean annual temperature (MAT) and precipitation (MAP)), vegetation types and the manipulation magnitudes of global change drivers was extracted from each of the 2,230 publications collected in our literature search to demonstrate the state-ofthe-art and analyse the future challenges of ecosystem manipulative experiments in global change research. Based on this information, the 2,230 publications were classified into two groups: (1) publications reporting single-factor experiments having a single control and at least one single-factor treatment (for example, eCO 2 or eN) and (2) publications reporting multifactor experiments using factorial designs having a single control, multiple single-factor treatments (for example, eCO 2 and eN), and combinations of factors (for example, eCO 2 + eN). In addition, we identified how many single-and multifactor experiments were documented in the 2,230 publications. If a publication reported several experiments conducted at different sites or in different plant communities, it would be considered as several independent experiments. If a single-factor experiment manipulated several different global change drivers with a same control at the same site, it would be broken up into several independent control-treatment pairs and thus considered as several independent experiments. Values of MAT and MAP were downloaded from Climate Model Intercomparison Project phase 5 (CMIP5; https://esgf-node.llnl. gov/projects/cmip5/) based on the site locations if individual studies did not report these values for their sites. Ecosystems subjected to manipulative experiments were classified into five typical types: forests (including mature forests and tree seedlings), grasslands (including grasslands, meadows, short-and tall-grass prairies, temperate/semi-arid steppes, shrublands, savannas, pastures and oldfields), tundra, wetlands (including peatlands, bogs, marshes and fens) and deserts. In this study, 86.5% and 59.6% of single-and multifactor forest experiments were conducted in mature forests, respectively.
We plotted the distributions of MAT and MAP against latitude, and calculated the De Martonne wetness index 85 of each experimental site using the data of MAT and MAP as In addition, the spatial distributions of single-and multifactor experiments were mapped in ArcMAP v.10.0. Moreover, to compare the differences between the manipulation magnitudes of global change drivers extracted from the 2,230 publications and realistic climate change scenarios, we calculated the changes in MAT, MAP and atmospheric CO 2 concentrations at each experimental site as the differences between the periods of 2081-2100 and 1986-2005 using the data reported in the IPCC RCP2.6, 4.5, 6.0 and 8.5. Maximum average N deposition over the world projected in RCP2.6 (1.55 g N m −2 yr −1 ), RCP4.5 (1.02 g N m −2 yr −1 ) and RCP8.5 (1.32 g N m −2 yr −1 ) by the end of the twenty-first century used in the synthesis were reported in Lamarque et al. 17 .

Effects of individual global change drivers on ecosystem carbon-cycle variables: meta-analysis.
This study focused on the responses of ecosystem carbon cycling to global change drivers. Data of mean values ( X I ), standard deviations (s.d.) or standard errors (s.e.m.) and sample sizes (number of plot replication) of ecosystem carbon-cycle variables (net primary productivity and its above-and belowground components, total, aboveground and root biomass, root-to-shoot ratios, litter mass, gross and net ecosystem productivity, and ecosystem and soil respiration) in control and treatment groups were extracted from each publication when possible. SigmaScan 5.0 (SPSS, Inc.) was used to digitize the figures and extract the numerical values when a publication presented the data graphically. In our meta-analysis, experimental data were excluded if the experiments were conducted over less than one year/growing season, except for tundra studies for which most of measurements were performed during the growing season (July-August). In total, 4,010 observations from single-factor experiments and 1,013 observations from multifactor experiments were collected (Supplementary Table 2). In addition, experimental data of the last year/growing season were used in our analyses when multiple measurements were taken at different years/growing seasons in a study to meet the statistical assumption of independence among samples in the metaanalysis 86 , based on the methods reported in several previous meta-analyses 1,64,69 . After checking for independence, that is, only keeping data of the last year/ growing season in the meta-analyses when multiple measurements were taken at different years/growing seasons in a study, 2,104 observations from single-factor experiments and 392 observations from multifactor experiments in total were used in our meta-analysis (Supplementary Table 3).
The natural log-transformed response ratio (RR) was used to estimate the effect size of global change treatments on ecosystem carbon-cycle variables as: with a variance of: where X I T , s:d: T I , and n T I represent the mean, standard deviation and sample size of treatment (T) group, respectively, X C I , s:d: C I and n C I are the mean, standard deviation and sample size of control (C) group, respectively. Based on the methods developed by a previous study 87 , the ln RR was calculated separately for each control-treatment pair and treated as independent data when data extracted from multifactor experiments which had multiple single-factor treatments and a single control.
The metafor 88 package in R software 89 was used to calculate the weighted response ratio (ln RR ++ ) and bias-corrected 95% bootstrap-confidence interval (CI) using inverse-variance weighted regressions and random-effects models. Moreover, the publication bias was also evaluated by funnel plots (that is, the scatter plots of the treatment effect sizes against their standard errors) 90 . If the studies were distributed symmetrically in a 'funnel' shape around a mean effect size, it would indicate the absence of publication bias. Therefore, we further assessed the potential asymmetry of the funnel plots using Egger's regression 90 . Results of funnel plots and Egger's regressions under the W, increased precipitation (IP), decreased precipitation (DP), eCO 2 and eN treatments are provided in Supplementary Figs. 7-11, respectively. We found that publication bias in our meta-analyses was rare. However, there was also some evidence for publication bias: for example, the studies examining the eCO 2 effects on GEP (Supplementary Fig. S10; z = 3.84, P < 0.01). Finally, percentage change (%) of treatment effects was evaluated as: If the 95% CIs of treatment effects did not cover zero, the responses of ecosystem carbon-cycle variables to global change drivers were statistically significant.
In addition, responses of ecosystem carbon-cycle variables with different ecosystem types were considered to be significantly different if their 95% CI values did not overlap. Furthermore, we standardized the effects of IP/DP across the studies that manipulated different levels of precipitation treatments to compare the differences in precipitation influences among ecosystems and wetness indices based on the methods developed by a previous study 23 . All of the levels of precipitation manipulations were converted into a percentage of MAP. Overall, the manipulations in the IP treatments ranged from 1.1% to 151.0% of the MAP, with a median of 35.1%, and the manipulations in the DP treatments ranged from 5.8% to 100% of the MAP, with a median of 32.3%. We pooled all the studies to determine the absolute values of their manipulation levels, and the median was 32.9% of the MAP for the whole database. The meta-regressions conducted in the metafor package indicated that the RRs of net primary productivity and its aboveand belowground components, total, aboveground and root biomass, root-to-shoot ratios, litter mass, gross and net ecosystem productivity, and ecosystem and soil respiration showed linear responses to precipitation manipulation levels across all of the precipitation studies ( Supplementary Fig. 12). Therefore, we performed the following linear transformation to normalize the RRs of all the carbon-cycle variables under the different treatment/manipulation levels to 32.9% of the MAP: where X NT I is the normalized value under 32.9% above or below the MAP, and P is the precipitation manipulation levels (% of the MAP), with positive values representing the IP treatments and negative values indicating the DP treatments. The normalized lnRR ++ and 95% CIs for changing precipitation studies were also calculated using the metafor package in R software. Results of funnel plots and Egger's regressions under the normalized IP and DP treatments are provided in Supplementary Figs. 13 and 14, respectively. In addition, meta-regressions were conducted to examine the correlations of the normalized lnRR ++ with wetness indices using the inverse-variance weighted regressions.
To explore which factors could determine the variations of the eCO 2 effects (RR) on ecosystem carbon-cycle variables among sites, meta-regressions were also used to examine the correlations of the eCO 2 effects with eCO 2 levels (ppm), MAP and wetness indices. In addition, we also investigated whether the depths of root biomass sampling affect the responses of root-to-shoot ratios to individual global change drivers using meta-regressions.

Interaction evaluations of multiple global change drivers on ecosystem carboncycle variables.
Given that there was not enough data available from three-and four-factor experiments to analyse interactions among three and four global change drivers, we only estimated interactive effects of two-driver pairs using 80 two-factor experiments (362 observations in total; Supplementary Table 4). If a two-factor experiment included four treatments (the control, treatment A, treatment B, and treatment A plus B), based on the techniques developed by two previous studies 62,87 , the interactions were calculated as: with a variance (the multiplicative property of variances) 87  and n C I represent the mean, standard deviation and sample size of the group of treatment A plus B, treatment B, treatment A and control, respectively. In the calculation of the interaction term, experimental data of the last year/growing season were used in the metafor package when multiple measurements were taken at different years/growing seasons. Results of funnel plots and Egger's regressions under the interactions between W and IP (W × IP), W and DP (W × DP), W and eCO 2 (W × eCO 2 ), W and eN (W × eN), IP and eN (IP × eN), DP and eCO 2 (DP × eCO 2 ) and eCO 2 and eN (eCO 2 × eN) are provided in Supplementary Figs. 15-20, respectively. However, caution should be taken when using our datasets, because we used fewer than three studies per group to test for publication bias for some two-driver interactions.
Based on a previous study that defined the calculation of interactions 62 , the interactions between two factors were classified into three types, including multiplication, positive interactions and negative interactions. If the 95% CI of the interaction term overlapped with zero, the interactive effect was considered to be multiplicative. If not, there were positive (the interaction term >0) or negative interactions (the interaction term <0).
Sensitivity of ecosystem carbon-cycle variables in response to global change drivers. Further calculations were performed to analyse sensitivity of ecosystem carbon-cycle variables to global change drivers as a standardized response magnitude per unit change in the global change driver. First, the main effects (E) of global change drivers on ecosystem carbon-cycle variables at each experiment were calculated using formula (4), which gives a percentage change. Second, the sensitivity (S) of ecosystem carbon-cycle variables to W (S W , percentage change per 1 °C −1 increase), IP (S IP , percentage change per 10 mm increase), DP (S DP , percentage change per 10 mm decrease), eCO 2 (SeCO 2 I , percentage change per 100 ppm increase) and eN (S eN , percentage change per 1 g N m −2 yr −1 addition) in each experiment was calculated as: The units of warming (1 °C increase), eCO 2 (100 ppm increase) and N addition (1 g N m −2 yr −1 addition) sensitivities were selected based on two previous studies 5, 44 . Moreover, given that effective rainfall is generally more than 5 mm, we calculated precipitation sensitivities as 10 mm increase or decrease in precipitation. Third, the weighted mean sensitivities of ecosystem carbon-cycle variables to global change drivers at global/ecosystem scales were calculated using the same method as the meta-analysis (bias-corrected 95% bootstrap-confidence interval using inversevariance weighted regressions and random-effects models).
To investigate whether the local climate conditions influence the sensitivity of ecosystem carbon-cycle variables in response to global change drivers, simple linear regressions with 95% CIs were conducted to examine the dependences of the warming sensitivity of ecosystem carbon-cycle variables on local MAT, and the relationships of the IP and DP sensitivity of ecosystem carbon-cycle variables on local MAP and wetness indices.
We also plotted the observed against their predicted effects of global change drivers on carbon-cycle variables at each experiment, to examine the predictive ability of our sensitivity database. The observed effects of global change drivers on carbon-cycle variables at an experiment are the influences reported directly in publications collected in this study. The predicted effects were calculated as the product of the magnitudes of experimental manipulations (for example, temperature increases in warming experiments) at the experiment and the sensitivity calculated from the other experiments (Supplementary Figure 5).
Reporting Summary. Further information on the research design is available in the Nature Research Reporting Summary linked to this article.

Statistics
For all statistical analyses, confirm that the following items are present in the figure legend, table legend, main text, or Methods section.

n/a Confirmed
The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement A statement on whether measurements were taken from distinct samples or whether the same sample was measured repeatedly The statistical test(s) used AND whether they are one-or two-sided Only common tests should be described solely by name; describe more complex techniques in the Methods section.
A description of all covariates tested A description of any assumptions or corrections, such as tests of normality and adjustment for multiple comparisons A full description of the statistical parameters including central tendency (e.g. means) or other basic estimates (e.g. regression coefficient) AND variation (e.g. standard deviation) or associated estimates of uncertainty (e.g. confidence intervals) For null hypothesis testing, the test statistic (e.g. F, t, r) with confidence intervals, effect sizes, degrees of freedom and P value noted

Data analysis
All meta-analyses and meta-regressions were performed using the package 'metafor' in R software (version 3.5.1). Details were reported in statistical analysis section of the Methods.
For manuscripts utilizing custom algorithms or software that are central to the research but not yet described in published literature, software must be made available to editors/reviewers. We strongly encourage code deposition in a community repository (e.g. GitHub). See the Nature Research guidelines for submitting code & software for further information.

Data
Policy information about availability of data All manuscripts must include a data availability statement. This statement should provide the following information, where applicable: -Accession codes, unique identifiers, or web links for publicly available datasets -A list of figures that have associated raw data -A description of any restrictions on data availability The data supporting the results of this study were archived in Figshare  All studies must disclose on these points even when the disclosure is negative.

Study description
Synthesis study.

Research sample
In total, we extracted 4010 observations from 850 single factor experiments and 1013 observations from 269 multi-factor experiments for ecosystem carbon-cycle variables including net primary productivity and its above-and below-ground components, total, aboveground, and root biomass, root-to-shoot ratios, litter mass, gross and net ecosystem productivity, and ecosystem and soil respiration (Supplementary Table 2), to perform meta-analysis calculations. We also extracted 362 observations from 80 two-factor experiments for ecosystem carbon-cycle variables (Supplementary Table 4), to estimate the interaction terms among multiple global change drivers.

Sampling strategy
The Timing and spatial scale Together, the studies published between 1973 to 2016, covered seven ecosystems including forests, grasslands, wetlands, tundra, and deserts. The mean annual precipitation (MAP) ranged from 52 mm to 4239 mm, and the mean annual temperature (MAT) ranged from less than -19°C to greater than 27°C.

Data exclusions
We excluded studies that were conducted in aquatic and agriculture ecosystems. In addition, we collected data from field studies, and excluded those from laboratory incubation studies.