Integrating Biophysical and Sociocultural Methods for Identifying the Relationships between Ecosystem Services and Land Use Change : Insights from an Oasis Area

Identifying the relationships between ecosystem services (ESs) and land use change is crucial for ES management and sustainable regional development. The Manas region in China has witnessed dramatic reclamation activities in its desert areas that resulted in ecological problems. The changes in eight ESs, including crop production (CP), livestock production (LP), soil conservation (SC), water yield (WY), sand fixation (SF), carbon sequestration (CS), habitat quality (HQ), and nature landscape recreation (NLR), were investigated by using biophysical and questionnaire methods. At the regional scale, provisioning services (i.e., CP and LP) showed some performance improvements, whereas most of the regulating services (i.e., WY, CS, and HQ) along with NLR showed a performance decline. Five ES bundles—Upper Mountain, Foothill, Oasis, Oasis–Desert Transition, and Desert bundle—were identified at the township scale via k-means clustering. From 2000 to 2015, the Oasis bundle sprawled as a result of oasisization, whereas the Oasis–Desert Transition and Foothill bundles decreased. We performed a questionnaire survey and a statistical analysis to identify the causes behind the performance improvement/decline of these ESs and found that the land use changes in the Manas region had a significant impact on these services. More than 50% of the survey respondents identified land use changes as the primary driver of the changes in some ESs (i.e., CP, CS, HQ, and NLR). In the correlation and partial correlation analyses, oasisization was significantly and positively correlated with CP but was negatively correlated with WY, CS, HQ, and NLR. We enhanced the reliability of our conclusions by integrating biophysical and sociocultural methods into our investigation of ES and land use change. In view of the huge losses in regulating and cultural services, the Manas region should limit its desert reclamation activities to control the expansion of its oasis and to improve the quality of its cropland. Our results can help formulate effective ES management and land use decisions in the Manas region or similar areas.


