P á ramo Ecosystems in Ecuador’s Southern Region: Conservation State and Restoration

: The p á ramo is home to a signiﬁcant proportion of global biodiversity and provides essential services for the development of life for millions of people in Ecuador. However, land use / land cover (LULC) changes threaten biodiversity and modify its functioning. The objectives of this study were: (1) to evaluate the conservation status of the herbaceous p á ramo (HP) ecosystem by analyzing its LULC in Ecuador’s southern region. (2) to identify possible regions where the native p á ramo ecosystem is being restored. We analyzed Landsat 8 images using Object-Based Image Analysis (OBIA) and a Classiﬁer Decision Tree (CDT) to achieve these objectives. The results show that the native herbaceous p á ramo (NHP) ecosystem is being transformed into an anthropogenic HP (AHP). The area covered by the NHP ecosystem (296,964 ha) has been reduced by 50% (149,834 ha). Nevertheless, we identiﬁed ﬁve regions where the NHP is upgrading. These regions are relevant for studying NHP regeneration in Ecuador’s southern region, where soils are mostly andosols. The LU of the p á ramo , with cycles of exploitation, abandonment, and regeneration in a secondary p á ramo , is transforming the NHP ecosystem. These exploitation practices, global climate change, and lack of knowledge about the NHP ecosystem’s regeneration and its soils’ recovery threaten to substantially reduce the NHP area, its functionality, and its ecosystem services.


Introduction
The páramos are humid tropical ecosystems located from the forest's upper limit with herbaceous and shrub vegetation [1,2]. The páramos is home to a significant proportion of the Ecuadorian biodiversity and provide essential services for the development of life for millions of people in Ecuador (e.g., supplying water and storing carbon) [2][3][4][5]. These ecosystemic services are possible due to a favorable relationship between biodiversity and ecosystem functioning [6].
When the native páramo ecosystem (PE) is destroyed, biodiversity recovery can take hundreds to thousands of years. Although the number of different species may increase over time, the plant community colonizing the secondary páramo remains distinct from that of the native páramo [7]. Therefore, recognizing that the native páramo is different from the secondary páramo will improve our understanding of the relationship between biodiversity, community assembly time, and ecological functioning. Knowing these differences shed light on understanding human activities' ecological consequences and global climate change on the native and secondary páramo ecosystems. A previous study assessed the páramos of central Ecuador's conservation state based on the differentiation between natural and secondary (anthropogenic) PE [8]. Given the global importance of the páramos, the páramo ecosystems' dynamics must be understood to advance our comprehension of the páramo ecosystems' response to environmental changes and anthropogenic disturbances such as farming and agriculture.
The Andean corridor hosts the Ecuadorian PEs [9]. Their cold and humid climate has allowed several plant species to develop unique adaptive strategies to cope with reasonably hostile environmental conditions [10]. These particular conditions have led to high endemism and a remarkable diversity of plants, the highest among all alpine landscapes in the world [11]. Another characteristic of the PE is its capacity to store water due to the high content of organic matter (OM) of its soils [12,13]. The high content of OM in its soils is favored by environmental conditions such as a low intensity but long duration of annual rainfall and high temperature during the day and low temperature during the night [14,15]. The alteration of native PEs caused by land use/land cover (LULC) change is the main threat to biodiversity and a predominant modifier of ecosystem functioning. These alterations continue to increase, driven primarily by economic and demographic imperatives associated with intensive exploitation of essential resources needed for human needs [7,16]. Therefore, the restoration of degraded PE is recognized as a global priority.
There are few reports about how the intensive exploitation of the páramo soils in Ecuador's southern region has affected the páramo ecosystems [3][4][5]17,18]. Very little is known about the regeneration of native PE. Sustainable development of the páramos will fail if we do not know how the exploitation of the páramos affects these ecosystems and what factors determine the regeneration of the native PE.
In previous studies, García et al. [8] reported a methodology based on the OBIA and CART algorithm to analyze the conservation state of the páramos in the central region of Ecuador. These authors found a decision tree with six predictor variables to classify the LULC at the PE. These authors reported that the methodology is simple, fast, and can be implemented in other Ecuador regions as an exploratory tool for evaluating and planning the páramos.
The objectives of this study were: (1) to evaluate the conservation state of the herbaceous paramo ecosystem by analyzing the LULC of the soils under the páramo ecosystems in Ecuador's southern region; and (2) to identify possible regions where the native herbaceous páramo (NHP) is being restored. To achieve these objectives, we analyzed Landsat 8 images using the object-based methodology and decision tree to classify objects. We generated an updated database of LULC in the páramos ecosystems in Ecuador's southern region.

