Exploring the relative importance of socio-ecological factors to ecosystem services clusters: a support to spatially targeted management

CA includes two mountain countries (Tajikistan and Kyrgyzstan) and three plain countries (Uzbekistan, Kazakhstan and Turkmenistan), covering an area of approximately 4 × 10⁶ square kilometers (figure 2(a)). The elevation gradually decreases from the Pamir Plateau, the Tianshan Mountains and the Altay Mountains to the lowlands in the Western Caspian Sea. Marked differences exist in the spatial pattern of temperature (TEM) and precipitation (PPT). The annual average TEM is below 0 °C in some areas of Kyrgyzstan and Tajikistan, while the annual average TEM is above 15 °C in Turkmenistan, and the overall annual average TEM is 8 °C in CA. The annual PPT in desert and plain area (<200 mm) is far less than that in mountain areas (>1000 mm). Nearly 75% of CA is covered by grassland, desert and semi-desert ecosystems distributed in some basins and remote mountains areas, as well as in lower foothills and slopes. During 1960–2016, the grain yields experienced an increase (+1000 kg ha⁻¹) in CA (Liu et al 2020). The urban population of CA also increased from 24.01 million in 1995 to 33.13 million in 2015, with an increase of 37.97% (Li et al 2019a). Considering that many administrative units could share common ES management problems and challenges, the 43 state-level administrative units were selected as a spatial analysis unit (Raudsepp-Hearne et al 2010) (figure 2(b)).

Zoom In Zoom Out Reset image size

The datasets used in this study included: climate data (2.2.1), land-use and normalized difference vegetation index (NDVI) data (2.2.2), topography and soil data (2.2.3) and socio-economic data (2.2.4).

2.2.1. Climate data

2.2.2. Land-use and NDVI data

2.2.3. Topography and soil data

2.2.4. Socio-economic data

The spatial distribution of water and soil resources in CA is uneven. The downstream countries (Kazakhstan, Turkmenistan and Uzbekistan) suffer frequent wind erosion, while two mountainous countries (Kyrgyzstan and Tajikistan) are prone to water erosion (Li et al 2020a). Both wind erosion and water erosion not only accelerate to soil degradation, but also threaten food security and ecological stability (Sun et al 2015). Improving land use and soil conservation in Tajikistan Project recently proposed by the international atomic energy agency has effectively promoted sustainable land management practices and improved the economic well-being and livelihood of local residents (Mirgharibov 2018). In addition, the Amu Darya and Syr River basins are considered to be historically rich irrigated agricultural areas, and Kazakhstan is the largest rainfed agricultural region in CA (Schierhorn et al 2020). These farmlands provide basic supply services for inhabitants. Moreover, the grassland, forestland and oasis in CA are important providers of carbon sources and high-quality habitats (Han et al 2016). With the aggravation of climate change, the negative impacts of unreasonable human agricultural practices have increased a series of critical ecological problems, such as food and water shortages, land desertification, soil erosion and habitat degradation. Therefore, we selected seven types of ESs (water yield (WY), habitat quality (HQ), food production (FP), sand fixation (SF), soil conservation (SR), CS and net primary productivity (NPP)) to explore the complex relationship between ESs and their driving factors. Table 1 gives the specific methods and algorithms used to quantify each ES.

Table 1. Estimation methods for quantifying ESs.

	ESs	Abbreviation	Unit	Models	Algorithms	Description
