Spatiotemporal Variation Characteristics of Ecosystem Service Losses in the Agro-Pastoral Ecotone of Northern China

Being subject to climate change and human intervention, the land-use pattern in the agro-pastoral ecotone of Northern China has undergone complex changes over the past few decades, which may jeopardize the provision of ecosystem services. Thus, for sustainable land management, ecosystem services should be evaluated and monitored. In this study, based on Landsat TM/ETM data, we quantitatively evaluated the losses of ecosystem service values (ESV) in three sections of the agro-pastoral ecotone from 1980–2015. The results were as follows: (1) the main characteristic of the land conversions was that a large area of grassland was converted into cultivated land in the agro-pastoral ecotone; (2) on the spatial scale, the ESV losses of the agro-pastoral ecotone can be called an “inclined surface” in the direction of the northeast to southwest, and the northeastern section of the agro-pastoral ecotone lost more ESV than the middle and northwest sections (p < 0.05), on the temporal scale, the order of losses was 1990–2000 > 1980–1990 > 2000–2015; (3) the agro-pastoral ecotone lost more ESV, which was mainly due to four kinds of land conversion, which were grassland that was transformed into cultivated land, grassland transformed into unused land, grassland transformed into built-up areas, and cultivated land transformed into built-up areas; (4) although these land conversions were curbed after the implementation of protection policies at the end of the 1990s, due to reduced precipitation and increasing temperatures, the agro-pastoral ecotone will face a more severe situation in the future; and, (5) during the period of 1990–2015, the overall dynamic processes of increasing population gradually expanded to the sparsely populated pastoral area. Therefore, we believe that human interventions are the main cause of ecological deterioration in the agro-pastoral ecotone. This study provides references for fully understanding the regional differences in the ecological and environmental effects of land use change and it helps to objectively evaluate ecological civilization construction in the agro-pastoral ecotone of Northern China.


