Grassland degradation by shrub encroachment: Mapping patterns and drivers of encroachment in Kyrgyzstan

Shrub encroachment is one of the major phenomena of grassland degradation throughout the world. One of the largest grassland areas is found in Central Asia. This area is playing an important role for biodiversity and food security, especially in a context of climate change. While sustainable land management is critical to avoid shrub encroachment, the dynamics of such encroachment are often unclear, especially in remote areas. We assessed the spatial and temporal dynamics of shrub encroachment in Suusamyr valley in Kyrgyzstan and analysed how these patterns relate to abiotic and social factors. To assess the extent of shrub encroachment we performed a random forest classification with Google Earth Engine to depict shrub encroachment between 2000 and 2019. By using a combination of vegetation indices and digital elevation derivatives in the classification model, the overall classification accuracy reaches above 90%. The mapping of temporal changes revealed an increase in shrub encroachment of 48% within the valley, adding up to an overall coverage of 2% of the grassland area. The spatial patterns of shrub encroachment are more strongly associated with anthropogenic factors, such as locations of roads and settlements, than with variations in precipitation and amount of organic carbon in soils. Around 45% of shrub encroachment occurs close to main roads and local paths, and within around 5 km distance from temporal settlements. This degradation is probably caused by overgrazing near settlements and the loss of grasslands in those areas puts additional pressure on local livestock farming, potentially leading to further degradation in neighboring areas. A shrub encroachment susceptibility map was created to identify areas vulnerable to further shrub expansion, affecting valuable grassland area. The identification of these areas at risk may inform sustainable grassland management by better allocating pastures according their capacity and reducing pressures on those places most susceptible to degradation.


