Interactive effects of global change factors on terrestrial net primary productivity are treatment length and intensity dependent

Individual effects of co‐occurring global change factors on net primary productivity (NPP) have been widely studied; however, their interactive effects remain highly debated. Here, we conducted a global meta‐analysis based on 919 multifactor observations from 120 published studies to examine the interactive effects on NPP of global change factors including elevated [CO2], warming, nitrogen addition, irrigation, drought and changes in species diversity. On average, of the factors studied, six pairs of factors had additive and two pairs had synergistic interactions. Importantly, some of those interaction types changed over time and with treatment intensity. The synergistic interaction between elevated [CO2] and nitrogen addition became additive at high nitrogen addition rates, whereas the synergistic interaction between irrigation and warming diminished at higher temperatures. Over time, the additive effect between elevated [CO2] and increased species richness switched to synergistic. Other global change factor pairs—including elevated [CO2] and warming, nitrogen addition and increased richness, irrigation and N addition, as well as drought and increased richness—remained additive regardless of their treatment intensity or experimental duration. Interaction types of those global change factor pairs did not vary with ecosystem types assessed in our study. Synthesis. Our results suggest that the assumptions of static effects through time or ignoring treatment intensity effects will provide inaccurate predictions of the interactive effects of global change factors on terrestrial NPP. Understanding the context‐dependent nature of interactive effects is crucial for validating Earth system models and predicting future NPP responses to co‐occurring global change drivers.

the combined effect of multiple factors is equal to or not differing from the sum of individual effects, while synergistic and antagonistic interactions occur when the combined effect is greater or weaker than the sum of the individual effects respectively Yue et al., 2017). Two previous syntheses demonstrated that ecosystem responses declined with an increasing number of manipulated factors (Dieleman et al., 2012;Leuzinger et al., 2011), indicating an antagonistic mechanism when global change factors are altered simultaneously. In contrast, recent meta-analyses suggest that additive, rather than synergistic or antagonistic, effects of two global change factors on terrestrial carbon/nutrient cycling are more common (Song et al., 2019;Yuan & Chen, 2015;Yue et al., 2017;Zhou et al., 2016). These divergent findings may be attributable to the selection of analytical methods, but more likely due to the potential change of interactive effects with experimental duration, for example, through changes in community composition (Komatsu et al., 2019), and treatment intensity (i.e. the level of the treatment applied).
Both the 'Liebig's law of the minimum' and the multiple resource co-limitation theory provide a conceptual basis for predicting the interactive effects of global change factors on NPP (Bloom, Chapin, & Mooney, 1985;Rastetter & Shaver, 1992). Usually, interactive effects occur when one global change factor's effect causes the limitation of a critical resource (i.e. carbon, nutrients and water), and other global change factors either alleviate or increase the severity of the limitation. However, the occurrence of interactive effects is highly dependent on the experimental duration and the level of the treatment applied. For example, the N status of ecosystems could be a critical determinant of how elevated [CO 2 ] and warming interact (Norby & Luo, 2004); their synergistic interactions may dissipate over time due to progressive N limitation Morgan et al., 2011).
Interactions between species diversity and climate changes are more likely to occur in the long term (e.g. several years after experiment initiation) because the ability of species diversity to alleviate the negative impact of climate change on ecosystem NPP through an insurance effect (e.g. a community with higher diversity of hydraulic traits has greater resilience to drought), or to enhance the positive impact through niche complementarity (e.g. plant species richness enhances the positive response of biomass production to elevated CO 2 and N addition), tends to increase over time (Anderegg et al., 2018;Hisano, Searle, & Chen, 2018;Niklaus, Leadley, Schmid, & Korner, 2001;Reich et al., 2001;Tilman, Reich, & Knops, 2006). Moreover, although irrigation may ameliorate the desiccating effect of moderate warming through increasing soil water availability , it is unlikely to reverse the physiological damage and mortality caused by extreme warming (Morgan et al., 2011). However, no consensus exists yet for how interaction types of global change factors on NPP may change with experimental duration and treatment intensity.
Over the last decade, an increasing number of experiments have been conducted to investigate the responses of terrestrial ecosystems to multiple global change factors with full factorial designs (Dieleman et al., 2012;O'Brien, Reynolds, Ong, & Hector, 2017;Zhou et al., 2016). As a result, substantial empirical data on productivity from field and laboratory experiments manipulating two or more global change factors are now available across levels of organization ranging from individual plants to whole ecosystems (Jagadish et al., 2014;Meyer-Grunefeldt, Friedrich, Klotz, Von Oheimb, & Hardtle, 2015). Some studies also included simultaneous manipulation of plant species richness in factorial combination with climate change factors (O'Brien et al., 2017;Smith, Lukac, Bambrick, Miglietta, & Godbold, 2013). Global change impacts on NPP involve a complex mixture of driving variables and will act on large spatial and long time-scales. Enough multifactor studies across a wide range of ecosystem types, experimental durations and treatment intensities now exist to allow for the predictive understanding of when and where various interaction types are expected to occur and whether these interaction types are context dependent.
Here, we conducted a meta-analysis based on 919 multifactor (≥2) observations from 120 published studies, which covered a broad range of terrestrial ecosystem types ( Figure 1)  (plant growth OR production OR productivity OR biomass production)'. Only primary studies that satisfied the following criteria were included in this meta-analysis: (a) at least a 2 × 2 full factorial design was used to examine the effects of global change factors including elevated [CO 2 ], warming, N deposition, irrigation and/or drought, and change in species diversity; (b) at least one component of production (e.g. above-ground, below-ground or total) was examined under all treatments and the control at the same temporal and spatial scales; (c) the plots used for all treatments were of the same ecosystem type and environmental condition as the control; and (d) the mean and sample size of the selected variables were available or could be calculated from the data provided in the publications. All original data were extracted from the text, tables, figures and appendices of the publications. When data were graphically presented, SigmaScanPro version 5 (Systat Software Inc.) was used to obtain numeric data. Our final dataset included 919 multifactor (≥2) observations from 120 published studies (Appendix S1), which covered a wide range of terrestrial ecosystem types ( Figure 1). Of all studies, half of the experiments were conducted in the field (64), with the others conducted in greenhouses (35), growth chambers (19) or mesocosms (2). The number of studies conducted at the community level (80) was greater than those conducted at the plant species level (40). For those conducted at the species level, NPP was measured on individual plant species. Across all studies, field-and laboratory-based measurements of individual and multiple global change effects on NPP did not differ significantly (Tables S1 and S2). Similarly, both individual and interactive effects of multiple global change factors on NPP did not differ between plant and community levels except for the individual effect of N addition (Tables S1 and S2).
Furthermore, data regarding experimental duration (years), treatment intensity (i.e. increased or decreased quantities applied in the treatment as compared with the control), ecosystem type (forest, grassland, cropland, tundra, wetland), production component (aboveground, below-ground or total) and growth form (woody or herbaceous), were also recorded from detailed descriptions in the Section 2 of papers included in the analysis. Diverse methods were used in applying the drought treatment, such as reduced precipitation (13), reduced field capacity (2), reduced water content (7) and the withholding water until plants exhibited symptoms of wilting (3). While data related to experimental duration, ecosystem type, production component, growth form, and treatment intensity of elevated [CO 2 ] and species richness were available in all studies, data on intensity were not available for some treatments (e.g. warming (1), N addition (4), irrigation (9) and drought (15)) and were thus not included the analysis.