Introduction
Land use has profound ecological and environmental effects [1]. Land use alters the natural surface landscape and it simultaneously causes changes in ecosystem function and structure [2], thus playing a decisive role in maintaining ecosystem services [3,4]. Changes in land use patterns inevitably lead to changes in ecosystem services [5]. The ecological footprint [6], eco-environmental quality index [7], ecosystem services value [5], etc can characterize the effects of land use change. Choosing the ecosystem service values (ESV) that can be expressed in the form of monetary value to quantify the effects of land use change will alert relevant managers and decision-makers [8]. Costanza et al. (1997) evaluated global ESV by using a meta-analysis, determining the estimation principle and the scientific significance of ecosystem services [9]. Since then, further studies have reported that change in land-use/land-cover (LULC) strongly influences ecosystem services, such as rapid urbanization on the ecosystem service value [10], soil carbon [3], the surface terrestrial biogeochemical cycle [11,12], nutrient cycling [13], and recreation and aesthetic value [14]. Xie et al. (2008Xie et al. ( , 2015 formulated a table of ESV per unit area of different terrestrial ecosystems in China, which was concise and reliable and it has been widely used [8,15,16]. The major ecosystem services vary with different land use types. The LULC influences the main ecological process of ecosystem, including energy exchange, hydrological cycle, soil erosion, biogeochemical cycle, and so on, which further change the provision of ecosystem services. For instance, for cultivated land, its supply capacity for agricultural products is higher than its service competence of regulation, culture, and support. On the other hand, for natural grassland, its service capacity for regulation and support is higher than its product supply. Based on the evaluation of farmland ecosystem services, Ma et al. (2015) indicated that farmland ecosystems not only provided services for human beings, but they also consumed them (consumed economic investment, caused soil inorganic erosion, gave rise to pesticide pollution, etc.) [17]. Changes in the land-use pattern have an obvious influence on ecosystem services. Su et al. (2012) indicated that the landscape fragmentation, which was generated by the expansion of artificial land-use in the process of urbanization, might have a negative influence on ecosystem services [18]. Due to the preference of ecosystem services requirement, when involved in land-use administration, people often focus on one or a few types of ecosystem services, and thus unintentionally or intentionally affect the provision of other types of ecosystem service. For example, due to the blind pursuit of food, people occupy natural grassland to substantially expand the cultivated land, in this way, with the increment of grain output, the services of gas regulation, climate regulation, hydrology regulation, and soil conservation are dramatically weakened, i.e., the ESV losses is the value corresponding to the weakened services.
The agro-pastoral ecotone of Northern China (hereafter called the agro-pastoral ecotone) is a transition zone between agricultural cultivation and animal husbandry [19]. The agro-pastoral ecotone is the largest in area and the longest in spatial scale of all the typical ecotones in China [20,21]. Due to high migration inwards from other regions, grassland reclamation has been frequently undertaken since the end of the Qing dynasty in the nineteenth century [22], which neglects the natural laws of landscape ecological evolution. Although the agro-pastoral ecotone has some attributes of a transition zone in ecological and geographical aspects, it does not have typical transition zone characteristics: in the context of the unstable soil matrix, the land use structure is extremely irrational; its species diversity is not rich but is rather is extremely poor [23]; drought, rainstorms, hail, sandstorms, freezing, and other meteorological disasters are frequent [24]; and, the regional economy has a weak foundation, with 56.3% of the region belonging to national poverty counties where the available land resources and ecological carrying capacity are in sharp decline [22].
The agro-pastoral ecotone is pivotal for ecological security, social stability, and economic development; it is the origin of major Chinese rivers, a key area in resource exploitation, a multiethnic area, and an economic link between an agricultural zone and a pastoral area. Since the implementation of Chinese reform and opening up policies, under the dual influence of human activities and climate change, the agro-pastoral ecotone has experienced drastic changes in land use patterns [25,26].
Many scholars have studied the change in ESV that is caused by land use change in the agro-pastoral ecotone, for example, the Chifeng in Inner Mongolia [27], Ordos City [28], west Jinlin [8], the whole agro-pastoral ecotone [29], and at the national scale [15]. All of these studies showed that human interference had greatly changed the LULC, simultaneously affecting the ecosystem service provision. However, land use change has strong regional characteristics with different driving factors. The agro-pastoral ecotone has a large span from northeast to southwest, and there may be significant differences in regional precipitation, temperature, and population. However, there is a lack of relevant studies (simultaneously holistic and segmented) to assess the spatial and quantitative variation in the ESV losses that are caused by land use change in the agro-pastoral ecotone.
Therefore, the main objectives of this study are to: (1) analyze the spatiotemporal processes of land-use conversion in the three sections of the agro-pastoral ecotone; (2) quantitatively assess the ESV losses under the dual influences of population growth and climate change; and, (3) evaluate the ecological consequences of LULC change and their practical implications.

Study Area
Zhao (1953) first proposed the concept of the agro-pastoral ecotone [19]. In recent years, the location of the agro-pastoral ecotone has been mainly redefined based on agro-meteorological factors [30][31][32]. In this study, the agro-pastoral ecotone is assumed to include nine provinces ( Figure 1A), and it is mainly divided by the 400-mm rainfall contour, the agricultural zone to the east and south, and the pastoral area to the west and north ( Figure 1B). The agro-pastoral ecotone has a total area of 654,000 km 2 , which mainly covers arid and semiarid regions, and the annual average precipitation is approximately 300-450 mm. The annual average temperature is approximately 2-8 • C, with a growing season from April to September, and there are 30-100 windy (>5 m/s) days of sand movement. Regarding the landscape, cultivated land and grassland are mutually staggered (Figure 2A,B). In this study, in order to fully understand the regional differences in the ecological and environmental effects of land use change, the agro-pastoral ecotone is divided into three parts: the northeastern section, middle section, and northwestern section ( Figure 1C). Jinlin [8], the whole agro-pastoral ecotone [29], and at the national scale [15]. All of these studies showed that human interference had greatly changed the LULC, simultaneously affecting the ecosystem service provision. However, land use change has strong regional characteristics with different driving factors. The agro-pastoral ecotone has a large span from northeast to southwest, and there may be significant differences in regional precipitation, temperature, and population. However, there is a lack of relevant studies (simultaneously holistic and segmented) to assess the spatial and quantitative variation in the ESV losses that are caused by land use change in the agro-pastoral ecotone. Therefore, the main objectives of this study are to: (1) analyze the spatiotemporal processes of land-use conversion in the three sections of the agro-pastoral ecotone; (2) quantitatively assess the ESV losses under the dual influences of population growth and climate change; and, (3) evaluate the ecological consequences of LULC change and their practical implications. Zhao (1953) first proposed the concept of the agro-pastoral ecotone [19]. In recent years, the location of the agro-pastoral ecotone has been mainly redefined based on agro-meteorological factors [30][31][32]. In this study, the agro-pastoral ecotone is assumed to include nine provinces ( Figure 1A), and it is mainly divided by the 400-mm rainfall contour, the agricultural zone to the east and south, and the pastoral area to the west and north ( Figure 1B). The agro-pastoral ecotone has a total area of 654,000 km 2 , which mainly covers arid and semiarid regions, and the annual average precipitation is approximately 300-450 mm. The annual average temperature is approximately 2-8 °C, with a growing season from April to September, and there are 30-100 windy (>5 m/s) days of sand movement. Regarding the landscape, cultivated land and grassland are mutually staggered ( Figure  2A and 2B). In this study, in order to fully understand the regional differences in the ecological and environmental effects of land use change, the agro-pastoral ecotone is divided into three parts: the northeastern section, middle section, and northwestern section ( Figure 1C).

Assignment of Ecosystem Service Values (ESV)
There are some issues that need to be addressed when applying the framework of Costanza et al. (1997) to China. For instance, the total ESV likely reflect the level of ecosystem services in western countries; thus, the ESV are vastly overestimated for developing countries, such as China. However, the ESV for cultivated lands are thought to be underestimated [15,33]. When considering the actual situation in China, Xie et al. (2015) further developed the framework of Costanza et al. (1997), emphasizing the vital function of wetlands [34,35]. Permanent ice and snow, beach and shore, bottomland, and swampland were extracted from the original primary land-use types as a new primary LULC type, namely wetland (Table S1). Therefore, there were seven final primary LULC types on the raster maps ( Figure S1).
In this study, we used the framework of Costanza et al. (1997), however we employed ESV that were developed for China by Xie et al. (2015) (Table S2). According to the theory of ecosystem service evaluation, there are 21 types of land-use conversions that result in ecosystems losing ESV (Table 1).

Assignment of Ecosystem Service Values (ESV)
There are some issues that need to be addressed when applying the framework of Costanza et al. (1997) to China. For instance, the total ESV likely reflect the level of ecosystem services in western countries; thus, the ESV are vastly overestimated for developing countries, such as China. However, the ESV for cultivated lands are thought to be underestimated [15,33]. When considering the actual situation in China, Xie et al. (2015) further developed the framework of Costanza et al. (1997), emphasizing the vital function of wetlands [34,35]. Permanent ice and snow, beach and shore, bottomland, and swampland were extracted from the original primary land-use types as a new primary LULC type, namely wetland (Table S1). Therefore, there were seven final primary LULC types on the raster maps ( Figure S1).
In this study, we used the framework of Costanza et al. (1997), however we employed ESV that were developed for China by Xie et al. (2015) (Table S2). According to the theory of ecosystem service evaluation, there are 21 types of land-use conversions that result in ecosystems losing ESV (Table 1).

Calculating ESV Losses
Based on the assessment model of Costanza et al. (1997), the equations used to evaluate ESV are as follows: where, V k is the ESV of the LULC type k, ESV ki is the i type of ESV s of the LULC type k, V j is the ESV of the LULC type j, and ESV ji is the i type of ESV of the LULC type j. If V k > V j , then: where ESV loss is the loss of ESV when the LULC type k transformed into the LULC type j; ∆A is the corresponding land-use transformation area, ESV t is the ESV losses per unit area for all of the land-use transformations in the research region, and A is the research region area. To make the agricultural zone, agro-pastoral ecotone, and pastoral area comparable, the ESV losses that are mentioned below refer to ESV t .

Data Sources
We utilized LULC maps with 30 m grids from 1980, 1990, 2000, and 2015. The Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences provided these maps. Based on Landsat TM/ETM data, the images were geometrically corrected and geo-referenced. Additionally, an outdoor survey and random sample check were performed, which confirmed that the average accuracy of the interpretation of LULC changes exceeded 90% [36,37]. The Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences provided the annual precipitation, average temperature, and population data with 1 km grids.
Furthermore, we used a dynamic ESV assessment method to validate our results. The net primary productivity (NPP) was derived while using the improved Carnegie-Ames-Stanford approach (CASA) model [35,38]. The spatiotemporal regulation factor of soil conservation was derived from the Revised Universal Soil Loss Equation (RUSLE), which is the most widely used equation in the world for estimating soil erosion.

Data Statistical Analysis
The trend analysis of the ESV losses was implemented by moving weighted trend surface analysis [39,40], which was calculated in MATLAB 2014a (MathWorks, Natick, MA, USA). Appendix A provides the details of the process. The Mann-Kendall test was used to analyze the shifting trends of annual precipitation, average temperature, and population. Appendix B provides the details of the method. The directional distributions of population patterns, as performed by standard deviation ellipse analysis [41,42], were analyzed in ArcGIS 10.2 (Esri, Redlands, CA, USA). More detailed information on the calculation process can be found in Appendix C.
Partial least squares-discriminant analysis (PLS-DA), which is a multivariate statistical method, is a combination of partial least squares regression and the discrimination analysis for classification tasks [43,44]. PLS-DA can reduce the influence of multicollinearity among indices, which makes the analysis of the relationship between the independent and dependent variables more reasonable [45]. Moreover, the variable importance in projection (VIP) score reflects the influence of the independent variable on the dependent variable [46,47], and its equation is: Rd(Y;t 1 ,··· ,t m ) , where VIPj represents the VIP value of the j-th index, p is the number of selected index variables, m is the extracted PLS principal component number, Rd(Y; t h ) represents the explanatory power of axis t h on Y, ω 2 hj is the j-th weight of axis ω h , and Rd(Y; t 1 , . . . , t m ) represents the cumulative explanatory power of axis t 1 , . . . , t m on Y. Herein, PLS-DA was used to differentiate the ESV losses at the spatiotemporal scale, which was performed in the SIMCA 14.1 software (Umetrics, Umeå, Västerbotten, Sweden).

Spatiotemporal Processes of Land-Use Conversions in the Agro-Pastoral Ecotone
For "land converted from" ( Correlation is the relationship between the values of two variables (or indicators) that are professionally related to each other to some extent, and the values of the variables either increase synchronously or one increases and the other decreases. The correlation matrix, which is also called the correlation coefficient matrix, is composed of the correlation coefficients among the columns of For "land converted to" ( Figure 3B,D,F), the overall characteristic was that most land types were transformed into cultivated land. During 1980During -1990During and 1990During -2000, the area change of cultivated land in the northwestern section was the greatest, presenting values of 15.30 and 15.96%, respectively. The next greatest change was in the cultivated land of the middle section in 1980-1990 and 1990-2000, with values of 11.37 and 11.53%, respectively. In addition to the cultivated land, the built-up areas and the unused land of the three sections showed a significant increasing trend in 2000-2015. The above results indicate opposing trends in grassland and cultivated land in the agro-pastoral ecotone, which may completely reflect the evolution process of the agro-pastoral pattern.
Correlation is the relationship between the values of two variables (or indicators) that are professionally related to each other to some extent, and the values of the variables either increase synchronously or one increases and the other decreases. The correlation matrix, which is also called the correlation coefficient matrix, is composed of the correlation coefficients among the columns of the matrix. That is, the element of the i-th row and j-th column in the correlation matrix is the correlation coefficient of the i-th and j-th column of the original matrix. In the northeastern section, the proportion of cultivated land was negatively correlated with grassland, water areas, and wetland, with statistically significant (p < 0.05) correlation coefficients that are lower than −0.8 ( Figure 4A), while the proportion of wetland was statistically significantly positively correlated with grassland and water areas. In the middle section of the agro-pastoral ecotone, the proportion of cultivated land was strongly negatively correlated with the built-up areas, and there was a statistically significant positive correlation between the proportion of grassland and wetland ( Figure 4B). In the northwestern section, the proportion of grassland was significantly negatively correlated (p < 0.05) with forestry areas and built-up areas, with correlation coefficients lower than −0.8. The built-up areas were significantly positively correlated with forestry areas and wetlands, with correlation coefficients that are greater than 0.8. The forestry areas and wetlands were significantly positively correlated ( Figure 4C).

Analysis of ESV Losses in the Agro-Pastoral Ecotone
The trend surface polynomial analysis method simulates the distribution of the geographical elements in space, displaying the variation trend of elements in geographical space. Generally, the higher the polynomial order, better the fitting effect of trend surface near data collection point will be. However, when the order is excessively high, the interpolation and extrapolation properties will significantly worsen, and therefore the order should not exceed five in the trend surface analysis [48]. Therefore, to select the optimal model, we simulated the ESV losses with five orders ( Figure 5A). The goodness of fit (R 2 ) of the first order trend surface model, which was the highest of all orders, was 0.58, 0.60, and 0.51 for the periods 1980-1990, 1990-2000, and 2000-2015, respectively. Thus, it was used to assess the trend of ESV losses in response to land-use conversions. The first order trend surface, which delineated the overall distribution of ESV losses, was similar to an "inclined surface". Figure S2 of Supplementary Files provides the schematic diagram of the first order trend surface. The losses in the northeastern section are the highest, followed by the middle section and the northwestern section ( Figure 5B-D). On the temporal scale, the order of ESV losses was 1990-2000 > 1980-1990 > 2000-2015 ( Figure 5E). The parameters of the first order trend surface model were statistically significant (p < 0.05), and the estimated regression equations fitted well with the measured values (R 2 > 0.5), which were high for the cross-section data (Table S3).
Additionally, we used another method to verify the reliability of our results, the spatiotemporal dynamic equivalence factor method, to assess the ESV losses in 2010-2015 (Appendix D and Figure  S3). In the spatiotemporal dynamic equivalence factor method, the equivalence factor table of ESV is revised and supplemented to establish the spatiotemporal dynamic evaluation methods for different ecosystem types. The result also showed that it was a first order trend surface model. Furthermore, when compared with the method of equivalence factor per unit area that was used in this study, the spatiotemporal dynamic equivalence factor method indicated that the three sections lost more ESV ( Figure S3). northwestern section, the proportion of grassland was significantly negatively correlated (p < 0.05) with forestry areas and built-up areas, with correlation coefficients lower than −0.8. The built-up areas were significantly positively correlated with forestry areas and wetlands, with correlation coefficients that are greater than 0.8. The forestry areas and wetlands were significantly positively correlated ( Figure 4C).

Analysis of ESV Losses in the Agro-Pastoral Ecotone
The trend surface polynomial analysis method simulates the distribution of the geographical elements in space, displaying the variation trend of elements in geographical space. Generally, the higher the polynomial order, better the fitting effect of trend surface near data collection point will be. However, when the order is excessively high, the interpolation and extrapolation properties will significantly worsen, and therefore the order should not exceed five in the trend surface analysis [48].   Figure 6B, Groups 1 and 2 lost more ESV, as they were close to more land-use conversions, while Groups 3 and 4 were far away from these land-use conversions and they thus did not lose too many ESV.
Additionally, an important discriminative variable usually has a VIP score that is greater than 1. As shown in Figure 6C, the VIPs of grassland transformed into cultivated land, grassland  Figure 6B, Groups 1 and 2 lost more ESV, as they were close to more land-use conversions, while Groups 3 and 4 were far away from these land-use conversions and they thus did not lose too many ESV. transformed into unused land, grassland transformed into built-up areas, and cultivated land transformed into built-up areas were greater than 1, which meant that the classifications of the three regions were distinct from each other on the spatiotemporal scale.    Additionally, an important discriminative variable usually has a VIP score that is greater than 1. As shown in Figure 6C, the VIPs of grassland transformed into cultivated land, grassland transformed into unused land, grassland transformed into built-up areas, and cultivated land transformed into built-up areas were greater than 1, which meant that the classifications of the three regions were distinct from each other on the spatiotemporal scale.
To show the areas of ESV change, the ESV losses have been visualized onto maps ( Figure S4)

Variance Analysis and Trend Analysis of Precipitation, Temperature and Population in the Agro-Pastoral Ecotone
It is necessary to analyze their regional differences and changing trends due to the significant influence of climate and human factors on LULC and ESV. Several analyses are made as follows: variance analysis (precipitation, temperature, and population) in the three sections of agro-pastoral ecotone, trend analysis of annual averages (precipitation, temperature, and population) in the three sections of agro-pastoral ecotone, interdecadal shift analysis of 400 mm rainfall contour, and standard deviation ellipse of population.
Although there were no significant differences in the average annual precipitation among the three regions ( Figure 7A), the 400 mm rainfall contour strongly fluctuated during 1980-2015, especially in the northeastern and middle section of the agro-pastoral ecotone ( Figure 7B). The average annual precipitation had decreasing trend lines in the three regions ( Figure 7C), i.e., the agro-pastoral ecotone tended towards aridness. However, the Mann-Kendall test showed that the trend of precipitation reduction was not significant (|UF| < 1.96) ( Figure 7D).

Variance Analysis and Trend Analysis of Precipitation, Temperature and Population in the Agro-Pastoral Ecotone
It is necessary to analyze their regional differences and changing trends due to the significant influence of climate and human factors on LULC and ESV. Several analyses are made as follows: variance analysis (precipitation, temperature, and population) in the three sections of agro-pastoral ecotone, trend analysis of annual averages (precipitation, temperature, and population) in the three sections of agro-pastoral ecotone, interdecadal shift analysis of 400 mm rainfall contour, and standard deviation ellipse of population.
Although there were no significant differences in the average annual precipitation among the three regions ( Figure 7A), the 400 mm rainfall contour strongly fluctuated during 1980-2015, especially in the northeastern and middle section of the agro-pastoral ecotone ( Figure 7B). The average annual precipitation had decreasing trend lines in the three regions ( Figure 7C), i.e., the agropastoral ecotone tended towards aridness. However, the Mann-Kendall test showed that the trend of precipitation reduction was not significant (|UF| < 1.96) ( Figure 7D). There was no difference in the average temperatures between the northeastern section and the middle section of the agro-pastoral ecotone, however the temperatures of these two sections were significantly different from those of the northwestern section ( Figure 8A). The average temperature had increasing trend lines in all three regions ( Figure 8B), i.e., the agro-pastoral ecotone tended to have a higher temperature. The Mann-Kendall test showed that this increasing trend was significant (|UF| > 1.96) ( Figure 8C). Furthermore, when compared with the overall agro-pastoral ecotone, this trend was more obvious in the middle section and the northwestern section ( Figure 8C). There was no difference in the average temperatures between the northeastern section and the middle section of the agro-pastoral ecotone, however the temperatures of these two sections were significantly different from those of the northwestern section ( Figure 8A). The average temperature had increasing trend lines in all three regions ( Figure 8B), i.e., the agro-pastoral ecotone tended to have a higher temperature. The Mann-Kendall test showed that this increasing trend was significant (|UF| > 1.96) ( Figure 8C). Furthermore, when compared with the overall agro-pastoral ecotone, this trend was more obvious in the middle section and the northwestern section ( Figure 8C).  As shown in Figure 9A, there were significant differences in the population among the three sections. Standard deviation ellipse analysis further analyzed the changes in population. During the period of 1990-2015, the overall dynamic processes of increasing population gradually expanded to the sparsely populated pastoral area ( Figure 9B). The Mann-Kendall test showed that the increasing trend of population was significant in the middle section and the northwestern section (|UF| > 1.96) ( Figure 9C). Furthermore, when compared with the trend in the overall agro-pastoral ecotone, the trend was more obvious in the middle section. As shown in Figure 9A, there were significant differences in the population among the three sections. Standard deviation ellipse analysis further analyzed the changes in population. During the period of 1990-2015, the overall dynamic processes of increasing population gradually expanded to the sparsely populated pastoral area ( Figure 9B). The Mann-Kendall test showed that the increasing trend of population was significant in the middle section and the northwestern section (|UF| > 1.96) ( Figure 9C). Furthermore, when compared with the trend in the overall agro-pastoral ecotone, the trend was more obvious in the middle section.

Human Interferences Led to ESV Losses in the Agro-Pastoral Ecotone Losing
In semi-arid areas, grassland husbandry is the main land use type. It is a common practice in the United States of America (USA), Australia, and other semi-arid areas. However, the agro-pastoral

Human Interferences Led to ESV Losses in the Agro-Pastoral Ecotone Losing
In semi-arid areas, grassland husbandry is the main land use type. It is a common practice in the United States of America (USA), Australia, and other semi-arid areas. However, the agro-pastoral ecotone of Northern China, which also contains semi-arid areas, is dominated by rain-fed agriculture. This agro-pastoral ecotone is a special ecological-socia-leconomic system that was formed by the large-scale arrival of immigrants and extensive reclamation. The three sections of the agro-pastoral ecotone have undergone drastic LULC change, with a sharp increase in cultivated land that is mainly caused by the overexploitation of traditional natural grassland, since the 1980s. The consequences of this were that the ecotone lost more ESV. Terrestrial ecosystems not only continuously provide food, fiber, fuel, and other products for humans, but also simultaneously supply public goods that cannot be delivered by the market, such as climate regulation, hydrology regulation, waste regulation, soil conservation, and aesthetic landscape provision, bringing substantial welfare to human [8,9,49]. This could cause ecosystems to have high environmental load and low sustainability, purely for the purpose of economic development and personal needs [23,50]. In this study, the results of the moving weighted trend surface analysis showed that the ESV losses pattern of the agro-pastoral ecotone was similar to an "inclined surface", i.e., the northeastern section of the ecotone lost more ESV than the other sections. This virtual "surface" is an expression of ecological vulnerability on a macro scale that is more sensitive to external disturbances [51,52]. The surface, which is a visible manifestation of strained human-land relations, was mainly formed by the long-term improper use of land resources and it provides a quantitative description of ESV losses.
The results show that the total population of the northeastern section of the agro-pastoral ecotone was significantly larger than that of the middle section and the northwestern section. Additionally, based on the results of the standard deviation ellipse, it was found that the overall population dynamic process gradually expanded to the sparsely populated pastoral area. Thus, rapid population growth may have been the main driver behind these land conversions. This finding is consistent with the analysis for Northern China and the Loess Plateau by Liu et al. (2014). It is noteworthy that the Mann-Kendall trend analysis indicated that the population showed a significantly increasing trend in the middle section and the northwestern section (|UF| > 1.96), while the population growth in the northeastern section was not significant (|UF| < 1.96).
By quantitative analysis of PLS-DA in this study, we have shown that land transformations (grassland transformed into cultivated land, grassland transformed into unused land, grassland transformed into built-up areas, and cultivated land transformed into built-up areas) were decisive land conversions that led to a vast reduction in ESV. Due to the restraint of water, soil fertility, and other factors, crop productivity in those areas was low, which led to a vicious circle of population growth-land reclamation-ecosystem deterioration-further reclamation. The grain yields were able to meet people's needs in years of plentiful rainfall; however, they dramatically decreased or fell to zero in drought years [36,53]. To acquire enough food, people cannot help but constantly reclaim new grassland, simultaneously abandoning infertile cultivated land [24,54], which gradually transforms into unused land (e.g., sandy land, salina, and bare soil). By analyzing the emergy differences among the terrestrial ecosystems, from the perspective of the environment and productivity, it has been proven that rain-fed artificial grassland has a higher development potential in the agro-pastoral ecotone [23]. Moreover, under urbanization and increasing population, urban buildings and rural settlements have encroached onto the grasslands of the agro-pastoral ecotone [55]. The ecotone is also a key area of resource exploitation, as it is rich in mineral resources [22]. Considerable amounts of grassland have been transformed into built-up areas, including large industries, oil fields, salt fields, quarries, roads, and other specialized land uses. However, most of the regions in the agro-pastoral ecotone were national poverty counties, where economic development came at the cost of the ecological environment [56]. Extensive management and long-term unrestricted mining have caused the agro-pastoral ecotone to lose substantial ESV [15,26].

Natural Factors that Have Led to ESV Losses in the Agro-Pastoral Ecotone
Furthermore, as an important ecological boundary, any small movement of the 400 mm rainfall contour will bring about great changes in an ecosystem [57]. In this study, we find that precipitation will be probably diminished in the future in the agro-pastoral ecotone. Additionally, the temperatures of the three regions of the ecotone have obviously increasing trends, which indicates that the ecotone will have to face a warm and dry climate, a more severe situation, in the future. Hydraulic erosion and wind erosion coexist in the agro-pastoral ecotone. These two types of soil erosion promote mutually, interlacing on the temporal scale and superimposing on the spatial scale. Therefore, the erosional energy of the agro-pastoral ecotone is significantly higher than that of the region of mainly hydraulic erosion (agricultural zone) and that of the region of mainly wind erosion (pastoral area) [22]. In the spring, when the land surface is thawed with sparse vegetation, sand movement is extremely serious, as there is faster soil water evaporation in dry years. It has been proven that the agro-pastoral ecotone is a source of sandstorms in many areas, especially in Beijing-Tianjin [25,58]. Based on the spatial distribution of the soil texture, the soil in the agro-pastoral ecotone is mainly sandy soil [20,59,60]. A low coverage, weak soil capacity and low water conservation capacity characterize the fragile surface structure, so that even minor disturbances can cause damage to the ecosystem. Moreover, the annual precipitation in the agro-pastoral ecotone has high variability. The annual precipitation can be as high as 500-600 mm, however it can also be less than 200 mm [61]. Precipitation is generally concentrated in the summer, with 60-70% of the yearly amount occurring in June-August, when the area is prone to rainstorms and it has extremely low water-use efficiency. Groundwater is continuously pumped for irrigation, creating funnel regions and shrinking wetland [62]. Natural conditions are inappropriate for farming in the agro-pastoral ecotone.
Based on the above, we conclude that the agro-pastoral ecotone is affected by a combination of natural factors and human interference, with the former determining its potential ecological fragility and the latter resulting in a truly fragile ecosystem.

Ecological Protection Policies Implemented in the Agro-Pastoral Ecotone
Long-term eco-restoration projects are critical in sustainable ecological and economic development [31,63,64]. Since the end of the 1990s, China has initiated specific ecological protection polices, such as the Grain for Green Project, the Beijing-Tianjin Dust Storms Sources Control Project, Grassland Ecological Protection Subsidy Incentives, and the Grain for Forage Project, with remarkable effects. Through the moving weighted trend surface analysis and PLS-DA that were conducted in this study, we found that the order for the ESV losses was 1990-2000 > 1980-1990 > 2000-2015 in the agro-pastoral ecotone. Based on these policies, it may be concluded that proper strategies can restore ecosystem functions and improve ecosystem services in limited resources.

Validity and Limitations of This Study
Due to the difficulty in estimating the value from the non-marketed components of ecosystem services, the existing socioeconomic statistical system does not fully account for ESV. Although many previous studies [9,35,[65][66][67][68] have evaluated the valuation of ecosystem services, it is still difficult to identify, quantify, and monetarize ESV. There is still no unified and complete set of scientific assessment methods for ESV. For most of the studies, the obvious barrier is a lack of relevant data in the assessment of ESV, which leads to a rough evaluation result being obtained [69][70][71]. In some small areas, the data that were collected from field observation and survey were able to meet their requirements; however, for large areas, if we still used the same methods to acquire data, except for longer periods and higher costs, large-scale models might abridge some indicators due to the shortage of data, which often caused the final results to only show the average level of the whole changing trend [70]. Reviewing the research literatures on ecosystem services, Martínez-Harms and Balvanera (2012) indicated that the studies of ESV mainly focused on the regional-and country-scale, and the service of regulation was the most studied type [71]. When compared to original data derived from field surveys and experimental observations, secondary data, such as remote sensing data and social economic data, were the more frequently used.
The method of equivalence factor per unit area used in this study is mainly applied to decision-making when estimating ESV at broad regional or national scales [15], as it allows for a quick assessment. More micro-level mechanisms were not studied due to a lack of reliable methods and relevant data. To verify reliability, we used another method, the spatiotemporal dynamic equivalence factor method, to assess the ESV losses in 2010-2015. The dynamic assessment method was consistent with our results. Since our method reflects the large-scale average level compared with the spatiotemporal dynamic equivalence factor method, the inclined first order trend surface from northeast to southwest is underestimated. Given the substantial uncertainties that are involved, we may never have an accurate evaluation of ESV losses [35]. However, this study may be a starting point in emphasizing the potential threat of continuing to squander ecosystem services.

Conclusions
Satellite-based data can substantially assist in the monitoring of large-area landscape changes and in analyzing the effects of landscape changes on ESV losses. In this study, the spatiotemporal process of land-use conversions and their impacts on ecosystem services were investigated in three sections of the agro-pastoral ecotone of Northern China. Significant expansion of cultivated land and reduction in grassland area characterized the trend of land conversion during 1980-2000, whereas the trend was effectively controlled during 2000-2015. When compared with the middle section and the northwestern section of the agro-pastoral ecotone, the irrational land use patterns caused more ESV losses in the northeastern section. In terms of ecosystem service provisions, the agro-pastoral ecotone can be called an "inclined surface" in the direction of northeast to southwest. We believe that human beings play a decisive role in the unique agro-pastoral ecotone of Northern China. Based on the assessment of ecosystem services in the northeastern section, middle section, and northwestern section, long-term restoration projects are crucial in sustainable development in the agro-pastoral ecotone. However, with the reduction in precipitation and increasing temperatures, it should be noted that the risk of ESV losses is still severe. The process of land conversion in the agro-pastoral ecotone requires sustained attention. This study provides some preliminary findings, which could be helpful in improving the understanding of ESV change in response to regional land-use conversions in the agro-pastoral ecotone. Due to the complexity of sustainable land management, it may be a starting point, and more needs to be done in the future.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A. Moving Weighted Trend Surface Analysis
The trend surface analysis, based on the least sum of squares, is a multivariate polynomial regression analysis that can be used to fit mathematical surfaces of increasing complexity to a variable distribution on a map. The expressions of the polynomial trend surfaces are Z i (x i , y i ) = a 0 + a 1 x i + a 2 y i + a 3 x i y i + a 4 x 2 i + a 5 y 2 i + · · · a p y n i (A2) where Z i (x i , y i ),Ẑ i (x i , y i ), and e i are the observed, trend, and residual values of variable Z at coordinate (x i , y i ), respectively, and a 0 to a p are the polynomial coefficients. The variables in Equation (A2) are written in matrix form: Based on the least-squares theory, the residual sum of squares (Q) between the observed value Z i and the trend valueẐ i should be minimized to achieve a "best fit". Thus: The matrix is . .
where p i is the weight, S i is the distance between two points, and R is the sampling circle radius.
For the terrain in this study, we chose Equation (A9). The matrix is: In this study, we used the goodness of fit (R 2 ) and F test with a 95% confidence level to evaluate different orders of polynomial regression equations to select the most suitable moving weighted trend surface model.

Appendix B. Mann-Kendall Trend Test
The Mann-Kendall test is used to detect trends that are monotonic but not necessarily linear. The null hypothesis is independent and randomly ordered data in the Mann-Kendall test. The Mann-Kendall test statistic S is computed, as follows: where x j and x i are the values in years j and i, j > i, respectively, and n is the number of data points. The value of sgn(x j − x i ) is computed, as follows: The test is conducted using a normal approximation, Z statistics, with the mean and the variance, as follows: E(S) = 0 (A14) Var(S) = 1 18 n(n − 1)(2n + 5) − ∑ q p=1 t p (t p − 1)(2t p + 5) where q is the number of tied (zero difference between compared values) groups, and t p is the number of data values in the p th group. The statistics S and Var(S) are used to compute Z, as follows: The Z value is used to evaluate the statistically significant trend, which is a positive value indicating an upward trend or a negative value indicating a downward trend, at the given significance levels.
of the factors over time. The major and minor axes of the ellipses indicate the directions and ranges of the data distributions. Azimuths convey main trend directions, and changes in azimuth can reflect spatial rotations of overall factors in a certain direction. A median center, the center of all the data, identifies the location that minimizes the overall Euclidean distance to the features in a dataset and coincides with the arithmetic mean center under small variability. Herein, we used this method to trace the changes in spatial patterns of populations. Due to data availability limitations, the population data are only from 1990-2015. The algorithm is as follows: The standard distances in the x and y directions are: where x i and y i are the coordinates for feature i, {x mc , y mc } represents the mean center for the features, and n is the number of features. The azimuth is calculated as follows: where xx i and yy i are the deviations of the x and y coordinates from the mean center. The standard deviations for the x-axis and y-axis are as follows: F nij = P ij * F n1 or F nij = R ij * F n2 or F nij = S ij * F n3 (A25) where F nij is the dynamic equivalence factor per unit area for ecosystem service n of an ecosystem in region i in month j; P ij is the spatiotemporal regulation factor of NPP in region i in month j; R ij is the spatiotemporal regulation factor of precipitation in region i in month j; S ij is the spatiotemporal regulation factor of soil conservation in region i in month j; F n1 is the equivalence coefficient for the services from food production, raw materials supply, gas regulation, climate regulation, waste regulation, biodiversity maintenance and aesthetic landscape provision; F n2 is the equivalence coefficient for the service from hydrology regulation; F n3 is the equivalence coefficient for the service from soil conservation; B ij is the NPP of the ecosystem in region i in month j; B is the national annual NPP of the ecosystem; W ij is the average precipitation per unit area in region i in month j; W is the annual average precipitation per unit area; E ij is the quantity of soil conservation of the ecosystem in region i in month j; and, E is the average annual quantity of soil conservation per unit area.