Enhancing Ecosystem Services in the Agro-Pastoral Transitional Zone Based on Local Sustainable Management: Insights from Duolun County in Northern China

.


Introduction
The world is facing various crises, such as climate change and loss of biodiversity, which seriously threaten the survival and development of human beings [1,2].In response to global challenges such as climate change and environmental degradation, the United Nations General Assembly declared 2021-2030 to be the Decade on Ecosystem restoration, dedicated to promoting and restoring ecosystem services (ESs) [3].Nature-based solutions (NbS) also refer to addressing a wide range of human challenges by protecting and sustainably using the vital ESs provided by natural ecosystems [4,5].However, rapid population growth and the overconsumption of natural resources have led to a deterioration in the capacity of ecosystems to provide ESs, and sustainable development is under serious threat [6][7][8].Therefore, it is imperative to monitor and evaluate the spatial and temporal changes in ESs and explore the driving mechanisms behind the observed changes to help policymakers formulate effective policies for managing ESs [9][10][11][12].
There is abundant evidence that the generation of ESs depends on socio-biophysical factors and ecological processes, and that ESs have scale dependence in space and time [1,13,14].It is necessary to quantify and map the ESs at broad temporal and spatial scales to help determine restoration priorities and sustainable development management of ecosystem services [1,15,16].Multiple methods have been developed and deployed to quantify ESs [17][18][19][20][21].With advances in modeling techniques, the use of ecological models to assess and spatially present ESs can better relate findings to ecosystem structures and functions, which can also provide more quantitative evidence for decision-making [8,20,21].Despite rapid progress in quantifying and mapping ESs, there are still some fundamental issues that have not been adequately addressed.One of these issues is that existing studies have mostly focused on quantifying ESs at a single point in time or over several time intervals [22,23].The results of these studies did not allow for general conclusions and are highly likely to mask some uncertainties due to extreme environmental disturbances or human activities during the year [5,24], especially in ecosystems where the landscape is characterized by rapid change and heterogeneity.To successfully manage natural resources and the related ESs, research aimed at long time series perspectives may provide deeper insights into ES changes and underlying ecological and anthropogenic drivers than time-point analyses [25,26].
The effectiveness of ES management is affected by many biophysical and anthropogenic drivers, such as terrain diversity, soil type, climate change, and land use/cover [6,13,24,[27][28][29].Different drivers affect the supply of ESs in time and space in different ways.For example, climate change affects the supply of ESs by influencing biomes through temperature and precipitation [28].Natural background conditions (e.g., topography, soils, and other drivers) directly influence regional ecosystem structures, resulting in the spatial heterogeneity of ESs [30].Whether the impact of these drivers is positive or negative depends on the ES, and the magnitude of these impacts varies from place to place [2,16,31,32].Identifying the impacts of the drivers on ESs on a time scale can help government departments optimize local response strategies in the context of climate change and intense human activity [6,33,34].Furthermore, identifying the drivers of spatial heterogeneity of the ESs has facilitated the implementation of regional ecosystem planning and decision-making [23,35].However, most studies have only analyzed the drivers of the temporal changes or spatial changes in ESs [28,32,33], and the results could not provide more detailed information for the formulation of ecological restoration policies.More efforts are needed to deeply explore the spatial and temporal changes and driving mechanisms of ESs in combination with time-series data to manage and improve multiple ESs by manipulating the drivers.
We focused on the agro-pastoral transitional zone of northern China (APTZNC), which divides agricultural areas from pastoral areas [36] and is situated in the transition from semiarid to arid regions [37].APTZNC includes highly diverse ecosystems (e.g., farmland ecosystems, grassland ecosystems, and woodland ecosystems) and supports the provision of rich ESs (e.g., food supply, grass production, carbon sequestration, and soil and water conservation) [38].To date, studies of the APTZNC have provided valuable information on ESs at large spatial scales, such as cities [39], urban agglomerations [36], and even the entire APTZNC [38,40].However, studies that have identified drivers of change in multiple ESs at local scales (e.g., county administrative districts) are scarce.There are significant differences in precipitation, temperature, biological communities, and human activities within the APTZNC [37], and ESs are prone to change at temporal and spatial scales [40].It is difficult to apply the results of watershed or urban-scale studies to ecological restoration management at county-level administrative units (the most basic administrative unit), which thus may not be effective in improving the supply capacity of ESs [41].Within the county administrative districts, combining the field research situation and ES assessment results can provide a more in-depth and accurate analysis of the ecological status of the region and detailed information for site-specific ecosystem management.On the other hand, analyzing the supply and drivers of ESs in different ecological production zones (e.g., the agricultural advantage zone, the pastoral advantage zone, the semiagricultural and semipastoral zone) is essential to balance the regional ecological protection and economic development.Therefore, it is necessary to carry out county-scale research regarding ESs and their drivers in the APTZNC to serve local policy decisions.
Combined with many years of field research, we focused on Duolun County, a typical county administrative district located in the APTZNC.The ultimate goal of our study was to identify the drivers of spatial and temporal changes in ESs at the local scale based on time-series data.Specifically, we aimed to understand the following: (1) What are the changing trends in ESs from 2000 to 2017?; (2) How do drivers influence the supply of ESs over time?; and (3) What are the main drivers influencing the spatial distribution of ESs in county administrative districts and different ecological production zones?Our intent was that our results could contribute to an improved understanding of the key anthropogenic and biophysical processes underlying the supply of the ESs in a county administrative district of the APTZNC and provide scientific support for promoting sustainable development management in each ecological production zone.

Study Area
Duolun County (41° 46′-42° 36′ N, 115° 51′-116° 54′ E), which covers 3 towns, 2 townships, and 65 administrative villages, is located in the middle of the APTZNC and covers a total area of 3.95 × 10 3 km 2 (Figure 1b).The topography of the study area is semi-annular basin-like, with an elevation range from 1149 to 1796 m (Figure 1c).The soil types can be divided into 7 soil types and 14 subtypes, and the main soil types are chestnut soil, aeolian sand soil, and meadow soil.Duolun County is in a typical agro-pastoral transitional zone.In August 2017, we visited Duolun County for field investigation, and reviewed and verified the vegetation types, land use types and soil types at the sampling sites (Figure 1c).Through investigating the ecological environment, we confirmed that the actual soil types and the land use types of the sampling sites were roughly the same as the soil data used in this study and the land use/cover data of 2015.We combined the administrative boundaries of townships and regional ecological production patterns to divide Duolun County into three regions (Figure 1c).There are large grasslands in the north and east of Duolun County, which are dominated by animal husbandry (Figure 1(cII)); the dominant plants in this grassland mainly include Leymus chinensis, Agropyron cristatum, Stipa krylovii, Cleistogenes squarrosa, Congsheng grass and Artemisia frigida.The south of the study area is an important grain-producing area; farmland mainly planting annual flax, oats, buckwheat, spring wheat and also some silage corn (Figure 1(cI)), The central part of the study area is an ecotone between agriculture and animal husbandry (Figure 1(cIII)), with a similar proportion of the development of agriculture and animal husbandry.In the 21st century, a series of ecological projects have been implemented in this area, such as the Grain for Green Project, grazing prohibition projects and Beijing-Tianjin sandstorm control engineering, to mitigate environmental pressure [41,42].The field research found that there was still land desertification around the county town and in the northern area.Duolun County has a temperate continental arid climate with annual precipitation of 344-399 mm and a range of mean annual temperatures of 3.1 to 4.5 °C.From 2000 to 2017, the climate conditions and vegetation cover in Duolun County underwent obvious changes (Figure 2).The annual precipitation and annual mean temperature increased insignificantly at a rate of 4.526 mm yr −1 and 0.013 C yr −1 , respectively, whereas the annual mean wind speed decreased insignificantly at a rate of 0.008 m s −1 yr −1 .This indicated that the climate in the study area is becoming warm and humid.In addition, the annual vegetation cover increased significantly at a rate of 0.614% yr −1 .The vegetation cover conditions in the study area have been improving.

Methodological Framework and Data Sources
After the field investigation in 2017, we selected and evaluated six key ESs of high relevance to stakeholders in the region for assessment and analysis, including food supply (FS), grass production (GP), water yield (WY), carbon sequestration (CS), soil water erosion (SLwater) and soil wind erosion (SLwind).Then, combined with the time-series data, the spatial and temporal changes in the ESs were analyzed in depth.Finally, the spatiotemporal change drivers of the ESs from the perspective of the zones were analyzed.We selected four major categories of eight drivers to research the causes of the spatial and temporal changes in the ESs.The selected drivers were soil type (ST), meteorological drivers (annual precipitation (Pre), annual mean temperature (Tem), and annual mean wind speed (WS), whereas topographic drivers were altitude (Alt), slope (Slo)), and the anthropogenic drivers were vegetation cover (FVC), land use/cover (LULC).Table 1 shows the datasets used in this study.All raster maps were converted to the UTM coordinate system at a spatial resolution of 100 m.

Food Supply
According to the positive NDVI, the yields of crops such as grains, oilseeds, and vegetables, as well as the yields of meat and milk in the statistical yearbook of the corresponding year, were allocated to the cultivated land and grassland grids, respectively.The formula can be expressed as follows [29]: where  represents the crop yield or meat and milk yield (kg/hm 2 ) of grid i,  is the crop yield or meat and milk yield (kg/hm 2 ) in the county,  is the NDVI value of farmland grid i or grassland grid i, and  is the sum NDVI of the farmland or grassland in the county.

Grass Production
In this study, an inversion model of regional grass yield was used to estimate the grass yield per unit area.The calculation formula is as follows [43]: where  represents the grass yield (kg/hm 2 ) of grid i and  is the NDVI value of grassland grid i.

Water Yield
In this study, water yield was simulated using the Integrated Valuation of ESs and Trade-offs model (InVEST, version 3.3.0),which was developed by Stanford University, The Nature Conservancy (TNC) and the World Wildlife Fund (WWF).The model is expressed as follows [17]: where Yx is the annual water yield (mm) for grid point x, AETxj is the annual actual evapotranspiration (mm) for pixel x, Px is the annual precipitation (mm) in pixel x, Wx is a nonphysical parameter that characterizes the natural climatic-soil properties, PETx is the potential evapotranspiration (mm), AWCx is the volumetric plant-available water content (mm), and Z is an empirical constant which is used to characterize the seasonal distribution of precipitation [41].

Carbon Sequestration
In this study, NPP was used as a proxy for carbon sequestration and estimated using the Carnegie-Ames-Stanford Approach (CASA) model [44].The model is expressed as follows: ,  =  ,  ×  ,  × 0.5 where NPP(x,t) is the net primary productivity of pixel x in month t (gC/m 2 ), APAR(x,t) is the photosynthetically active radiation absorbed by pixel x in month t (MJ/m 2 ), ε(x,t) is the light utility efficiency of pixel x in month t (gC/MJ), SOL(x,t) is the total solar radiation on pixel x in month t (MJ/m 2 ), and FPAR(x,t) is the ratio of the absorption of the incoming photosynthetically active radiation by the vegetation layer (dimensionless) [41].The constant 0.5 reflects the proportion at which the effective solar radiation accounted for the total solar radiation, εmax is the maximum light use efficiency of the vegetation (gC/MJ), whereas Tε1(x,t), Tε2(x,t), and Wε(x,t) refer to parameters describing the stress coefficients at the highest temperature, the stress coefficients at the lowest temperature and the water stress coefficient in cell x in month t, respectively [45].

Soil Water Erosion
The soil water erosion was calculated using the revised universal soil loss equation (RUSLE) [46].The formula is expressed as follows: where SLwater represents soil erosion water (t/(hm 2 )), R is the rainfall erosion factor (MJ is the vegetation cover index (dimensionless), P is the soil erosion control practice factor (dimensionless), and LS is the slope length and slope gradient factor (dimensionless).Detailed information on the parameter localization of RUSLE is presented in Supplementary Table S1.

Soil Wind Erosion
The soil wind erosion was calculated using the revised wind erosion equation (RWEQ) model [47].The basic equations are as follows: where SLwind is soil wind erosion under the conditions of the ground cover vegetation (kg/m 2 ), z is the distance from the upwind edge of a field (m), Qmax is the maximum transport (kg/m), S is the critical field length (m), K′ is the surface roughness, WF is the climatic factor (kg/m), EF is the soil erodibility factor (%), SCF is the soil crust factor, and COG is the combined vegetation factor.Detailed information on the parameter localization of the RWEQ is presented in Supplementary Table S2.

Trend Analysis of the Time Series Data
The Theil-Sen median (Sen) estimation and Mann-Kendall nonparametric test [48,49] were used to detect the statistically significant changing trends and trend slopes in the long time series ESs at a significance level of 0.05.Sen's slope indicates the direction and amplitude of the variables' changes with time, whereas the positive and negative slopes indicate increasing and decreasing trends, respectively.The Mann-Kendall test does not require the samples to have normal distributions and is less sensitive to outliers, which is extensively employed to analyze long time series data [33].The calculation was performed with MATLAB R2020b programming.

Drivers of the Temporal Changes in ESs
Correlation analysis and partial correlation analysis were used to detect the relationship between ESs and drivers (Pre, Tem, WS and FVC).When multiple drivers are related to ESs at the same time, the use of partial correlation analysis allows removal of the influence of the remaining drivers, and the relationship between a single driver and ESs is analyzed separately.The calculation was performed with R statistical software and OriginPro 2022 software [50].In addition, the annual average values of each ES in individual land cover categories (only unchanged land cover) in Duolun County and different ecological production zones were calculated to compare the ES supply capacity of the major land cover (farmland, woodland, and grassland).The calculations were performed with ArcGIS 10.5 software and OriginPro 2022 software [51].

Drivers of the Temporal Changes in ESs
A geographical detector is a statistical method used to detect the spatially stratified heterogeneity of geographic phenomena and reveal nonlinear associations between potential drivers and geographic phenomena [52].A geographical detector assumes that if there is spatial consistency between independent variable X and dependent variable Y, then a statistical association is present between them.The advantages of this method are that it does not need a linear hypothesis, and its physical meaning is clear.It contains four formulas: a factor detector, an interaction detector, a risk detector and an ecological detector.In this study, the factor detector module was selected to evaluate the explanatory power of the independent variable X to the dependent variable Y.The independent variable X was assigned to the eight drivers, and the dependent variable Y was assigned to the seven types of ES supply changes.The formula to measure the q value is as follows: where q is the explanatory power of variable X to the spatial variation in variable Y.The q value ranges from 0 to 1, and the value means that X explains q × 100% of Y.The larger (or smaller) the q-statistic is, the stronger (or weaker) the explanatory power of the independent variable to explain the dependent variable.L is the number of classifications or partitions of X; N and Nh are the numbers of units in the entire study area and subregion h, respectively. 2 and  ℎ 2 are the variances of variable Y over the entire study area and subregion h, respectively.The input data of geographic detectors must be in the form of categorical layers (such as soil type and land use type); therefore, continuous datasets (such as precipitation and temperature) must be categorized.In this study, the precipitation, temperature, wind speed, vegetation coverage, elevation, and slope were divided into six strata by the natural break method.

Temporal Changes in Ecosystem Services
During 2000-2015, the six ESs showed similar trends in annual mean values in the county administrative district, agro-zone, pastoral zone and agro-pastoral zone (Figure 3).The amount of FS exhibited an insignificant upward trend (p > 0.05), with the largest slope in the agro-zone (slope = 7.1 kg m −2 yr −1 ) and the smallest slope in the agro-pastoral zone (slope = 5.025 kg m −2 yr −1 ).There was a trend of significant increase (p < 0.05) in GP both in the county administrative district and different ecological production zones.Notably, there are strong interannual fluctuations in CS and WY with a non-significant increasing trend (p > 0.05), with the largest increase rate in WY in the agro-zone (slope = 2.007 mm yr −1 ) and CS in the pastoral zone (slope = 1.623 gC m −2 yr −1 ).Generally, SLwater and SLwind showed a downward trend, especially in the agro-pastoral zone, and the rate of decrease reached −0.3 t hm −2 yr −1 (p = 0.068) and -0.016 kg m −2 yr −1 (p < 0.05), respectively.Combined with the results of the pixel-by-pixel trend analysis (Figure 4), we found that the trends in the six ESs were spatially different.Except for the two negative ESs (SLwater and SLwind), all four ESs showed a good trend in that the gain areas were larger than the loss areas.Interestingly, the decrease in WY (6.2%) was concentrated in the central part of this study area, where GP and CS increased significantly.The increase in CS was mostly located in the pastoral zone, whereas the losses were mainly located in the agro-zone and the agro-pastoral zone.SLwater and SLwind showed good trends in that the areas of improvement were much larger than the areas of their gain (Figure 4), with the areas of increases in SLwater scattered throughout the county.

Spatial Changes in Ecosystem Services
The spatial patterns of the ESs in Duolun County exhibited regional heterogeneity (Figure 5), and the spatial distribution patterns of each ES remained stable (Figures S1-S6).Taking 2017 as an example (Figure 5), FS shows a decreasing trend from the agrozone to the pastoral zone (Figure 5a).In contrast, GP and CS show an opposite pattern compared with FS, i.e., the high supply was concentrated in the northern and eastern areas of the pastoral zone (Figure 5b,d).The spatial distribution of WY patterns changed slightly in different years, but most of the high-value zones were distributed in the agro-pastoral zone, with a few in the northern and eastern parts of the pastoral zone (Figure 5c).SLwater showed a decreasing characteristic from southwest to northeast, with most of the higher areas concentrated in the agro-zone and the eastern part of the pastoral zone, whereas SLwater was less concentrated in the agro-pastoral zone (Figure 5e).Although the SLwind has decreased significantly over the past 18 years (Figure S6), the relatively higher areas can still be found mainly in the northern pastoral zone of Duolun County, and the lower areas are more stably distributed in the southern and eastern regions (Figure 5f).The results showed that (1) Pre, Tem and FVC exhibited a positive correlation with FS and GP, whereas WS exhibited a negative correlation with FS and GP; (2) FVC and Pre were the two drivers that positively affected the changes in WY, whereas Tem and WS showed insignificant inhibitory effects on WY; (3) Pre, Tem and FVC had a positive effect on CS, whereas WS had an inhibitory effect on CS; (4) Pre, WS and FVC mainly showed an insignificant promotion effect on SLwater, whereas Tem showed an inhibitory effect on SLwater; (5) FVC and Pre had suppressive effects on SLwind, whereas Tem and WS had facilitating effects on SLwind.Combined with the results of partial correlation analysis (Figure 7): (1) FVC was significantly and positively correlated with FS and GP, whereas WS showed a negative correlation with FS and GP; (2) the correlation between Pre and WY was still significant and positive, but the correlation between FVC and WY became significantly negative from significantly positive; (3) the relationship between CS and drivers did not change, but the driver with the greatest correlation between CS in agro-zone, pastoral zone, and agropastoral zone was different, namely, Pre, WS, and WS; (4) FVC and SLwater changed from a positive correlation to a negative correlation, whereas Pre and SLwater were still positively correlated; (5) Tem and WS were significantly and positively correlated with SLwind, whereas Pre and FVC were still negatively correlated with soil wind erosion.From 2000 to 2017, the sum of farmland, grassland and woodland in Duolun County exceeded 80% of the total area of the study area; thus, this paper analyzed the changes in the annual average ES supplies of the three key land use/cover types (only unchanged land cover) over the past 18 years (Figure 8).The results indicate that ES supplies differ by land-use type.Specifically: (1) The FS of farmland fluctuated widely from year to year, with a multiyear average of approximately 491.35 kg/hm 2 in Duolun County.In contrast, the FS of grassland was more stable, with a multiyear average of approximately 69.1 kg/hm 2 ; (2) All the GP in this study was provided by grassland.The multiyear average value of GP in Duolun County was 1633.4 kg/hm 2 and the highest GP (1672 kg/hm 2 ) was in the pastoral zone; (3) There were significant differences in the effects of land-use types on the WY, specifically farmland > grassland > woodland.The differences in the supply of WY by similar land-use types in different zones were not significant; (4) The interannual fluctuations of the CS in the three land-use types were smaller and show the characteristics of woodland > grassland > farmland, with multiyear average values of 502.6 gC/m 2 , 392.6 gC/m 2 , and 347.5 gC/m 2 in Duolun County, respectively.The supply of CS on woodland in the agro-zone (415.62 gC/m 2 ) was significantly lower than in other zones; (5) The characteristics of SLwater in the three land-use types are the same as those of CS, i.e., the interannual fluctuations were smaller and woodland > grassland > farmland.The highest SLwater was found in woodland in the agro-zone.(6) In terms of the multiyear average values in the county, the highest SLwind erosion (0.1 kg/m 2 ) was found in grassland, followed by farmland and woodland.Among the three zones, the highest SLwind was found in the agro-pastoral zone, followed by the pastoral zone.and the semiagricultural and semipastoral zone (agro-pastoral zone), respectively.The first to sixth rows indicate the distribution of annual mean values of food supply (FS), grass production (GP), water yield (WY), carbon sequestration (CS), soil water erosion (SLwater), and soil wind erosion (SLwind) across different land use cover, respectively.

Drivers of Spatial Changes in Ecosystem Services
The dominant drivers affecting the spatial heterogeneity of the six ESs and their qvalues differed significantly within the county administrative district (Figure 9a).The top three drivers affecting the distributions of FS and GP were the same (i.e., FVC > ST >Tem), with FVC having the largest explanatory power for FS and GP, indicating that the distributions of the two ESs were more influenced by FVC.The largest explanatory power for WY was LULC (q = 0.73), and the following drivers were FVC, WS, ST, Alt, Pre, Tem, and Slo, in order of q.The important drivers for the distribution of CS were FVC and LULC, which explained more than 49.77% and more than 23.2% of the distribution, respectively, whereas the other drivers had weaker explanatory power.The slope was the dominant driver determining the distribution of SLwater (q = 0.36), followed by ST (q = 0.136) and Alt (q = 0.129), indicating that soil type and topographic drivers are important environmental drivers affecting SLwater.The FVC had the greatest effect on SLwind (q = 0.15), followed by ST (e.g., sandy soil distribution areas with high SLwind).The explanatory power of meteorological drivers, LULC, and topographic drivers on SLwind did not exceed 5%.In a comprehensive view, the explanatory powers of the four major categories of drivers on the distribution of SLwater are ranked as follows: topography drivers > soil type > anthropogenic drivers > meteorological drivers.The explanatory powers of the four major drivers on the spatial distribution of the other five ESs are as follows: anthropogenic drivers > soil type > meteorological drivers > topographic drivers.The q values of drivers affecting the spatial distributions of the ESs (FS, GP, WY, CS, SLwater, and SLwind) in the agricultural advantage zone (agro-zone).(c) The q values of drivers affecting the spatial distributions of the ESs (FS, GP, WY, CS, SLwater, and SLwind) in the pastoral advantage zone (pastoral zone).(d) The q values of drivers affecting the spatial distributions of the ESs (FS, GP, WY, CS, SLwater, and SLwind) in the semiagricultural and semipastoral zone (agro-pastoral zone).FS: food supply, GP: grass production, WY: water yield, CS: carbon sequestration, SLwater: soil water erosion, SLwind: soil wind erosion.
The dominant drivers for the distributions of ESs in the three zones remained consistent with those of the county, and the importance of the other drivers changed with the ecological production zone.FVC was the dominant driver in the distributions of FS and GP in each ecological production zone.Furthermore, in the agro-zone, Slo (q = 0.136) and Pre (q = 0.104) contributed more to the distributions of FS, and Pre (q = 0.221) and Alt (q = 0.175) contributed more to the distributions of GP.In the pastoral zone, Tem (q = 0.097) and Alt (q = 0.094) had greater explanatory powers for FS, and Tem (q = 0.091) and WS (q = 0.071) had greater explanatory power for GP.The Slo had a greater explanatory power for the distributions of FS and GP in the agro-pastoral zone.LULC was still the dominant driver of WY in different ecological production zones, and the other important influencing drivers were meteorological drivers.FVC and LULC were the two most important drivers affecting the distributions of CS in all three ecological production zones, indicating that anthropogenic drivers had an important role in the distribution of CS.Comparing the three ecological production zones, the explanatory power of the slope on the distribution of SLwater exceeded 30% in all three zones.The explanatory power of Alt for SLwater was higher in the agro-zone and the agro-pastoral zone, 26.7%, and 34.1%, respectively, but only 1.5% for SLwater in the pastoral zone.Different from SLwater, SLwind in the different ecological production zones was influenced less by topographic drivers.The dominant drivers of SLwind in the agro-zone and pastoral zone were FVC (q = 0.388) and FVC (q = 0.21), respectively, whereas the dominant driver of SLwind in the agro-pastoral zone was ST (q = 0.246), followed by FVC (q = 0.175).In addition, ST and WS had greater explanatory power for SLwind in the pastoral zone and agro-pastoral zone, whereas Pre and Tem had greater explanatory power in the agro-zone

The Impacts of Anthropogenic and Meteorological Drivers on the Temporal Variations in the ESs
As a result of the Grain for Green Project and the Beijing-Tianjin Sand Source Control Project (initiated in 2000), vegetation coverage has increased in Duolun County [42], effectively increasing grass production.Although there was a large amount of farmland conversion in the implementation of the ecological project [53], the food supply in Duolun County gradually increased through the gradual development of unused land and the construction of agricultural mechanization demonstration parks.In addition, our results showed that the increase in FVC reduced the probability of SLwater and SLwind (Figure 10).On the one hand, precipitation intercepted by the vegetation canopy can reduce the direct erosion of soil by rainfall, which, in turn, reduces the probability of SLwind [33].On the other hand, the vegetation canopy can reduce wind speed, and the vegetation root system has the function of consolidating soil, which effectively enhances the soil resistance to erosion [32,54].However, it is worth noting that WY is negatively influenced by FVC.Studies on the Loess Plateau [55] also pointed out that extensive vegetation restoration threatens the water supply required for human survival and vegetation growth.The rapid growth of vegetation in Duolun County since 2000 has consumed a large amount of soil water and increased evapotranspiration, thus significantly reducing the WY [32,34].In particular, the study area is located in an arid and semi-arid region, and the actual precipitation and available water resources in recent years need to be taken into consideration when carrying out vegetation restoration projects to minimize the pressure on water resources caused by vegetation restoration.Although vegetation cover is considered to be the main factor in conserving soil from water erosion [54], this protection may be less than the higher SLwater potential due to increased precipitation in the ecologically fragile APTZNC.We noted that increased precipitation can significantly contribute to the occurrence of SLwater.Ecosystem structure is fragile within the APTZNC, and despite revegetation efforts on slope farmland, the damaged soil structure has not been fully restored due to the long history of farming [8,34], and soil erodibility is still high.Despite SLwater in the study area having improved, it is necessary to continue to strengthen the water and fertility retention capacity of the soil through ecological protection in the future, especially on sloping farmland.In addition to precipitation, changes in temperature and wind speed lead to changes in the ecological environment and ES supplies.A warmer climate can increase evaporation from the surface, with drier soil surfaces and greater susceptibility to wind erosion [56].As our results show, there is a significant positive correlation between temperature and SLwind.However, the mean annual temperature in Duolun County is low, and warmer temperatures can offset some of the wind erosion by promoting plant growth and increasing surface roughness [6,57].In line with other studies [58], our results show that wind speed is an important driver of the occurrence of SLwind.The increase in wind speed not only tends to disperse the soil, but also accelerates the evaporation of soil water, which inhibits vegetation growth [33] and increases the risk of SLwater.The decrease in wind speed in the study area directly reduces the potential for SLwater and SLwind.
Land use/cover changes are also the greatest pressures affecting the provision of ESs [34,59].Over the last 20 years, the growth of plantation forests in Duolun County has gradually entered semimature and mature stages, and the strong transpiration of the canopy consumes a large amount of water [34,54], resulting in a lower WY of woodland than that of farmland and grassland.In contrast, farmland has a higher capacity of WY due to less evaporation from vegetation, although its CS is lower than that of woodland and grassland [29,60].Grassland has a lower water demand than woodland, and its CS capacity is greater than that of farmland.Increasing the area of grassland in arid and semiarid areas may be a compromise in terms of increasing CS and WY at the same time.SLwater and SLwind erosion were higher in the woodland and grassland in the study area than in the farmland, which is similar to our hypothesis.This is related to the spatial geography of this study area, and our numerous field studies have revealed that most plantation forests and grasses in the agro-zone have been planted in erodible landscapes characterized by relatively steep slopes and poor soils.Simultaneously, previous studies have shown that if landscape patches are disturbed, a patch may have difficulty blocking any erosive action [8].Tree planting in the study area (especially in the agro-zone) is dispersed, leading to the fragmentation of woodland and grassland landscape, which, to some extent, weakens the soil and water conservation capacity of the two land use types.Although SLwater and SLwind were higher in woodland and grassland, the decreasing trends in SLwater and SLwind were higher in both land types than in farmland, which indicates that woodland and grassland are more capable of erosion control, especially woodland.Appropriately increasing woodland area and aggregation degree can improve the benefit of soil and water conservation.

The Impacts of the Eight Drivers on the Spatial Changes in the ESs
Ecological control measures and approaches can be explored in a targeted manner by clarifying the driving characteristics of ES spatial changes [61].Our study shows that anthropogenic drivers (FVC and LULC) have a stronger influence on the spatial distributions of multiple ESs than other environmental factors in county administrative districts and different ecological production zones (Figure 11), indicating that the supplies of the ESs in the APTZNC depend to a large extent on the degree of restoration of ecological engineering and the rational allocation of land resources [62].The regional heterogeneity of ecological production zones means that other drivers (meteorological drivers, topographic drivers, ST) affecting the spatial distributions of the six ESs show significant differences, which are closely related to the human and physical geography of the different regions.For example, topographic drivers (Alt and Slo) have a greater influence on the spatial distribution of CS in the agro-zone.However, in the pastoral zone and agro-pastoral zone, the influence of topographic drivers decreased, and the influence of the meteorological drivers increased.This pattern is formed because the spatial heterogeneity of Slo and Alt in the pastoral zone and agro-pastoral zone is small, and the water and heat conditions required for vegetation growth are almost entirely dependent on nature.In contrast, in the agro-zone, humans intervene and regulate moisture and temperature according to vegetation growth needs, and meteorological drivers such as natural precipitation have relatively little influence on the ecological environment.At the same time, the topographies of the agro-zone are complex (high altitudes and steep slopes) and crops are grown with significant spatial heterogeneity, which can have a significant impact on the FVC and ESs [62].Therefore, meteorological drivers have stronger explanatory power on the spatial distribution of CS in the pastoral zone and agro-pastoral zone, and topographic drivers have more influence on the spatial distributions of the ESs in the agro-zone.In addition, ST has a greater explanatory power than other drivers for SLwind in the agropastoral zone.Although soil types are abundant in Duolun County, they are mostly reflected in the agro-pastoral zone.As the foundation of terrestrial ecosystems, soil provides many important benefits and ESs to society, and changes in the physical properties of soils can affect ecosystem processes, and thus, the ESs [63,64].Therefore, the spatial heterogeneity of soil type should be considered in regional ecosystem management to ensure a precise fit between management measures and soil type, thus contributing to the improvement of the various ESs.
Although the drivers of ES spatial heterogeneity change with the region, the effects of some environmental factors on ESs still show consistent characteristics, and this effectively reduces information redundancy needed for managing multiple ESs simultaneously to some extent.For example, LULC and meteorological drivers are always the dominant drivers affecting the distributions of WY, FVC is the dominant driver for FS, GP, CS and SLwind.More importantly, Slo has greater influences on SLwater than FVC.Soil erosion is more likely to occur in areas with steep slopes, and regional topographic features should be considered when developing management and soil erosion prevention measures.

Implications for Integration of Ecosystem Services in APTZNC Management
As a bridge between natural ecosystems and socio-economic systems, ESs are critical to conserving biodiversity, maintaining ecological security, and meeting human livelihood needs [37,38,40].This study focused on the APTZNC, a transition zone where agricultural and pastoral areas are connected, which plays a critical role in sustaining the stability of natural ecosystems in northern China and provides an important guarantee for the livelihood of local and surrounding residents [37,41].The sustainable management of ESs in the APTZNC is conducive to restoring the ecological environment and enhancing human well-being in the region [38,40].Based on widely used biophysical models and remote sensing data, this study analyzed the drivers of spatiotemporal changes in ESs in the most basic administrative unit (i.e., county), which can provide some clues for ES management studies in other administrative units in the APTZNC.Our research suggests that as compared with the use of time node variations, the continuous change in ESs over a long period provides more meaningful results.Attempts to use discontinuous years as a study period to assess ESs for management purposes must be taken cautiously, and researchers should consider modeling ESs over continuous time periods and finer spatial scales [1,32].In addition, ES management should be tailored to local conditions and zoning.
In the face of ecological engineering and climate change, ecosystems in the agro-zone still exhibit vulnerability.Although FS was highest in the agro-zone, CS in farmland, woodland, and grassland was the lowest of the three zones, and SLwater was significantly higher than in the other zones.Our field survey found that although the implementation of ecological projects has enhanced the restoration of vegetation on slopes, the damage caused by long-term tillage to the ecosystem has not been fully remediated.The study area should continue to strengthen the ecological protection of sloping land and reclassify areas to avoid the fragmentation of ecological woodlands and grasslands [8,65] to give full play to the water and fertility retention capacity of both land types.To promote the synergistic development of FS and the other ESs in the agro-zone, it is recommended to ensure the irrigated area of arable land through water conservancy engineering measures, to improve the water and fertility retention capacity of the soil through deep plowing and deep loosening techniques and to increase the application of organic fertilizers to maintain and enhance the ES supply capacities.

Driving factors
On the other hand, the supply of the six ESs in the agro-pastoral zone was at a low level in the whole county.There are bare soils and abandoned farmlands in the dryland regions of the agro-pastoral zone, which is the county seat, and it is necessary to actively promote vegetation restoration in unused land to improve the various ESs.It is important to pay attention to the characteristics of diversified soil types and restore vegetation in this area.At this time, attention needs to be paid to the recovery of sandy vegetation, which may result in a shortage of water resources in the local and surrounding areas [34].To address this phenomenon, in addition to relying on natural precipitation, there is a need to reasonably exploit groundwater and take advantage of the county's large impervious layer to properly collect surface runoff as reserve water.
Finally, the distributions of the ESs in the non-intensively managed pastoral zone are susceptible to the influence of soil type and meteorological drivers, in addition to human activities.The northern part of the pastoral zone is distributed with large areas of grassland sandy soil, which still needs attention to prevent SLwind.We believe that the northern region should be well protected from the wind and that the expansion of farmland should be prohibited to avoid negative impacts on the other ESs due to deteriorating soil conditions [64].In addition, we recommend that the northern part should increase the degree of grassland aggregation by planting some drought-tolerant or less water-consuming grass species [33], and attention should be given to the negative impact of revegetation on terrestrial water storage (especially in arid and dry areas).The eastern part of the pastoral zone is rich in water and has larger areas of natural forests and grasslands, which are tourist destinations.The region should protect the existing vegetation, improve the quality and stability of the forest ecosystems, and avoid overexploitation.

Limitations and Future Research Directions
Some constraints in our analysis should be considered.First, due to the lack of annual high-resolution land use data, we had to use land use/cover data from one point in time to assume a gradual change in land use over a five-year period, which may have had some modest effect on the WY and some modest effect of the soil P-factor on SLwater.In the future, the use of consecutive years of land use/cover data for interannual WY and SLwater assessments deserves further study.Second, the correlation analysis only characterized the numerical relationships between the ESs and the drivers, and cannot characterize the threshold values of these relationships [54].Combining analytical methods such as multiple regression and constraint lines to investigate the nonlinear relationships and interaction thresholds between driving factors and ecosystem services will be further explored in future research.

Conclusions
This study simulated spatiotemporal changes in six ESs over the last 18 years (2000-2017) and explored the drivers of changes in ESs of the county administrative district.This provides more detailed information for local ecosystem service studies and has the potential to assist decision-making processes.FS, GP, WY, and CS increased over time, and the capacity of the landscape to control water erosion and wind erosion was enhanced.The drivers of temporal changes in a single ES over time showed similarities across zones.The significant increase in FVC since 2000 has improved the FS, GP, CS, SLwater, and SLwind, but at the same time, has put pressure on water resources.Precipitation contributed to the improvement of WY and CS, but increased SLwater.The reduction in wind speed improved the six ESs, and the temperature had a significant promoting effect on SLwind.Exploring the drivers influencing the spatial distribution of the six ESs underscored the importance of anthropogenic drivers for the spatial distribution of ES over other environmental factors.Our findings also suggest the importance of integrating ES management with the ecosystem characteristics of each zone, because the influence of meteorological drivers, topographic drivers, and ST on the spatial distribution of ESs varied in different zones.For example, it is important in the agro-zone and agro-pastoral zone to account for Alt having a greater influence on SLwater, and in the agro-pastoral zone to account for ST being the main driver of SLwind.The study area should continue to strengthen the soil's ability to retain water and fertilizer in sloping fallow areas.The agro-pastoral zone needs to actively promote the revegetation of unused land and strive to improve the six ESs.The overall level of ESs in the pastoral zone is high, but attention is still needed to prevent and control SLwind.

Figure 1 .
Figure 1.Locations in the study area and general description of geographical information: (a) the location of the APTZNC in China; (b) the location of Duolun County in the APTZNC; (c) a digital elevation model (DEM) of the study area; (cI) survey photographs of the agricultural advantage zone (agro-zone); (cII) survey photographs of the pastoral advantage zone (pastoral zone); and (cIII) survey photos of the semiagricultural and semipastoral zone (agro-pastoral zone).

Figure 2 .
Figure 2. Temporal variation in the annual (a) precipitation, (b) mean temperature, (c) mean wind speed and (d) vegetation cover in Duolun County during 2000-2017.The gray areas represent the 95% confidence intervals.

Figure 3 .
Figure 3.The temporal dynamics of the six ESs during 2000-2017.We determined the average annual values of each ES and plotted their overall trends.The first column (a), second column (b), third column (c) and fourth column (d) show the trends in the six ESs in Duolun County, the agricultural advantage zone (agro-zone), the pastoral advantage zone (pastoral zone) and the semiagricultural and semipastoral zone (agro-pastoral zone), respectively.The first to sixth rows indicate the trends of food supply (FS), grass production (GP), water yield (WY), carbon sequestration (CS), soil water erosion (SLwater), and soil wind erosion (SLwind), respectively.The gray areas represent the 95% confidence intervals.

Figure 4 .
Figure 4. Spatial patterns of the six key ESs change trends and significance levels (p < 0.05) from 2000 to 2017.Zone I: the agricultural advantage zone (agro-zone), Zone II: the pastoral advantage zone (pastoral zone), Zone III: the semiagricultural and semipastoral zone (agro-pastoral zone).FS: food supply, GP: grass production, WY: water yield, CS: carbon sequestration, SLwater: soil water erosion, and SLwind: soil wind erosion.

Figure 6
Figure 6 represents the correlation coefficients based on the time series data between ESs and drivers (including Pre, Tem, WS, and precipitation) for Duolun County and different regions.The results of the study showed that the important factors affecting changes in ESs over time were essentially the same in Duolun County and different zones.The results showed that (1) Pre, Tem and FVC exhibited a positive correlation with FS and GP, whereas WS exhibited a negative correlation with FS and GP; (2) FVC and Pre were the two drivers that positively affected the changes in WY, whereas Tem and WS showed insignificant inhibitory effects on WY; (3) Pre, Tem and FVC had a positive effect on CS, whereas WS had an inhibitory effect on CS; (4) Pre, WS and FVC mainly showed an insignificant promotion effect on SLwater, whereas Tem showed an inhibitory effect on SLwater; (5) FVC and Pre had suppressive effects on SLwind, whereas Tem and WS had facilitating effects on SLwind.

Figure 6 .
Figure 6.Correlation analysis between the ESs and the four drivers (vegetation cover (FVC), precipitation (Pre), temperature (Tem), and wind speed (WS)) in different regions from 2000 to 2017.The blue and red circles above the diagonal indicate positive and negative correlations, respectively.The asterisks in the circles show the significance degree ( * for p < 0.05).The numbers below the diagonal indicate Spearman's correlation coefficients, with their color matching those of the corresponding circles.Darker colors demonstrate stronger correlations.FS: food supply, GP: grass production, WY: water yield, CS: carbon sequestration, SLwater: soil water erosion, SLwind: soil wind erosion.agrozone: the agricultural advantage zone, pastoral zone: the pastoral advantage zone, agro-pastoral zone: the semiagricultural and semipastoral zone.

Figure 8 .
Figure 8.The supply of the six ESs in a different land-use cover (farmland, grassland, and woodland) from 2000 to 2017.The first column (a), second column (b), third column (c) and fourth column (d) indicate the distribution of annual mean values of ESs on different land use cover in Duolun County, the agricultural advantage zone (agro-zone), the pastoral advantage zone (pastoral zone)

Figure 9 .
Figure9.The q values of drivers affect the spatial distributions of ESs in Duolun County and different zones.The drivers include precipitation (Pre), temperature (Tem), wind speed (WS), soil type (ST), altitude (Alt), slope (Slo), vegetation cover (FVC), and land use/cover (LULC).(a) The q values of drivers affecting the spatial distributions of the ESs (FS, GP, WY, CS, SLwater, and SLwind) in Duolun County.(b) The q values of drivers affecting the spatial distributions of the ESs (FS, GP, WY, CS, SLwater, and SLwind) in the agricultural advantage zone (agro-zone).(c) The q values of drivers affecting the spatial distributions of the ESs (FS, GP, WY, CS, SLwater, and SLwind) in the pastoral advantage zone (pastoral zone).(d) The q values of drivers affecting the spatial distributions of the ESs (FS, GP, WY, CS, SLwater, and SLwind) in the semiagricultural and semipastoral zone (agro-pastoral zone).FS: food supply, GP: grass production, WY: water yield, CS: carbon sequestration, SLwater: soil water erosion, SLwind: soil wind erosion.

Figure 10 .
Figure 10.Visualization of the correlation between temporal changes in ecosystem services and drivers (precipitation, temperature, wind speed, vegetation cover) from 2000 to 2017.

Table 1 .
Descriptions and sources of the study data.