1	Water yield	WY	(mm km−2)	InVEST	${Y_{x,y}} = \left( {1 - \frac{{{\text{AE}}{{\text{T}}_{x,j}}}}{{{P_x}}}} \right) \times {P_x}$	${P_x}$ refers to the annual precipitation on grid x; ${\text{AE}}{{\text{T}}_{x,j}}$ refers to the annual actual evapotranspiration for land use type j on grid x; $\frac{{{\text{AE}}{{\text{T}}_{x,j}}}}{{{P_x}}}$ refers to the approximation of the Budyko curve.
2	Habitat quality	HQ	Dimensionless index (0–1)	InVEST	${Q_{xj}} = {H_j}\left( {1 - \frac{{D_{xj}^2}}{{{k^2} + D_{xj}^2}}} \right)$	Dxj is the total threat level for land use type j on grid x; The k refers to half of the maximum value of Dxj. Hj refers to the habitat suitability of land use type j.
3	Food production	FP	t		Statistical yearbook	FP mainly includes wheat, barley, oats, maize and rice in this study.
4	Sand fixation	SF	t km−2	RWEQ	${\text{SL}} = \frac{{2z}}{{{s^2}}}{Q_{{\text{max}}}}{{\text{e}}^{ - {{\left( {\frac{z}{s}} \right)}^2}}}$	Qmax is the maximum transport capacity (kg m−1 ); s refers to the critical field length (m); z is the calculated distance of downward wind (m).
5	Soil conservation	SR	t km−2	RULSE	$A = R \times K \times LS \times \left( {1 - C \times {\text{P}}} \right)$	R refers to the rainfall erosivity factor (MJ·mm ha−1·h·a); K refers to the soil erodibility factor (t·ha·h/ha·MJ·mm); LS refers to the slope length factor and steepness factor; C is the vegetation management factor; and P is the soil practice factor.
6	Carbon storage	CS	g C m−2	InVEST	${C_{{\text{tot}}}} = {C_{{\text{above}}}} + {C_{{\text{below}}}} + {C_{{\text{soil}}}} + {C_{{\text{dead}}}}$	Cabove is aboveground biomass; Cbelow is belowground biomass; Csoil is soil carbon storage; Cdead is dead organic matter.
7	Net primary productivity	NPP	g C m−2	CASA	${\text{NPP}}\left( {x,t} \right) = {\text{APAR}}\left( {x,t} \right) \cdot \varepsilon \left( {x,t} \right)$	t refers to time; x refers to a grid location; APAR (x, t) represents the photosynthetically active radiation (MJ m−2); (x, t) is the light use efficiency (g C/MJ).

We selected evaluation models for the following reasons (table 1): (a) the InVEST (Integrated Valuation of ESs and Trade-offs) module operates quickly based on readily available data and is sensitive to changes in drivers such as climate or human activities (Sharps et al 2017). (b) The revised wind erosion equation (RWEQ) performs well in assessing global wind erosion, especially in semi-arid areas where climate information is scarce (Li et al 2020a). (c) The revised universal soil loss equation model has been widely applied to measure water erosion because of its advantages of less data requirements and integration with GIS databases (Kayet et al 2018). (d) The Carnegie–Ames–Stanford approach model, combined with remote sensing observations, successfully processed the spatiotemporal dynamics of NPP at regional scales (Li et al 2021a). We used the zonal statistics tool in ArcGIS10.2 to integrate the calculation results at grid scale for each ES.

ES hotspots (coldspots) can be considered the regions with relatively strong (weak) capabilities of providing one or more ESs (Spanò et al 2017). In this study, the top 25% of state units for each service were defined as ES hotspots, while the bottom 25% of state units for each service were defined as ES coldspots (Qiu and Turner 2013). Maps of multiple ES hotspots and coldspots were obtained by overlaying the single spatial distribution of hotspots and coldspots for each ES.

In this study, the ES clusters were explored by applying the following three steps. First, principal component analysis (PCA) was used to evaluate whether the ES co-occurs in spatial clusters (Marsboom et al 2018). Second, clustering algorithms (k-means cluster) were applied to determine ES groups that were significantly associated with the relevant PCA axes (Spake et al 2017). The ideal number of clusters was defined by evaluating the dominance of various ESs in each cluster. Third, the spatial distribution maps were used to visualize ES clusters, and radar diagrams were used to illustrate the relative ES value (Zhao et al 2018). K-means clustering can help managers to determine the regions that display similar relationships among ESs and quantitatively explain them through connecting with broad socio-ecological variables (Turner et al 2014).

To eliminate unit differences between ESs, the ES values were normalized to 0–1 using the minimum–maximum normalization method. The create random points and extract multi values to Points tools in ArcGIS 10.2 were used to sample the ES values of 10 000 random points in each cluster. Spearman's correlation analysis (R v3.5 statistical software) was adopted to determine trade-offs and synergies among ESs in the different ES clusters (Qiu et al 2018).

