Analysis of spatial and temporal patterns of aboveground net primary productivity in the Eurasian steppe region from 1982 to 2013

Abstract To explore the importance of the Eurasian steppe region (EASR) in global carbon cycling, we analyzed the spatiotemporal dynamics of the aboveground net primary productivity (ANPP) of the entire EASR from 1982 to 2013. The ANPP in the EASR was estimated from the Integrated ANPPNDVI model, which is an empirical model developed based on field‐observed ANPP and long‐term normalized difference vegetation index (NDVI) data. The optimal composite period of NDVI data was identified by considering spatial heterogeneities across the study area in the Integrated ANPPNDVI model. EASR's ANPP had apparent zonal patterns along hydrothermal gradients, and the mean annual value was 43.78 g C m−2 yr−1, which was lower than the global grasslands average. Compared to other important natural grasslands, EASR's ANPP was lower than the North American, South American, and African grasslands. The total aboveground net primary productivity (TANPP) was found to be 378.97 Tg C yr−1, which accounted for 8.18%–36.03% of the TANPP for all grasslands. In addition, EASR's TANPP was higher than that of the grasslands in North America, South America, and Africa. The EASR's TANPP increased in a fluctuating manner throughout the entire period of 1982–2013. The increasing trend was greater than that for North American and South American and was lower than that for African grasslands over the same period. The years 1995 and 2007 were two turning points at which trends in EASR's TANPP significantly changed. Our analysis demonstrated that the EASR has been playing a substantial and progressively more important role in global carbon sequestration. In addition, in the development of empirical NDVI‐based ANPP models, the early–middle growing season averaged NDVI, the middle–late growing season averaged NDVI and the annual maximum NDVI are recommended for use for semi‐humid regions, semi‐arid regions, and desert vegetation in semi‐arid regions, respectively.


| INTRODUCTION
Aboveground net primary productivity (ANPP) represents the major input of nutrition and energy into ecosystems, and it is an integral indicator of ecosystem functions (McNaughton, Oesterheld, Frank, & Williams, 1989). ANPP is one of the main components of the carbon cycle and one of the most important and fundamental fluxes that reflect carbon sinks/sources of ecosystems (Scurlock, Cramer, Olson, Parton, & Prince, 1999). ANPP can indicate the growth status of vegetation, for which variations over time reflect the response of ecosystems to climate change (Roy, Mooney, & Saugier, 2001). In addition, ANPP is a good index of potential economic production (food, fuel, fiber) (Scurlock et al., 1999). The spatiotemporal dynamics of ANPP have been a key research topic of Global Change and Terrestrial Ecosystems (GCTE) (Fang et al., 2003;Steffen et al., 1998).
Globally, grasslands are one of the most widespread biomes, and they are known as the prairie in North America, the pampas in South America, the veld in South Africa, the steppe in Eurasia, and the savanna in Africa and Australia (Woodward, 2008). These grasslands account for approximately 20% of the world's land surface (Lieth, 1978). From a perspective of carbon cycles of ecosystems, grasslands likely contribute an annual carbon sink of up to ~0.5 Pg C (Scurlock & Hall, 1998), and these systems amount to approximately 18% of the total current global terrestrial carbon sink (Canadell et al., 2007), playing a key role in balancing the concentration of global atmospheric greenhouses gases through carbon sequestration (Lauenroth, 1979;O'Mara, 2012). From an applied perspective, grasslands significantly contribute to resources needed for human activities by providing food for herbivores (O'Mara, 2012;Scurlock et al., 1999). Therefore, examining spatiotemporal ANPP patterns in grasslands is critical to understanding terrestrial carbon sequestration and fundamental for determining the sustainable use of grassland resources (O'Mara, 2012;White, Murray, & Rohweder, 2000).
The Eurasian steppe region (EASR), which is the largest continuous grassland biome worldwide, plays an important role in global grasslands, as do North American, South American, and African grasslands (Woodward, 2008). The EASR, which is located in northern midlatitudes, is influenced by monsoon, continental, and Mediterranean climates, and it is sensitive to global environmental change ( Figure 1).

