Materials and Methods

Field Experiment

The field experiment was conducted in Duolun County (42°02’N, 116°70’E, 1324 m a.s.l.), a semiarid temperate steppe in Inner Mongolia, northern China. The multiple-year mean annual precipitation and mean annual temperature are approximately 385 mm and 2.1°C, respectively. The topography is characterized by low foothills at elevations of 1150-1800 m. Soil is Calcis-orthic Aridisol (the US Soil Taxonomy classification), with 62.75 ± 0.04% sand, 20.30 ± 0.01% silt and 16.95 ± 0.01% clay (Wu et al. 2010). The dominant plant species include Stipa krylovii , Artimesia frigida , Potentilla acaulis ,Cleistogenes squarrosa , Allium bidentatum andAgropyron cristattum .
To assess the nitrogen (N) and phosphorus (P) impacts on CH4 flux, a block design experiment with different combinations of N and P additions was established in early 2013 and run to 2016. A complete random design, with three block replicates, was adopted to address the high spatial heterogeneity. There were four experimental treatments (i.e., control (CK), N addition (N), P addition (P), both N and P additions (N+P)), which were randomly assigned to four 6 m × 6 m plots in each block. Two chambers for greenhouse gas measurement were set up in each plot. All blocks were separated with a 3-m buffer zone. The N and P were applied twice per month from May to September during 2013-2016; we sprayed the fertilizer solution to ensure that the application was evenly distributed in the plots. The dose of N addition was 100 kg N ha-1 y-1 as NH4NO3 solution; this dose was selected because it meets the N required to sustain the local maximum vegetation productivity (Bai et al. 2010). We monitored the atmospheric N deposition and estimated it to be 20.4 kg N ha-1y-1 in our experimental site in 2012. The dose of P was 100 kg P ha-1 y-1, equivalent to previous nutrient addition experiments in grasslands (Phoenix et al. 2003). The plots with both N and P additions received the same amounts of N and P as in the N-only or P-only addition treatments. Control plots were not fertilized, but rather watered with the same amount of water as used in the fertilizer solutions; the water used to dissolve N and P was approximately 800 ml, which is equivalent to 0.8 mm for a 1 m2 plot, a pre-treatment field experiment found that this small amount of water addition did not cause any significant changes in water supply, thus unlikely to have altered the ecosystem functions.
Over the study period, soil and plant properties were measured (Supplementary Online Materials). Soil samples were collected to a depth of 10 cm once per year during 2013-2016. Six soil cores were randomly taken in each plot and mixed completely. The measured soil variables included ammonium (NH4+) and nitrate (NO3-), soil pH, soil organic carbon (SOC), soil total nitrogen (STN), soil total phosphorus (STP), soil microbial biomass carbon (MBC) and soil microbial biomass nitrogen (MBN). It should be noted that no nitrite was detected at our field site. The measured vegetation variables included aboveground biomass, plant total carbon (PTC) and total nitrogen (PTN).

Field Measurements of CH4Flux

A static chamber technique was used to measure the CH4flux (Song et al. 2009; Zhang et al. 2017). Stainless steel permanent bases (50 × 50 × 12 cm) with a 3-cm-deep groove for water sealing were inserted into soil down to a depth of 12 cm in the plots approximately one month before the experiment started in the first year. The chamber base was left in the field for one month before any flux measurement; this was to let the system evolve into an equilibrium state, avoiding any potential disturbances to methanogens and methanotrophs and thus the CH4 flux. An opaque stainless-steel top chamber (50 cm height) with heat-isolation was installed over the base, with bottom in the groove. The grooves were filled with water for sealing the chamber. Two electric fans were installed inside the chamber to mix the air in the headspace continuously and thoroughly during the measurements (Fan et al.2007). Sixty-ml gas samples were collected at 10-min intervals for 30 min using sixty-ml syringes with airtight stopcocks. Simultaneously, the air temperature and pressure in the chamber were measured, and the soil temperature and moisture were measured at a 5 to 10-cm depth using a thermometer (Spectrum Technologies, Inc. East - Plainfield, Illinois) and a portable soil moisture measuring kit ML2x (ThetaKit, Delta-T Devices, Cambridge, UK), respectively. The CH4concentrations of the gas samples were analyzed using a gas chromatograph (Agilent 7890A, Agilent Technologies Limited Co., USA) equipped with flame ionization detector (FID) using 60–80 mesh 13 XMS column (2 mm inner diameter and 2 m long), with an oven temperature of 55°C. Pure nitrogen was used as carrier gas at a flow speed of 30 ml min-1, and the CH4 flux was determined from the linear slope of the mixing ratio changes in four samples taken at 0, 10, 20 and 30 min after chamber closure. The CH4flux was measured weekly from May to September in 2013, 2014, 2015 and 2016, respectively.
The CH4 fluxes were linearly interpolated and sequentially cumulated to estimate the total flux over growing seasons (Zhang et al. 2017). The linear interpolation was adopted due to two reasons: (1) it has been widely used in previous studies and has been suggested to be effective for seasonal interpolation of ecological variables (Song et al. 2009; Nikiema et al. 2011); (2) the auxiliary data to assist annually interpolation was lacking. Throughout this analysis, positive fluxes represent gas uptake by the grassland ecosystems.

