The selected area is the district of Marangani (Canchis, Cusco), comprised within the Raya Mountain range (Figure 1). It is the third most populated district of Canchis, and the population at the end of 2017 was 10,554 (National Institute of Statistics and Informatics, 2022). It covers an area of 440.32 km², comprising 29 watersheds that range from 70 to 0.5 km² with an average of 15 km². This landscape is dominated by an expansion of livestock breeding in the upper lands and farming in the fertile lowlands. It is typical of many mountain agroecosystems across the world. In that sense, Marangani has natural grasslands, sparsely vegetated areas, glaciers upstream, bofedales in the mid-stream, and agricultural land in the downstream reaches. The most important river is Vilcanota. The district has a semi-frigid rainy season with a dry winter climate (58%). Almost 47.63% of the population mainly works in the agricultural sector (National Institute of Statistics and Informatics, 2022), characterized by annual crops such as alfalfa, bean, potato, and wheat (National Institute of Statistics and Informatics, 2012). The natural pastures present in the territory make an optimal fodder for camelids, sheep, and cattle.

To estimate the ESV, we used the benefit transfer method, transferring available data from previous studies in various locations and from different ecosystems globally that have similitudes with the biomes of the study area. The search was carried out on the ESV database (Foundation for Sustainable Development, 2021) and published research (Costanza et al., 2014; Xie et al., 2017; de Groot et al., 2020) for 38 ecosystem services and ten biomes. After the search, only values for 23 ecosystem services were found. According to Costanza et al., (1997), the ESV coefficient (unit values in USD ha⁻¹ yr⁻¹) of a given LULC was counted as all the ESV for that LULC. To achieve the total ecosystem service value of each LULC at the municipality and the watershed level, we multiply the area of each LULC type at the corresponding level by the ESV coefficient. Finally, the ESV coefficients were adjusted using the Consumer Price Index (CPI) obtained from the website of the World Bank (2022) and applying the inflation rate formula (Table 3).

The spatiotemporal distribution of LULC of the Marangani district among the study years revealed some variations (Figure 2). The precise area of each LULC is shown in Table 4. Natural grasslands coincided with being the most abundant class each year, spatially dispersed covering more than 65% of the landscape. Sparsely vegetated areas delimited more than 15% of the landscape over time. This LULC is spatially associated with glaciers in the highland. The third LULC type was peat bogs and high-Andean wetlands (bofedales) with an average of 8%, mainly located in the upland depressions. Agricultural areas (>4%) were evident in the lowland valley around the main river (Vilcanota river) related to the urban center, the farmhouses, and even topography. However, shrublands mainly drew a transitional zone in the lowland that separated agriculture from the natural grassland. Last, built-up areas, forest plantations, and water surfaces enclosed less than 1% of the territory.

The details of LULC changes for 1986–1995, 1995–2007, and 2007–2019 are presented in the transition matrices (Supplementary Tables S3–S5). In terms of the total area, 4,217 ha (9.6%), 3,999 ha (9.1%), and 3,914 ha (9%) were transformed in each period, respectively. Two transitions (>2,900 ha) involved the deterioration and the recovery of natural grasslands and bofedales during each period. A third important transition was related to glacier retreat, being more intense during the first period (543 ha). Another significant land-use change over time (strongest during the first period, >300 ha) refers to agricultural expansion, which implicated the conversion of natural grasslands, bofedales, and shrublands. Last, the farming land abandonment represented an increase in shrublands, mainly shown during the third period (>200 ha).

The search of 38 ecosystem service values for the ten study area biomes resulted in 23 ESVs found (Supplementary Table S7). The values (Supplementary Table S8) focused on provisioning (5), regulating (12), and cultural services (6) that are classified in the MEA (Millennium Ecosystem Assessment, 2005). This last group of services provides the highest total economic value (74,260.98 USD ha⁻¹ yr⁻¹), mainly due to the value of recreation (57,248.23 USD ha⁻¹ yr⁻¹), while regulating and provisioning afford lower values (35,356.34 and 33,124.41 USD ha⁻¹ yr⁻¹), being the value of regulation of extreme events and the provision of habitat for wildlife species (9,634.75 and 18,444.80 USD ha⁻¹ yr⁻¹) those who contribute most in each category. Bofedales and water surfaces are the LULCs that supply the highest ESV (∼33,000 USD ha⁻¹ yr⁻¹ each), mainly because of the value of the habitat for wildlife species, erosion prevention, regulation of extreme events, and recreation. However, human-related areas retain values that range between 3,211.69 and 6,687.02 USD ha⁻¹ yr⁻¹ due to the provision of food production, recreation, and air quality regulation.