F I G U R E 1
The geographic extent of the Eurasian steppe region and the spatial distribution of ANPP field sites. ANPP denotes the aboveground net primary productivity per year and per square meter. EASR denotes the Eurasian steppe region, BKSSR denotes the Black Sea-Kazakhstan steppe subregion, MPSSR denotes the Mongolian Plateau steppe subregion, and TPSSR denotes the Tibetan Plateau alpine steppe subregion Studies on the spatiotemporal dynamics of the ANPP of the EASR mainly focus on typical geographical units within it, such as on the Inner Mongolian temperate grasslands Ma, Liu et al., 2010;Ma, Yang, He, Hui, & Fang, 2008) and the Tibetan Plateau alpine grasslands (Zhang et al., 2014). In recent years, some studies have also analyzed the spatiotemporal dynamics of the ANPP in the grassland areas of Mongolia and Kazakhstan at a national scale (Bao et al., 2016;Eisfelder et al., 2014). However, studies on the spatiotemporal dynamics of ANPP have not yet been specifically reported for the entire EASR. There is thus a knowledge gap regarding the spatiotemporal dynamics of the ANPP for the entire EASR and the roles that the EASR plays in the global carbon budget. Therefore, studying the spatiotemporal patterns of the EASR's ANPP will further our understanding of carbon cycling mechanisms in grassland ecosystems and will prove central to the assessment of global carbon budgets.
The Normalized Difference Vegetation Index (NDVI) is the normalized reflectance difference between the satellite near-infrared band and the visible red band (Rouse, Haas, Schell, & Deering, 1974;Tucker, 1979). The NDVI represents the photosynthetic potential of the vegetation canopy and is extensively used in ecosystem monitoring (An, Price, & Blair, 2013;Box, Holben, & Kalb, 1989;Gu, Wylie, & Bliss, 2013;Gu, Wylie, & Howard, 2015;Hobbs, 1995;Paruelo, Epstein, Lauenroth, & Burke, 1997;Rouse et al., 1974;Tucker, 1979). Previous studies have shown strongly positive relationships between the NDVI and ANPP in grasslands. Developing an NDVI-based empirical remote sensing inversion model for ANPP estimation involves obtaining NDVI values over a specified time period (composite periods of NDVI data) (Reed et al., 1996). Composite periods of NDVI data may affect the accuracy of empirical NDVI-based annual ANPP estimation models.
The optimal period of the NDVI composite varies according to regional climatic conditions Mkhabela, Bullock, Raj, Wang, & Yang, 2011) and vegetation type (Jin et al., 2014;Mao, Wang, Li, & Ma, 2014). Therefore, in this study, we attempted to achieve the following objectives based on field-observed ANPP and long-term NDVI timeseries data: (1) to identify the best composite period of NDVI data for developing a robust annual ANPP estimation model designed for the entire EASR, (2) to evaluate the ANPP of the entire EASR, and (3) to explore the temporal dynamics of the EASR's TANPP and further discuss the role of the EASR in the global carbon budget.

