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

Abstract

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.

1. 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 [16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42]. 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.2 trillion 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.

2. Materials and Methods

2.1. 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² [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.

2.2. 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 Section 2.3 and Section 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.

2.3. Ecosystem Services Supply: Biophysical Methods

2.3.1. 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.

2.3.2. Data Requirements and Preparation

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

Data Types	Data Contents	Data Purposes	Data Sources
Statistical data	Crop production and livestock inventory	Evaluating CP and LP	Statistical yearbooks [52]
Meteorological and hydrological data	Daily solar radiation, temperature, precipitation, and wind speed	Input parameters of SC, WY, SF, and CS evaluation	National Climatic Bureau [53]
	Annual runoffs	Evaluating rationality of WY change
Land use data	Cropland, woodland, grassland, water, urban, and desert	Evaluating land use change and NLR; Input parameters of SC, SF, CS, and HQ evaluation	Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences [54]
Basic geographic data	Road network and water system	Input parameters of HQ evaluation
	Elevation	Input parameters of SC evaluation
Remote sensing data	MOD13 NDVI and MOD15 LAI	Input parameters of SC, WY, SF, and CS evaluation	NASA’s Earth data [55]
	MOD16 ET and MOD17 NPP	Verification of estimation accuracy of ET and NPP
Soil attribute data	Content of sand, silt, clay, and organic carbon	Input parameters of SC and SF evaluation	Cold and Arid Regions Sciences Data Center at Lanzhou [56]

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; NDVI—normalized difference vegetation index; LAI—leaf area index; ET—evapotranspiration; NPP—net primary production.

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. 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.

2.3.3. Mapping of Individual Ecosystem Services Supply

Each ES was assessed by using the indicators and approaches described in

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

Ecosystem Services (ESs)	Indicators	Units	Approaches
Provisioning
Crop production (CP)	Annual crop output	t/km2	Statistical data in township-level units
Livestock Production (LP)	Livestock inventory	number/km2	Statistical data in township-level units
Regulating
Soil conservation (SC)	Soil retention amount	t/km2	Revised universal soil loss equation (RUSLE) [57]
Water yield (WY)	Water output	mm	Water balance equation [58]
Sand fixation (SF)	Binding sand quantity	t/km2	Revised wind erosion equation (RWEQ) [59]
Carbon sequestration (CS)	Carbon sequestration of noncrop plants [33]	t/km2	Carnegie Ames Stanford Approach (CASA) [60]
Habitat quality (HQ)	Habitat quality index	Dimensionless index (0–1)	InVEST (Integrate valuation of ES and tradeoffs) habitat quality model [61]
Cultural
Nature landscape recreation * (NLR)	Proportion of forest and desert cover	%	Remote sensing interpretation (Landsat)

. 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.

$$A_p=R\times K_{sl}\times LS$$

(1)

$$A_r=R\times K_{sl}\times LS\times C_{sl}\times P$$

(2)

$$A_c=A_p-A_r$$

(3)

(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)

$$ET=c+b\ast PPT\ast PET+a\ast PPT\ast LAI$$

(5)

$$PET=0.0162\times \left(\raisebox{1ex}{$RS$}\!\left/ \!\raisebox{-1ex}{$58.5$}\right.\right)\times \left(T+17.8\right)$$

(6)

(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].

$$S{F}_c=S{F}_p-S{F}_r$$

(7)

$$S{F}_p=\frac{2z}{S_p}\times Q_{\max p}\times e^{-{(z/S_p)}^2}$$

(8)

$$S{F}_r=\frac{2z}{S_r}\times Q_{\max }\times e^{-{(z/S_r)}^2}$$

(9)

$$S_p=150.71\times {(WF\times EF\times SCF\times K)}^{-0.3711}$$

(10)

$$Q_{\max p}=109.8\times (WF\times EF\times SCF\times K)$$

(11)

$$S_r=150.71\times {(WF\times EF\times SCF\times K\times C_s)}^{-0.3711}$$

(12)

$$Q_{\max }=109.8\times (WF\times EF\times SCF\times K\times C_s)$$

(13)

where SF_c is the amount of SF (kg/m²); SF_p is the potential amount of sand erosion (kg/m²); SF_r is the reality amount of sand erosion (kg/m²); 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₁), water-limiting factors (f₂), 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

$$NPP=0.5\times R_s\times FPAR\times f_1\times f_2\times W\times e_{\max }$$

(14)

(7) Habitat quality

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

$$Q_{xj}=H_j\left[1-D_{xj}^z/(D_{xj}^z+k^z)\right]$$

(15)

$$D_{xj}={\displaystyle \sum _{r=1}^R{\displaystyle \sum _{p=1}^{P_r}\left(w_r/{\displaystyle \sum _{r=1}^R{w}_r}\right)}}{r}_p{i}_{rxp}S{E}_{jr}$$