Introduction
Grasslands are one of the most important global ecosystems that are not only providing carbon storage, but also habitat for pollinators, and forage for livestock and wild herbivores (Scurlock and Hall, 1998). Grassland land-use change and degradation of natural grasslands are amongst the largest threads to biodiversity decline worldwide, but receive less attention than forest loss and degradation. It is estimated that 23% of global grasslands are degraded or poorly managed (Brierley, 2013;Conant et al., 2001;Oldeman, 1992). While many studies address losses or degradation of forests, grasslands dynamics are not frequently assessed and monitored. This is partly due to their heterogeneity and rapidly evolving structure.
Studies that have addressed changes in grasslands have mostly focused on grassland degradation, which can be assessed through changes in biomass and productivity (Fassnacht et al., 2019;Xu et al., 2020;Zhou et al., 2017), changes in the area covered by grass Lyu et al., 2020;Zha et al., 2003) or shrubs (Guido et al., 2017;Knapp et al., 2008;Peng et al., 2013), or changes in soil conditions (Žížala et al., 2019), for example through the responses of soil bacterial and fungal community to degradation (Fu et al., 2012;Zhang and Fu, 2021;Zong and Fu, 2021), and biodiversity dynamic (Plantureux et al., 2005). One of the largest grassland ecosystems is located in Central Asia. Countries from the Central Asia region that are dominantly covered by grassland, such as Kazakhstan, Kyrgyzstan, Tajikistan, Turkmenistan, and Uzbekistan, are considered as a fragile ecological area. The preservation of grasslands in this area gets attention from non-governmental organizations (NGOs), governments and other environmental agencies (Lipovsky, 1995). Abiotic factors in conjunction with socio-economic conditions limit agricultural and livestock production in the region, and impede sustainable usage of the natural resources (Berkum, 2015). Natural conditions coupled with anthropogenic pressures on grasslands have caused a transformation of pastures into deserts or woodlands (van Auken, 2000). Climatic variation causing droughts, changes in soil nutrients and water supply put additional strain on grasslands, especially in arid and semi-arid areas (Eldridge et al., 2011;Wang et al., 2017).
Overgrazing, prescribed and out-of-control fire and mismanagement change the structural characteristics of grasslands vegetation, reduce biodiversity, and cause soil erosion (Brierley, 2013). In this context, the conventional approach on grassland area usage and allocation should be altered based on an understanding of the ongoing degradation and grazing processes, especially for those locations most sensitive to degradation.
This paper focuses on the grasslands of Kyrgyzstan and in particular on one of the largest natural pasture areas: Suusamyr. Grasslands in Kyrgyzstan are represented mainly by pastures and cover around 45% of the country's area (NSC, 2021). They are an essential land asset for well-being of rural inhabitants in the region. According to the National Statistics Committee report (NSC, 2019), the livestock population has reached about 15 million animals in the past 10 years. Some surveys (Isakov and Thorsson, 2015) suggest that these numbers are underestimated (for around 30%) because farmers tend to avoid paying for the rights to use pastures. According to the Ministry of Agriculture and Land Reclamation, about 70% of the pastures are estimated to be affected by degradation (ELD, 2017). Such assessment of grassland degradation, however, is carried out through inventory information from local authorities, which does not provide spatial patterns of changes in grassland conditions, which are essential to identify focus areas for improved management.
Pastures in Kyrgyzstan are under pressure due to population growth and an increasing dependence on grassland biomass resources as result of reduced investments after the collapse of the Soviet Union. Intensive use of the pastures without rotation led to overgrazing of the most accessible pastures and underuse of remote parcels (Yu and Kasymov, 2020). Diminishing of the grassland productivity due to the degradation processes directly influence well-being of the rural inhabitants (Zhaparov, 2015). Therefore, the loss of these pastures is not only an ecological concern, but also threatens an important livelihood option for local population. Several projects supported by United Nations Development Programme (Abdylakykova and Usubaliev, 2014), World Bank (World Bank, 2019), and The Deutsche Gesellschaft für Internationale Zusammenarbeit (GIZ) (Quitmann et al., 2019) aimed to improve pasture monitoring and livestock management in Kyrgyzstan. The results of these projects indicate an ongoing degradation in grasslands, including in Suusamyr. Several aspects of grassland degradation were detected: shrub encroachment, soil erosion and productivity decrease. Replacement of palatable species with shrubs or other types of unpalatable vegetation has been a critical issue for decades, and some authors report that up to 36% of pastures have been have encroached by shrubs or other unpalatable species (Semenova et al., 2020;Zhumanova et al., 2018). The major problem identified in these studies is the encroachment of an unpalatable shrub species Caragana (Caragana arborescens) (ELD, 2021;Kilyazova et al., 2019). Caragana is a widespread shrub species within grassland ecosystems, especially in arid and semi-arid areas. The stem of this shrub consists of compound leaves and thorns (Zhang et al., 2006). On the one side, it plays a significant role in preventing soil erosion, particularly in mountainous regions, it contributes to the restoration of degraded soil by fixing nitrogen  and it can serve as a fuel for local population. On the other side, it poses an additional threat to cattle since its thorns and prickles can harm animals (Fig. 1). Furthermore, its uncontrolled encroachment reduces natural grass cover and replaces palatable vegetation. Some studies (Kilyazova et al., 2019) demonstrate that preventive measures can impede the encroachment, while elimination of existing shrubs requires usage of chemical agents, which is financially demanding and can lead to soil deterioration.
It is hypothesized that during the last decades, the area covered by this shrub in the Suusamyr pasture has been constantly increasing, impairing quality and sustainability of the pastures. However, there is no assessment of the current spatial distribution of shrub encroachment, and how it has changed over the time. This information is essential for developing grassland management strategies that ensure food security in the region, sustainable land management and pasture parcel allocation (Shamshiev et al., 2017).
While other studies focus on aspects of combatting existing shrub encroachment in Suusamyr pasture (Kilyazova et al., 2019;Semenova et al., 2020) or assess effects of Caragana shrub on grassland dynamics (Umanova, 2021), this paper aims to understand and explain the spatial and temporal patterns of shrub encroachment to help target these management interventions and better understand the drivers of the degradation process.
Our research aims to assess shrub cover changes in the natural pasture of Suusamyr in Kyrgyzstan from 2000 until 2019. We used Landsat-5/7/8, Sentinel-2 and a digital elevation model in Google Earth Engine (GEE) to identify the spatial distribution of shrub cover changes in relation to potential natural and social determinants of this encroachment. We hypothesize that spatial patterns of shrub encroachment can be explained by a combination of social and biophysical variables, similar to earlier findings of Pazur et al. (2021) for a grassland area in northern Eurasia.
Based on the research results we aim to make inferences about potential areas of future shrub encroachment that have similar conditions as currently encroached areas.