| Individual and combined effects
We used the natural log response ratio (lnRR) to quantify the effects of individual and combined global change factors on plant growth and ecosystem NPP (Hedges, Gurevitch, & Curtis, 1999). The natural lnRR was calculated as: where X t and X c are mean values of productivity under treatment (i.e. with factor A, and with factor B, and with both factors, A + B) and primary control respectively. Within each study, lnRR was calculated separately for each treatment level and experimental duration.
Effect size estimates and subsequent inferences in the metaanalyses may be dependent on how individual observations are weighed (Ma & Chen, 2016). In our dataset, sampling variance was not reported in 6 of the 120 studies, and more importantly, weightings based on sampling variance had a very wide range (Table S7) which could assign extreme importance to a few individual observations, and consequently, the average lnRR would be determined predominantly by a small number of studies (Ma & Chen, 2016).
Similar to previous studies (Pittelkow et al., 2015), we instead employed the number of replications for weighting.
where W r is the weight associated with each lnRR observation, n c and n t are the number of replications in the control and treatment respectively.
We tested whether the average lnRR differed from zero, and whether the lnRR was affected by experimental duration and treatment intensity using the following models: where β 0-3 , π study and ε are model coefficient, the random effect of 'study' and sampling error respectively. For the combined effect, treatment intensities of both factors were included in the model. The random effect explicitly accounts for interdependence among observations within each 'study'. Akaike information criterion (AIC) was used to compare the models with and without the interaction term of Dur × Int (Equations 3 vs. 4). Since Equation (3) had a lower AIC value and thus was used in our analysis to avoid overfitting (Johnson & Omland, 2004).
We conducted the analysis using maximum likelihood estimation with the lme4 package in r (Bates, Maechler, Bolker, & Walker, 2015). Weight equals to W r in Equation (2). When continuous predictors, that is, Dur and Int in Equation (3) We also examined whether lnRR changed with ecosystem type, production component, growth form, response level and experimental facility by adding the individual term to Equation (3). We further examined the AIC values of the models with and without the interaction term of Dur × another factor (i.e. ecosystem type, production component, growth form, response level and experimental facility) and Int × another factor (same as above). The models without the interaction term had lower or similar AIC values (Table S3) and considered as the most parsimonious models (Johnson & Omland, 2004).
To investigate whether the experimental facility influenced the sensitivity of productivity in response to global change drivers, we conducted a sensitivity analysis (Table S5)

