2.1. The Study Area
The Inner Mongolia Plateau (IMP) is chosen as the study area (
Figure 1). IMP boasts the most important grassland resource in northern China [
25]. IMP extends from 37.40° N to 53.38° N and from 97.2° E to 126.06° E, covering an area of 1.18 million square kilometers [
3]. IMP has an annual precipitation 50–450 mm and an annual average temperature 1–10 °C. Both precipitation and temperature display a regional transition from semi-humid areas in east to semiarid and arid areas in west. Grassland degradation is a serious threat to sustainability of these ecosystems [
26,
27,
28]. Northern China and, especially, Inner Mongolia Plateau have suffered severe ecological degradation in recent decades [
29,
30,
31]. It is important for researchers, policy makers and resource managers to understand how each vegetation community responds to climate change. There are an increasing number of studies reported in current literature. Several examine status change of vegetation evolution [
2,
32]. Some investigate coupled impact of socioeconomic transformation and climate change on vegetation growth [
9,
13,
33]. Few are about composition structure and its relationship with climate factors [
15,
34]. Most studies are focused on regional scale and interested in identifying regional patterns [
2,
35,
36]. Little effort has been exerted to explore vegetation growth response to climate change among different plant communities over various precipitation zones [
37,
38].
For the above reasons, we focused on revealing spatial and temporal change patterns of grassland communities over an east-to-west transection and how the changes are related to the fluctuations of precipitation and temperature at both regional and local scales. Therefore, we designed a horizontal transect extending from east to west across the entire study area. Due to the unique shape of the study area, we bisected the transect i at 113° E, which is an ideal position to separate the transect into two sections that cover the majorities of vegetation communities in the study area (
Figure 1). In addition, we examined four major land covers and uses in the study area, typical steppe, bare steppe, desert and farmland.
2.2. The Data and Pre-Processing
The data used in this paper include three datasets: Modis EVI (enhanced vegetation index), meteorological observations (of precipitation and temperature), and IMP land cover map 2010. EVI was adopted as a proxy of the vegetation primary productivity instead of NDVI [
35,
39,
40], since EVI overcomes the saturation and soil noise problems of NDVI by decreasing the atmospheric effect and reducing the canopy background [
41]. The 16-day EVI data with 250 m spatial resolution for 2000–2014 were collected from MODIS 13Q1, which is available at
https://lpdaac.usgs.gov/data_access/data_pool. Nine 16-day periods during the growing season each year (8 May to 28 September) were included in the analysis [
42] totally 135 periods (15 years × 9 periods). Further, the transects were clipped out in ArcGIS by two buffers with 1750 m (a length of 7 pixels) along 41° N and 43° N. This clipping process was also applied to the following two datasets.
Meteorological data were extracted from the China Meteorological Data Service Center (CMDC, 2013). There are 50 meteorological stations in the study area. Two meteorological factors, precipitation and temperature, were selected for the analysis. The original climate data were the daily record from 2000 to 2014. To be consistent with the EVI data, the accumulation of precipitation and the mean temperature for each of the nine 16-day periods during the growing season were calculated. These values at 50 stations were spatially interpolated as precipitation and temperature surfaces at 250 m resolution, respectively. The interpolations were processed using inverse distance weighting (IDW) method. Cross-validation was used to calibrate the best pair of interpolation parameters (power and radius) over the data of the nine periods of 2014 [
43]. Forty stations were chosen as the training sample to execute IDW interpolation and ten were reserved as the validation sample. The root mean square error (RMSE) was used to judge the validation outcomes. The validation result showed that the RMSE of temperature between the estimated value and validation value was around 6% and the RMSE of precipitation was about 18%. Due to the dramatic spatial variation of precipitation in IMP [
21,
44,
45], the RMSE of precipitation was much higher than the RMSE of temperature but was reasonably small.
The 2010 Inner Mongolia land cover classification map was provided by Institute of Remote Sensing and Digital Earth, China [
9]. The map was made through interpreting the Landsat images along with ground reference points. The original spatial resolution of the product was 30 m. The vegetation map was resampled to 250 m so that it had the same resolution as EVI. In brief, EVI, precipitation, and temperature were pre-processed with the same spatial and temporal resolutions and each of them had a 135-period time series (2000–2014, with 9 periods per year in growing season).
2.3. The Analysis Method
We synthesized the integrated research method depicted in
Figure 2. First, we applied Empirical Model Decomposition (EMD) to extract long-term trends of EVI, precipitation, and temperature to investigate the spatial and temporal patterns of primary productivity and climate change. Secondly, we designed two east-to-west transects across the entire study area to examine spatial variations of vegetation response to climate change over different precipitation zones. Thirdly, we analyzed different responses of four types of plant communities over the two east-to-west transects, respectively. Finally, we applied redundancy analysis to explore partial and joint effects of precipitation and temperature on primary productivity and to determine whether the effects of precipitation or temperature, or both, on primary productivity varied regionally or over different types of vegetation.
Fluctuations of vegetation growth and climate change exhibit scale-dependent cycles and a trend [
37,
46,
47]. To understand how these cyclic phenomena and the trend aligning with climate change is a key to reveal how vegetation growth respond to climate change. Algorithms in signal decomposition have been studied to extract the trend of change from non-linear and non-stationary time series data. Martínez and Gilabert [
24] successfully detected the dynamics of land cover changes by applying wavelet transformation through normalized deviation vegetation index (NDVI) in inter-annual scale. The basic function, decomposition scale and threshold selection for wavelet transformation are manually set. In addition, these manual settings are multivariate based on specific applications [
24]. As a self-adaptive approach, EMD decomposes the non-stationary signal based on its own peculiarity [
48]. Kong et al. [
39] presented a rigorous decomposition approach for MODIS NDVI with ensemble EMD. In his study, the trend is generated by the combination of lower frequency components decomposed from the original data and the residual which is an overall trend of time series [
49]. Wu and Huang [
50] claimed that the overall trend that is derived as an intrinsically fitted monotonic function where there can be at most one extremum within a given span. Therefore, in this study, we adopted the residual, the trend of overall time series [
39], to analyze vegetation response to climate change.
The algorithm of EMD was originally developed to decompose non-stationary and non-linear signals with several zero mean components [
48], which was improved by adding new stopping criteria to make EMD implementation more effective [
51]. Rilling’s EMD algorithm code was adapted in this paper. The following is a brief description of the EMD approach. Signal x(t) is the given time series of EVI, or precipitation or temperature (for details, see [
51]):
Create the top and bottom envelope with the cubic splines fitted by maximums and minimums extras of x(t).
Calculate the mean of envelopes: mean(t) = 0.5 × (emax(t) + emin(t)) (where e stands for envelop)
Extract the mean value out: x’(t) = x(t) − mean(t)
Repeat previous steps until stopping criteria are met and then extract out the intrinsic mode function (IMF): IMF = x(t) − ∑(mean(t))
Extract IMF out of x(t): x′ = x(t) − IMF
Repeat Steps 1–5 until the criterion that it cannot extract IMFs is met.
The residual is what is left: Residual = x(t) − ∑(IMF)
Our research question is whether precipitation or temperature has varied impact on vegetation growth. Pearson’s correlation analysis is often used to test the relationship between one vegetation series and one corresponding environment variable [
2]. However, it is long known that climate factors themselves show very interactive and complex correlation. Redundancy analysis (RDA), a special case of Canonical Correspondence Analysis (CCA) that integrates correspondence analysis and multi-variable regression, is currently adopted to test the partial effect of each climate variable on vegetation growth [
35,
40]. CCA is often used for an analysis where the interested variation is embedded in compositional data and is based on the unimodal model. RDA is typically applied to non-complex data based on the linear model, where the gradient detected by CCA is less than 3 [
52]. Han [
35] confirmed that the precipitation was a major climate factor affecting vegetation communities in Inner Mongolia through CCA. Sun et al. [
40] found that the precipitation was the primary factor in governing the dynamics of vegetation via RDA in Tibetan Plateau. Considering that our data are not composite data, RDA was chosen and applied with the software package CANOCO 5. In other words, we adopted RDA to interpret the relationships between types of vegetation and two climate factors.
RDA is a kind of constrained ordination, first demonstrated by Rao [
53], which aims to detect the explanation proportion of the variation of response variable contributed by each explanatory variable [
52] (pp. 630). Tsai et al. [
54] used RDA to analyze the explanation proportion of climate factors on vegetation change in Taiwan. Moreover, Li et al. [
34] found that temperature and rainfall noticeably affect vegetation composition structure through RDA. The explanation proportion is calculated through the following functions
where
is the same as
R square in multi-variable linear regression (the square of the Pearson’s coefficient),
Y is the response matrix with n cases and p response variables, and
X is the explanatory matrix with n cases and m explanatory variables.
is the sum of square of
(estimated
Y through linear regression) and
is the sum of square of
Y (observed values). Notice that
or
is sometimes described as sum of squared deviation from mean [
52] (pp. 633).
The following is the adjusted
R square equation [
55], which is also considered as explanatory proportion of
X to variance of
Y [
35].
where
n is the number of cases in study, and
m is the number of explanatory variables. To clarify the statics effectiveness of the adjust
R square, also explanation proportion, the F-statistic to test the significant should be conducted as follows [
56]:
where
p is the number of response variables. Notice that
Y here should be the standardized
Y (for details, see [
52] (pp. 634 [
56,
57]).
We conducted RDA to analyze the reactions of vegetation to its climate factors in CANOCO, which is effective software for analyzing the relation between a set of environment factors and a set of species in ecology [
35,
58]. The CANOCO v5 software package is available at
http://www.canoco5.com/. Before running the RDA, which is a multivariate regression analysis, the statistical assumptions concerning a regression analysis need to be verified. In CANOCO, partial RDA [
59], which is the RDA integrated with variation partitioning analysis, was used to evaluate the marginal explanation rate of precipitation and temperature on EVI variation.