Study area
The major natural pasture of Kyrgyzstan, Suusamyr is located in the central part of Kyrgyzstan and lies between 41 • 60 ′ N and 42 • 40 ′ N and 73 • 10 ′ E and 74 • 60' E ( Fig. 2). The study area covers an area of approximately 471.000 ha measuring about 150 km from west to east and 60 km from north to south. The territory is represented by various landforms including ridges, hills, rocks, river valleys and flats. The central part of the valley forms a plateau at an altitude between 2200 and 2600 m a.s.l. There is a main river Suusamyr and about 200 rivers and streams fed by glaciers are flowing into this river. The climate of the valley is influenced by the mountains. The area includes different natural landscapes such as steppes, hills, mountains, rocks, valleys and grasslands, as well as agricultural areas. The average annual temperature is below 0 • C, precipitation averages around 345 mm. Snow covers the valley from November and lasts in some areas until the end of June.
During the Soviet period, state authorities were responsible for the regular monitoring of pastures. The livestock movements between seasonal pastures were an important component of grazing. However, Kyrgyzstan was subsidized with winter nutrition for cattle, which meant that migration could be less distant and allowed the livestock numbers to reach a historical maximum (Robinson and Fabian, 2015). After the collapse of the Soviet Union subsidies for winter nutrition along with other pasture improvement programs such as fertilizing, sowing, tillage and pest control have been terminated. As a result, grazing became very dependent on natural pastures.
New legislation on pasture use and management was introduced in 2009 and since then local pasture committees became organizations responsible for the monitoring of pasture state. However, due to insufficient administrative capacity, lack of clarity regarding their responsibilities and issues of intersectional collaboration, their operation is considered as not very effective (Liechti, 2012). The absence of reliable information on livestock numbers in each committee complicates the estimation of the carrying capacity of pastures and establishing recommendations about the sustainable number of livestock. Statistical reports (NSC, 2019) show that the number of livestock has been increasing in the past 10 years reaching the maximum recorded during Soviet Union period. At the same time pasture degradation in the country is reported to reach 40% (Kyrgyzsgiprozem, 2015).