Based on the published literature and CA's actual socio-ecological condition (Spake et al 2017, Chen et al 2020, Dou et al 2020), 16 predictors related to climate, topography, vegetation, soil, population, economy and landscape structure were selected in this study (table 2). Key climate parameters such as TEM and PPT can directly affect ecological processes and indirectly affect human well-being (Pedrono et al 2016). The choice of topographic factors (DEM and slope) supports the assumption that land accessibility and isolation limit human activities and their impact on the local ES distribution. PD and GDP variables were chosen because of the impact of socio-economic and urban development on ecological degradation (Chen et al 2020). In addition, the landscape (e.g. cropland and urban) pattern can directly alter the ES supply by influencing ecosystem types, patterns and ecological processes (Wu 2013).

The multi-commonality among driving factors was examined by using the variance inflation factor (VIF) (Murray et al 2012). The VIF of all 16 explanatory variables was less than 10, indicating that the selected driving factors were representative and could explain the formation mechanism of ES cluster. RDA, an appropriate quantitative approach, was applied to detect the relationships between ES distribution and a set of environmental factors (Grimaldi et al 2014). The advantage of RDA is that it can visualize the order of response variables constrained by explanatory variables (Spake et al 2017). Monte Carlo permutation tests were used to verify the significance of 16 driving factors in the interpretation of ES distribution (Chen et al 2019a).

Figure 3 shows the spatio-temporal distribution (figures 3(a1)–(g1)) and changes (figures 3(A)–(G)) for the seven ESs from 1995 to 2015 in the 43 states of CA. WY (figures 3(a1)–(a3)), HQ (figures 3(b1)–(b3)), SR (figures 3(e1)–(e3)), CS (figures 3(f1)–(f3)) and NPP (figures 3(g1)–(g3)) all follow a similar spatial distribution pattern, being prominent mostly in the northern and eastern states of Kazakhstan and in some states of Tajikistan and Kyrgyzstan. This would indicate that there are intrinsic interrelationships between these ESs. In contrast, low values for these ESs are found in the states of Uzbekistan and Turkmenistan. However, we observed that high SF services (>200 t km⁻²) were provided in Navoiy, Dashgovuz and Lebap where there are large areas of desert (the Karakum and Kyzylkum Deserts) (figures 3(d1)–(d3)). Unlike other services, there are obvious spatial differences in FP (figures 3(c1)–(c3)). High FP values are found in North Kazakhstan (2678.50 t–5076.38 t), Akmola (1957.25 t–4529.73 t), and Kostanay (1769.418 t–3774 581.43 t) from 1995 to 2015.

Zoom In Zoom Out Reset image size

Based on spatial location of high and low values of the seven ESs from 1995 to 2015 (figures 3(A)–(G)), the results show that 93% and 88% of the states experienced increases in WY (+2.32 mm km⁻² 379 to +217.26 mm km⁻²) and SR (+21.61 t km⁻² to +287.34 mm km⁻²) (figures 3(A) and (E)), while 74% and 84% of the states experienced decreases in CS (−13.65 g C m⁻² 381 to −244.85 g C m⁻²) and NPP (−7.35 g C m⁻² to −37.95 g C m⁻²) (figures 3(F) and (G)). SF services increased significantly (>50 t km⁻²) in Dashgovuz and Akhal (figure 3(d)). Seventy six percent the states experienced obvious increases in FP (+2.32 mm km⁻² to +217.26 mm km⁻²) (figure 3(c)). Notably, there was a sharp downward trend in HQ in Karakalpakistan (figure 3(b)).

As shown in figure 4, an average of 37 of the 43 states contained at least one ES hotspot and an average of 27 states contained at least one ES coldspot over the period 1995–2015 (figures 4(a)–(c) and (e)–(g)). There is a concentration of ES hotspots (>4) in East Kazakhstan and Chüy (figures 4(a)–(c)). ES coldspots (>4) are mainly located in the Manghystau, Dashgovuz, Lebap, Balkan, Navoiy and Bukhara (figures 4(e)–(g)). Over the period of study, 55.81% of the states had no changes in hotspots and 46.51% had no changes in coldspots (figures 4(d) and (h)). Notably, there has been a decrease in ES hotspots in Pavlodar, West Kazakhstan, Akhal, Surkhandarya, Zhambyl, Jalalabad and Andijan (figure 4(d)). Meanwhile, the number of ES coldspots increased in three states in Eastern Kazakhstan, Karakalpakistan, Dashgovuz, Akhal, Lebap, Chüy, Ysyk-Köl and Gorno-Badakhshan (figure 4(h)).