Meta-analysis

A meta-analysis was conducted to investigate the N and P impacts on CH4 across global grasslands. We collected publications by searching for “nitrogen and phosphorus”, “methane”, “grassland”, and “upland” in Google Scholar in July 2016, and later updated in September 2018. The search returned 6230 publications. A few criteria were then used to screen these publications for our purpose. The criteria applied to determine whether or not to use the data were: (1) it must be manipulation experiments with either external N and/or P addition; (2) the field measurements must cover at least one full growing season, which makes the estimation of annual budget more reliable; (3) the studies report clear information for the field site that is useful when extracting edaphic and meteorological data from global datasets. For sites without clear latitude and longitude, we Googled and found their geographic coordinates with country and site names. Finally, 35 papers were selected for data extraction. When the data were presented in figures, we extracted mean values and standard errors using GraphClick (http://www.arizona-software.ch/graphclick/). For studies with measurements from different N and P addition levels in one paper, they were extracted and treated as independent data points.
Finally, we compiled a comprehensive dataset of 382 data points for meta-analysis (Figs. S4 and S5 ). The dataset covers the major grassland types across the globe and all experiments were carried out between 1980-2017. For all field experiments, the atmospheric depositions of N and P from the global dataset were treated as background rate and were added to the reported N and/or P addition rates. Since a majority of the data points from Yu et al. did not contain N or P deposition (Yu et al. 2017), we used the extracted contemporary nutrient deposition (Mahowald et al. 2008), based on latitudes and longitudes, as nutrient inputs. Every CH4flux rate corresponds to one site with auxiliary information including latitude, longitude, and factors such as N deposition, P deposition, soil temperature (ST), soil moisture (SM), soil pH, soil organic carbon (SOC), bulk density (BD), and clay content (CL). Across the studies in our data set, the N and P treatments fell within the ranges of 0~200 kg N ha-1 and 0~50 kg P ha-1, respectively. The non-growing season CH4 flux is normally not available for most field experiments, therefore, the annual rate of CH4 uptake (kg C-CH4ha-1 yr-1) thereafter was calculated by the ratio of the growing season to non-growing season CH4 uptake as reported by a few studies (Li et al. 2012; Yue et al. 2016).
The soil properties for each point data site were extracted from global datasets according to their latitude and longitude. Specifically, the soil pH, SOC, BD and clay were retrieved from the Re-gridded Harmonized World Soil Database v1.2 in the Oak Ridge National Laboratory Distributed Active Archive Center for Biogeochemical Dynamics (available online: https://daac.ornl.gov/SOILS/guides/HWSD.html); soil temperature and moisture for the top 10 cm were from the NCEP/DOE AMIP-II Reanalysis (Reanalysis-2): http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.derived.surfaceflux.html).

Global Extrapolation

The data for the grassland distribution were derived from the Global Mosaics of the standard MODIS land cover type data product (MCD12Q1) from http://glcf.umd.edu/data/lc/index.shtml. The spatial coverage is -180.0° ~ 180.0° in longitude, and -64.0° ~ 84.0° in latitude. The global land cover data sets are available as ESRI ASCII Grid format files and was re-projected to be consistent with soil and climate datasets (Fig. S4 ). The global simulations were carried out at a spatial resolution of 0.5° × 0.5°. According to the definition of Food and Agriculture Organization of the United Nations (FAO), any geographic area dominated by natural herbaceous plants including grasslands, prairies, steppes, and savannahs with coverage of at least 10% was designated as grassland (Lathamet al. 2014). The annual CH4 uptake was determined by summing up all grassland grids (equation 1). For the model simulation, we set N and P depositions to pre-industrial level to represent ambient condition, N set to pre-industrial level and P set to contemporary level represent P addition, P set to pre-industrial and N set to contemporary represent N addition, both N and P were set to contemporary level to represent N and P concurrent treatment.
Both climate and soil data were resampled and re-projected to be consistent through NCL (NCAR Command Language, current version 6.4.0). The data with higher resolution were uniformly transformed to the lowest resolution at 0.5° × 0.5°. Based on the calculated CH4uptake rate in every grid and grid area of grassland, we scaled up the results from this analysis by multiplying them for target CH4 uptake with the corresponding grid areas currently summarized:
\(F(CH4)=\sum_{k=0}^{n}{\left(\frac{n}{k}\right){\text{CH}_{4}}_{k}A_{k}}\)(1)
where F(CH4) is the sum of the global grassland CH4 uptake expressed as Tg C-CH4year-1, and CH4k is the CH4 uptake rate in kth grid as kg C-CH4ha-1 y-1, and Ak is the area of the kth grid (Fig. S8 ).

Statistical Analysis and Uncertainty Analysis

We used a one-way ANOVA analysis (ANOVA; R package) followed by Duncan’s post-hoc test to examine the N and P impacts on annual CH4 uptake rate, soil NH4+ and NO3-, and plant N content. A stepwise multivariate analysis was adopted to evaluate the response of the soil plant and microbial communities to environmental factors as well as their correlations with each other. A number of key variables, including N input rate, P input rate, soil temperature, soil pH, soil organic carbon density, bulk density, and clay content, are kept in the multiple linear equation to quantify their impact on CH4 flux (Equation 2).
A structure equation model (SEM) was used to explore how nutrient (N, P and N + P)-induced changes in CH4 uptake considering meteorological, plant, edaphic factors, and microbial biomass. The “lavaan” package (https://cran.r-project.org/web/packages/lavaan/lavaan.pdf) for R program (version 3.4) was used for the structural equation modeling. Data used in the SEM and linear regression analyses were calculated as the mean of every year during the four years (Figs. S2 and 3 ). Prior to the SEM procedure, we conducted principal component analysis (PCA) to remove the redundant variables for meteorological factors, soil edaphic parameters, plant and microbial parameters. Meanwhile, due to the strong multicollinearity among other variables, the variance inflation factor (VIF) was used to quantitatively select the right factor for SEM model (VIF < 5). One SEM model (Fig. S3) was used to disentangle the direct and indirect environmental impacts on the CH4 uptake, and the other (Fig. 3) was used to analyze the N and P impacts on CH4 uptake. A conceptual model was created to represent the key linkages and connecters between different parameters. The fitted conceptual model was further modified with “mod_ind” (within the R package “lavaan”) to yield a reliable multivariate causal network. The SEM results were evaluated with the comparative fit index (CFI), the normed fit index (NFI), and the chi-square test (χ2). The chi-square value is the traditional measure for evaluating the overall model fit and accessing the magnitude of discrepancy between the sample and fitted covariance matrices; NFI is an incremental fit index that assesses the model by comparing the chi-square value of the model to the null model that assumes that all variables are uncorrelated; CFI is a revised form of the NFI that takes into account sample size; the model fits with NFI > 0.95 and CFI > 0.95, indicating a good SEM (Hooper et al. 2008).
The Latin Hypercube Sampling based Monte Carlo method was adopted to quantify the uncertainties of the global grassland CH4uptake under different treatments (Control, N, P and N+P conditions) estimated by the empirical models. The LHS approach has been used in our previous study (Xu 2010); detailed procedure can be found inSupplementary Online Material .