Mapping shrub encroachment
Following similar studies on grassland degradation we define grassland degradation as a process in which grassland productivity, biodiversity and ecosystem services decline due to the combination of climatic factors and anthropogenic activities (Chen et al., 2019;Zhang et al., 2018;Zhou et al., 2014). Grassland degradation is normally studied by combining field investigation and RS data. A significant number of studies introduces vegetation indices to identify grassland degradation (Fassnacht et al., 2015;Kong et al., 2019;Peng et al., 2013;Tanser and Palmer, 1999). For long term analysis researchers often use imagery from Landsat sensors (Fassnacht et al., 2015;Saltz et al., 1999;Zha et al., 2003) or MODIS (De Leeuw et al., 2019;Liu et al., 2019) to assess different aspect of degradation. For mapping grassland conditions in recent years, Sentinel-2 imagery (launched in 2015) is now more frequently used (Asam et al., 2019;Bekkema and Eleveld, 2018;Vermeulen et al., 2021). In-situ observations usually include ground data collected by spectrometers Lyu et al., 2020;Xu et al., 2020). However, it is challenging to combine in-situ samplings with low and medium resolution data due to the difference of scale, especially in heterogeneous landscapes. Remoteness and study area size pose further challenges in terms of accessibility and affordability of collecting field data.
Our research is based on long term time series of MODIS and Landsat data. Sentinel-2 was used only for the 2019 mosaic, because there were no Landsat data available for that year. Considering variation in spectral reflectance over different seasons, it is important to select the best period during the year for observing the phenomenon of interest and use that period consistently (Lawley et al., 2016). The reviewed studies indicate that data acquired during the growing season provide a higher likelihood of achieving reliable results (Alcantara et al., 2013;Ali et al., 2016;Dusseux et al., 2014;Price et al., 2002). Previous papers have shown the highest accuracy for distinguishing grass from shrubs during the spring season due to high shrub visibility (Filippa et al., 2015;Soubry and Guo, 2021). However, a major part of the case study area is covered by snow till the end of June. Therefore, we have used vegetation indices to identify the most appropriate month for distinguishing shrubs in the study area (Prishchepov et al., 2012;Seong et al., 2020). We manually digitized areas covered with shrub and grasses (100 polygons for each class) for which we calculated average NDVI using a MODIS dataset with 250 m spatial resolution (MOD13Q1.006 Terra Vegetation Indices). The NDVI values for grass and shrub areas differ less in the months of May, June and September as compared to other growing season months. Therefore, July and August were chosen as the most suitable months for differentiating shrubs from grass in this study area (Fig. 3).
While MODIS datasets are deemed highly applicable for generating time series, the coarse spatial resolution constrains its application in heterogeneous areas such as grasslands (Ali et al., 2016). A classification of grass and shrub areas based on MODIS NDVI data resulted in only 60% overall accuracy. For dry years with average July precipitation of less than 2.5 mm (2005,2010,2013,2015) the accuracy was even lower (between 30 and 50%). Therefore, datasets with finer spatial resolution (images from Landsat 5/7/8 missions (with 30 m spatial resolution) or Sentinel -2 (with 15 m spatial resolution)) were considered to be required for the accurate identification of locations of shrub encroachment. A complicating factor is the location of the case study area in a cloud-dominated region (based on review of available data) impeding the creation of cloudless mosaics. Therefore, only the following years were selected for analysis, optimizing the availability of cloudless mosaics: 2000 (July-August), 2005, 2016 and 2019 (July).
Image classification was conducted on the Google Earth Engine (GEE) platform using mosaics with a spatial resolution of 30 m. GEE is a cloud-based solution providing an access to repositories of georeferenced datasets (satellite images, weather and climate datasets, land cover maps, digital elevation models etc.) as well as powerful computation for analysis of large datasets (Gorelick et al., 2017). Classification processes in mountainous areas are more challenging and require more sophisticated methods than in flatlands because of the variability in elevation affecting the reflectance of remote sensing imagery (Shi et al., 2009). To classify the shrubs within the area of interest our samples should be exclusive and exhaustive, representing each category of shrub and non-shrub cover (Gartzia et al., 2013). To fulfil this requirement, we manually digitalized and homogeneously allocated 150 non-overlapping sample polygons (10708 training pixels for 2000; 10680 -for 2005; 11464 -for 2016; 11693 -for 2019) for shrub and non-shrub classes for each year (Ramezan et al., 2019).
We choose to use a pixel-based supervised Random Forest (RF) algorithm for its ability to prevent overfitting and its accuracy of capturing the class of interest by creating an ensemble of decision trees (Belgiu and Drȃguţ, 2016;Na et al., 2010;Shelestov et al., 2017;Teluguntla et al., 2018). In order to assess the accuracy of the RF we applied an internal cross validation technique by setting aside a third of the samples (so called out-of-the bag samples) as validation data (Breiman, 2001;Janitza and Hornung, 2018). As variables for the RF classification we selected different vegetation indices (VI) and a digital elevation model with derivatives. We calculated the following VIs: Land Surface Water Index (LSWI), Normalized Difference Vegetation Index (NDVI), Enhanced vegetation index (EVI), Modified Soil-Adjusted Vegetation Index (MSAVI), Bare Soil Index (BSI), Index-based Built-up Index (IBI). VIs provides valuable indicators of vegetation activity and its conditions (Bounoua et al., 2000;Shan and Xu, 2013;Usman et al., 2015) based on the ability of healthy vegetation to absorb most part of the visible spectrum and reflect a large amount of the near-infrared lights. Whereas the reflection of unhealthy or sparse vegetation mostly consists of visible light, the reflection of healthy vegetation is affected mostly by wavelength in near infrared part of the spectrum (Tucker et al., 1994;Zhu et al., 2016). LSWI uses the shortwave infrared (SWIR) and the NIR wavelength allowing to estimate the water content in vegetation. Chandrasekar et al. (2010) analysed the response of LSWI to rainfall showing that this index is useful for monitoring the increase of liquid water content in vegetation and soil background. MSAVI and BSI allow to compensate the soil background influence and thus minimize the problems with classification of bare soil (Ren and Zhou, 2014;Rikimaru et al., 2002). Usage of IBI enhanced differentiation of built-up areas and suppresses background noise (Xu, 2008). Additionally, we integrated a number of Digital Elevation Model (DEM) derivatives: aspect, eastness, northness and slope. The DEM used in this research was taken from the Shuttle Radar Topography Mission (SRTM) with a resolution of 30 m (Farr et al., 2007).
For all variables we computed in GEE the importance using the mean decrease in impurity (or gini importance) mechanism (Carvajal et al., 2018). IBI has the lowest importance among the indices which is explained by the only sparse and scattered settlements locations in the study area. On the other hand, the importance of slope and absolute elevation is high due to the mountainous and rugged terrain. As a next step we analysed the 4 layers depicting shrub coverage for 2000, 2005, 2016 and 2019 using ArcGIS Pro 2.1 and identified areas where and when changes in coverage had occurred.