Introduction
Ecosystems provide humans with a diverse range of goods and services that support their survival and development [1,2].To meet the requirements of economic development amid a growing population, humans have changed and disturbed natural ecosystems more rapidly in the past 50 years compared with any other period in human history [3,4].Land use change can directly or indirectly affect ecosystem types, patterns, and processes and subsequently change the available ecosystem services (ESs) [5][6][7].Several studies have identified land use change as a crucial factor that contributes to ES change [3,6,[8][9][10].Specifically, these studies have identified the impact of land use changes on pollination, water yield (WY), soil conservation (SC), carbon sequestration (CS), and cultural services [11][12][13][14][15].However, given the differences among ES types (e.g., provisioning and regulating services) and land use change modes (e.g., grassland cultivation, urbanization, and ecoengineering afforestation), the specific effects (i.e., positive and negative) of land use change on ESs remain inconclusive, and require a careful and systematic analysis [8,10].
Multiple research methods have been applied to analyze the relationships between land use changes and ESs .Many studies have estimated ES responses to land use changes using the approach of economic value.These studies first interpreted land use changes via remote sensing and then assessed the changes in the economic value of ESs [16][17][18][19][20].For example, Costanza et al. [17] found that the global ESs value was reduced by $4.3-$20.2trillion due to land use changes.Experimental methods have been applied to uncover the land use change processes affecting ESs [21][22][23][24][25].For example, Allan et al. [21] evaluated the biodiversity, functional traits, and 14 ESs in 150 plots of grass with different land use intensities and revealed the land use change mechanism that affects ESs via changes in biodiversity.Recently, some studies have used biophysical models to quantify the physical quantity of ESs.These studies then estimated the impact of land use changes to ESs via scenario simulation or correlation analysis [26][27][28][29][30]. Wei et al. [28] quantified the change in CS using the Carnegie-Ames-Stanford approach (CASA) model and found that although afforestation decreased the cropland area in the Northern Shaanxi Loess Plateau, the CS increased because of investments in natural capital (i.e., afforestation, chemical fertilizer, and agricultural machinery power).Except individual ESs, a few studies have found that ES bundles are linked to land use [30,31].A set of consistently associated ESs forms an ES bundle by repeatedly appearing together across space [32,33].Because the formation of ES bundles is associated with natural and humanistic features of landscape, ES bundles can be interpreted as specific landscape units [32].In the past few years, ES bundles have attracted considerable attention [30,31,[33][34][35][36], and have been studied at various scales, including urban-rural ecotones [34,37], metropolitan areas [33], watersheds [30,31,38,39], nations [40,41], and continents [42], which provided valuable information or recommendation for ES and land use management, such as regarding ecological conservation [30] or land sharing strategies [34].Identifying the ES bundles and their dynamics in a specific area can show the comprehensive influence of land use changes on ESs, which could aid to understand the relationships between ES and land use change, especially the impact of land use change on the distribution of ESs.
Biophysical methods are usually applied to analyze the relationships between land use changes and ESs, but the results of most studies in the field do not meet the information requirements of decision makers and have not yet been incorporated into the relevant policy decisions given that these results do not take into account the perspectives of local stakeholders [16,18,20,26,28,30].Given that biophysical methods are unable to capture the perceptions of local residents toward ESs and land use change, multidisciplinary approaches, including biophysical and sociocultural methods, must be applied to understand the relationships between ESs and land use change.Given their capability to quickly collect information on the perceptions of stakeholders, questionnaire surveys have great potential in addressing research gaps [35,43].In sum, we can fully understand the supply and demand sides of ES and obtain sufficient information about the relationship between ESs and land use change by integrating biophysical and sociocultural methods into ES research.
Although recent studies have investigated different land use change modes in various areas, such as urbanization in metropolitan areas [27,29], wetland reclamation in marshes and coastal regions [44,45], ecoengineering afforestation in watersheds [7,30], and agricultural expansion in agricultural country [46], they have largely ignored the impact of land use change on ESs in oases.Oasisization refers to the process of transforming a natural desert into a man-made oasis to meet the requirements of social and economic development [47].Land use changes associated with oasisization mainly involve the conversion of deserts and low-coverage grasslands into croplands.Many studies examined the ecological problems caused by oasisization [28,47,48], such as salinization, shrinking of lakes, and degradation of vegetation, while few works have identified long-term variations in the land use and ESs available in oases [48][49][50].However, very few systematic studies have estimated the impacts of oasisization on ESs.Oases are often vulnerable and easily disturbed by human activities.Therefore, identifying the relationships between ES and land use change in oases can provide valuable information for oasis management.
According to the above literature review, sociocultural methods have been rarely applied in identifying the relationships between ESs and land use change and oases have attracted limited research attention as a study area.We selected the Manas region, which exhibits a typical oasis ecosystem, as our study area.By integrating biophysical and sociocultural methods, we attempt to identify the relationships between ESs and land use change for policy decisions in the Manas area.In recent years, population growth and economic development have resulted in dramatic land use changes in this area, and a large part of its desert has been reclaimed in the Manas area.Such conflict has been observed between ESs and land use also present has implications in the sustainable development of the region and threatens the well-being of local inhabitants.To address such challenges we (i) identified the patterns and dynamics in individual ESs and ES bundles from 2000 to 2015 using biophysical methods, (ii) investigated the perceptions of local residents to ES and land use change by questionnaire, and (iii) analyzed the relationships between ESs and land use.We also recommend some important land use management and ecology policies to promote ES management in oases.

Study Area
The Manas region, which is located in the northern part of Xinjiang Uygur Autonomous Region, China (43 • 27 -45 • 21 N, 85 • 01 -86 • 32 E; Figure 1), is a critical area of the economic belt on the northern slope of the Tianshan Mountains.The Manas region spans a total area of 22,900 km 2 [50].Our study area exhibits a complex terrain, and the elevation ranges between 277 m and 5146 m with an average value of 1148 m.The distribution pattern of mountains and basins has resulted in a distinct regional climate in the study area.The temperature (T) exhibits higher values in the north than in the south, whereas the precipitation (PPT) exhibits a reverse pattern (Figure A1).The annual average T and PPT of the Manas region in 2015 were 4.02 • C and 230.29 mm, respectively.The Manas region comprises four main rivers-Manasi, Bayingou, Jingou, and Taxi Rivers (Figure 1).Therefore, the water resources are relatively rich.The study area includes 47 township-level units that belong to three county-level administrative regions, namely, the Shihezi reclamation area, the Shawan County, and the Manas County.Our study area includes multiracial populations, and the main ethnic groups comprises Han, Uyghur, Kazak, Hui, and Mongolian.The population is approximately one million, and the per capita gross domestic product is approximately CNY 68,100 (i.e., USD 10,150).The urbanization rate in 2015 was less than 50%, and local inhabitants mainly included three categories, citizens who lived in urban areas and most work in enterprise or public institution, farmers who were engaged in agricultural production, and herders who were engaged in animal husbandry.Regional human activities, including agricultural and industrial activities, are mainly concentrated in the oasis area, which belongs to the piedmont plain.A series of ecological and environmental problems, including landscape fragmentation between desert and oasis, salinization inside the oasis, and grassland degradation, have arisen with the agricultural and industrial activities in Manas.All these activities threaten the sustainable development of the oasis.engaged in agricultural production, and herders who were engaged in animal husbandry.Regional human activities, including agricultural and industrial activities, are mainly concentrated in the oasis area, which belongs to the piedmont plain.A series of ecological and environmental problems, including landscape fragmentation between desert and oasis, salinization inside the oasis, and grassland degradation, have arisen with the agricultural and industrial activities in Manas.All these activities threaten the sustainable development of the oasis.

Research Route
By integrating biophysical and sociocultural methods, our study attempts to identify the relationships between ES and land use change and offer policy decisions based on the proposed framework (Figure 2).We selected eight ESs-crop production (CP), livestock production (LP), soil conservation (SC), water yield (WY), sand fixation (SF), carbon sequestration (CS), habitat quality (HQ), and nature landscape recreation (NLR)-following the criteria presented in Section 2.3.2.Based on the research framework, we first evaluated the supply of these eight ESs by applying biophysical methods and then highlighted the changes in each ESs and land use.Second, we identified several ES bundles and their corresponding changes by conducting a clustering analysis.Third, we conduct a questionnaire survey to evaluate the perceptions of stakeholders toward land use change, the importance of different ESs, the changes in these services, and the causes of such changes.These approaches are described in detail in Sections 2.3 and 2.4.Fourth, we perform comprehensive and comparative analyses of our findings to identify the relationships between ESs and land use change.The results of our analysis can aid in the formulation of ES management and land use decisions for the study area.We also perform correlation and partial correlation analyses to further examine the impact of land use change on ESs.

Research Route
By integrating biophysical and sociocultural methods, our study attempts to identify the relationships between ES and land use change and offer policy decisions based on the proposed framework (Figure 2).We selected eight ESs-crop production (CP), livestock production (LP), soil conservation (SC), water yield (WY), sand fixation (SF), carbon sequestration (CS), habitat quality (HQ), and nature landscape recreation (NLR)-following the criteria presented in Section 2.3.2.Based on the research framework, we first evaluated the supply of these eight ESs by applying biophysical methods and then highlighted the changes in each ESs and land use.Second, we identified several ES bundles and their corresponding changes by conducting a clustering analysis.Third, we conduct a questionnaire survey to evaluate the perceptions of stakeholders toward land use change, the importance of different ESs, the changes in these services, and the causes of such changes.These approaches are described in detail in Sections 2.3 and 2.4.Fourth, we perform comprehensive and comparative analyses of our findings to identify the relationships between ESs and land use change.The results of our analysis can aid in the formulation of ES management and land use decisions for the study area.We also perform correlation and partial correlation analyses to further examine the impact of land use change on ESs.

Selection of Key Ecosystem Services
The Millennium Ecosystem Assessment [51] has identified more than 20 types of ES (e.g., food production, climate regulation, soil formation, primary production, and aesthetics) and classified them into provisioning, regulating, supporting, and cultural services.Land use change influences the agriculture, vegetation, water resources, and inhabitants in the study area.The eight ESs presented in Figure 2 reflect the agriculture, water, vegetation, and recreational uses of the natural ecosystem in the Manas region and indicate that these uses are all associated with land use changes.Based on the categorization proposed by the Millennium Ecosystem Assessment, the eight ESs selected for this study can be categorized into provisioning (i.e., CP and LP), regulating (i.e., SC, WY, SF, CS, and HQ), and cultural services (i.e., NLR).We selected these ESs based on the feasibility of their quantification, their importance, and the availability of related data.Crop and livestock production are selected given their importance to the agricultural and animal husbandry development as well as the wellbeing of the local inhabitants of the study area.Regulating services reflect the physical geographical characteristics (e.g., elevation), climate features (e.g., drought), and ecological risks (e.g., soil and wind erosion) of the study area.Specifically, SC and SF were selected given that the altitude difference between mountains and oases can lead to considerable soil erosion and desert winds can lead to severe sand erosion.WY and CS were selected as they can potentially mitigate the risks brought by drought and climate change, both of which are observed in the study area.HQ was chosen as the key ES mainly because Manas is a key area for biodiversity conservation in China.For cultural services, we selected NLR (i.e., forest recreation (FR) and desert recreation (DR)) given that unique natural scenery, including forests and deserts, can provide local residents with recreational services and can be easily influenced by oasisization.2.3.Ecosystem Services Supply: Biophysical Methods

Selection of Key Ecosystem Services
The Millennium Ecosystem Assessment [51] has identified more than 20 types of ES (e.g., food production, climate regulation, soil formation, primary production, and aesthetics) and classified them into provisioning, regulating, supporting, and cultural services.Land use change influences the agriculture, vegetation, water resources, and inhabitants in the study area.The eight ESs presented in Figure 2 reflect the agriculture, water, vegetation, and recreational uses of the natural ecosystem in the Manas region and indicate that these uses are all associated with land use changes.Based on the categorization proposed by the Millennium Ecosystem Assessment, the eight ESs selected for this study can be categorized into provisioning (i.e., CP and LP), regulating (i.e., SC, WY, SF, CS, and HQ), and cultural services (i.e., NLR).We selected these ESs based on the feasibility of their quantification, their importance, and the availability of related data.Crop and livestock production are selected given their importance to the agricultural and animal husbandry development as well as the well-being of the local inhabitants of the study area.Regulating services reflect the physical geographical characteristics (e.g., elevation), climate features (e.g., drought), and ecological risks (e.g., soil and wind erosion) of the study area.Specifically, SC and SF were selected given that the altitude difference between mountains and oases can lead to considerable soil erosion and desert winds can lead to severe sand erosion.WY and CS were selected as they can potentially mitigate the risks brought by drought and climate change, both of which are observed in the study area.HQ was chosen as the key ES mainly because Manas is a key area for biodiversity conservation in China.For cultural services, we selected NLR (i.e., forest recreation (FR) and desert recreation (DR)) given that unique natural scenery, including forests and deserts, can provide local residents with recreational services and can be easily influenced by oasisization.Statistical data on CP and LP were acquired from the relevant township-level statistical yearbooks for 2011 and 2016 [52], which included the production of food crops, cash crops, fruits and livestock inventory.Meteorological data on daily solar radiation (RS), T, PPT, and wind speed were obtained from the National Climatic Bureau [53] and used for the input parameters of ESs evaluation.Ordinary kriging interpolation and multiple regression analysis were combined to generate the raster maps of the meteorological data.Land use datasets for 2000 and 2015 were obtained from the Data Center for Resources and Environmental Sciences (RESDC), Chinese Academy of Sciences [54] and they have a spatial resolution of 30 × 30 m. Land use types used in our study were classified into seven categories: cropland, woodland, grassland, water, urban, desert, and others.Basic geographic data on road network, water system and elevation also were derived from RESDC, Chinese Academy of Sciences [54] and used for the input parameters of HQ or SC evaluation.A series of Moderate Resolution Imaging Spectroradiometer (MODIS) datasets were obtained from NASA's Earth data [55], including MOD13 normalized difference vegetation index (NDVI), MOD15 leaf area index (LAI), MOD16 evapotranspiration (ET), and MOD17 net primary production (NPP).The NDVI and LAI was used to estimate SC, WY, SF, and CS and the annual ET and NPP products were used for comparison with the simulated ET and NPP results.The soil attribute data (i.e., the content of sand, silt, clay, and organic carbon) were obtained from the Cold and Arid Regions Sciences Data Center at Lanzhou [56] and used for the input parameters of SC and SF evaluation.

Mapping of Individual Ecosystem Services Supply
Each ES was assessed by using the indicators and approaches described in Table 2. CP, LP, and NLR were evaluated in township-level units because of limited data, whereas SC, WY, SF, CS, and HQ were evaluated in cell units based on grid data.To maintain consistency at the spatial scale, the supply level of SC, WY, SF, CS, and HQ in township-level units was calculated by using the spatial analysis techniques in ArcGIS.The same software was employed to estimate the supply of SC, WY, SF, CS, and HQ and to obtain the mappings of each ES.(1) Crop production We initially obtained the annual crop output (i.e., grain, cotton, and fruit) at the township level for two points in 2000 and 2015 to map the CP.Then, we calculated the crop output per unit area.CP services improve as the crop output number per unit area increase.
(2) Livestock production LP number at the township level were initially obtained for two points in 2000 and 2015 to map the LP.Then, we calculated the LP number per unit area.LP services improve as the LP number per unit area increase.
(3) Soil conservation SC service reflects the capacity of the ecosystem to avoid soil erosion caused by PPT.The quantification of soil retention level (A c ) has three steps and uses the RUSLE method [57].The first step was to obtain the potential amount of soil erosion (A p ), which was the product of rainfall erosivity (R), soil erodibility (K sl ), and topography (LS).The second step was to obtain the actual amount of soil erosion, which was the product of rainfall erosivity, soil erodibility, topography, vegetation cover (C sl ), and SC measures (P).The last step was to calculate the soil retention level via the difference between the potential and actual amount of soil erosion (A r ).The calculation processes for rainfall erosivity, soil erodibility, topographic, and SC measures can be found in Wischmeier et al. [62], Sharpley and Williams [63], McCool et al. [64,65], and Lufafa et al. [66].The key formulas were as follows.
(4) Water yield WY service reflects the capacity of the ecosystem for hydrological regulation.The water balance equation [58] could be applied in arid areas to quantify the water output.PPT and ET were crucial parameters in the equation.The spatial distribution of PPT was obtained via the interpolation of meteorology stations.ET could be estimated by PET, PPT, and LAI [67].PET in arid areas could be estimated by combining the T and RS of the earth surface [68].The key equations are presented as follows.
WY = PPT − ET (4) (5) Sand fixation SF service reflects the capacity of the ecosystem to avoid sand erosion caused by wind.The calculation processes for binding sand quantity included three steps and used RWEQ [59].The first step was to obtain the potential amount of sand erosion cause by wind, which could be estimated via weather, soil erodibility, soil crust, and surface roughness factors.The second step was to calculate the actual level of sand erosion caused by wind, which could be estimated via weather, soil erodibility, soil crust, surface roughness, and vegetation cover factors.The last step was to obtain the binding sand quantity via the difference between the potential and actual amounts of sand erosion caused by wind.The weather factor included wind speed, PPT, T, snow coverage, and RS.The key equations for the SF amount of were as follows [59].
where SF c is the amount of SF (kg/m 2 ); SF p is the potential amount of sand erosion (kg/m 2 ); SF r is the reality amount of sand erosion (kg/m 2 ); Z is the distance from the upwind edge of the field (m); S p and S r is the critical field length without and with vegetation cover, respectively (m); Q maxp and Q max is the maximum transport capacity without and with vegetation cover, respectively (kg/m); WF reflects the influence of climate condition on sand erosion by wind; EF is soil erodible factor; SCF is soil crusting factor; and K is the surface roughness factor; C s is the vegetation factor.The calculation processes for soil erodibility, soil crust, and surface roughness factor can be found in Ouyang et al. [69].
(6) Carbon sequestration CS service reflects the capacity of the ecosystem to balance carbon and oxygen, as well as regulate global temperature.Only the CS levels of noncrop plants were quantified using the CASA model because the carbon sequestered by crop plants ultimately returns to the atmosphere [60,70].In this approach, solar radiation (Rs), fraction of photosynthetically active radiation (FPAR) absorbed by the vegetation canopy, temperature limiting factors (f 1 ), water-limiting factors (f 2 ), and maximum light use efficiency of vegetation (e max ) were pivotal parameters.FPAR could be calculated by using NDVI, and the detailed calculation processes for the other parameters can be found in Zhu et al. [71].The key equation was as follows Many areas of the Manas region are considered hotspots for biodiversity conservation, and the HQ is critical for the species.InVEST [61] was applied to estimated HQ.This method considers the effects of land use and external disturbance on HQ.The key equations for the calculation of HQ are presented as follows where Q xj is the HQ for pixel x in land use type j, H j is the habitat suitability of land use type j, D xj is the habitat degradation degree for pixel x in land use type j, k is the half-saturation constant and is set equal to half of the maximum habitat degradation degree, z is the default parameter, R is the number of habitat threat factors, P r is the pixel number of threat factor r, w r is the weight of threat factor r, r p is the number of threat factors in each pixel, i rxp is the impact of threat factor r in pixel p to pixel x, and SE jr is the sensitivity of land use type j to threat factor r.
In view of the actual conditions of the Manas region, settlement, farmland, river, roads, and population were selected as the threat factors.The habitat suitability of each land use type and the sensitivity of a land use type to a threat factor were set (Table A1) [61].
(8) Nature landscape recreation The Manas region ecosystem is complex and distinctive, and the cultural services of the ecosystem are mainly manifested as NLR.The biophysical indicators (i.e., coverage rate of forest and desert) were applied to quantify the supply level of NLR [10].At the township level, the mapping process for NLR was conducted in ArcGIS.The main step was t was to calculate the coverage rate of forest and desert area in each township in 2000 and 2015.NLR services increase with the coverage rates of forest and desert.

Identification and Mapping of Ecosystem Service Bundles
After quantification and mapping of individual ESs supply, we identified ES bundles at the township scale using the supply level of the eight ESs.In this study, we used administrative boundaries as clustering unit and the explanations were as follows.Firstly, the supply of some ESs (i.e., CP, LP, and NLR) was only available at township level.Secondly, social processes shaped the production and consumption of ESs and the use of administrative boundaries allowed us to gather ecosystem and social system and identify different types of social-ecological systems on a landscape [33].Thirdly, the ES bundles of administrative unit could provide information at the decision-making level.A few studies have used k-mean clustering [30,31,33], principal component analysis [72,73], self-organizing map [74], and multivariate regression tree [75] to identify ES bundles.The selection of method in bundle studies were mainly based on the structure of ESs data [32].Considering the discrete characteristic of eight ESs data in our study, we mainly employed k-means clustering in the SPSS 20 statistical software to identify the ES bundles at the township level in 2000 and 2015 in the study area.The detailed steps of identification and mapping of ES bundles were as follows.Firstly, principal component analysis was used in identifying the ES variations of each town in 2000 and 2015 to determine the appropriate number of clusters for the eight ESs in 47 townships.The first three principal components explained more than 75% of the ES variation, whereas the first five principal components explained more than 90% of the ES variation.In view of the local characteristics of the Manas region, we classified the 47 townships into five categories using k-means clustering.The mapping of ES bundles was conducted in ArcGIS software.We also analyzed the change of each ES bundle in area and township numbers to identify the impact of land use change to ES.

Questionnaire Survey
Questionnaire survey was applied to capture the perception of local residents to ES and land use change in the study area.In August and September 2017, 900 questionnaires were distributed across the different landscape areas in 19 points (Table S3).The numbers of surveys distributed in the study area were in accordance with population distributions.Questionnaire responses received were examined for completeness based on respondents' information, including age, education, sex, residency period, and occupation, and identified 815 as valid (Table A2).We categorized the respondents into farmers (farming), herders (grazing), and citizens (enterprise or public institution) on the basis of their occupation.The questionnaire contents were presented in the Supplementary Material Text S1 and could be divided into two parts: (a) the awareness of the interviewees on land use change (i.e., oasisization and urbanization), benefits of nature, and importance of ESs, and (b) the cognitions of the interviewees on the improvement or degradation of ESs and the reasons (e.g., land use change and climate change).We explained oasisization, urbanization, and eight ESs to the interviewees in the survey to improve the understanding of the interviewees about the questionnaires.All results of the valid questionnaires could be used in analyzing the awareness of the benefits of nature and the importance of ESs.However, considering age and residency period, only part of the questionnaire results could be applied in analyzing the perceptions of land use change, ES change, and the causes of ES change.Only 733 questionnaires were valid in analyzing the cognitions of oasisization or urbanization and the improvement or degradation of ESs.These respondents knew the historic situation of the study area.

Data Analysis
A correlation analysis was performed to analyze the impact of land use change on each ES.The expansion of croplands and buildings (i.e., urban areas) was identified as the main driving mechanism of land use change in the townships of Manas.Cropland and building coverage rates were used to quantify the relationships between each ES and land use change.Climate change is another important driving factor of changes in ESs [56].Accordingly, we used T and PPT as the relevant variables and performed a partial correlation analysis to analyze their correlation and to exclude the impact of other variables [28].For example, when analyzing the correlation between cropland and ESs, partial correlation analysis excludes the impact of building and climate factors.Correlation analysis was also performed to evaluate the credibility of estimating the ET and NPP.In analyzing ES bundles, spider diagrams were created to compare the supply level of ESs across different bundles [31].
In the questionnaire survey, the perceptions of respondents toward land use change, nature benefits, changes in ESs, and the reasons behind ESs changes were explored.To evaluate the importance of each ES, the answers of the respondents were coded from 0.1 (lowest score) to 1.0 (highest score) as can be seen in the Supplementary Material Text S1.A higher score indicates the greater importance of each ES to these respondents.The importance scores of each ES were calculated by taking the average of the scores given by all respondents.The importance of ESs in different ES bundles was quantified based on the mean importance score assigned by the respondents to each bundle.

Land Use Change
Croplands, grasslands, and deserts were the main land use types in the study area on the basis of the land use maps for 2000 and 2015 (Figure 3).The cropland coverage in 2000 was 19.18% and increased by 3.94% by 2015 (Table A3).The increased cropland areas came mostly from conversion of grasslands and desert areas.Grasslands and deserts dropped by 407 km 2 (1.77% of the total area) and 264 km 2 (1.15% of the total area), respectively.Moreover, the areas of woodlands and water bodies decreased by 0.29% (68 km 2 ), and 0.03% (6 km 2 ), respectively, whereas urban coverage increased by 0.26% (60 km 2 ) (Table A3).The increased urban areas came mostly from conversion of grasslands and croplands.and 264 km 2 (1.15% of the total area), respectively.Moreover, the areas of woodlands and water bodies decreased by 0.29% (68 km 2 ), and 0.03% (6 km 2 ), respectively, whereas urban coverage increased by 0.26% (60 km 2 ) (Table A3).The increased urban areas came mostly from conversion of grasslands and croplands.

Spatiotemporal Changes of Individual Ecosystem Services
Total CP in the Manas region substantially changed and increased by 1.85 times the original value, ranging from 6.27 × 10 5 t in 2000 to 1.79 × 10 6 t in 2015 (Figure A2).The increase in CP was detected in nearly all the townships, except for some townships proximate the central city area (Figure 4a).The growth in CP in most towns ranged between 1.00 × 10 4 t and 6.00 × 10 4 t.The total level of LP increased by 28.60%, ranging from 1.59 × 10 6 t in 2000 to 2.04 × 10 6 t in 2015.Decreases in LP were detected in the south and north areas of the Manas region, whereas remarkable increase in LP was observed in the central region (Figure 4b).
The average value of SC modulus changed from 3253 t/km 2 in 2000 to 1838 t/km 2 in 2015, representing a decline of 76.99%.The spatial change patterns (Figure 4c) showed that a total of 86.13% of the area experienced a decrease in SC, and these areas were distributed in the south and north of the Manas region.The areas where the decline exceeded 1000 t/km 2 were mainly found in the south.The average WY in the study area decreased by 93.17% from 2000 to 2015, and the average value of WY was less than 10 mm in 2015.Spatially, a decrease was observed in 92.77% of the total area, and the decline in WY for most areas ranged from 30 mm to 150 mm (Figure 4d).The average value of SF modulus decreased by 82.72%, changing from 3877 t/km 2 in 2000 to 670 t/km 2 in 2015.The spatial change patterns showed that more than 95.00% of the total study area experienced a decline in SF.The areas where the declines exceeded 3000 t/km 2 were mainly concentrated in the northern (Figure 4e).The average value of CS changed from 61 t C/km 2 in 2000 to 52 t C/km 2 in 2015, representing a decline of 14.75%.Spatially, a decline was detected in 84.21% of the total area (Figure 4f).The average value of the HQ index changed from 0.36 in 2000 to 0.34 in 2015 representing a decline of 5.56% due to increased habitat threats.The spatial change patterns indicated that 83.18% of the total area experienced a decrease in HQ (Figure 4g).NLR decreased by 21.54% due to the decline in forest and desert coverage from 2000 to 2015, respectively.The decline in NLR was mainly distributed in the northern portion of Manas region (Figure 4h).

Spatiotemporal Changes of Individual Ecosystem Services
Total CP in the Manas region substantially changed and increased by 1.85 times the original value, ranging from 6.27 × 10 5 t in 2000 to 1.79 × 10 6 t in 2015 (Figure A2).The increase in CP was detected in nearly all the townships, except for some townships proximate the central city area (Figure 4a).The growth in CP in most towns ranged between 1.00 × 10 4 t and 6.00 × 10 4 t.The total level of LP increased by 28.60%, ranging from 1.59 × 10 6 t in 2000 to 2.04 × 10 6 t in 2015.Decreases in LP were detected in the south and north areas of the Manas region, whereas remarkable increase in LP was observed in the central region (Figure 4b).
The average value of SC modulus changed from 3253 t/km 2 in 2000 to 1838 t/km 2 in 2015, representing a decline of 76.99%.The spatial change patterns (Figure 4c) showed that a total of 86.13% of the area experienced a decrease in SC, and these areas were distributed in the south and north of the Manas region.The areas where the decline exceeded 1000 t/km 2 were mainly found in the south.The average WY in the study area decreased by 93.17% from 2000 to 2015, and the average value of WY was less than 10 mm in 2015.Spatially, a decrease was observed in 92.77% of the total area, and the decline in WY for most areas ranged from 30 mm to 150 mm (Figure 4d).The average value of SF modulus decreased by 82.72%, changing from 3877 t/km 2 in 2000 to 670 t/km 2 in 2015.The spatial change patterns showed that more than 95.00% of the total study area experienced a decline in SF.The areas where the declines exceeded 3000 t/km 2 were mainly concentrated in the northern (Figure 4e).The average value of CS changed from 61 t C/km 2 in 2000 to 52 t C/km 2 in 2015, representing a decline of 14.75%.Spatially, a decline was detected in 84.21% of the total area (Figure 4f).The average value of the HQ index changed from 0.36 in 2000 to 0.34 in 2015 representing a decline of 5.56% due to increased habitat threats.The spatial change patterns indicated that 83.18% of the total area experienced a decrease in HQ (Figure 4g).NLR decreased by 21.54% due to the decline in forest and desert coverage from 2000 to 2015, respectively.The decline in NLR was mainly distributed in the northern portion of Manas region (Figure 4h).The Upper Mountain bundle covered 5909.52 km 2 and accounted for 25.75% of the total area.The area covered by Upper Mountain bundle did not change from 2000 to 2015.The proportion of cultivated land in this area is low.Therefore, the CP service exhibited low levels.We also found relatively strong levels for the SC and WY services (Figure 5c,d) due to high altitude and forest coverage.The Upper Mountain bundle covered 5909.52 km 2 and accounted for 25.75% of the total area.The area covered by Upper Mountain bundle did not change from 2000 to 2015.The proportion of cultivated land in this area is low.Therefore, the CP service exhibited low levels.We also found relatively strong levels for the SC and WY services (Figure 5c,d) due to high altitude and forest coverage.
The area covered by Foothill bundle decreased by 183.66 km 2 , changing from 3627.48 km 2 in 2000 to 3443.82 km 2 in 2015.In response, the proportion of the total study area changed from 15.80% in 2000 to 15.00% in 2015.The ES levels were relatively high for CS, HQ, and LP services because of the high amount of vegetation (e.g., grass) coverage in these types of area (Figure 5c,d).
The Oasis bundle covered 2246.57km 2 in 2000 and 3717.47 km 2 in 2015, thereby increasing by 1470.90 km 2 .In response, the proportion of the total study area changed from 15.80% in 2000 to 15.00% in 2015.The land use types in these areas were mainly croplands, followed by grasslands and buildings.The proportion of forest was low.The provisioning services (i.e., CP and LP) in this bundle were relatively strong (Figure 5c,d).
The area covered by Oasis-Desert Transition bundle decreased by 835.54 km 2 , changing from 5229.00 km 2 in 2000 to 4393.46 km 2 in 2015.Their proportion of the total study area decreased from 22.78% in 2000 to 19.14% in 2015.Cultivated and sandy lands were the main land use types, and windbreaks were planted at the edge of the desert.Therefore, the SF service was high (Figure 5c,d).
The area covered by Desert bundle decreased by 451.71 km 2 , changing from 5930.31 km 2 in 2000 to 5478.60 km 2 in 2015.As a result, the ratio of these areas to the total study area changed from 25.84% in 2000 to 23.86% in 2015.These areas were nearly completely covered by sandy land, and the proportion of cropland was low.Moreover, the annual average PPT was only approximately 120 mm.Given the aforementioned characteristics, the HQ and DR in this bundle were high, whereas the provisioning services (e.g., CP and LP) were low (Figure 5c,d).

Perceptions of Ecosystem Services and Land Use Change
On the basis of the questionnaire results of land use change, 93.59% and 90.86% of the respondents believed that the study area experienced oasisization and urbanization, respectively.Only 32.73% and 25.45% of the herders perceived oasisization and urbanization, respectively, because the herders were mainly distributed in the upstream areas, and land use change was relatively weak.However, the majority of farmers and citizens (more than 90%) recognized cropland and urban expansion.
ESs refer to the benefits to human well-being that originate from nature, so survey results about perceptions of nature benefits were important to the analysis of ESs.Based on the questionnaire results of perceptions of nature benefits, the majority of respondents recognized the different benefits of nature (e.g., CP, LP, WY, SC, SF, and NLR) (Figure A3), which indicated that most of them admitted exists of ESs.Only 1.64% of the interviewees admitted not knowing any benefits of nature, whereas 41.20% believed that the benefits of nature were highly significant.
The scores of the importance of different ESs were calculated by the average value of all the respondents (Figures A4 and 6A).The average score of the importance of provision services was the highest (0.71), followed by regulating services (0.58) and NLR (0.51; Figure A4).For different stakeholders (i.e., farmers, herders and citizens), the difference of income level, place of residence, occupation and education level caused their demand for ESs varied.Therefore, under the influence of stakeholders' characteristics, the scores of the importance of ESs from the different stakeholders varied.Provisioning services were generally viewed as more important by farmers and herders than citizens, whereas NLR was viewed as more important by citizens than farmers and herders (Figure 6A).The average score of the importance of regulating services for the farmers was the highest (0.60), followed by citizens (0.57), and herders (0.46).Farmers in the study area viewed CP, WY, and SF as more important than other services.Herders viewed LP as important, whereas citizens viewed CS and NLR as important (Figure 6A). Figure 6B showed the importance of each ES in different ES bundles.The scores of the importance of CP were notable high (more than 0.65) and was viewed a relatively important ES in all ES bundles.LP and SC was viewed as more important in the Upper Mountain and Foothill bundles than Oasis and Oasis-Desert Transition bundle.The scores of the importance of WY and CS were notable high (more than 0.65) in the Oasis and Oasis-Desert Transition bundle.SF was viewed as more important in the Oasis-Desert Transition bundle, while HQ and NLR were viewed as most important in the Oasis bundle relative to other bundles.Overall, the scores of the importance of provisioning services were higher in the Upper Mountain and Foothill bundles than other bundles, whereas the scores of the importance of most regulating and cultural services were higher in the Oasis and Oasis-Desert transition bundles than other bundles.
The perceptions of the improvement or degradation of ESs showed that most respondents believed that provisioning services had been improved; however, most regulating services, except for SF, had declined in our study area (Figure 7).The majority of the respondents considered CP to have improved.In the study area, agricultural technology and social level developed rapidly; thus, CP was believed to have improved by most respondents.Although 64.53% of the respondents believed that LP improved, 25.38% considered that degradation occurred in LP, which was associated with the ecological policy, such as closing hillsides to grazing.The majority or 74.35% of the respondents believed that HQ declined, followed by WY (71.08%),CS (58.53%),SC (53.21%), and NLR (49.93%).Moreover, 47.89% of the respondents considered NLR to have improved, which may be associated with residential location.Figure 6B showed the importance of each ES in different ES bundles.The scores of the importance of CP were notable high (more than 0.65) and was viewed a relatively important ES in all ES bundles.LP and SC was viewed as more important in the Upper Mountain and Foothill bundles than Oasis and Oasis-Desert Transition bundle.The scores of the importance of WY and CS were notable high (more than 0.65) in the Oasis and Oasis-Desert Transition bundle.SF was viewed as more important in the Oasis-Desert Transition bundle, while HQ and NLR were viewed as most important in the Oasis bundle relative to other bundles.Overall, the scores of the importance of provisioning services were higher in the Upper Mountain and Foothill bundles than other bundles, whereas the scores of the importance of most regulating and cultural services were higher in the Oasis and Oasis-Desert transition bundles than other bundles.
The perceptions of the improvement or degradation of ESs showed that most respondents believed that provisioning services had been improved; however, most regulating services, except for SF, had declined in our study area (Figure 7).The majority of the respondents considered CP to have improved.In the study area, agricultural technology and social level developed rapidly; thus, CP was believed to have improved by most respondents.Although 64.53% of the respondents believed that LP improved, 25.38% considered that degradation occurred in LP, which was associated with the ecological policy, such as closing hillsides to grazing.The majority or 74.35% of the respondents believed that HQ declined, followed by WY (71.08%),CS (58.53%),SC (53.21%), and NLR (49.93%).Moreover, 47.89% of the respondents considered NLR to have improved, which may be associated with residential location.

Impact of Land Use Change on Ecosystem Services
On the basis of the correlation and partial correlation coefficients between cropland coverage change and ES change in each township (Table 3), changes in croplands had significant influence on most ESs.These changes were significant positive with CP (provisioning service) and negative with WY, CS, HQ, FR, and DR (regulating and cultural services).However, the influences of cropland changes on LP, SC, and SF service were weak and the significance levels of correlation coefficients were less than 0.05.Changes in the LP service were mainly influenced by ecological policy, such as closing hillsides to grazing and returning grazing land to pasture, whereas changes in the SC and SF services were mainly caused by climatic factors.SC and SF services were significantly correlated with PPT (p < 0.01).Analysis of building coverage and ESs in each township showed that building coverage exhibited a weak influence on most ESs and the significance levels of correlation and partial correlation coefficients were less than 0.05.Only WY and CS services exhibited significant correlations with building coverage (Table 3) (p < 0.05).Overall, changes in ES within the study area were mainly affected by the expansion of croplands.This result agrees with the conclusions of Wei et al. [10], who identified the significant relationships between cropland expansion and ESs via temporal serial data.The expansion of croplands occurred in natural areas, such as grasslands and deserts, which are providers of regulating and cultural services.Therefore, most regulating and cultural services (e.g., WY, CS, HQ, FR, and DR) declined.Some ES bundles changed with the land use changes from 2000 to 2015.This transformation illustrates that land use changes (e.g., expansion of croplands) in the Manas region have caused impacts on the ES bundles.During the study period, areas included in the Foothill, Oasis-Desert Transition, and Desert bundles decreased.However, the Oasis bundle increased in area due to the

Impact of Land Use Change on Ecosystem Services
On the basis of the correlation and partial correlation coefficients between cropland coverage change and ES change in each township (Table 3), changes in croplands had significant influence on most ESs.These changes were significant positive with CP (provisioning service) and negative with WY, CS, HQ, FR, and DR (regulating and cultural services).However, the influences of cropland changes on LP, SC, and SF service were weak and the significance levels of correlation coefficients were less than 0.05.Changes in the LP service were mainly influenced by ecological policy, such as closing hillsides to grazing and returning grazing land to pasture, whereas changes in the SC and SF services were mainly caused by climatic factors.SC and SF services were significantly correlated with PPT (p < 0.01).Analysis of building coverage and ESs in each township showed that building coverage exhibited a weak influence on most ESs and the significance levels of correlation and partial correlation coefficients were less than 0.05.Only WY and CS services exhibited significant correlations with building coverage (Table 3) (p < 0.05).Overall, changes in ES within the study area were mainly affected by the expansion of croplands.This result agrees with the conclusions of Wei et al. [10], who identified the significant relationships between cropland expansion and ESs via temporal serial data.The expansion of croplands occurred in natural areas, such as grasslands and deserts, which are providers of regulating and cultural services.Therefore, most regulating and cultural services (e.g., WY, CS, HQ, FR, and DR) declined.Some ES bundles changed with the land use changes from 2000 to 2015.This transformation illustrates that land use changes (e.g., expansion of croplands) in the Manas region have caused impacts on the ES bundles.During the study period, areas included in the Foothill, Oasis-Desert Transition, and Desert bundles decreased.However, the Oasis bundle increased in area due to the expansion of croplands (Figure 5).The ESs in the Oasis bundle were dominated by provisioning services, whereas the those in the Foothill, Oasis-Desert Transition, and Desert bundles were dominated by regulating and cultural services.The changes in the area of ES bundles were in agreement with the changes in the ESs, that is, increases in the provisioning services and decreases in the regulating services and cultural services.
On the basis of the questionnaire results about the reasons for ES change, more than 50% of the respondents believed that land use change was the main cause of the changes in CP, CS, HQ, and NLR (Figure 8).However, less than 40% of the respondents believed that land use change was the main cause of the changes in WY, SF, and LP.The changes in WY, SC, and SF were considered by most respondents to be caused by climate change.However, on the basis of the questionnaire results, the cause of the changes in LP was unclear.Overall, for each ES, at least 30% of the respondents believed that land use change was the main cause of ES change (Figure 8).expansion of croplands (Figure 5).The ESs in the Oasis bundle were dominated by provisioning services, whereas the those in the Foothill, Oasis-Desert Transition, and Desert bundles were dominated by regulating and cultural services.The changes in the area of ES bundles were in agreement with the changes in the ESs, that is, increases in the provisioning services and decreases in the regulating services and cultural services.
On the basis of the questionnaire results about the reasons for ES change, more than 50% of the respondents believed that land use change was the main cause of the changes in CP, CS, HQ, and NLR (Figure 8).However, less than 40% of the respondents believed that land use change was the main cause of the changes in WY, SF, and LP.The changes in WY, SC, and SF were considered by most respondents to be caused by climate change.However, on the basis of the questionnaire results, the cause of the changes in LP was unclear.Overall, for each ES, at least 30% of the respondents believed that land use change was the main cause of ES change (Figure 8).

Integrating Biophysical and Sociocultural Methods into Ecosystem Services and Land Use Change Research
Identifying the relationships between ESs and land use change and incorporating effective results into management practices present a key challenge in achieving sustainability [43,56].The requirements for assessing the different aspects of ESs, engaging with stakeholders, and meeting the demands of environmental policy makers motivate the integration of biophysical and sociocultural methods into ESs and land use change research.A few conceptual models have been built to facilitate such integration, including the framework for integrating ES supply and demand for land use management or landscape planning [76] and the nitrogen elemental framework for linking land use change, ES, and human well-being [77].Some comprehensive approaches, such as the specialist marking method based on the biophysical characteristics of land use types [78,79], were applied to assess the impact of land use change on ESs supply, demand, and budget.Our proposed framework (Figure 2) integrates multiple approaches, including biophysical methods (e.g., ecological and hydrological model), a sociocultural method (questionnaire survey), and statistical analyses (e.g., kmeans clustering and partial correlation analysis), to identify the relationships between ESs and land use change.By using these approaches, this study shows great improvements over previous works that examine the relationships between ESs and land use change [16][17][18][19][20][26][27][28][29][30].First, we integrated biophysical and social perception indicators in assessing ESs and land use change, and then obtain sufficient information of the supply and demand sides, thereby facilitating the development of efficient policies for ecological conservation and land use management.Second, we apply k-means clustering to identify the pattern and dynamics of ES bundles and to highlight the comprehensive influences of land use change on ESs.
Integrating biophysical and sociocultural methods into ES research presents several advantages.Benefits of nature (i.e., ESs) could be perceived by most respondents in the Manas region, which made

Integrating Biophysical and Sociocultural Methods into Ecosystem Services and Land Use Change Research
Identifying the relationships between ESs and land use change and incorporating effective results into management practices present a key challenge in achieving sustainability [43,56].The requirements for assessing the different aspects of ESs, engaging with stakeholders, and meeting the demands of environmental policy makers motivate the integration of biophysical and sociocultural methods into ESs and land use change research.A few conceptual models have been built to facilitate such integration, including the framework for integrating ES supply and demand for land use management or landscape planning [76] and the nitrogen elemental framework for linking land use change, ES, and human well-being [77].Some comprehensive approaches, such as the specialist marking method based on the biophysical characteristics of land use types [78,79], were applied to assess the impact of land use change on ESs supply, demand, and budget.Our proposed framework (Figure 2) integrates multiple approaches, including biophysical methods (e.g., ecological and hydrological model), a sociocultural method (questionnaire survey), and statistical analyses (e.g., k-means clustering and partial correlation analysis), to identify the relationships between ESs and land use change.By using these approaches, this study shows great improvements over previous works that examine the relationships between ESs and land use change [16][17][18][19][20][26][27][28][29][30].First, we integrated biophysical and social perception indicators in assessing ESs and land use change, and then obtain sufficient information of the supply and demand sides, thereby facilitating the development of efficient policies for ecological conservation and land use management.Second, we apply k-means clustering to identify the pattern and dynamics of ES bundles and to highlight the comprehensive influences of land use change on ESs.
Integrating biophysical and sociocultural methods into ES research presents several advantages.Benefits of nature (i.e., ESs) could be perceived by most respondents in the Manas region, which made the quantification of ESs supply meaningful.The results of these two methods for each ES can be compared.Both of these methods highlighted a performance improvement for the provisioning services and a performance decline for most of the regulating and cultural services.Specifically, the results for CP and LP revealed that the supply levels of provisioning services have increased, which was noticed by most of the survey respondents.The statistical analysis results identified land use change as the main driver of changes in few ESs (i.e., CP, CS, HQ, and NLR).This finding was supported by more than 50% of the respondents in the questionnaire survey.However, the quantification of RWEQ demonstrated a performance decline for SF, which was not reported in the survey results.Such inconsistency can be ascribed to two reasons.First, most of the respondents were living away from the desert (Table S3) and were thus less exposed to sandstorms.Second, the wind speed significantly decreased [10] in the study area, which might have led most respondents to believe that SF has improved.
By integrating biophysical and sociocultural methods into ESs, the perceptions of the importance of ESs and the ES supply level could be linked.Across diverse ES bundles, biophysical supply of ES was not completely accordant with the importance of ES.For example, in the Upper Mountain and Foothill bundles, the supply levels of CP and LP were less than 0.4 and relatively low (Figure 5c,d); however, the scores of the importance of the two services were more than 0.7 (Figure 6b).In the Oasis bundle, NLR had high importance and the score was more than 0.6; however, the supply level of the two services were less than 0.1.Although WY was viewed highly important in the Oasis-Desert Transition bundle and the score was more than 0.8, the supply level was less than 0.2.The factors of supply-side and demand-side both could cause the mismatches.For example, weak function of CP in Upper Mountain and Foothill bundles could cause the mismatches between the supply level and the importance of CP.In addition, high social demand for WY in the Oasis-Desert Transition bundle could also cause the mismatches between the supply level and the importance of WY.The identification of ES mismatch also found by other studies [36,76,[78][79][80][81] was crucial for land use management, but the systematic and detailed research still need to be explored by integrating biophysical and sociocultural methods into ESs and land use change.

Methodological Concerns
Although our study used multidisciplinary methods and incorporated local perspectives on ESs and land use change, some limitations should be considered when interpreting the results.In the biophysical methods, five models-RUSLE, water balance equation, RWEQ, CASA, and InVEST HQ model-were employed to estimate the supply of SC, WY, SF, CS, and HQ; the credibility of estimating results should also be evaluated.Because the experimental monitoring data on the regulating services are not existed in the Manas, we only can evaluate the credibility of estimating results on some ESs by indirect approaches.Based on the availability of data and the feasibility of approaches, we only evaluate the relevant results of WY and CS.Water balance equation and CASA were applied to estimate the WY and CS (Table 2).In the processes of estimation, ET and NPP are crucial variables and affected by vegetation physiological factors and complex environment factors, and the credibility of them should be considered [28].We proposed two possible approaches to illustrate the rationality of the results.The first approach was to compare the monitoring data from hydrological stations with the WY value (Figure A5).The second approach was to compare existing remote sensing products (i.e., MOD16 ET and MOD17 NPP) with the estimated ET and NPP.
WY simulated by water balance equation decreased from 105.76 mm in 2000 to 7.22 mm in 2015.The records of hydrological stations indicate that the runoffs showed a downward trend since 2000 (Figure A5), which enhances the credibility of the water balance equation in the study area.ET is a critical and uncertain parameter in the process of simulating water yield.It has strong variability and is affected by physiological vegetation factors and complex environmental factors.Therefore, the rationality of the water balance equation in the study area should be evaluated.At present, the period of MOD16 products only covers 2000-2014, and the simulated ET is compared with MOD16 ET in 2000 (Figure 9A).The significant correlation between MOD16 ET and simulated ET (R = 0.66; p < 0.01) indicates the rationality of the water balance equation in the study area.MOD17 NPP data were compared with the simulated NPP in this study.MOD17 NPP values are significantly correlated with the simulated NPP in 2000 and 2015 (R > 0.65; p < 0.01; Figure 9 B,C), which indicates the rationality of the CASA model in the study area.There seems to be nonlinear relationships between MODIS products and the simulated values (Figure 9A,C).The reason may be that the values of MODIS products (i.e., NPP and ET) in forests were underestimated in arid region of China, which also found by other research [28].A face-to-face questionnaire survey was used to collect information about the perceptions of stakeholders toward ES and land use change.However, performing such questionnaire survey entailed high costs given the complex topography and inconvenient transportation in Manas.The mean via post or e-mail also can be applied to some interviewees (e.g., specialist, professor, and government agents), except farmers and herders in the next steps.Moreover, during the interviews, the structured basis of the questionnaire does not allow much flexibility, whereas the motivations for the perceptions of the respondents toward ESs and land use change cannot be easily captured.In the future, socioeconomic data and advanced approaches can be used to understand the mechanisms behind the perceptions of respondents toward ESs and land use change.
We used correlation and partial correlation coefficients to identify the impacts of land use change on ESs.Given that some of these services were mapped by using indicators related to land use characteristics (i.e., land use types, NDVI, and LAI), the correlations between land use change and some ESs seems predictable.Therefore, the effective approaches for eliminating the influences of land use types, NDVI, and LAI on these correlations need to be explored.By combining biophysical and sociocultural methods, the results for ESs and land use change can be enhanced.For example, most respondents believed that CP and LS increased whereas SC, WY, CS, HQ, and NLR decreased, and the results agreed with the ES supply changes derived from the biophysical method.Several problems were also encountered in this study that need to be explored in the future.First, the data on the perceptions of the respondents toward ESs and land use change cannot be well linked to biophysical data, probably because of the mismatch in the spatial scale and the varying dimensions of the evaluation results.Specifically, the results of the questionnaire survey and biophysical methods for ES did not match, thereby influencing the depth integration development of ES biophysical and social attributes.In the future, questionnaire survey and objective data can be combined to explore the raster expression of ES questionnaire results.

Policy Implications
According to the evaluation results of ESs supply, provisioning services has been improved but regulating services are decreasing in the study area.Agriculture and animal husbandry experienced huge progress with the cost of declining regulating services.The decline in SC, WY, CS, and SF will threaten the sustainable production of CP and LP, for example, the drought and soil erosion caused by the decreasing in WY and SC may threaten the agriculture production.Ecological conservation policies are urgently needed in the Manas region especially some key area.Five ES bundles were identified at the township level and the patterns and dynamics of ES bundles could provide insights for land use management.In comparison with the other areas, the Upper Mountain bundle exhibited A face-to-face questionnaire survey was used to collect information about the perceptions of stakeholders toward ES and land use change.However, performing such questionnaire survey entailed high costs given the complex topography and inconvenient transportation in Manas.The mean via post or e-mail also can be applied to some interviewees (e.g., specialist, professor, and government agents), except farmers and herders in the next steps.Moreover, during the interviews, the structured basis of the questionnaire does not allow much flexibility, whereas the motivations for the perceptions of the respondents toward ESs and land use change cannot be easily captured.In the future, socioeconomic data and advanced approaches can be used to understand the mechanisms behind the perceptions of respondents toward ESs and land use change.
We used correlation and partial correlation coefficients to identify the impacts of land use change on ESs.Given that some of these services were mapped by using indicators related to land use characteristics (i.e., land use types, NDVI, and LAI), the correlations between land use change and some ESs seems predictable.Therefore, the effective approaches for eliminating the influences of land use types, NDVI, and LAI on these correlations need to be explored.By combining biophysical and sociocultural methods, the results for ESs and land use change can be enhanced.For example, most respondents believed that CP and LS increased whereas SC, WY, CS, HQ, and NLR decreased, and the results agreed with the ES supply changes derived from the biophysical method.Several problems were also encountered in this study that need to be explored in the future.First, the data on the perceptions of the respondents toward ESs and land use change cannot be well linked to biophysical data, probably because of the mismatch in the spatial scale and the varying dimensions of the evaluation results.Specifically, the results of the questionnaire survey and biophysical methods for ES did not match, thereby influencing the depth integration development of ES biophysical and social attributes.In the future, questionnaire survey and objective data can be combined to explore the raster expression of ES questionnaire results.

Policy Implications
According to the evaluation results of ESs supply, provisioning services has been improved but regulating services are decreasing in the study area.Agriculture and animal husbandry experienced huge progress with the cost of declining regulating services.The decline in SC, WY, CS, and SF will threaten the sustainable production of CP and LP, for example, the drought and soil erosion caused by the decreasing in WY and SC may threaten the agriculture production.Ecological conservation policies are urgently needed in the Manas region especially some key area.Five ES bundles were identified at the township level and the patterns and dynamics of ES bundles could provide insights for land use management.In comparison with the other areas, the Upper Mountain bundle exhibited the lowest population density and economic development level but the highest values for regulating services (e.g., WY and SC).The water resources conveyed from the Upper Mountain bundle are the cornerstone of socioeconomic development in the Oasis bundle.As a result, ecological conservation in the Upper Mountain bundle is a crucial link in the ES of this region, and deforestation should be banned in this bundle.Closing hillsides to ban grazing should be implemented in this bundle.In the Desert bundle, the habitat service maintains biodiversity, which underlies the sustainability of other ESs.Therefore, strengthening the habitat service is crucial for ecological conservation.Establishing a nature reserve in the key area of the Desert bundle can help maintain biodiversity and raise the ecological awareness of the local population.The increase in the Oasis bundle reduced the regulating services in the Manas region.The direct social benefits to human well-being produced by regulating services were weak.However, the maintenance of regulating services could not be ignored because they provide the foundation for provisioning services [82] and could exert an impact on human well-being in the long term [43].The increased amount in the Oasis bundle arose mainly from the conversion of the Oasis-Desert Transition bundle.Therefore, restricting reclamation in the bundle could alleviate decreases in the regulating services.Moreover, a large amount of the land in the Oasis bundle was salinized.Therefore, improving the quality of cultivated land via restricting reclamation may improve provisioning services.
According to the results of the questionnaire survey, the importance of each ES differed across various stakeholders and ES bundles.The herders mainly live in the Upper Mountain bundle, which is the key area of ES maintenance.The supply level of provisioning services (i.e., CP and LP) was viewed as highly important by these herders but was given low importance in the Upper Mountain bundle.Oasis bundle provided a high level of provisioning services, and this finding can provide important insights into the management of ESs.For example, a high supply of provisioning services can be transferred from the Oasis to the Upper Mountain bundle to improve the well-being of herders.Moreover, a high supply of regulating services in the Upper Mountain bundle should be maintained, and those enterprises that heavily pollute the area should be banned.The high regulating services in the Upper Mountain bundle could benefit the Oasis bundle, where most regulating services (e.g., WY) were viewed as highly important.However, the value of regulating services could not be reflected in economic markets.Therefore, implementing the ecological compensation system from the Oasis to the Upper Mountain bundle (e.g., exporting key provisioning service to the Upper Mountain bundle for compensation) can improve the well-being of residents in the Upper Mountain bundle.Apart from material compensation and monetary compensation, some alternative ways for ecological compensation should be further improved.For example, the governments of Oasis bundle can provide comfortable housing to herdsmen to improve their residence status or build highways to enhance the communication between Upper Mountain bundle and Oasis bundle.

Conclusions
In the Manas region, spatial-temporal changes in CP, LP, SC, WY, SF, CS, HQ, and NLR services from 2000 to 2015 were mapped using statistical data (i.e., annual crop output and livestock number), RUSLE, water balance equation, RWEQ, CASA, InVEST HQ model, and remote sensing interpretation data.Cluster analysis was applied to identify the ES bundle types, which are closely linked with land use.We then used a questionnaire survey to capture the perception of local residents toward ES and land use change.At the regional scale, the results indicated that provisioning services increased, whereas regulating and cultural services declined.Under the combined influence of natural and socioeconomic factors, multiple ESs clustered in five different locations of the Manas region-Upper Mountain, Foothill, Oasis, Oasis-Desert Transition, and Desert-and we grouped these into five ES bundles.From 2000 to 2015, the Oasis bundle expanded as a result of oasisization.However, the Oasis-Desert Transition and Foothill bundles declined accompanied by a decline in regulating services.We applied

Figure 1 .
Figure 1.(a) Location of the Manas region and (b) a digital elevation model of the study area.

Figure 1 .
Figure 1.(a) Location of the Manas region and (b) a digital elevation model of the study area.

Figure 2 .
Figure 2. Framework for integrating biophysical and sociocultural methods into ecosystem services (ESs) and land use change.CP-crop production, LP-livestock production, SC-soil conservation, WY-water yield, SF-sand fixation, CS-carbon fixation, HQ-habitat quality, and NLR-nature landscape recreation.

Figure 2 .
Figure 2. Framework for integrating biophysical and sociocultural methods into ecosystem services (ESs) and land use change.CP-crop production, LP-livestock production, SC-soil conservation, WY-water yield, SF-sand fixation, CS-carbon fixation, HQ-habitat quality, and NLR-nature landscape recreation.

Figure 3 .
Figure 3. Land use of the study area in 2000 (a) and 2015 (b).The land use datasets were interpreted by Landsat TM and was provided by the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences [53].The accuracies of the land use maps were more than 95%.

Figure 3 .
Figure 3. Land use of the study area in 2000 (a) and 2015 (b).The land use datasets were interpreted by Landsat TM and was provided by the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences [53].The accuracies of the land use maps were more than 95%.

Figure 4 .
Figure 4. Spatial distribution of ecosystem service (ESs) changes in the township level from 2000 to 2015.The figure shows (a) crop production change, (b) livestock production change, (c) soil conservation change, (d) water yield change, (e) sand fixation change, (f) carbon sequestration change, (g) habitat quality change, and (h) nature landscape recreation change.

Figure 4 .
Figure 4. Spatial distribution of ecosystem service (ESs) changes in the township level from 2000 to 2015.The figure shows (a) crop production change, (b) livestock production change, (c) soil conservation change, (d) water yield change, (e) sand fixation change, (f) carbon sequestration change, (g) habitat quality change, and (h) nature landscape recreation change.3.2.2.Changes in Ecosystem Service Bundles Based on the value of eight ESs in each township, k-means clustering was applied to classify the 47 townships into five categories.These categories represent five different ES bundle types.The number of townships that belong to each ES bundle ranged from three to 21 (Figure 5a,b).We utilized spider diagrams to visualize the level of different ESs in the five ES bundles using the normalized ES values (Figure 5c,d).ES bundle types were associated with landforms and land cover types, and the details of each ES bundle are shown in Figure 5a,b.The pattern and dynamic change of each ES bundle is as follows.

Figure 6 .
Figure 6.Survey results related to the importance of ESs across diverse stakeholders (A) and ES bundles (B).CP-crop production; LP-livestock production; SC-soil conservation; WY-water yield; SF-sand fixation; CS-carbon sequestration; HQ-habitat quality; NLR-nature landscape recreation.Note: The scores of the importance of ESs are not presented because no relevant questionnaire survey was performed on the Desert bundle.ES bundles in 2015 were applied to calculate the scores of the importance of ESs.

Figure 6 .
Figure 6.Survey results related to the importance of ESs across diverse stakeholders (A) and ES bundles (B).CP-crop production; LP-livestock production; SC-soil conservation; WY-water yield; SF-sand fixation; CS-carbon sequestration; HQ-habitat quality; NLR-nature landscape recreation.Note: The scores of the importance of ESs are not presented because no relevant questionnaire survey was performed on the Desert bundle.ES bundles in 2015 were applied to calculate the scores of the importance of ESs.

Figure 7 .
Figure 7. Survey results related to the cognitions of the improvement or degradation of ESs (CPcrop production; LP-livestock production; SC-soil conservation; WY-water yield; SF-sand fixation; CS-carbon sequestration; HQ-habitat quality; NLR-natural landscape recreation).
048 0.008 −0.643 ** −0.189 −0.711 ** −0.439 ** −0.545 ** PR2 −0.235 −0.042 −0.050 0.271 * −0.066 −0.430 ** −0.089 −0.133 CP-crop production; LP-livestock production; SC-soil conservation; WY-water yield; SF-sand fixation; CS-carbon fixation; HQ-habitat quality; NLR-Nature landscape recreation; R1correlation coefficient between cropland coverage change and ES change in the 47 towns of our study area; R2-correlation coefficient between building coverage change and ES change in the 47 towns of our study area; PR1-partial correlation coefficient between cropland coverage change and ES change (excluding the impact of building and climate factors) in the 47 towns of our study area; PR2-partial correlation coefficient between building coverage change and ES change (excluding the impact of cropland and climate factor) in the 47 towns of our study area.** correlation or partial correlation is significant at the 0.01 level; * correlation or partial correlation is significant at the 0.05 level.

Figure 7 .
Figure 7. Survey results related to the cognitions of the improvement or degradation of ESs (CP-crop production; LP-livestock production; SC-soil conservation; WY-water yield; SF-sand fixation; CS-carbon sequestration; HQ-habitat quality; NLR-natural landscape recreation).
CP-crop production; LP-livestock production; SC-soil conservation; WY-water yield; SF-sand fixation; CS-carbon fixation; HQ-habitat quality; NLR-Nature landscape recreation; R1-correlation coefficient between cropland coverage change and ES change in the 47 towns of our study area; R2-correlation coefficient between building coverage change and ES change in the 47 towns of our study area; PR1-partial correlation coefficient between cropland coverage change and ES change (excluding the impact of building and climate factors) in the 47 towns of our study area; PR2-partial correlation coefficient between building coverage change and ES change (excluding the impact of cropland and climate factor) in the 47 towns of our study area.** correlation or partial correlation is significant at the 0.01 level; * correlation or partial correlation is significant at the 0.05 level.

Figure 8 .
Figure 8. Survey results related to the factors of the improvement or degradation of ESs (CP-crop production; LP-livestock production; SC-soil conservation; WY-water yield; SF-sand fixation; CS-carbon sequestration; HQ-habitat quality; NLR-nature landscape recreation).

Figure 8 .
Figure 8. Survey results related to the factors of the improvement or degradation of ESs (CP-crop production; LP-livestock production; SC-soil conservation; WY-water yield; SF-sand fixation; CS-carbon sequestration; HQ-habitat quality; NLR-nature landscape recreation).

Sustainability 2019 ,
11, x FOR PEER REVIEW 19 of 27compared with the simulated NPP in this study.MOD17 NPP values are significantly correlated with the simulated NPP in 2000 and 2015 (R > 0.65; p < 0.01; Figure9 B,C), which indicates the rationality of the CASA model in the study area.There seems to be nonlinear relationships between MODIS products and the simulated values (Figure9A,C).The reason may be that the values of MODIS products (i.e., NPP and ET) in forests were underestimated in arid region of China, which also found by other research[28].

Figure 9 .
Figure 9. Validation of simulated ET and NPP through comparisons with MOD16 ET (A) and MOD17 NPP (B,C) data at cell scale.

Figure 9 .
Figure 9. Validation of simulated ET and NPP through comparisons with MOD16 ET (A) and MOD17 NPP (B,C) data at cell scale.

Figure A1 .
Figure A1.Spatial distribution of T (a) and PPT (b) in 2015 in the Manas region.

Figure A2 .Figure A2 .
Figure A2.Spatial distribution of ESs at the township level in 2000 and 2015.The figure shows (A) crop production, (B) livestock production, (C) soil conservation, (D) water yield, (E) sand fixation, (F) carbon sequestration, (G) habitat quality, and (H) nature landscape recreation.Sustainability 2019, 11, x FOR PEER REVIEW 23 of 27 Figure A2.Spatial distribution of ESs at the township level in 2000 and 2015.The figure shows (A) crop production, (B) livestock production, (C) soil conservation, (D) water yield, (E) sand fixation, (F) carbon sequestration, (G) habitat quality, and (H) nature landscape recreation.

Figure A3 .
Figure A3.Survey results related to awareness of the benefits of nature.

Figure A3 .
Figure A3.Survey results related to awareness of the benefits of nature.

Figure A3 .
Figure A3.Survey results related to awareness of the benefits of nature.

Figure A4 .
Figure A4.Survey results related to the importance of different ecosystem services (CP-crop production; LP-livestock production; SC-soil conservation; WY-water yield; SF-sand fixation; CS-carbon sequestration; HQ-habitat quality; NLR-nature landscape recreation).

Figure A5 .
Figure A5.Historical changes of runoffs in the study area.

Figure A4 .
Figure A4.Survey results related to the importance of different ecosystem services (CP-crop production; LP-livestock production; SC-soil conservation; WY-water yield; SF-sand fixation; CS-carbon sequestration; HQ-habitat quality; NLR-nature landscape recreation).

Figure A3 .
Figure A3.Survey results related to awareness of the benefits of nature.

Figure A4 .
Figure A4.Survey results related to the importance of different ecosystem services (CP-crop production; LP-livestock production; SC-soil conservation; WY-water yield; SF-sand fixation; CS-carbon sequestration; HQ-habitat quality; NLR-nature landscape recreation).

Figure A5 .
Figure A5.Historical changes of runoffs in the study area.

Table 1
summarized the data information required in the quantification of ESs, including data types, contents, purposes, and sources.The research data represent the following years, 2000 and

Table 1
summarized the data information required in the quantification of ESs, including data types, contents, purposes, and sources.The research data represent the following years, 2000 and 2015.

Table 1 .
Data information, including data types, contents, purpose, and sources.

Table 2 .
Indicator selection and quantification approaches for ecosystem services (ESs).

Table 3 .
Relationships between land use change and ES change.

Table 3 .
Relationships between land use change and ES change.

Table A1 .
Habitat suitability of different land use types and sensitivity of different land use types to threatening factors.

Table A1 .
Habitat suitability of different land use types and sensitivity of different land use types to threatening factors.

use Types Habitat Suitability Threatening Factors Settlement Cultivation Population River Road
Figure A5.Historical changes of runoffs in the study area.

Table A1 .
Habitat suitability of different land use types and sensitivity of different land use types to threatening factors.

Table A2 .
Main characteristics of the respondents.

Table A3 .
Land use change from 2000 to 2015 in the study area (10 2 km 2 ).The land use datasets in 2000 and 2015 were interpreted by Landsat TM and was provided by the RESDC, Chinese Academy of Sciences [53].*: The value is proportion of land use change area accounting for the total area.