| Interactive effects
The interaction between two factors could be inferred by comparing additive individual effects with their combined effect (Folt, Chen, Moore, & Burnaford, 1999). In our study, however, direct comparisons of the estimated average effect sizes for individual and multiple global change factors were hampered by their unequal sample sizes.
Additionally, because most multiple factor models and individual studies based on ANOVA assume additivity, our objective was to specifically test how often multifactor effects deviate from additive. Therefore, we favoured the use of the approach in Crain et al. (2008) and performed a pair-specific meta-analysis (taking into account only those studies that simultaneously tested individual and combined effects for each global change factor pair) to identify interaction types based on the Hedge's d ( Figure S1). We note that Hedge's d and lnRR provided consistent results for the individual effects of global change factors in our study.
The individual Hedge's d of a factor was calculated using Equation (5).
where X t and X c are mean values of productivity under a treatment and the control, respectively, and s and J(m) are the pooled standard deviation and correction term for small sample bias (Hedges & Olkin, 1985), which were estimated by Equations (6 and 7) respectively.
where n t , n c , s t and s c are the sample sizes (n t , n c ), and standard deviations (s t , s c ) in the treatment and control groups. m is the degree of freedom (m = n t + n c − 2). The interaction (d I ) of factor A and B was calculated by Equation (8).
where X C , X A , X B and X AB , are means of a variable in the control and treatment groups of A, B and their combination (A + B) respectively.
The standard deviation (s) and degree of freedom (m) were estimated by Equations (9 and 10), respectively, for the interactive effect.
Within global change pairs where a sufficient number of studies were available, a meta-analysis was performed. We tested whether average Hedge's d differed from zero, and whether the interactive effects (Hedge's d I ) were affected by experimental duration (Dur, year) and treatment intensity (Int) using the following models: where β 0-3 , π study and ε are model coefficient, the random effect of 'study' and sampling error respectively. For the interactive effect d I , treatment intensity of both factors was included in the model. The random effect explicitly accounts for interdependence among observations within each 'study'. Based on the AIC values, Equation (11) was selected for our analysis. We conducted the analysis using maximum likelihood estimation with the lme4 package in r (Bates et al., 2015).
Weight equals W r in Equation (2). When continuous predictors, that is, Dur and Int in Equation (11), were scaled, β 0 was the overall mean Hedge's d at the mean Dur and mean Int (Cohen et al., 2014).
According to the methods of Crain et al. (2008), we used the individual and interaction effect size and 95% CI of Hedge's d to classify each global change factor pair as additive, synergistic or (11) d, d I = 0 + 1 ⋅ Dur + 2 ⋅ Int + study + , (12) d, d I = 0 + 1 ⋅ Dur + 2 ⋅ Int + 3 ⋅ Dur × Int + study + , F I G U R E 2 Individual and combined effects of different global change factors on net primary productivity. Values are mean ± 95% confidence intervals of the percentage effects between treatments and control. Individual and combined effects are represented by red and green respectively. The number of observations is noted beside each factor without parentheses, and the number of studies is listed in parentheses [Colour figure can be viewed at wileyonlinelibrary.com] antagonistic ( Figure S1). Additive effects were those with their 95% CI of interaction effect size overlapping zero. For factor pairs whose individual effects were both positive, interaction effect sizes greater than zero were synergistic, and those less than zero were antagonistic. In cases where the individual effects were both negative, the interaction was interpreted in the opposite manner (interaction effect sizes >0 is antagonistic and <0 is synergistic).
We also examined whether Hedge's d I changed with ecosystem type, production component, growth form and response level by adding the individual term to Equation (11). We further examined the AIC values of the model with and without the interaction term of Dur × another factor and Int × another factor. The models without the interaction term had lower or similar AIC values (Table S4) and considered as the most parsimonious models (Johnson & Omland, 2004). While low sample size within global change factor pairs limited the power of our pair-wise meta-analyses, we did not find that pairs with more replicate studies always had significant interactions while those with few studies were additive; therefore, we chose a minimum of five replicate studies as a reasonable cut-off despite the limited availability of data.