| Study region
The Eurasian steppe region (EASR) (Figure 1) in the northern midlatitudes extends over approximately 110 longitudinal units from the grassy plains at the mouth of the Danube River in the west; across Russia, Kazakhstan, and Mongolia to the Songliao Plain in China to the east; and to the Himalayas in China to the southwest (Woodward, 2008). The EASR is the largest continuous grassland biome in the world, covering an area of 8.65 million km 2 , and the region is preserved relatively well (Woodward, 2008). The EASR is influenced by the Mediterranean climate, by the southwest monsoon of the Indian Ocean, by the East Asian monsoon and by the westerlies. Multiyear mean annual precipitation levels in the region vary from 60 mm to 1,100 mm (Appendix S1), and the mean annual temperatures (MATs) range from −9 to 20°C (Appendix S1) across the entire EASR, with rainy and hot conditions characterizing this period (Appendix S2).
Throughout the EASR, natural vegetation mainly includes meadows, meadow steppes, typical steppes, desert steppes, alpine steppes, and alpine meadows (Appendix S3, Olson et al., 2001;Editorial Committee of Vegetation Map of China Chinese Academy of Sciences, 2007;Woodward, 2008;Bao et al., 2016). The dominant grasses include perennial bunchgrasses, and constructive species belong to the Stipa genus, which refers to species controlling the structure and function of the ecosystem while their amount being not always maximum in the ecosystem (Woodward, 2008;Zhou, 1980). Chestnut soil is the main soil type (Woodward, 2008). Phytogeographically, the EASR can be divided into three subregions ( Figure 1): the Black Sea-Kazakhstan steppe subregion (BKSSR), the Mongolian Plateau steppe subregion (MPSSR), and the Tibetan Plateau alpine steppe subregion (TPSSR) (Li, 1979;Wu, 1979;Лaвpeнкo, 1959).
It is important to note that the geographical extent of the EASR has not yet been clearly defined, so we identify the regional extents of the EASR and of its three subregions based on Moderate Resolution Imaging Spectroradiometer (MODIS) Land cover data according to descriptions of the EASR provided in Лaвpeнкo (1959), Wu (1979), Li (1979 and Woodward (2008).

| Field-observed ANPP data
In this study, field-observed ANPP (ANPP obs ) was estimated as the peak aboveground biomass for harvesting during the growing season (from April to October of the year). These data were primarily obtained from three sources: (1) 1015 ANPP field observations from 206 publications (a list of the data sources can be found in Appendix S11) cited in the Web of Science (www.Webofknowledge.com) and China National Knowledge Infrastructure (http://epub.cnki.net), (2) 7 ANPP field observations from the Class A dataset of the global ANPP database of the Oak Ridge National Laboratory Distributed Active Archive Center (https://daac.ornl.gov/cgi-bin/dataset_lister. pl?p=13#grassland_anchor) and (3) 809 ANPP field observations provided by the researchers of this study. In total, 1990 site-year ANPP observations from 1831 field sites were initially collected for the entire EASR over the past three decades .
Before carrying out our analysis, we completed four tasks to eliminate unsuitable field-observed ANPP data: (1) we excluded observations missing site-description metadata (e.g., latitude or longitude), (2) we excluded observations without specific sampling times, (3) we excluded observations in ecotones of grasslands and other systems according to the Moderate Resolution Imaging Spectroradiometer (MODIS) Land cover product (MCD12C1) (https://lpdaac.usgs.gov/ dataset_discovery/modis/Modisproducts_table/mcd12c1), and (4) we excluded observations with ANPP outliers (falling outside a range of mean ± 3 standard deviations).
Thus, field ANPP observations of 1717 site years from 1539 field sites were examined in this study. These ANPP field sites spanned 28°N to 53°N in latitude, 36°E to 125°E in longitude, and 20 m to 5600 m in elevation ( Figure 1). In addition, the sampling time span ranged from 1982 to 2013 (Appendix S4).
It is worth noting that ANPP (g C m −2 yr −1 ) in this study denotes the aboveground net primary productivity per year and per square meter, and TANPP (Tg C yr −1 ) denotes the regional total aboveground net primary productivity level per year. Field data reported in dry matter form (g m −2 yr −1 ) from previous studies were converted into units of C using a conversion factor of 0.45 (Lieth & Whittaker, 1975).

Long-term NDVI time-series data
The biweekly NDVI data for 1982 to 2013 used in this study were obtained from third-generation NDVI (NDVI3 g) datasets produced through the Global Inventory Modeling and Mapping Studies (GIMMS). GIMMS NDVI3 g data (http://ecocast.arc.nasa.gov/data/ pub/gimms/3 g.v0/) were obtained at a spatial resolution of .083° by applying the 15-day maximum-value composition (MVC) technique (Holben, 1986) to observations generated by Advanced Very High Resolution Radiometers (AVHRRs) aboard National Oceanic and Atmospheric Administration (NOAA) satellites (Tucker et al., 2005;Zhu et al., 2013). GIMMS NDVI3 g data were corrected for sensor degradation, intersensor differences, cloud cover, solar zenith angle, and viewing angle effects resulting from satellite drift as well as the presence of volcanic aerosols. These data were widely used to monitor long-term vegetation activation trends Wu & Liu, 2013).

Land cover data
The land cover data used in this study were collected from the Land  (Friedl et al., 2010). The land cover data used in this study were based on the IGBP global vegetation classification scheme.

| ANPP estimation model development
The best ANPP estimation model for the entire EASR was accomplished through following three steps:

Data preprocessing
Before developing the ANPP estimation model, we completed the following three data preprocessing tasks.
where i is the month from month 1 to month 12, MNDVI i is the maximum of two NDVI images available for month i, and NDVI ia and NDVI ib are NDVI images of the first and second halves of month i, respectively.

2.
The optimal composite period of monthly NDVI data for generating annual NDVI values for ANPP estimation varies with the climate and vegetation conditions in the study area Mkhabela et al., 2011). Therefore, 13 different annual NDVI values were obtained by calculating the annual maximum NDVI and 12 different averaged NDVI values of various composite periods for April to October. Theses 12 different periods were identified by changing the start and end date of the monthly NDVI composite periods (Appendices S5 and S6).

3.
According to the geographical locations and the corresponding sampling year of field-observed ANPP, we extracted 13 annual NDVI values for each ANPP field site. In turn, we generated a dataset in which every record included field-observed data and 13 corresponding annual NDVI values. Approximately 75% of the field-observed ANPP data were randomly selected to develop the ANPP estimation model, and the remaining 25% were used to validate the model. Given that NDVI images of sparsely vegetated areas can be affected by the spectral characteristics of soils, we only analyzed grasslands with an annual maximum NDVI value of > 0.1 in this study.

Model development
Based on the preprocessed data listed above, we developed ANPP estimation models for the entire EASR through two different schemes: the Entirety Overall Scheme and the Subregions Integrated Scheme.
Then, the best ANPP estimation model specific to the entire EASR was obtained by selecting the model with higher validation accuracy. In the Entirety Overall Scheme, spatial heterogeneities across the EASR were not considered in the development of the ANPP estimation model. Thus, NDVI-based variables for ANPP estimation were considered universal for the entire EASR in this scheme. By contrast, the ANPP estimation model developed in the Subregions Integrated Scheme considered spatial heterogeneities between three subregions of the EASR.

Model validation and optimization
The best ANPP estimation model (the Overall ANPP NDVI model) for the entire EASR was identified through two phases in the Entirety Overall Scheme: (1) 52 regression models, including linear, exponential, power, and logarithmic models, were developed for field-observed ANPP data and 13 annual NDVI values for the entire EASR, and (2) the performance of these models was assessed using the coefficient of determination (R 2 , Equation 2) and the root mean square error (RMSE, Equation 3). The model that met the maximum R 2 and minimum RMSE criteria of the 52 regression models was regarded as the optimal ANPP estimation model for the entire EASR. Finally, the corresponding period was considered the optimal composite period of the NDVI for ANPP estimation in the EASR.
where R 2 is the coefficient of determination between field-observed ANPP data (ANPP obs ) and modeled ANPP data (ANPP mod ), which denotes a similar pattern between ANPP obs and ANPP mod and the fraction of ANPP obs variation that can be explained by the model. RMSE is the root-mean-square error between ANPP obs and ANPP mod , which represents biases that cause modeled ANPP data to differ from fieldobserved ANPP data. n is the number of field ANPP observations included in the dataset for model validation.
As stated above, the model with the highest R 2 and lowest RMSE of the 52 regression models developed was identified as the best model for estimating the ANPP of the entire EASR (Taylor, 2001). When no model satisfies both maximum R 2 and minimum RMSE, the existing 52 models need to be optimized (Taylor, 2001). In this study, the optimal model was obtained by generating the averaged result of ANPP estimations using the model with the highest R 2 and the model with the lowest RMSE at the pixel scale.
The optimal ANPP estimation model (the Integrated ANPP NDVI model) specific to the entire EASR in the Subregions Integrated Scheme was also determined through two steps: (1) optimal ANPP estimation models for three subregions in the EASR were obtained by developing model approaches that reflect the Entirety Overall Scheme, and (2) the best ANPP estimation models of the three subregions were integrated to generate the optimal model for the entire EASR. Finally, after comparing the Overall ANPP NDVI model with the Integrated ANPP NDVI model, the one presenting greater validation accuracy was used to estimate the ANPP of the entire EASR.

| Analysis of temporal dynamics of TANPP
TANPP time series trends for 1982-2013 were analyzed using a simple regression model. The turning points (TPs) at which trends in the TANPP time series significantly changed were identified using a piecewise linear regression (PLR) model (Toms & Lesperance, 2003) (Equation (4)).
where t is the year from 1982 to 2013 and TANPP is the regional total aboveground net primary productivity per year. α is the estimated TP of the time series, which denotes the timing of a trend change. β 0 , β 1 , and β 2 are regression coefficients, and ξ is the residual of the fit. β 1 and β 1+ β 2 denote linear TANPP trends before and after the turning point, respectively. The three regression coefficients were determined through least-square linear regression. In addition, a t-test was conducted to evaluate the necessity of introducing a TP based on the following null hypothesis: "β 2 is not significantly different from zero." A p value of <.05 was considered significant. The TANPP trends for each subperiod defined by the TP were also analyzed.

| Development and validation of the ANPP estimation model
Model accuracy levels directly affect the reliability of research conclusions. Accuracy assessments based on the remaining 25% of fieldobserved ANPP data (Appendices S5 and S6) indicated that models with both maximum R 2 and minimum RMSE were absent from the 52 regression models fitted in the Entirety Overall Scheme or in the Subregions Integrated Scheme. According to the approaches for developing an ANPP estimation model, the Overall ANPP NDVI model and the Integrated ANPP NDVI model were composed of the fitted regressions meeting the maximum R 2 criteria and satisfying the minimum RMSE criteria among the 52 regression models in both the Entirety Overall Scheme and the Subregions Integrated Scheme, respectively.
In the Entirety Overall Scheme, the regression relationship between field-observed ANPP and the annual maximum NDVI met the maximum R 2 criteria (Figure 2a). In addition, the regression correlation of field-observed ANPP with the average NDVI of the period running from July to September satisfied the minimum RMSE criteria ( Figure 2b). Fitted regressions with the maximum R 2 criteria (Figure 2a) and the minimum RMSE criteria (Figure 2b) between field-observed ANPP and NDVI values together formed the Overall ANPP NDVI model (Equation (5)). The validation of the Overall ANPP NDVI model ( Figure   4a) shows that the R 2 and RMSE between the estimated ANPP and field-observed ANPP were 0.58 and 17.24 g C m −2 yr −1 , respectively, for the entire EASR (Figure 4a).
(2) In each subregion, the fitted regressions meeting the maximum R 2 and minimum RMSE criteria formed the best ANPP estimation model for the subregion. The best ANPP estimation models for the three subregions integrated together formed the Integrated ANPP NDVI model (Equation (6)), which was optimal for ANPP estimations of the entire EASR. The validation of the Integrated ANPP NDVI model showed that the R 2 and RMSE between the estimated ANPP and field-observed ANPP were 0.65 and 11.20 g C m −2 yr −1 for the BKSSR (Figure 3c for year t. NDVI GS0710 (t), NDVI GS0609 (t), NDVI GS0408 (t), and NDVI GS0508 (t) denote, respectively, the averaged NDVI values for July to October, June to September, April to August, and May to August for year t.
Comparisons between accuracy assessments of the Overall ANPP NDVI model and the Integrated ANPP NDVI model ( Figure 4)

| Geographic ANPP patterns
The relationship meeting the maximum R 2 criteria (a) and minimum RMSE (b) criteria between field-observed ANPP data and the NDVI values of the corresponding composite period in the Entirety Overall Scheme. R 2 and RMSE denote the coefficient of determination and the root mean error, respectively. ANPP denotes the aboveground net primary productivity per year and per square meter. ANPP obs denotes field-observed ANPP data. NDVI max denotes the annual maximum NDVI. NDVI GS0709 denotes the averaged NDVI of the period running from July to September. *** indicates that a regression equation was significant at the .001 level The multiyear average ANPP exhibited pronounced spatial variations in the EASR ( Figure 5) which corresponded to different grassland types reflecting variations in hydrothermal conditions (Appendices S3 and S7). Regarding climatic conditions, the BKSSR and the MPSSR are located in temperate semi-arid regions in which water levels typically limit vegetation growth (Appendices S1 and S2).
In the BKSSR and MPSSR, desert steppes, typical steppes, meadow steppes, and meadows were found along an increasing precipitation F I G U R E 3 Fitted regressions that meet the maximum R 2 criteria and minimum RMSE criteria between field-observed ANPP data and NDVI values of the corresponding composite period for the BKSSR (a, b), MPSSR (d, e), and TPSSR(g, h) in the Subregions Integrated Scheme, and the validation of the Integrated ANPP NDVI model in the BKSSR (c), MPSSR (f), and TPSSR (i). R 2 and RMSE denote the coefficient of determination and the root mean error, respectively. ANPP denotes the aboveground net primary productivity per year and per square meter. ANPP obs denotes field-observed ANPP data. ANPP mod denotes the modeled ANPP of the Integrated ANPP NDVI model. NDVI max denotes the annual maximum NDVI. NDVI GS0710 , NDVI GS0609 , NDVI GS0408 , and NDVI GS0508 denote the averaged NDVI values for July to October, June to September, April to August, and May to August, respectively. *** indicates that a regression equation was significant at the .001 level F I G U R E 4 Validation of the Overall ANPP NDVI model (a) and Integrated ANPP NDVI model (b) for the entire Eurasian steppe region. ANPP denotes the aboveground net primary productivity per year and per square meter. ANPP obs denotes field-observed ANPP data. ANPP mod denotes modeled ANPP data. *** indicates that a regression equation was significant at the .001 level gradient from the center to the boundaries (Appendices S1 and S3), and the increasing ANPP values reached 27.69 g C m −2 yr −1 , 47.61 g C m −2 yr −1 , 56.30 g C m −2 yr −1 , and 75.74 g C m −2 yr −1 (Appendix S7). In addition, the TPSSR was subjected to a unique cold alpine climate. From the northwest to the southeast of the TPSSR with increasing precipitation and temperatures, representative vegetation types included alpine steppes and alpine meadows (Appendices S1 and S3) where ANPP values reached 51.11 g C m −2 yr −1 and 33.93 g C m −2 yr −1 , respectively (Appendix S7).

| Temporal TANPP dynamics
Interannual variations in TANPP for 1982-2013 were analyzed using a simple regression model ( Figure 6). For the entire EASR, TANPP values significantly increased from 1982 to 2013 at an annual rate of 0.84 Tg C yr −1 or 0.49% (Figure 6a). TANPP changes in the BKSSR and the MPSSR were similar to those for the entire EASR and exhibited an obvious increase (Figure 6b, c). However, the TANPP of the TPSSR showed no significant changes (Figure 6d).

EASR's TANPP variations were not continuous over the 32-year
period. The piecewise linear regression between EASR's TANPP and time (year) indicated that 1995 and 2007 were two turning points at which the EASR's TANPP time series significantly changed (Figure 6a).
The EASR's TANPP experienced a significant increase from 1982 to 1995, followed by a marked decrease that occurred from 1996 to 2007. After 2008, the EASR's TANPP increased slightly (p = 0.10).
Changes in the EASR's TANPP over the three subperiods constituted a superposition of variation trends in TANPP within the three subregions (Figure 6b, c,

| ANPP estimation model development theories
The NDVI could be used to model the spatiotemporal dynamics of ANPP for the entire EASR. Additionally, the best composite period of monthly NDVI data for estimating annual ANPP varied according to regional climatic conditions and vegetation types. The optimal composite period of monthly NDVI data for annual ANPP estimations included the middle to late growing season (June to September or July T A B L E 1 Mean (ANPP) and total values (TANPP) of annual ANPP in the Eurasian steppe region and three subregions. ANPP denotes the aboveground net primary productivity per year and per square meter. TANPP denotes the regional total aboveground net primary productivity per year. EASR denotes the Eurasian steppe region, BKSSR denotes the Black Sea-Kazakhstan steppe subregion, MPSSR denotes the Mongolian Plateau steppe subregion, and TPSSR denotes the Tibetan Plateau alpine steppe subregion. Mean denotes the mean value of annual ANPP for a region. SD denotes the standard deviation of mean annual ANPP proposed by Köppen (1923) (Appendix S8). Therefore, water is a major factor limiting vegetation growth in this subregion. Therefore, the phenological period is affected not only by temperature but also by precipitation in this subregion. In addition, biomass normally peaks later in the growing season (Barnes, Tieszen, & Ode, 1983;Reed et al., 1994). Therefore, the averaged NDVI from the middle-late growing season can be used to effectively assess annual ANPP variations in the MPSSR.
For the BKSSR, the climate, similar to that of the MPSSR, is semiarid (Table 2). Therefore, the averaged NDVI of the middle-late growing season is a key variable for assessing annual ANPP variations.
However, according to the climate index (CI Köppen ), vegetation in the BKSSR is subjected to more severe drought stress compared to the MPSSR. More specifically, desert steppes account for over 50% of the BKSSR area (Appendix S3). Desert steppes follow less predictable phenological patterns (Barnes et al., 1983;Reed et al., 1994) because they depend on less reliable precipitation events (Rauzi & Dobrenz, 1970). For desert steppes, the annual maximum NDVI can effectively reflect annual productivity levels. Therefore, ANPP variations of the entire BKSSR should be assessed by combing the annual maximum NDVI with the averaged NDVI for the middle-late growing season.
For the TPSSR, climatic patterns in the TPSSR are semi-humid as a whole (Table 2). Therefore, the phenological period in the TPSSR is relatively stable owing to weak precipitation limitations affecting vegetation growth. Moreover, biomass typically reaches its maximum value at the beginning or middle of the summer season (Barnes et al., 1983;Reed et al., 1994). As a result, the averaged NDVI of the early-middle growing season reflects annual ANPP variations in the TPSSR.
However, in the ANPP estimation models, most previous studies used NDVI values from a predefined time period based on the "normal," or mean, growing season or a subjective time period such as calendar months (i.e., from April to October, from May to September, etc.) (Gu et al., 2013;Guo et al., 2012;Paruelo et al., 1997). This method may have affected the robustness of models (Reed et al., 1996). Only a few studies conducted in the arid rangelands located in Senegal (Fuller, 1998) and central Australia (Hobbs, 1995) explored the important composite periods of the NDVI in ANPP estimation models for the specific regions.
This study suggested that the averaged NDVI of the middle-late growing season or maximum NDVI was strongly related to ANPP in the MPSSR and the BKSSR, in which the climate was semi-arid. This result agreed with conclusions obtained from research conducted in arid regions located in Africa (Fuller, 1998;Rasmussen, 1992) and arid rangeland located in central Australia (Hobbs, 1995 (Gao et al., 2013;Jiang et al., 2015;Ma, Liu, et al., 2008;Yang, Fang, Pan, & Ji, 2009)

for different grassland types in Inner
Mongolian temperate and Tibetan Plateau alpine grassland areas.
Moreover, the mean annual ANPP values of alpine steppes and alpine meadows were comparable to values reported in previous studies (Jiang et al., 2015;Ma, Liu et al., 2010). In addition, the spatial distributions of ANPP for the Tibetan Plateau alpine grasslands and the Inner Mongolia temperate grasslands were consistent with the conclusions of previous studies (Guo et al., 2012;Yang et al., 2009).
It should be noted that some uncertainties may exist in the ANPP estimation model. Similar to most traditional studies based on field investigation, our study design did not allow for quantitative assessment of the sampling quality. The field-observed ANPP data in this study combined multiple field surveys and datasets from previous studies (a list of the data sources can be found in the Appendix S11) without a consistent sampling design. Moreover, the previous studies disproportionally focused on the Mongolian Plateau steppe subregion and the Tibetan Plateau steppe subregion, while there were only a few samples from the Black Sea-Kazakhstan steppe subregion owing to either its remote location or political restrictions (Li et al., 2015). The bias of the spatial distribution of ANPP field sites might generate some uncertainties in the ANPP estimation model. Unfortunately, we are unable to evaluate the uncertainties generated from this spatial bias.
However, in our opinion, the spatial bias of ANPP field sites would have little effect on the estimation of ANPP. The theoretical connection between ANPP and the NDVI is the Monteith's (1981) equation (Paruelo et al., 1997). In addition, the NDVI is a linear indicator of absorbed photosynthetic active radiation (Sellers, 1985(Sellers, , 1987Sellers, Berry, Collatz, Field, & Hall, 1992). Therefore, the equation was applicable for the entire EASR, which expressed the relationships between ANPP and the NDVI based on few field-observed ANPP data of the BKSSR.

| The importance of the EASR in global carbon cycling
The EASR has been playing a significant role in global carbon sequestration. EASR's TANPP was found to be 378.97 Tg C yr −1 (Table 1), which represented 8.18%-36.03% of that of all grasslands. According to previous studies (Bazilevich et al., 1971;Parton et al., 1995;Whittaker & Likens, 1975;Xia et al., 2014), the TANPP of all grasslands is 1423 Tg C yr −1 -4635 Tg C yr −1 (Appendix S10). EASR's TANPP was higher compared to the TANPP for grasslands in North America, South America, and Africa (Appendix S10). The mean value of annual ANPP for the entire EASR was recorded as 43.78 ± 22.77 g C m −2 yr −1 (Table 1), which was lower than that for the global grasslands average (Bazilevich et al., 1971;Parton et al., 1995;Whittaker & Likens, 1975;Yang et al., 2008;Xia et al., 2014) and North American (Bazilevich et al., 1971;Xia et al., 2014), South American (Xia et al., 2014), and African grasslands (Xia et al., 2014)  the MAP being comparable between them (Appendix S1). The other reason was that ANPP was affected by not only annual precipitation but also the precipitation seasonal distribution (Guo et al., 2012).
Rain and heat occurred over the same period both in the MPSSR and in the TPSSR (Appendix S2). However, the rainy season did not coincide with the hot season in the BKSSR (Appendix S2), which was not good for annual herbaceous growth. The mean annual ANPP of the TPSSR was lower that of the MPSSR owing to the limitations of low temperature for vegetation growth under the special cold alpine environment of the Tibetan Plateau Kato et al., 2006).

| TANPP variation trends
The TANPP for the entire EASR showed an obvious increase from 1982 to 2013. Compared to other important natural grassland regions, the increasing TANPP of EASR estimated in this study (0.53  was lower than that of the global grasslands average (2.43 Tg C yr −1 ) and African grasslands (1.21 Tg C yr −1 ) and higher than that of North (0.33 Tg C yr −1 ) and South (-0.44 Tg C yr −1 ) American grasslands (Xia et al., 2014). Despite a statistically significant uptrend in EASR's TANPP for the past three decades , change trends were not temporally homogeneous over the whole period ( Figure 6). The EASR's TANPP significantly increased from 1982 to 1995, followed by a marked decrease from 1996 to 2007 and a weakening uptrend from 2008 to 2013.
Vegetation growth was influenced by precipitation in the BKSSR and the MPSSR (Bao et al., 2014(Bao et al., , 2015Jiao et al., 2017), and by temperature in the TPSSR (Zhang et al., 2014). During 1982-1995, the TANPP of EASR increased significantly because TANPP increases in the three subregions occurred owing to increasing precipitation in the BKSSR and MPSSR (Bao et al., 2015) and warming in the TPSSR (Piao et al., 2011;Zhou et al., 2001). In the period of 1996-2007, EASR's TANPP decreased apparently because of a decrease in TANPP in the MPSSR and TPSSR, which was caused by decreasing precipitation, especially summer precipitation in the MPSSR (Bao et al., 2014(Bao et al., , 2015Piao et al., 2011) and by decreasing temperature in the TPSSR (Piao et al., 2011).
The fact that the TANPP trend in EASR reversed from positive during 1982-1995 to negative during 1996-2007 was similar to TANPP variations in other grasslands (e.g., North America, South America and Africa). Xia et al. (2014) showed that the increase in the TANPP of other grasslands globally during the period of 1982-2006 reversed around 1994. In addition, a similar prominent reversal occurred in vegetation NDVI trends by the mid or late-1990s for the temperate ecosystems of Eurasia, with a pronounced increase occurring before the mid or late 1990s and a decline (or a weakened increase) occurring afterward (Mohammat et al., 2013;Piao et al., 2011). This conclusion was further confirmed by our analysis results.
In addition, our results based on an extended period ( It is important to note that the interpretations of the causes of variation trends in EASR's TANPP mentioned above were generated by reviewing current studies conducted in the MPSSR, TPSSR, and Kazakhstan. The scientific validity of these interpretations needs to be further confirmed in the future. The mechanisms of the temporal dynamics of TANPP are complex in the EASR because of its vast area, complex topography, and diverse climate regimes, which will be discussed in greater detail and depth in our future studies.

| CONCLUSIONS
To the best of our knowledge, this study was the first to assess the role of the entire EASR in the global carbon cycle. According to our analysis, although EASR's ANPP is lower than that of North American, South American, and African grasslands, EASR's TANPP is higher than that of grasslands in North America, South America, and Africa, accounting for 8.18%-36.03% of that of all grasslands.
The EASR's TANPP displayed an obvious uptrend over the past three decades, for which the increasing rate was higher than that in North and South American grasslands over the same period. This result indicates the indispensable and ever-increasing role of the EASR in global carbon sequestration. Moreover, there were several important turning points of the EASR's TANPP trend in the past 30 years.
In addition, our analysis also demonstrates that the best composite period of NDVI data for annual ANPP estimation varies with climate and vegetation in the study region. More specifically, the early-middle growing season averaged NDVI, the middle-late growing season averaged NDVI, and the annual maximum NDVI should be, respectively, applied to semi-humid regions, semi-arid regions, and desert vegetation in semi-arid regions.

ACKNOWLEDGMENTS
We are grateful to Oak Ridge National Laboratory (ORNL) Active Archive Center for offering field-observed data in global ANPP da-

CONFLICT OF INTEREST
None declared.

DATA ACCESSIBILITY
Field-observed ANPP data from published papers in this study can be available by contacting with Guirui Yu (yugr@igsnrr.ac.cn) and Cuicui Jiao (jiaocuicui1987@sina.cn).