Zoom In Zoom Out Reset image size

Five different ES clusters were explored by a principal components analysis and k-means clustering in this study (figure 5), including 'ESC 1: agricultural cluster', 'ESC 2: carbon cluster', 'ESC 3: sand fixation cluster', 'ESC 4: habitat cluster' and 'ESC 5: Soil and water cluster'.

Zoom In Zoom Out Reset image size

Agricultural cluster (ESC 1): four states accounted for 15.96% of the study area and were mainly concentrated in Northern Kazakhstan (figure 5(a)), including North Kazakhstan, Akmola, Kostanay and Pavlodar. This cluster contains a high percentage (60.83%) of cropland (table 3), which provided higher provisioning services (FP) (figure 5(b)).

Carbon cluster (ESC 2): seven states accounted for 38.35% of the study area and were mainly distributed in eastern Kazakhstan and northwestern Kyrgyzstan, providing high CS (22.40%) and WY (15.70%) (figures 5(a) and (b)). The main vegetation types included grassland (60.70%) and farmland (20.40%) (table 3).

Sand fixation cluster (ESC 3): 11 states accounted for 28.35% of the study area and were mainly concentrated in southwest CA (figure 5(a)), such as Karakalpakistan, Kyzylorda, Dashgovuz, Lebap, and Navoiy. This cluster mainly comprised 66.06% desert or bare land with low vegetation coverage (table 3), providing high SF (32.01%) services (figure 5(b)).

Habitat cluster (ESC 4): six states accounted for 32.05% of the study area and mainly provided HQ (36.28%) and CS (18.57%) (figures 5(a) and (b)). This cluster was mainly concentrated in the central grassland of Kazakhstan, including Atyrau, Aktobe, West Kazakhstan, Karaganda, South Kaz and Zhambyl (figure 5(a)). Main vegetation types included grassland (77.49%) (table 3).

Soil and water cluster (ESC 5): 15 states accounted for 8.86% of the area of the study region and primarily included the southeastern mountain states (figure 5(a)). The land cover mainly comprised grassland (47.25%) and cropland (28.63%) (table 3). More rainfall and vegetation coverage could substantially promote SR (17.68%), WY (17.44%) and healthy ecosystem habitats (HQ: 26.38%) (figure 5(b)).

A total of 21 potential trade-offs and synergies relationships between seven ESs were analyzed by Spearman's correlation coefficients at the ESC scales (figures 6(a)–(e)). For five different ESCs (figures 6(a)–(e)), SF displayed negative correlations with other ESs, especially in ESC1 (r = −0.34 to −0.66; p < 0.01) and ESC3 (r = −0.38 to −0.54; p < 0.01). FP was positively correlated (r = −0.2; p < 0.05) with NPP and SR in all ESCs (figures 6(a)–(e)). Good synergistic relationships existed pairwise between WY, SR, CS and NPP except for ESC3. Significant negative correlations (r = −0.21 to −0.66; p < 0.01 and p < 0.05) were found between FP and HQ in ESC1, ESC2 and ESC5 (figures 6(a), (b) and (e)), while a positive correlation (r = 0.39; p < 0.01) existed in ESC3 (figure 6(c)).

Zoom In Zoom Out Reset image size

Sixteen driving factors affecting the spatial distribution of seven ESs were identified by a RDA in this study (figure 7). The explanatory variables of RDA account for 56.32%, 58.46%, 54.37%, 57.69% and 60.37% of the variation in ESC1, ESC2, ESC3, ESC4 and ESC5, respectively, indicating that the selected factors represented each ESC distribution well. Among all ESCs (figures 7(a)–(e)), PPT was positively correlated with WY, SR, CS and NPP; VC were positively correlated with SR, NPP and CS; S_d and TEM were positively correlated with SF; DEM, SL, PF and PG were closely related to HQ; and FP was explained more by PD, GDP, GD, PU and PC. Notably, trade-offs between socio-economic factors (e.g. PD, GDP, GD and PC) and HQ were found in figures 7(a), (b), (d) and (e).