| RE SULTS
Across all studies, the experimental duration ranged from 0.3 to 18 years. Elevated [CO 2 ] treatment ranged from 110 to 450 ppm above ambient. Warming treatment caused an increase in air temperature (49 studies) ranging from 0.9 to 20°C, and in soil temperature (8 studies) ranging from 0.6 to 10°C. The N application rate ranged from 1 to 90 g m −2 year −1 , increased precipitation treatments ranged from 23.7 to 967 mm above ambient, reduced precipitation ranged from 9 to 735 mm below ambient, and species richness ranged from 1 to 32.

| Individual, combined and interactive effects
Under the mean experimental duration and treatment intensity,

F I G U R E 3
Results of pair-specific meta-analysis using weighted Hedge's d for individual and interactive effects of given global change factor pairs. Values are mean ± 95% confidence intervals of effect size. Colours indicate interaction types: red = synergistic; green = additive. The numbers outside and inside the parentheses represent the number of observations and studies respectively [Colour figure can be viewed at wileyonlinelibrary.com]

| Regulation of experimental duration, treatment intensity, ecosystem type, production component and growth form on global change effects on NPP
With increasing experimental duration, the positive effects of elevated [CO 2 ] and the negative effects of drought on NPP on average decreased marginally (p = 0.06 and 0.08 respectively), while the positive effects of increased species richness on average increased marginally (p = 0.10) (Figure 4a). The interaction effect size (d I ) between elevated [CO 2 ] and increased species richness increased significantly with experimental duration (Figure 5a). Their interaction types switched from additive to synergistic with the increase in experimental duration (Figure 6b).
The positive effects of warming on NPP decreased (p < 0.001) with the increasing magnitude of temperature treatment, while those of increased species richness became significantly stronger with increasing species richness levels ( Figure 4b). The interaction effect size between elevated [CO 2 ] and N addition increased significantly with rising [CO 2 ] but decreased with increasing N addition (Figure 5b), and the observed synergistic interactions became additive at high N application rates (Figure 6a). The interaction effect size between N addition and warming increased significantly with the N application rates (Figure 5b) but their interaction type remained additive (Figure 6c). With increasing temperatures, the interaction effect size between irrigation and warming decreased significantly (Figure 5b), where their synergistic interactions became additive at higher temperatures ( Figure 6d).
Ecosystem type did not affect the response of NPP to any individual or combined global change factors (Table S1).
The individual effects of warming, N addition, increased species richness, as well as the combined effects of elevated [CO 2 ]-N addition, elevated [CO 2 ]-increased species richness, N addition-increased species richness on NPP, were greater for the above-ground component than below-ground (Table S1; Figure S2a). Both individual effects of elevated [CO 2 ] and warming, as well as their combined effects were greater for woody F I G U R E 4 Effects of (a) experimental duration (years) and (b) treatment intensity on the natural log response ratio (lnRR) of productivity to individual and combined global change factors. Values are mean ± 95% confidence intervals of the coefficient. Individual and multiple factors are represented by red and green respectively. Grey circles represent the effects of treatment intensity of the second factor of each global change factor pair. The numbers outside and inside the parentheses represent the number of observations and studies respectively [Colour figure can be viewed at wileyonlinelibrary.com] species over herbaceous species (Table S1; Figure S2b). The interaction effect size between N addition and increased species richness differed significantly between production components (Table S2; Figure S3), while the interaction effect sizes of all eight global change factor pairs were not affected by ecosystem type nor growth form (Table S2).