(16)

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. Habitat suitability of different land use types and sensitivity of different land use types to threatening factors.

Land use Types	Habitat Suitability	Threatening Factors
		Settlement	Cultivation	Population	River	Road
Cropland	0.2	0.7	0	0.2	0.1	0.5
Woodland	1	0.4	0.1	0.2	0.2	0.3
Grassland	0.6	0.2	0.5	0.4	0.1	0.3
Water	0.5	0	0.1	0.2	0	0.5
Desert	0.6	0.1	0.5	0.2	0.2	0.7
Others	0.1	0.4	0.1	0.2	0.2	0.3

) [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.

2.3.4. 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.

2.4. 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. Main characteristics of the respondents.

Characteristics	Category/Number (Proportion)
Sex	Male/462 (56.69%); Female/353 (43.31%)
Age	<20 years/26 (3.19%); 20–39 years/254 (31.17%); 40–59 years/505 (61.96%); ≥60 years/30 (3.68%)
Education	None/10 (1.23%); Primary/272 (33.37%); Secondary/269 (33.01%); University/264 (32.39%)
Residency period	<5 years/11 (1.35%); 5–9 years/57 (6.99%); 10–19 years/54 (6.63%); 20–39 years/354 (43.44%); ≥40 years/339 (41.60%)
Occupation	Farming/347 (42.58%); Grazing/69 (8.47%); Enterprise or public institution/312 (38.28); Others/87 (10.67%)

). 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.

2.5. 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.

3. Results

3.1. 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. Land use change from 2000 to 2015 in the study area (10² km²).

Land Use Types	Croplands	Woodlands	Grasslands	Urban	Water	Desert	Others
2000	44.07	9.45	79.57	3.17	14.88	53.33	25.34
2015	53.12	8.77	75.5	3.77	14.82	50.69	23.14
2000–2015	9.05	−0.68	−4.07	0.60	−0.06	−2.64	−2.20
Proportion of the total area (%) *	3.94	−0.30	−1.77	0.26	−0.03	−1.15	−0.96

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.

). The increased cropland areas came mostly from conversion of grasslands and desert areas. Grasslands and deserts dropped by 407 km² (1.77% of the total area) and 264 km² (1.15% of the total area), respectively. Moreover, the areas of woodlands and water bodies decreased by 0.29% (68 km²), and 0.03% (6 km²), respectively, whereas urban coverage increased by 0.26% (60 km²) (Table A3). The increased urban areas came mostly from conversion of grasslands and croplands.

3.2. Changes in Ecosystem Services Supply

3.2.1. 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⁵ t in 2000 to 1.79 × 10⁶ 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⁴ t and 6.00 × 10⁴ t. The total level of LP increased by 28.60%, ranging from 1.59 × 10⁶ t in 2000 to 2.04 × 10⁶ 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² in 2000 to 1838 t/km² 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² 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² in 2000 to 670 t/km² 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² were mainly concentrated in the northern (Figure 4e). The average value of CS changed from 61 t C/km² in 2000 to 52 t C/km² 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).

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.

The Upper Mountain bundle covered 5909.52 km² 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², changing from 3627.48 km² in 2000 to 3443.82 km² 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.57 km² in 2000 and 3717.47 km² in 2015, thereby increasing by 1470.90 km². 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², changing from 5229.00 km² in 2000 to 4393.46 km² 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², changing from 5930.31 km² in 2000 to 5478.60 km² 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).

3.3. 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 (Figure A4 and Figure 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.

3.4. 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. Relationships between land use change and ES change.

ES	CP	LP	SC	WY	SF	CS	HQ	NLR
R1	0.504 **	0.119	0.190	−0.679 **	−0.255	−0.537 **	−0.432 **	−0.526 **
R2	−0.149	0.245	0.157	0.349 *	−0.144	−0.231	0.021	0.017
PR1	0.428 **	−0.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; 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.

), 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).

4. Discussion

4.1. 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.

4.2. 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.

4.3. 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.

5. 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 the questionnaire survey and statistical analyses to identify reasons for improvement or degradation of ESs and found that land use changes in the Manas region had a significant impact on ESs. Results from the questionnaires for some ES such as CP, CS, HQ, and NLR indicated that more than half of the respondents cited land use changes as the main cause of change. Correlation and partial correlation analyses showed that oasisization was significantly and positively correlated with CP (a provisioning service) but was negatively correlated with WY, CS, HQ, and NLR (regulating and cultural services). By integrating the results of biophysical and sociocultural methods into ES research, the reliability of our conclusions was enhanced. In view of the large reductions in regulating and cultural services, the Manas region should limit reclamation in the desert to control the oasis expansion. This measure will improve the quality of cropland. Furthermore, a biodiversity conservation reserve should be established in the key area within the desert area to limit further reclamation.