Zoom In Zoom Out Reset image size

Monte Carlo permutation tests were used to test whether the explanation of the ES distribution provided by the 16 environmental factors was significant. The factor contributions of the 16 factors presented differences throughout the ESCs (figures 8(a)–(e)). The spatial results of ESC1 (agricultural cluster) was dominated by provisioning services (FP) (figures 7(a) and 8(a)) and mainly determined by PPT (17.33%), PC (15.26%), TEM (13.61%) and PG (11.26%). ESC 2 (Carbon Cluster) (figures 7(b) and 8(b)), for which the importance of regulating services (CS) was high, was mainly affected by large values for PPT (21.14%) and VC (13.94%). PPT (20.22%) was still the most important driver of ESC3 (sand fixation cluster) distribution (figures 7(c) and 8(c)), followed by cropland (10.51%) and soil variables (9.60%–10.08%). In ESC4 (Habitat Cluster) (figures 7(d) and 8(d)), PPT (22.60%), PC (13.69%), SL (12.57%) and PG (7.44%) strongly impacted the multi-ESs. Figures 7(e) and 8(e) show that PPT, VC, TEM and DEM were the most important drivers affecting ESs for ESC5 (Soil and water cluster) with higher regulating services (SR and WY), accounting for 19.38%, 16.71%, 10.70% and 10.65% of the total explanatory variables, respectively. Notably, the explanatory power of PD, GDP and GD was low and not significant in 5 ESCs (figures 7(a)–(e) and 8(a)–(e)). By comparison, the contribution of the natural variables was much higher than that of the socio-economic variables, indicating that natural factors play a major role in the spatial distribution of ESCs (figures 7(a)–(e) and 8(a)–(e)).

Zoom In Zoom Out Reset image size

This study demonstrated that ES clusters in CA were much more sensitive to natural variables (e.g. climate and vegetation) than to socio-economic variables (e.g. population and economy) (figures 7(a)–(e) and 8(a)–(e)). Previous studies have indicated that stressed ecosystems like dryland are mainly controlled by natural environment factors (Chen et al 2019a, Zhu et al 2019, Li et al 2021). Among the environmental variables, we found that PPT and VC were the main factors influencing the associations among ESs (figures 8(a)–(e)). PPT (17.33%–22.60%) was the most important driver for promoting positive correlation between the WY, SR, CS and NPP services (figures 7(a)–(e) and 8(a)–(e)). The increase in rainfall directly changes surface runoff and detaches soil particles, indirectly leading to the increase in water production and soil water erosion (Ran et al 2012, Fang et al 2015). Furthermore, PPT has been proved to be the dominant factor regulating carbon stocks of dryland in CA (Zhu et al 2019). The ephemeral plants determine approximately 75% of the water/carbon and biomass dynamics in the temperate desert communities (Gang and Li 2014), and their biomass is closely correlated with the rainfall fluctuations during the growing season (Wang et al 2006). With a 0.4 °C decade⁻¹ warming rate during the past 30 years, CA is a hotspot of climate change (Hu et al 2014). Considering the high sensitivities of ES clusters to climate change and the high uncertainties of future climate trajectories in this region (Thevs et al 2013, Han et al 2016, Li et al 2020a), it is advisable for policymakers to closely monitor climate change patterns and prepare adaptive ecological management strategies (e.g. limiting irrigation water usage and reducing grazing stress to protect VC in the ESC2 and ESC5 which were most sensitive to PPT) in advance based on future climate scenarios (figures 8(b) and (e), Saidi and Spray 2018, Schierhorn et al 2020). When PPT was not a limitation factor, we found that the synergies among SR, NPP and CS among ES clusters were closed related to VC (figures 7(a)–(e)), especially in the mountain cluster (ESC5) (figure 8(e)). This finding is consistent with a previous research (Zhang et al 2020a), and is reasonable because VC is a good indicator of soil erosion, carbon sequestration, ecosystem productivity, etc (Chaitra et al 2018, Deng et al 2021).