Area of Study
The study area is in Ecuador's southern region. It extends from the north of Cañar province, between the coordinates UTM 18 S, Datum WGS84 (X = 742,040; Y = 973,6413), and UTM 17 South (X = 684,253; Y = 9,463,053) to the south of Loja province ( Figure 1). The study area was segmented into 62 regions of 20 × 20 km (40,000 ha) each region was processed individually.
The study area has a diverse topography where the elevations vary from 2600 to 4480 m above sea level (m.a.s.l.), which gives rise to various ecological interactions and plant formations [19,20]. Cañar and Azuay's region has a cold and humid climate influenced by the continental air masses of the Amazon basin and the west's dry and cold air masses. The annual average precipitation is 1210 ± 101 mm per year [17]. In the provinces of Loja and Zamora, the average annual temperature varies from 7 to 25 • C, with a rainfall of 500 to 800 mm per year. The rain predominates in the afternoon, Agronomy 2020, 10, 1922 3 of 15 and part of it is of low intensity, known as a drizzle. The Andean mountain range determines the precipitation regime of the entire study area. Precipitation is low but constant throughout the year, so paramos' ecosystems play the role of water regulators by continually providing the water produced by rain. This water flows down the slopes and riparian areas [21,22]. The study area has a diverse topography where the elevations vary from 2600 to 4480 m above sea level (m.a.s.l.), which gives rise to various ecological interactions and plant formations [19,20]. Cañar and Azuay's region has a cold and humid climate influenced by the continental air masses of the Amazon basin and the west's dry and cold air masses. The annual average precipitation is 1210 ± 101 mm per year [17]. In the provinces of Loja and Zamora, the average annual temperature varies from 7 to 25 °C, with a rainfall of 500 to 800 mm per year. The rain predominates in the afternoon, and part of it is of low intensity, known as a drizzle. The Andean mountain range determines the precipitation regime of the entire study area. Precipitation is low but constant throughout the year, so paramos' ecosystems play the role of water regulators by continually providing the water produced by rain. This water flows down the slopes and riparian areas [21,22].
The geology of the study area is dominated by the Tarqui Formation, occupying 51% of the area (Agglomerates, tuff agglomerates, and tufas of intermediate to acid composition), while 10% is occupied by the Zamora series (Sandstones, siltstones, and shales), 8% by the Saraguro Formation (Pyroclasts: tufas and thick agglomerates with lava blocks and lavas: andesites with alternating porphyritic texture), 4% metamorphic rocks, 5% glacial deposits, and the remaining 22% is dominated by a variety of geological units, each with less than 3% representation in the paramo ecosystems in the area of study [18].
Most of the soils under páramo located in the south are classified as Andosols [3]. Regardless of their parent material, páramo's soils have very dark epipedons consisting of Ah and A horizons with little differentiation, light overlying colored subsoil (C horizon) rich in minerals (dominated by clay) [17]. Among the most important soil properties are high organic C content (170-200 g/kg), very low bulk density (400 kg/m 3 ), and absence of allophane. The soils have developed strong hydric properties, with water retention at −1500 kPa of 1.4 g/g [18]. The geology of the study area is dominated by the Tarqui Formation, occupying 51% of the area (Agglomerates, tuff agglomerates, and tufas of intermediate to acid composition), while 10% is occupied by the Zamora series (Sandstones, siltstones, and shales), 8% by the Saraguro Formation (Pyroclasts: tufas and thick agglomerates with lava blocks and lavas: andesites with alternating porphyritic texture), 4% metamorphic rocks, 5% glacial deposits, and the remaining 22% is dominated by a variety of geological units, each with less than 3% representation in the paramo ecosystems in the area of study [18].