At the municipality level, almost all the LULCs faced variations over the study period of 33 years (Table 4). It should be noted that glaciers—as fragile ecosystems perturbed by climate change—had the highest change (accumulated ratio of −85.51%), which agreed with studies in the Himalayas (e.g., −23.8% and −35%, Rai et al., 2018; Shrestha et al., 2019). However, Andean glaciers seem to be under critical conditions that accelerate their disappearance (Juřicová and Fratianni, 2018; Vicente-Serrano et al., 2018). In the Marangani district, this situation reached its peak between the first period (1986–1995, ratio of −60.65%). It can be explained by the very strong ENSO warm events between 1987 and 1992 (Gergis and Fowler, 2009), which are related to the annual glacier mass balance (Wagnon et al., 2001). Consequently, the current ice loss in the study area signifies a variability of water availability on downstream socio-ecosystems that may impact agricultural irrigation of permanent crops and human-domestic use. In that sense, the cropland and the bofedales have shown a moderate increasing trend over the first period that may be driven by the increase of the population (National Institute of Statistics and Informatics) and the water supply from glaciers, respectively, while for the second and third periods, agricultural frontier expansion showed a deceleration influenced by a decrease in the population. This decline overlaps with the displacement by domestic terrorism in the second period (Transfer Commission, 2014) and internal migration by work opportunities in the third (Bergmann et al., 2021). However, existing farmland abandonment is very low, but if this land-change dynamic continues, it can lead to risks associated with land degradation (Narendra Raj and Teiji, 2006). Furthermore, during the last two periods, the bofedales faced a reduction coexisting with a decline of the melting of glaciers, causing augmentation of natural grassland, sparsely vegetated areas, and the conversion by humans at work activities as agricultural land and forest plantation.

At the watershed level, the study observes a higher tendency of land changes in the territory without glaciers, and those with glacier areas count on permanent larger bofedales. Then, glaciers provide a water flow that stabilizes the fluctuations (appearance balance) of the bofedales. This biome is considered sensitive to changes in extreme rainfall and glacial meltwater supply (Dangles et al., 2017). The main reason is that perennial bofedales are connected to the hydrological cycle through runoff generated by melting glaciers, while temporal wetlands depend directly on precipitation (Otto et al., 2011). It suggests that the total disappearance of the glaciers in the study area would cause the main water supply contribution by snow melting and rainfall. In this case, the mountain depositional landforms—glacial rocks immersed in the biome of the sparsely vegetated area—will constitute effective water distribution and storage. This function has an operation in the short and intermediate-term (Reato et al., 2021). It may cause a predominance of temporal bofedales. This condition could have a possible variation in the benefits provided by ecosystem services in the same biome for the different compositions of species. In general, temporal bofedales (mainly species of the Poaceae family) provide food for raising animals and wildlife during the flood season, whereas perennial bofedales (mainly cushion-forming species) accumulate organic soil, forage, and water storage also in the dry season (Cooper et al., 2010; Ruthsatz, 2012). In this last service, the hydrological importance enables change in the water flow paths and discharge velocity (Otto et al., 2011). In addition, perennial wetlands provide vital habitats for wildlife and extensive grazing of camelids (Loza Herrera et al., 2015; Dangles et al., 2017).

For the first time, the results presented here (Table 5) calculated the spatiotemporal changes of ESV in an environment of tropical glacier retreat. Concretely, the territory is on a trajectory of losing ESV appreciated in the scenario without inflation adjustment. This degradation is mainly due to reduction of provisioning and regulating services provided by the bofedales and the glaciers. The disappearance of these services not only has implications for a direct economic loss but also, in the case of the reduction of the flood protection provided by the glaciers, can pose the risk of dangerous glacial outburst floods on downstream communities that may lead to catastrophic economic damages (Chen et al., 2022). At the same time, the disappearance of bofedales endangers the regulation of extreme events that can increase these impacts. Likewise, the augmentation of water erosion, induced by the decline of bofedales, may affect negatively food security and natural resource management practices, principally in the poorest societies (Sartori et al., 2019), such as the study area. In the same way, the drop in water supply amplified by the reduction of both fragile ecosystems may affect livestock production, irrigated agriculture, availability and quality of water for human consumption, and generation of hydroelectricity (Carey et al., 2016), increasing the economic losses in the territory. Biophysical processes and human activities are the highest threats and causes of this situation. The ESV supplied by the human-related classes is in a crescendo, whereas the natural and semi-natural classes are losing their capital. This may result in a critical situation since the growth of provisioning services provided by agricultural development is at the expense of regulating services like water flow regulation and erosion prevention. This trajectory of ecosystem services is characteristic in mountain environments (Locatelli et al., 2017; Madrigal-Martínez and Miralles i García, 2019a).

