One of the major uncertainties for carbon-climate feedback predictions is an inadequate understanding of the mechanisms governing variations in ecosystem productivity response to warming. Temperature and water availability are regarded as the primary controls over the direction and magnitude of warming effects, but some unexplained results signal that our understanding is incomplete. Using two complementary meta-analyses, we present evidence that soil nitrogen (N) availability drives the warming effects on ecosystem productivity more strongly than thermal and hydrological factors over a broad geographical scale. First, by synthesizing temperature manipulation experiments, meta-regression model analysis showed that the warming effect on productivity is mainly driven by its effect on soil N availability. Sites with higher warming-induced increase in N availability were characterized by stronger productivity enhancement and vice versa, suggesting that N is a limiting factor across sites. Second, a synthesis of full-factorial warming×N addition experiments demonstrated that N addition significantly weakened the positive warming effect, because the additional N induced by warming may not further benefit plant growth when N limitation is relieved, providing experimental evidence that N regulates the warming effect. Further, we demonstrated that warming effects on soil N availability were modulated by changes in dissolved organic N and soil microbes. Overall, our findings enrich a new mechanistic understanding of the varying magnitudes of observed productivity response to warming, and the N scaling of warming effects may help constrain climate projections.

Literature search and screening protocol

We searched peer-reviewed publications using the Web of Science (http:// apps.webofknowledge.com) and the China National Knowledge Infrastructure Database (http://www.cnki.net/). Additional studies were also identified by screening the cited references in the obtained literature. The following criteria were adopted to identify eligible studies: (1) The temperature manipulation experiments should be performed in the field condition, and greenhouse and incubation experiments were excluded. (2) In Database 1, the candidate studies should report at least one variable representing ecosystem productivity (i.e., GEP, ANPP or NEP) and soil N availability (i.e., SIN) simultaneously in control and warming treatments under the same condition. The BNPP and total NPP are also important parameters representing vegetation C gains, but they were not included here due to the limited data in the collected studies. (3) In Database 2, the selected studies should report at least one productivity variable in control, warming, N addition and warming×N addition treatments simultaneously. Finally, our selection process yielded 105 warming experiments from 83 papers in the first database, and 53 full-factorial warming×N addition experiments from 38 papers in the second database. Among them, 31 papers in Database 2 are overlapped with those in Database 1. Raw data were either obtained from tables directly or extracted from figures using Origin Pro 2021 (OriginLab, Northampton, MA, USA). In the individual experiments of the two databases, GEP and NEP were measured using a static chamber with infrared gas exchange analyzer technique. For ANPP estimation, various methods (harvest of aboveground live plants or estimated by percent cover or normalized difference of vegetation index (NDVI) in grasslands and tundra, and estimated by tree size increase for forests) were accepted. The SIN was generally determined by extracting NH₄⁺-N and NO₃^–-N from soil solutions, resin bags or pore water. Five out of the 105 experiments used two methods for SIN determination (i.e., soil solutions & resin bags). In such a case, we calculated the mean SIN of the two measures. No significant methodological influence was observed on the warming effects on SIN (P=0.63; Appendix S1: Fig. S1). Sometimes the productivity and N availability measured in the same temporal and spatial scale were reported in different papers. For example, in a warming experiment conducted in a moist acidic tundra, the SIN (2012-2013) and ANPP (2009-2013) were shown in Salmon et al. (2016), while GEP and NEP (2009-2015) were reported in Mauritz et al. (2017). In such a case, we extracted the three productivity variables and SIN measured in the most recent year from the intersection of different measurement periods (that is, 2013) from the two papers, and combined them into one data item.

In addition to ecosystem productivity and soil N availability, we also retrieved other factors that potentially influence the warming effects, including ecosystem type, mycorrhizal type, warming method and experimental duration, mean annual temperature (MAT), mean annual precipitation (MAP), canopy temperature (CT), soil temperature (ST), soil moisture (SM), soil pH in both databases. Ecosystem types included forests, grasslands and tundra. Mycorrhizal type of the dominant species for the individual experiments was classified as ectomycorrhizal (ECM) and arbuscular mycorrhizal (AM) fungi according to Wang & Qiu (2006) and Maherali et al. (2016). Warming methods included active warming (i.e., infrared radiator and heating cable) and passive warming (i.e., open-top chamber, greenhouse, curtain and snow fence). Experimental duration ranged from 1 to 21 years. The MAT and MAP varied from -19.0 to 19.1°C and 60.0 to 1750 mm, respectively. Some studies did not report climatic characteristics; in such a case, we extracted the multi-year (1970-2000) average MAT and/or MAP data from WorldClim database (http://www.worldclim.org). The N fertilizer type was also extracted for the second database, including NH₄NO₃ and other fertilizer types (i.e., urea, Ca(NO₃)₂ and (NH₄)₂SO₄). In addition, to investigate the mechanisms of the warming-induced changes in soil N availability, we collected the dissolved organic N (DON) (reflecting substrate quantity) and microbial biomass C (MBC) (reflecting soil microbial biomass) which potentially explain the direction and magnitude of SIN response to warming (Kou et al., 2018; Li et al., 2019; Risch et al., 2019).

Meta-analysis

We adopted the response ratio (RR) to quantify the response of target variables to experimental warming. The RR is defined as the natural logarithm of the ratio between treatment (X_t) and control groups (X_c) (Hedges et al., 1999). The individual RRs are usually weighted by the inverse of the sampling variance (e.g., the standard deviation) or the number of replicates (Adams et al., 1997). Since several individual studies did not report the sampling variance, and weighting based on sampling variance may render the mean effect size depending largely on a small number of observations with extremely low variance, we used the number of replicates from the collected studies as the weighting factor (Adams et al., 1997):

w = n_cn_t / (n_c+n_t)

where w is the weighting factor, n_c and n_t are the number of replicates for the control and warming treatments, respectively. Case studies with more replicates were considered to have greater contribution to the grand mean effect size (RR₊₊). We calculated the overall effects in a weighted, mixed-effects model using the “rma.mv” function in “metafor” package. Several studies examined the interactions of warming with other manipulations (e.g., grazing, water addition or drought, and elevated CO₂, etc.). When an effect with additional resources was given, we took this effect as the control, and the combined treatment (e.g., warming plus grazing) as the warming treatment (Janssens et al., 2010; Wang et al., 2014). Also, some studies used different plant species in the same experiment. In such cases, the data were treated as non-independent estimates of the effect (Nakagawa & Santos, 2012; Noble et al., 2017). The non-independency of observations from the same study site was accounted for by including observations nested in ‘study’ as a random factor (Chen et al., 2018; Terrer et al., 2021). The effects of warming were considered significant if the 95% bootstrapped confidence interval (CI) did not overlap with zero. Percentage changes of variables defined here as (e^RR++ – 1) × 100% were calculated to ease interpretation. The warming effects on CT, ST and SM were also evaluated. Because about a half of the collected studies only reported the changes in CT (ΔCT), ST (ΔST) and SM (ΔSM) rather than the respective values for warming and control treatments, we then adopted these actual changes but not the RRs in our analyses.

Statistical analysis

Data were processed in the following three steps. First, to determine the relative importance of experimental and environmental factors on the variations in the RRs of productivity in the first database, we analyzed all possible factors in a mixed-effects meta-regression model using the “glmulti” package with observation nested in study as a random effect (Calcagno & de Mazancourt, 2010). The importance of each predictor was expressed as the sum of Akaike weights for models that included this factor, which can be considered as the overall support for each variable across all models. A cutoff of 0.8 was set to differentiate between essential and nonessential predictors (Chen et al., 2018). Because the model selection process requires complete moderator values, we filled the missing values of ΔCT, ΔST, ΔSM and pH-RR before running the model selection using the multivariate imputation by the chained equations with random forests method (Hou et al., 2021) in the ‘miceranger’ package (Wilson, 2020). Results showed that the interpolation did not alter the relationships between dependent and independent variables (Appendix S1: Fig. S2).

Second, to compare the warming effects on ecosystem productivity between no N and N addition treatments in the second database, we conducted a meta-regression analysis using the “rma.mv” function in “metafor” package with observations nested in ‘study’ as a random effect (Terrer et al., 2021). In the meta-regressions, the RR of ANPP was categorized into non-N and N-added treatments and treatment was included as a categorical moderator to test the between-group heterogeneity (Gurevitch et al., 2018). Further, we explored the relationships between ANPP-RR and SIN-RR for no N and N addition treatments separately using meta-regressions that included SIN-RR as a continuous predictor (Gurevitch et al., 2018). The relationships of GEP-RR and NEP-RR with SIN-RR were not investigated because of fewer data points reported on these variables. 

Third, to investigate the regulations of the warming-induced changes in soil N availability, the associations of SIN-RR with various biotic and abiotic factors were explored. We applied meta-regressions to fit the above relationships using the “rma.mv” function from the “metafor” package. The explanatory variables were included as fixed effects and observation nested in study was treated as a random effect. 

All statistical analyses were conducted in software package R 3.6.3 (R Development Core Team, 2020).