The Herbaceous Páramo
Most of the soils under páramo located in the south are classified as Andosols [3]. Regardless of their parent material, páramo's soils have very dark epipedons consisting of Ah and A horizons with little differentiation, light overlying colored subsoil (C horizon) rich in minerals (dominated by clay) [17]. Among the most important soil properties are high organic C content (170-200 g/kg), very low bulk density (400 kg/m 3 ), and absence of allophane. The soils have developed strong hydric properties, with water retention at −1500 kPa of 1.4 g/g [18].

The Herbaceous Páramo
The HP is characterized by species of the genera Calamagrostis, Agrostis, Festuca, Cortaderia, Stipa, and an abundant diversity of lifeforms [20]. In south-central Ecuador, in the Cajas páramo, the plant community varies as humidity conditions drop, and associations are created between Calamagrostis sp. and Viola humboldtii. In the southern part of Ecuador, the paramo ecosystem descends to 2600 m.a.s.l. and is composed of Agrostis breviculmis, Calamagrostis spp., Festuca asplundii, and Stipa ichu; in areas with Agronomy 2020, 10, 1922 4 of 15 steep slopes greater than 70% [2]. In lower altitude areas, the HP ecosystem is usually in association with shrub species.
The Ministry of Environment of Ecuador (MEE) defines the HP according to vegetation types, precipitation, and soils. The HP covers 73% of the study area, where the HP ecosystem (61%) and the HP ecosystem with evergreen shrubs (12%). A semi-deciduous shrub ecosystem covers 15% of the area. Different ecosystems cover the remaining 12% of the area [20]. For the sake of comparison, Figure 5a shows the HP páramo ecosystem distribution reported by the MEE [20], which served as a baseline. The methodology implemented by the MEE [23] includes six environmental variables: bioclimate, ecological floors, flood regime, relief, phenology, and biogeography. The bioclimatic model was built from the variables bioclimate, ambrotypes, thermotypes, and ecological floors. The flooding and fluvial origins are combined with the general relief map, macro, and mesorelief to give the geoforms model. The ecosystems' map was built from the bioclimatic model, geoforms model, phenology, biogeography, and land cover (natural vegetation) map. Landsat and Rapideye images acquired in 2007 by the MEE were used to generate the land cover map (natural vegetation). These images were pre-processed and corrected. A supervised classification was performed; the resulting environmental units were verified in the field. The vegetation map used the classes defined by the ecosystem type developed by the ecosystem classification system for continental Ecuador [20]. Finally, the ecosystem map was subjected to a process of (1) spatial debugging to eliminate spurious or unrepresentative patterns, considering a minimum mapping unit of 25 ha., and (2) verification of thematic consistency with floristic data and expert criteria.

Workflow
The methodology comprised four main stages: (1) image pre-processing, implementation of OBIA, (2) attribute extraction, (3) validation and performance evaluation of the decision tree, and (4) land cover and land use classification in the study area, map generation, database, and analysis ( Figure 2).

Landsat-8 Images
The Landsat 8 images were downloaded from the official NASA website, "earthexplorer.usgs.gov." Cloudiness of less than 30% was considered acceptable to make the visible field as complete as possible. The date of the downloaded images was 20 November 2016, and the grid number on the map (path = 10/row = 62 for the South Central zone and path = 10/row = 63 for the South zone).
The study area contemplated the use of two satellite images with the following labels: (1) Image pre-processing included downloading Landsat 8 images, panoramic sharpness, and brightness, and reflectance correction. The implementation of OBIA included image segmentation Agronomy 2020, 10, 1922 5 of 15 and object generation. Attribute extraction included extracting the maximum, mean, and minimum reflectance values, spectral indices, and ancillary data for each object. Decision tree validation involved sampling 610 random objects categorized by the decision tree and validating them using the Google Earth platform. Classification and analysis included a Python code that allowed us to classify each object according to the category and generate a map with the different LULC categories.