At the watershed level, the variation of ESV is mainly influenced by the changes of bofedales into farming and grasslands. These changes harm the value of regulation services. However, the value of provisioning services like habitat for wildlife species, genetic libraries, livestock production, and raw materials increased due to expansion of pastures. In general, the economic value of these provisioning services could be underestimated by the global sources researched by Costanza et al. (2014) because the principal cattle in the study area are llamas, alpacas, and vicuñas (National Institute of Statistics and Informatics, 2012). These two last species are considered of great interest for their finest fibers worldwide and are used for luxury costumes (Vilá and Arzamendia, 2020; Zarrin et al., 2020). In the same way, the production of meat from llamas and alpacas has interesting value for their use as a principal protein source for Andean rural communities (Pérez et al., 2000); and their healthy characteristics (low in fat and cholesterol) make them attractive for local and international markets, signifying an important income for small- and medium-scale local producers (Mamani-Linares and Gallo, 2014). Moreover, llamas have smaller and simpler nanobodies beneficial for human immunology and therapeutic applications (Wesolowski et al., 2009). Furthermore, although the density of alpacas in the study area is high (123 alpacas km⁻²), according to the levels established in other locations (Muñoz et al., 2015), biodiversity conservation is guaranteed with appropriate management practices (Alkemade et al., 2013).

In the watersheds with glacier retreat, the significant quantity of bofedales and natural grasslands controls the fluctuations of ESV. It is a consequence of the extension of grasslands that varies between 48 and 72.2% and the bofedales between 5.5 and 12%. Then, these watersheds still have an advantage from the water flow of glaciers and supply multiple provisioning and regulating services such as water yield, habitat for wildlife species, food production, carbon sequestration, erosion control, soil fertility, and regulation of extreme events. Therefore, these high-mountain watersheds play an essential role as hotspots of ecosystem services, giving benefits and value to local communities. The identification of ES hotspots guides the prioritization of areas for conservation (Li et al., 2017). In that sense, ES hotspot maps have potential benefits that can be transmitted to stakeholders to improve management planning (Hauck et al., 2013). Also, the implication of knowing this information provides a scientific basis for making relevant decisions in a territory (Li et al., 2022). Consequently, these ES hotspot watersheds should get the priority of being strategically managed for maintaining their essential ecosystems. So, understanding the relationship between nature and people in mountain socio-ecological systems is significant for planning sustainable development strategies (Madrigal-Martínez and Miralles i García, 2019b; Payne et al., 2020), including conservation strategies at different spatial levels, such as prioritizing the small wetlands of these watersheds, according to Otto and Gibbons (2017).

In this study, the analyses presented should be understood as using the best existing data of acceptable quality to admit robust evidence. Even so, the procedures used in processing satellite images may have a situation of some biome coverage. It is mainly related to four factors: 1) the definition of biome classes, 2) reliability of Landsat images with their limitations of the sensor used, 3) auxiliary variables in the classification, and 4) validated data and no clarity of training (Calderón-Loor et al., 2021). The selection of the training points was based on the national vegetation coverage map (Ministry of Environment, 2015b). Also, to best interpret the satellite image, the authors made an expert judgment. In that sense, subjectivity at the decision point regarding the pixel classification could cause statistical noise in the training and validation data. However, this noise in the final categorization could have a low effect on the final product (Zhu et al., 2016). In addition, the research applied the transfer method of benefits for valuation ecosystem services. Despite the benefits, it has several limitations that imply transfer errors from generalization (Schmidt et al., 2016). However, it has been taken into account that the extrapolation of the valuation is in similar geographies.