Exploring shrub encroachment
The selection of possible drivers of shrub encroachment was based on the conceptual framework of proximate drivers and underlying causes of land change suggested by Geist and Lambin (2001). We analysed the role of differences in location conditions and potential drivers by comparing the patterns of shrub encroachment to the variation in a range of geographic factors that either represent the variations in environment or reflect the influence of the potential drivers. This approach and research concept is frequently used in studies applicable to the region and land cover change studies in general (Li et al., 2012;Overmars and Verburg, 2005;Veldkamp et al., 2001).
Shrub encroachment into grasslands is occurring due to increases of anthropogenic disturbances and climatic fluctuations (Knapp et al., 2008;Urbina et al., 2020). As indicators of human activity, we used variables indicating the distance from settlements (permanent and temporal housing), the distance from the road network (main roads and local paths) and the accessibility time to the nearest populated area. Datasets containing information about settlements and roads were created by manually digitizing objects. As a basemap for identification of small-scale settlements, roads and paths we used the ESRI high-resolution aerial imagery (ESRI, 2019). Current grazing pressure was assessed through existing grassland parcels allocation (parcels management was introduced in 2009). This map was georeferenced and manually digitized. There are 156 parcels related to 88 committee, 26 parcels are not classified. The average area of parcels is 50 sq.km. The largest area of grasslands (25%) is owned by 3 committees: Suusamyr, Leskhos and Governmental (GZZ). As abiotic factors we analysed precipitation, temperature for the growing season as well as soil organic carbon content. Table 1 lists the details and sources of the different datasets used..
Multivariate logistic regression analysis allows to estimate relationship between a dependent variable (in our case shrub encroachment) and independent variables (in our case the variables in Table 1). For the analysis we took a balanced, random, sample of 24299 pixels using a minimum distance of three pixels to reduce spatial autocorrelation. A balanced sample ensures that the same number of shrub and non-shrub locations are selected. Multicollinearity between the explanatory variables was avoided by removing those variables with a Variance inflation factor (VIF) value > 5 and Tolerance value < 0.2, as these suggest multicollinearity issues within a dataset (Bera et al., 2020;Menard, 1995; Supplementary Materials 1). The correlation matrix (using Pearson's correlation coefficient) was calculated in IBM SPSS Statistics 23 for showing the relationship between individual variables (Supplementary Materials 3). In this study temperature was excluded from the analysis due to high VIF values and strong correlation with the precipitation variable. The area under the curve (AUC) of the Receiver Operating Characteristic (ROC) was calculated to evaluate the diagnostic ability of the model (Mandrekar, 2010).
Information about potential shrub encroachment is a crucial component for appropriate grassland management. To create a map of locations of potential shrub encroachment we used the predicted probabilities of the logistic regression model for all locations using the explanatory spatial datasets (distance to the roads and paths; distance to temporal settlements; distance to permanent settlements; accessibility time; precipitation).

Spatial patterns of shrub encroachment
The RF model prepared for the classification shows good overall validation accuracy -92-94%, the kappa coefficient is 0.89-0.91 for all years, the average omission error is 7% and the commission error is 10% (confusion matrix in Supplementary Materials 5). The inclusion of slope and absolute elevation in the model increased the overall accuracy by 7-10%. The analysis shows that the overall shrub area has increased by 48% since 2000 with a rather abrupt increase in the period up to 2005 followed by a gradual increase afterwards (Fig. 5). Most of the changes occur along the edges of existing shrub patches. The area of shrub encroachment for three periods (Fig. 6) shows that the changes in shrub cover in 2005 (blue) occurred mainly in the center part of the valley, within close proximity to roads and settlements while shrub encroachment in 2016 (red) and 2019 (yellow) gravitates to more remote territories.
The distribution of shrub encroachment in relation to elevation shows a pattern of increasing altitude of areas subjected to degradation (Supplementary Materials 6). A simple logistic model between altitude and shrub encroachment shows a fit of 48.2% and an odds ratio of 1.124 for altitude. The central part of the valley (between 2200 and 2600 m) has been more exposed to shrub encroachment (Fig. 7). Shrub encroachment between 2000 and 2005 occurred mainly (68%) within the flood plain of the main river Suusamyr and its tributaries (1800-2600 m). Whilst in 2016 and 2019 it extended towards higher altitudes.
Spatial analysis demonstrates that shrub encroachment in 2005 was occurring within the grassland parcels allocated to the largest committees: Suusamyr (16.2%), in 2016 within Leskhos (17.6%) and in 2019 in Governmental (14.3%). Low-angle slopes (<3 • ) were most exposed to shrub encroachment: more than 40% of changes occurred within floodplain and adjacent lands. Spatial correspondence reflects a tendency of shrub encroachment occurring within 1 km distance from roads and within 1 h accessibility to larger residential areas and within a short distance from temporal settlements (around 45% of all changes occurred within 5-7 km distance from temporal settlements) (Supplementary Materials 4).
The logistic regression model demonstrates that the spatial distribution of shrub encroachment has a strong (negative) association with anthropogenic factors such as increasing distance from roads and temporal settlements and increasing precipitation quantity. Associations with the distance from permanent settlements, accessibility time and SOC are less pronounced in the model (Supplementary Materials 2). The fit of the model is good with a significance value of the Hosmer-Lemeshow test larger than 0.05 for all time periods, and relatively high values of the Cox and Snell R 2 and Nagelkerke R 2 (Supplementary Materials 2). The AUC of the ROC is larger than 0.80 for all years,  indicating good model fit.

Susceptibility of shrub encroachment
The probability map (Fig. 8) shows the potential areas within the valley that might be affected by shrub encroachment. The probability ranges between 0.012 and 0.994. Within areas with a probability of more than 0.5, 80% of the pixels have medium (less than 0.75) probability for shrub encroachment. High probability can be observed only within 8% of the potential shrub expansion area. The main areas susceptible to future shrub encroachment are located within low altitude 2200-2500 m (65%) and flat slopes (<3 • ) (48%), within 1 km distance from roads and temporal settlements. Only 7% of probable shrub encroachment is in remote areas and on steep slopes (>30 • ).

Discussion
The resulting spatial patterns of shrub encroachment over a period of almost twenty years allow us to better understand ongoing grassland degradation in Kyrgyzstan. This understanding can inform improved pasture management towards more sustainability by the identification of locations susceptible to future degradation. Reduction of grass cover and expansion of unpalatable grass species accompanied with shrub encroachment have been associated with grassland degradation in many studies (Andrade et al., 2015;Duarte et al., 2006;Knapp et al., 2008;Urbina et al., 2020). However, only a few studies address shrub encroachment at the local level . Most of the studies on local scale are based on long series in-situ observations (Alvarez et al., 2011;Kilyazova et al., 2019;Zhang et al., 2006), which may hinder the analysis of isolated and remote areas. Our research combines a variety of remote sensing resources, vegetation indices and DEM derivatives, allowing to assess ongoing shrub cover changes within a remote mountainous region. Our research demonstrates that 36% of shrub encroachment occurred within 500 m distance from roads and local paths and within 60 min from larger residential centers. More than 40% of shrub encroachment took place within 7 km distance from temporal settlements, affecting grass coverage in more accessible areas.
A substantial number of studies shows that grassland degradation is caused by a combination of factors: climate variations and climate change (through temperature and precipitation variations) and anthropogenic activity (Urbina et al., 2020;van Auken, 2000;Zhou et al., 2017). Some studies illustrate that human activity sometimes only plays a subordinate role, aggravating or mitigating degradation processes . This is not the case in our study, as is indicated by the importance of the distance from roads and temporal settlements in shrub appearance, reflecting the influence of grazing on grassland condition. Our analysis reveals that anthropogenic disturbances may contribute more to shrub encroachment in this area as in other grasslands around the world (Costello et al., 2000;Peng et al., 2013). Among abiotic factors, precipitation shows stronger relationship with shrub encroachment than organic carbon content of the soil (SOC). While in this area grassland dynamics may be more dependent on variations in precipitation (Craine et al., 2012;Gang et al., 2014), others have indicated that on a global level shrub encroachment might be attributed to SOC as well (Knapp et al., 2008).
In the past, Kyrgyzstan was provided with pasture management programs, offering fertilization, pest control and winter nutrition for livestock. These programs led to unbalanced pressure on the grasslands with remote pastures not used. After the collapse of the Soviet Union, the local community was forced to exploit more distant areas. However, the rugged roads and absence of necessary infrastructure impeded the usage of remote natural grasslands, imposing additional strain on grasslands adjacent to settlements, explaining some of the patterns observed in our study. Over 60% of the five million residents of Kyrgyzstan live in rural areas and are directly dependent on natural grasslands as a source of survival. Therefore, the degradation tremendously affects the local community. The population in the valley is highly reliant on grassland productivity and availability: more than 86% of the inhabitants are dependent on livestock production (UNDP, 2005). Local herdsman have used parts of the grasslands in the valley for decades, which resulted in reluctance to adopt changes in pasture management. The pastures around settlements within the main pasture committees are used more intensively due to uncontrolled grazing. Our results demonstrate that main changes in shrub encroachment have occurred in the most productive (Kyrgyzsgiprozem, 2008) parts of the valley (Suusamyr and Leskhos committees). The pastures within those committees are used actively and intensively, whilst the area of palatable plants has been declining due to shrub encroachment, imposing additional pressure on the remaining grasslands. Hence, the decreasing grassland area will rise more concerns in local community on pasture management and the ability to sustain livelihoods in the region. The estimation of future shrub encroachment probability demonstrates that new shrubs tend to clump with a high likelihood around existing shrub areas, intensifying grassland degradation. Our findings indicate that more than 60% of potential areas exposed to shrub invasion are located within the largest pasture committee which has already been subject to shrub encroachment.
Shrub encroachment is generally seen as a process of land degradation, however, at the same time some benefits can be received from shrubs. The most represented shrub is Caragana. Caragana shrub can serve as a unique habitat for pollinators and provide fuel for local communities. The characteristics of this shrub play especially an important role in mountainous area, protecting it from soil erosion and fixing nitrogen. However, these benefits mostly refer to steep slopes (more than 30% of steepness) (Kilyazova et al., 2019), whilst our findings demonstrate that only around 11% of shrub encroachment occurred on steep areas. Therefore, the larger part of shrub encroachment in our study area should be regarded as a degradation of the pastures, which may induce extra pressure on remaining areas due to a decrease in grazing area, leading to a further degradation feedback.
Overgrazing was not assessed directly in this research due to the lack of the data on livestock numbers and the use of particular parcels in each season. The current findings enriched with statistical pasture usage can provide a more comprehensive overview of ongoing degradation at the parcel level and help finding causal relations between overgrazing and shrub encroachment. Taking into account the projected increase (doubling by 2040) in cattle number (FAO, 2018) and expansion of shrubs around settlements on the most valuable pastures, it can be inferred that grazing pressure is likely to further increase. Therefore, future studies should consider including up-to-date information on livestock numbers to analyse spatial patterns of shrub encroachment in terms of grazing activity. Moreover, the results of our study may have missed some shrub patches, especially in heterogeneous areas, due to the resolution of remote sensing data. Therefore, a hot-spot assessment based on high resolution data (e.g WorldView, Pléiades constellation) can be considered in future research when evaluating shrub encroachment.
Sustainable pasture management, which takes into account spatial patterns of shrub encroachment, can mitigate ongoing degradation. Such change in management is especially needed within the buffer zones of settlements. Redistribution of income sources for people at risk of poverty can reduce the dependence on natural grasslands, alleviating grazing pressure. Therefore, the information on ongoing and potential shrub encroachment should be incorporated in the process of pastures allocation, allowing to better distribute parcels according to their productivity and capacity, and leading to more sustainable pasture management in the region.

Declaration of competing interest
The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Khazieva Elizaveta reports a relationship with VU Amsterdam Faculty of Sciences that includes: non-financial support.

Data availability
Data will be made available on request.