Landsat-8 Images
The Landsat 8 images were downloaded from the official NASA website, "earthexplorer.usgs.gov." Cloudiness of less than 30% was considered acceptable to make the visible field as complete as possible. The date of the downloaded images was 20 November 2016, and the grid number on the map (path = 10/row = 62 for the South Central zone and path = 10/row = 63 for the South zone).
The study area contemplated the use of two satellite images with the following labels:

Categorization of Land Use and Land Cover
Land use and the cover were categorized into five classes: NHP, anthropogenic herbaceous páramo (AHP), forest (FRS), water bodies (WTR), and snow (SNW). Figure 3 shows the classes of LULC based on visual observation of the landscape and imagery: Water-watercourses and water bodies, including rivers, small lakes, and reservoirs. The water class includes the lake system and adjacent areas with high susceptibility to flooding. Forest-forest, and shrubs. The forests in the province of Cañar y Azuay include forests and shrubs located in low windy and humid areas. They are dominated by the genus Polylepis with heights ranging from 5 to 7 m. The forests in the provinces of Loja y Zamora are mainly evergreen shrubs that reach heights up to 3 m and found from 2800 to 3000 m.a.s.l.; the NHP is dominated by Calamagrostis spp. In regions well to the south, the NHP is associated with shrub species, extending from 2800 to 3900 m.a.s.l. The AHP is the páramo that has been transformed by humans. It includes farmland, pasture, bare soil, built/infrastructure areas, roads, burned areas, and regenerated grasslands with other species' association and exhibits alterations in its phenology. This class is in areas with altitudes lower than 3800 m.a.s.l. and slopes lower than 11%. They are buffer zones and human settlements that allow the transition from Andean forests to crops, plantations, or infrastructure. The snow class corresponds to the areas with an altitude greater than 3827 m.a.s.l., with small deposits of ice crystals giving the highest values of reflectance.
been transformed by humans. It includes farmland, pasture, bare soil, built/infrastructure areas, roads, burned areas, and regenerated grasslands with other species' association and exhibits alterations in its phenology. This class is in areas with altitudes lower than 3800 m.a.s.l. and slopes lower than 11%. They are buffer zones and human settlements that allow the transition from Andean forests to crops, plantations, or infrastructure. The snow class corresponds to the areas with an altitude greater than 3827 m.a.s.l., with small deposits of ice crystals giving the highest values of reflectance.