Overall, our results are similar to other studies in arid regions (Lyu et al 2019, Dou et al 2020). Although the explanatory power of some socio-economic factors (e.g. PD, GDP and GD) was low in CA (figures 8(a)–(e)), the antagonistic effects of these factors on HQ cannot be ignored in cross-regional ecological management (figures 7(a)–(e), Sallustio et al 2017). Good habitat environment was provided by large areas of grassland (77.49%) to maintain the biodiversity in habitat cluster (ESC4) (figure 5(b) and table 3). This area has been listed as a key biodiversity protection area by international organizations (Wagner et al 2020). Nature reserves and ecological corridors should be established to effectively reduce intensive grazing activities and raise the species protection awareness of the local population (Liu et al 2009, Negev et al 2019).

An important recommendation provided by this study is that ES clusters (rather than administrative or land-cover/land-use boundaries) should be treated as the basic ecological management unit in CA (Spake et al 2017), and different management strategies should be designed in accordance to the major interactions among the ESs in each ES cluster (Chen et al 2020). For example, an obvious negative relationship between high provisioning services (FP) and regulating and maintenance services (HQ) was found in ESC1 (figures 5(b) and 6(a)), which indicates that it is necessary to balance the trade-off between FP and HQ to promote the overall ESs in this region (Bernués et al 2019). Land management activities (like converting ecological land to cropland) that enhance certain ES (e.g. FP) at large costs to other ESs should be prevented (Searchinger et al 2015). In contrast, good synergies among CS and other regulating services were found in ESC2 (figures 5(b) and 6(b)). Measures such as increasing vegetation coverage and forestland area should be taken to stabilize the synergies between these ESs (figures 7(b) and 8(b), Mouchet et al 2017, Zhang et al 2020a).

However, effective ecological management policies based on ES clusters would require cross-country coordination (Kumar 2002, Batsaikhan and Dabrowski 2017). Since 2000, with the assistance of the United Nations and other regional organizations (Mangalagiu et al 2019), policies and laws on ecological protection, such as the Implementation of integrated water resource management (UNDP 2004), the Implementation of Integrated Pollution Prevention and Control (IPPC) (UNECE 2010), Law on Environmental Protection (Hao et al 2019), etc, have been continuously stipulating and integrating into international practices. Under the coordination of UNDP and UNECE, the five CA countries have made great progress in ecological protection in recent years (Suleimenov et al 2012, Turayeva 2012, Thevs et al 2013), but different attitudes towards ES management have created challenges to the Sustainable Development Goals. One of the challenges is how to balance the short-term and local ES benefits with the long-term and global ES values (Milner-Gulland 2012, McPhearson et al 2013). For example, promoting FP through farmland expansion is attractive to CA countries because it brings quick and materialized benefits to local societies (Schierhorn et al 2020). However, farmland expansion would also lead to degradation of other continental or globally important ESs (Li et al 2015, 2020a), such as SF (CA is the major sand-dust source to the northeastern Asia) and HQ (CA accounts for over 80% of global temperate desert ecosystems).

Therefore, cross-ES clusters coordination is important for understanding the associations among ESs and managing regional ES as an integrated system (Malmborg et al 2021). For example, although the ESC5 was mostly distributed in the mountainous areas (figures 1(a) and 5(a)) with low PD and economic development level, its water resources were the cornerstone of the normal supply of ESs in the sand fixation cluster (ESC3) (figure 5(b), Hao et al 2019). Also, the regulating services (e.g. WY and SR) from the upstream cluster (ESC5) had regional importance to other ES clusters (figures 1(a) and 5(a), Li et al 2021). To avoid losing the erosion regulating service and promote sustainable soil management (Jin et al 2008), deforestation should be banned and cropland on steep slopes should be converted into forestland in these areas (Decocq et al 2016, Eguiguren et al 2019). However, the value of regulating services not reflected in economic markets (Raudsepp-Hearne et al 2010). For all clusters, policy makers should prioritize implementing a series of compensation for ES (PES) policies (Fu et al 2018, Lin et al 2020), including 'Conversion of Farmland to Forests Compensation' (Yang et al 2019) and 'Returning Grazing Land to Grassland (Shao et al 2016),' to reduce the economic losses of residents through government subsidies (McElwee et al 2020).