| D ISCUSS I ON
We show that global change factors can interactively affect plant growth and ecosystem NPP, rendering the results from assessing single global change factor effects less useful for predicting the real effect of global change that is inherently multifactorial in its nature.
More importantly, our study highlights that those interaction types are unstable and could change over time and with treatment intensity, assumptions of static effects through time or ignoring treatment intensity effects will provide inaccurate predictions of the interactive effects of global change factors on terrestrial NPP.
The responses of productivity and plant growth to individual global change factors in this meta-analysis are qualitatively consistent with those of manipulative experiments, model simulations and meta-analysis syntheses, where the effect of experimental duration and treatment intensity was not accounted for in their analyses (Dieleman et al., 2012;Luo et al., 2008;Rustad et al., 2001;Shaw et al., 2002;Song et al., 2019;Wu, Dijkstra, Koch, Penuelas, & Hungate, 2011;Xia & Wan, 2008;Yue et al., 2017). Among the global change factors, N addition and increased species richness resulted in the largest enhancement in NPP, indicating that N is often a limiting factor for ecosystem NPP and the average effect of diversity loss is comparable to and has even exceeded in magnitude to that of numerous other global change factors that have already mobilized major international concerns and remediation efforts (Duffy et al., 2017;Hooper et al., 2012). The combined effect of elevated [CO 2 ] and N addition, and drought and increased richness were greater than the sum of their individual effects (Figure 2), indicating synergistic interactions. However, direct comparisons of the estimated average effect sizes for individual and combined global change factors were hampered by their unequal sample sizes (Figure 2). We favoured the results from pair-specific meta-analysis for the identification of interaction types (Crain et al., 2008).
Our pair-specific meta-analysis reveals that, under mean experimental duration and treatment intensity, the effects of many global change factor pairs were additive, lending support to cumulative models that assume additivity of global change effects (Hutchison & Henry, 2010;Song et al., 2019;Yue et al., 2017;Zavaleta, Shaw, Chiariello, Mooney, & Field, 2003), but also highlights that certain global change factor pairs have overall interactive effects worth incorporating in future modelling efforts. Synergistic interactions on NPP occurred primarily via a combination of elevated [CO 2 ] and N addition, as well as irrigation and warming (Figure 3). The availability of mineral nutrients, particularly N, can greatly modify the responses of plant growth to elevated [CO 2 ] (Idso & Idso, 1994;Kimball, 1983). For example, N addition can enhance an ecosystem's potential to take advantage of elevated [CO 2 ] by increasing Rubisco activity and production Reich, Hungate, & Luo, 2006). The synergism between irrigation and warming observed in our study is consistent with a modelling analysis showing a general positive interactive effect between warming and the doubling in precipitation . The most likely explanation is that irrigation strongly increased soil water availability and enhanced the response of water-limited processes to warming . Nevertheless, previous meta-analyses show additive effects of elevated [CO 2 ] and N addition, irrigation and warming on plant C pools (Song et al., 2019;Yue et al., 2017). We further demonstrate that factors such as experimental duration and treatment intensity could explain the discrepancy between our results and those from Yue et al. (2017) and Song et al. (2019).
Our analyses showed that interaction types of global change factors on NPP vary with experimental duration. Both the individual effect of richness and its interaction with elevated [CO 2 ] on NPP increased with increasing experimental duration (Figures 4a and 5a).
Their interaction types switched from additive to synergistic with the increase in experimental duration (Figure 6b), indicating that a community with diverse species could benefit more from the elevated [CO 2 ] in the long term. Over time, the response of NPP to elevated [CO 2 ] could be constrained by nutrient availability  as indicated by the decreased individual effect of elevated [CO 2 ] (Figure 4a), however, increased niche complementarity over time might alleviate this limitation, and thus result in a synergistic interaction (Reich et al., 2012). Global change factors can alter community structure (e.g. species composition and species diversity) and consequently, biomass production (Hisano et al., 2018).
Therefore, changing plant community composition over long timescales may also lead to the alteration of interactive effects among multiple factors (Komatsu et al., 2019). For example, N addition can reduce plant species richness (Harpole et al., 2016), which may consequently limit the response of NPP to elevated [CO 2 ] (Niklaus et al., 2001). Unfortunately, for studies conducted at the community level, information regarding species composition changes over time and species-specific responses were not reported in most studies. Thus, more studies need to be conducted to explore the potential role of changes in community composition in long-term productivity-global change feedbacks.
Our results further indicate that interaction types of global change factors on NPP can change with treatment intensity.
Although individual effects of elevated [CO 2 ] and N addition on NPP did not change with treatment intensity (Figure 4b), their interaction effect size increased with rising [CO 2 ] but decreased with increasing N addition (Figure 5b), and the observed synergistic interactions became additive at high N application rates (Figure 6a).
One possible explanation is that the greater positive response of above-ground NPP to N addition ( Figure S2a) may enhance standlevel evapotranspiration and reduce soil water reserves, which in turn dampens the benefit of both elevated [CO 2 ] and N addition (Litton, Raich, & Ryan, 2007). The limitation of other nutrients (e.g. phosphorus) might be another explanation for the additive effect of elevated [CO 2 ] and N addition at high N availability (Li, Niu, & Yu, 2016). Additionally, we found that both the individual effect of warming on NPP and its interaction with irrigation decreased with increasing temperatures (Figures 4b and 5b), where their synergistic interactions disappeared at higher temperatures ( Figure 6d). While mild warming may increase NPP, extreme warming could reverse this trend and lead to substantial declines exacerbated by reductions in water availability (D'Orangeville et al., 2018). Although irrigation may ameliorate the desiccating effect of moderate warming through increasing soil water availability , it is unlikely to reverse the physiological damage and mortality caused by extreme warming (Morgan et al., 2011). Together, these results suggest that interactive effects of multiple global change factors on NPP change with both experimental duration and treatment intensity most likely through depleting or alleviating limiting resources over time, altering species composition, or by causing irreversible damage to plants.
Both individual and interactive effects of global change factors on NPP may also change with other moderator variables such as ecosystem type, production component and growth form.
Different ecosystems are generally limited by the availability of different resources which may interact with climate change to regulate biomass production (Huxman et al., 2004;Reich et al., 2018).
For example, warming may negatively affect NPP in ecosystems in dry climates but positively affect NPP in those in wet climates (Song et al., 2019). Elevated [CO 2 ] and warming are more likely to exhibit synergistic interactions on biomass production in waterand nutrient-limited ecosystems since the CO 2 -induced increase in water use efficiency and the warming-induced increase in nutrient mineralization will allow the full effects of another factor (Dieleman et al., 2012). In our study, however, we did not observe differences among ecosystem types assessed for both individual global change factors and their interactions (Tables S1 and S2).
This may be attributable to the fact that most of the multifactor field experiments included in our study were conducted in grasslands, while other ecosystem types such as forests, croplands, tundra and wetlands were under-represented (Table S6), which limits our ability to determine the potential differences among ecosystem types. However, the enhanced response of aboveground NPP to individual global change factors (i.e. warming, N addition, and increased species richness), as well as some of their combinations ( Figure S2a) indicate an increased biomass allocation to above-ground components (Song et al., 2019). This may have been caused by increased soil nutrient availability directly through N addition, or indirectly through increased nutrient mineralization induced by warming or complementary effects caused by increased species richness (Litton et al., 2007;Rustad et al., 2001;Tilman, Lehman, & Thomson, 1997). Above-ground biomass and productivity are commonly used to assess the responses of plant growth and ecosystem NPP to global change; however, our results suggest that below-ground productivity should also be measured to better quantify ecosystem-level responses, given the potential shift of production allocation in response to global change (Song et al., 2019). Furthermore, both the individual and combined effects of elevated [CO 2 ] and warming were higher for the NPP of woody species than for herbaceous species ( Figure S2b). This finding is consistent with a meta-analysis showing that the enhancement of photosynthesis by elevated [CO 2 ] was greater for woody species than for herbaceous species (Wang, Heckathorn, Wang, & Philpott, 2012). The low productivity of herbaceous plants may have limited their capacity to take up mineralized N induced by warming, or the herbaceous plants become more shaded out by the overstorey plant species as they take advantage of the greater resource availability, thereby reducing the strength of this positive feedback (Wu, Dijkstra, Koch, & Hungate, 2012). Thus, woody plant productivity will benefit more than the herbaceous vegetation under elevated [CO 2 ] and warming. The interactive effect size of N addition and increased species richness differed between production components (Table S2); however, their interaction types were still additive ( Figure S3).

| CON CLUS IONS
The last several decades have witnessed an enormous research effort Importantly, the majority of existing multifactor studies focused on two-factor pairs in certain ecosystem types (i.e. grasslands). Thus, more multifactors studies, manipulating three or more global change factors, in under-represented ecosystems such as forest, tundra, cropland and wetland ecosystems, over long time-scales are needed to improve our ability to quantify terrestrial productivity-global change feedbacks.

ACK N OWLED G EM ENTS
We thank Eric B. Searle for editorial comments and Si Chen for editing the graphs. We thank the Natural Science and Engineering Research Council of Canada (NSERC) for partially supporting this research.

DATA AVA I L A B I L I T Y S TAT E M E N T
Data deposited in the Dryad Digital Repository: https://doi. org/10.5061/dryad.8cz8w 9gkz .