Object-Based Image Analysis
The images of each scene were processed in ArcGIS software 10.3, Redlands, California 92373 8100, USA (http://www.esri.com/software/arcgis). We followed the workflow suggested by Urbanski [25] to convert Landsat-8 OLI images into a land cover map where the fundamental units are geographical objects. In total, 1,876,324 geographical objects with common attributes were obtained for the 62 regions under study.
A map of LULC in the study area was generated using the CDT shown in Figure 4. This CDT served as an exploratory tool for data-based discovery and prediction and gaining new knowledge about LULC under the paramo ecosystems in Ecuador's central region [8]. The CDT requires only six predictor variables (reflectance, VARIg, altitude, EVI2, NBR2, and slope). In this study, we used this CDT to predict the LULC in the páramo ecosystems in Ecuador's southern region.

Object-Based Image Analysis
The images of each scene were processed in ArcGIS software 10.3, Redlands, California 92373 8100, USA (http://www.esri.com/software/arcgis). We followed the workflow suggested by Urbanski [25] to convert Landsat-8 OLI images into a land cover map where the fundamental units are geographical objects. In total, 1,876,324 geographical objects with common attributes were obtained for the 62 regions under study.
A map of LULC in the study area was generated using the CDT shown in Figure 4. This CDT served as an exploratory tool for data-based discovery and prediction and gaining new knowledge about LULC under the paramo ecosystems in Ecuador's central region [8]. The CDT requires only six predictor variables (reflectance, VARIg, altitude, EVI2, NBR2, and slope). In this study, we used this CDT to predict the LULC in the páramo ecosystems in Ecuador's southern region.

Predictor Variables
The predictor variables were recovered from the Landsat 8 bands and the DEM. A total of 6 predictor variables were used, including the necessary spectral information on reflectance, three spectral vegetation indices (SVI), and two indicators derived from DEM (Table 1). Band algebra was used through the ENVI Math Band tool to generate the detailed indices according to their

Predictor Variables
The predictor variables were recovered from the Landsat 8 bands and the DEM. A total of 6 predictor variables were used, including the necessary spectral information on reflectance, three spectral vegetation indices (SVI), and two indicators derived from DEM (Table 1). Band algebra was used through the ENVI Math Band tool to generate the detailed indices according to their corresponding Landsat 8 bands. The SVI were used as predictor variables because they are the primary source of information for monitoring and evaluating land cover [8,12,26]. Table 1. Features used for classification [8].
Lineally sensitive to vegetation fraction, it exhibits a good correlation with nitrogen contents [27].
Provide greater sensitivity in regions with high biomass while minimizing the soil's influence and atmosphere [28].

Decision Tree
Although other classifiers could achieve better classification results, the CDT represents a clear set of rules with values that define thresholds used in similar contexts without training [30]. Cao et al. [31] reported that decision trees could provide stable performance and reliable machine learning and data mining research results. Thus, decision trees have been used to identify spectral bands with the highest class discrimination capabilities and low misclassification rates [31]. Decision trees have also been used to identify greenhouse crops from remote sensing data [30]. However, in general, CDT rules are restricted by the root variable. In this study, only rules that are indicated by "Reflectance" were extracted. The CDT in Figure 4 exhibited excellent performance in the learning and validation phase for the five categories used in the characterization of LULC in the páramo ecosystem in Ecuador's central region [8]. In the CDT's performance, the predictor variables' relevance was 100% for the reflectance, 63% for the VARIg index, 56% for the altitude, 46% for the EVI2 index, 34% for the NBR2 index, and 23% for the slope.

Validation of the Classification Decision Tree
The use of the CDT in Figure 4 was validated in Ecuador's southern region using 610 random objects. Among the 610 objects, 50 belong to the category WTR, 100 to FRS, 225 to NHP, 225 to AHP, and 10 to SNW. The allocation of the number of samples per category was proportional to the total number of objects in each category. We choose 610 objects for the sake of practicality and statistical reasons. In reality, it could be any number larger than 385, the number of samples statistically representative of a population of 1,876,324 objects, with a confidence level of 95% and a confidence interval of 5%. The allocation of samples per category was proportional to the total number of objects in each category. The classification results by the CDT were verified with the Google Earth platform's information. Table 2 shows the table of confusion obtained in the validation of the CDT in the study area. The CDT performance was evaluated by a binary crosstabulation derived from the confusion matrix (Table 3).  The CDT's overall accuracy was above 94%. The misclassification rate was below 6%. The CDT's specificity (how well the CDT can recognize negative samples) was above 94%. The probability of an informed classification for each class was above 86%. How much each class is marked as a predictor or possible cause of another was above 70%. Matthew's correlation coefficient between observed and predicted classification was above 84%. Therefore, the percentage of correct categorizations of the 610 known ground objects was 91% (Table 2). This value is acceptable if we consider the heterogeneity of the study area. In addition, there is little information on land use and land cover under paramo ecosystems in the southern region of Ecuador. Therefore, the CDT can be used as an exploratory tool to classify and differentiate LULC objects under the páramo ecosystem in Ecuador's southern region. The CDT allowed the mapping of geographical objects' input space into a predefined class with good performance. There are 14 regions (4,9,11,13,24,25,27,31,38,42, 51, 57,58, and 59) where the MEE map and our analysis agree that there is no area covered by the NHP ecosystem. Meanwhile, in regions 1, 3, 6, 10, and 16, the MEE map reports no NHP, and our results report the existence of 42, 837, 133, 517, and 539 ha of NHP, respectively. This result suggests that the CDT identifies that the HP ecosystem is upgrading in these regions. Similarly, in regions 12 and 26, the MEE reports 1021 and 108 ha of NHP, and our results show that in these regions, there is no area covered by NHP. On the other hand, the CDT identified 23 regions (18,22,23,28,29,32,33,34,35,36,40,41,43,44,46,48,49,50,53,54,55, 60, and 61) where the surface area covered by the NHP ecosystem reported by the MEE has decreased by more than 2000 ha. This result suggests intensive land use that has led to NHP degradation; for example, regions 29 and 18 have lost 15,035 and 14,616 ha, respectively. In regions 2, 3, 5, 10, 16, and 17, the CDT identifies an increase of more than 500 ha in the surface area covered by the NHP ecosystem concerning the values reported by the MEE. Particularly in region 17, the CDT finds an increase of 3042 ha under the NHP. This result suggests that in these regions, the NHP is regenerating.

Results
In summary, the sum of all positive differences between the NHP surface's value achieved with our analysis and the MEE value suggests a total gain of 7607 ha. The sum of all negative differences between the total area given by OBIA and the reported by MEE in each region suggests a total loss of 154,618 ha. In contrast, regarding the MEE report, the results suggest a total net loss of 147,129 ha (296,963-149,833) of NHP. However, the MEE only reported the existence of the NHP ecosystem. The MEE reported in 2013 [20] the existence of 296,964 ha of NHP. There is no report for the AHP surface area. In our study, we found that the surface area covered by NHP, AHP, FRS, WTR, and SNW was 149,834 ha (39%), 146,829 ha (38%), 61,404 ha (16%), 27,916 ha (7%) and 726 ha (0.2%), respectively. Table 4 shows the distribution of LULC for each region in the study area. The regions with the largest area of NHP, AHP, FRS, WTR, and SNW are region 47, with 27,377 ha, region 40 (12,431 ha), region 29 (5152 ha), region 22 (3210 ha), and region 46 (313 ha), respectively (Table 4). Each region has an area of 40,000 ha. Therefore, 74% of the surface of region 47 is covered by NHP. Region 40 exhibits the largest anthropogenic activity area, with 35% of its surface area of AHP. Forests cover 33% of the surface area of Region 29. Water covers 19% of Region 22, and 3% of Region 46 is covered by snow.
There are 14 regions (4,9,11,13,24,25,27,31,38,42, 51, 57,58, and 59) where the MEE map and our analysis agree that there is no area covered by the NHP ecosystem. Meanwhile, in regions 1, 3, 6, 10, and 16, the MEE map reports no NHP, and our results report the existence of 42, 837, 133, 517, and 539 ha of NHP, respectively. This result suggests that the CDT identifies that the HP ecosystem is upgrading in these regions. Similarly, in regions 12 and 26, the MEE reports 1021 and 108 ha of NHP, and our results show that in these regions, there is no area covered by NHP. On the other hand, the CDT identified 23 regions (18,22,23,28,29,32,33,34,35,36,40,41,43,44,46,48,49, 50, 53, 54, 55, 60, and 61) where the surface area covered by the NHP ecosystem reported by the MEE has decreased by more than 2000 ha. This result suggests intensive land use that has led to NHP degradation; for example, regions 29 and 18 have lost 15,035 and 14,616 ha, respectively. In regions 2, 3, 5, 10, 16, and 17, the CDT identifies an increase of more than 500 ha in the surface area covered by the NHP ecosystem concerning the values reported by the MEE. Particularly in region 17, the CDT finds an increase of 3042 ha under the NHP. This result suggests that in these regions, the NHP is regenerating. In summary, the sum of all positive differences between the NHP surface's value achieved with our analysis and the MEE value suggests a total gain of 7607 ha. The sum of all negative differences between the total area given by OBIA and the reported by MEE in each region suggests a total loss of 154,618 ha. In contrast, regarding the MEE report, the results suggest a total net loss of 147,129 ha (296,963-149,833) of NHP. However, the MEE only reported the existence of the NHP ecosystem.

Discussion
This study found that the NHP ecosystem is being transformed into AHP. The results show that the NHP ecosystem (296,964 ha) has been reduced by 50% (149,834 ha). However, the rate of changes in the extent of NHP and AHP cannot be clearly defined. There is a strong possibility that the MEE classified all paramos (NHP and AHP) together as one unit and reported it as NHP. Similarly, we do not rule out the possibility that the 50% loss in NHP occurred over a period that started before the baseline MEE map. Nevertheless, this result agrees with what several authors have suspected [2,32]. Among the factors that may be contributing to the transformation of NHP are (1) intensive land use or expansion of anthropogenic activity, (2) global climate change, and (3) limited knowledge about NHP regeneration, among others.

Discussion
This study found that the NHP ecosystem is being transformed into AHP. The results show that the NHP ecosystem (296,964 ha) has been reduced by 50% (149,834 ha). However, the rate of changes in the extent of NHP and AHP cannot be clearly defined. There is a strong possibility that the MEE classified all paramos (NHP and AHP) together as one unit and reported it as NHP. Similarly, we do not rule out the possibility that the 50% loss in NHP occurred over a period that started before the baseline MEE map. Nevertheless, this result agrees with what several authors have suspected [2,32]. Among the factors that may be contributing to the transformation of NHP are (1) intensive land use or expansion of anthropogenic activity, (2) global climate change, and (3) limited knowledge about NHP regeneration, among others.
Intensive land use has modified a significant part of the páramo ecosystems in the study area. In 23 of the regions studied, the NHP ecosystem area decreased by more than 2000 ha. This decrease is because the paramos' inhabitants, to exploit the soils of the páramo, introduce replacement systems of the NHP ecosystem ( Figure 6). In the páramos, there is a wide variety of production systems with tremendous social and environmental impacts on soil, water quality, and social-environmental dynamics [34]. Humaninduced transformations, such as fire and clearing activities, have decreased the area covered by páramo ecosystems and increased their fragmentation [35]. In the páramos, there is a wide variety of production systems with tremendous social and environmental impacts on soil, water quality, and social-environmental dynamics [34]. Human-induced transformations, such as fire and clearing activities, have decreased the area covered by páramo ecosystems and increased their fragmentation [35].
The transformation of the natural paramo due to LULC changes affects soil carbon retention. When soil is cultivated, the A horizon is degraded, and the layer of organic remains disappears due to erosion [36]. When the soil is unproductive, it is abandoned and left to "recover," leading to a secondary herbaceous páramo ( Figure 6) with unpredictable results. However, very little is known about the recovery of the NHP and the páramo's soils.
Although there may be a wide range of soil rehabilitation recommendations, it is essential to conduct a complete ecosystem analysis before site rehabilitation. The analysis must integrate the description of the developed soils and the slow pedogenic processes [16]. Thus, the soils' study under páramo is indispensable to understand the páramo ecosystem's functioning and optimize its ecological rehabilitation efficiently.
The páramos are among the most vulnerable ecosystems to global climate change. It is expected that the páramo temperature change will be lower than in other lower altitude regions. However, a 10% reduction in annual rainfall is expected in the páramo by 2040. The change in temperature can cause an altitudinal displacement of a vast number of species and, with them, of the bioclimatic band. This temperature change may reduce the surface covered initially by páramo ecosystems by 15 to 25% [34]. The growing transformation of the NHP, the increase in temperature, and the decrease in humidity reduce taxonomic and functional diversity, and the functioning and supply of ecosystem services that the páramo provides.
Soils in the southern Ecuadorian páramos are burned after three to four years of exploitation and left to recover. The burning of vegetation to grow grass is a widespread management practice in the region [4]. When burned and abandoned, the secondary HP ecosystem (AHP) regenerates through secondary succession, with limited functionality and reduced taxonomic diversity ( Figure 6) [2,33,35]. There is little information on how soil properties under paramo ecosystems impact herbaceous and shrub species' natural regeneration. However, this secondary HP may or may not lead to an NHP. The difficulty in the NHP regeneration is favored by the alteration of the pedogenesis of its andosol soil [12,13,18]. In some cases, the exploitation has been so intensive that only bare soil remains. In this condition, the ecosystem demands an active restoration with direct intervention on the degraded páramo ecosystem's structure and characteristics to overcome stressors that prevent regeneration and ensure recovery [2,7,16,18,32,34]. However, few know how to or even if it is possible to regenerate the NHP, the resilience of NHP, and the resilience of andosols to continue supporting the NHP ecosystem. There are few studies on the páramos' soil regarding the impact resulting from exploiting the páramos or their restoration after exploitation [2,4,18,32].
In our study, the CDT identified five regions where the NHP ecosystem is upgrading. In particular, region 17 shows an increase in the surface covered by HPN in 3000 ha. The CDT used spectral indices of vegetation as biophysical indicators of the LULC. This outcome suggests that in region 17, the NHP is regenerating. Few studies on the recovery of paramo ecosystems after the human intervention have ceased, limiting the implementation of restoration actions [5]. In páramo ecosystems, the geomorphological position regarding LU and management practices conditions the availability of carbon in the soil. It affects pedogenetic processes that compromise the regeneration of NHP. However, pedogenetic processes are prolonged and even blocked by specific environmental parameters [16]. The processes involved in the evolution of intensively exploited soils in the páramos have not yet been addressed. The nature of pedogenetic processes in the páramo remains unexplored. All of the NHP ecosystem's compartments must be considered, particularly the soil, to know whether to intervene after LU ceases and restore the ecological communities' functioning in the long term.

Conclusions
The intensive LU of the páramo with cycles of exploitation, abandonment, and regeneration in a secondary páramo transforms the NHP ecosystem. These exploitation practices, global climate change, and the lack of knowledge about the regeneration of the ecosystem and recovery of its soils threaten to substantially reduce the ecosystem's area, functionality, and provision of the NHP ecosystem services. Analysis of the LULC database in the páramo ecosystem of Ecuador's southern region shows that 39% (149,834/386,709) of the total study area remains as a native HP (NHP). However, 38% (146,829/386,709) have been transformed into anthropogenic HP (AHP), and 23% ((61,404 + 27,916 + 726)/386,709) remains FRS, WTR, and SNW. Overall, 386,709 ha is the sum of all classes' area determined in our study (see Table 4). Regarding the MEE report, 50% (149,834 ha) of the surface area of the páramo ecosystem (296,964 ha) remain as a native HP. Among the factors that may be contributing to the transformation of NHP are (1) intensive land use or expansion of anthropogenic activity, (2) global climate change, and (3) limited knowledge about NHP regeneration, among others.
However, we found evidence of regeneration of the NHP, indicating its resilience. Additional studies to integrate new technological tools based on artificial intelligence algorithms, evolutionary processes, functional diversity, community dynamics, and ecosystem services of the NHP and AHP are necessary for exploring NHP's response to global climate change and all páramo communities' adaptability to such change. The future of NHP depends on the study of the factors that favor NHP regeneration. Thus, the páramos' conservation depends on the new knowledge we can acquire by studying areas where NHP regeneration has been identified. This study will be primarily useful to scientists and decision-makers interested in páramo ecosystems in Ecuador's southern region.