Drivers and Spatial Patterns of Burned Area in High-Andean Páramos

1. Introduction

Páramos are high-mountain tropical ecosystems located primarily in the Andes. Occupying the transition zone above dense tropical mountain cloud forest until reaching glaciers, ranging from altitudes between 3000 and 4800 m [1]. Páramos cover approximately 35,000 km², spanning Venezuela, Colombia, Ecuador, northern Peru, and small regions in Costa Rica and Panama [2]. Páramos are highly valued for their ecosystem services, particularly as "water towers" of South America [3]. They store vast amounts of water that supply major Andean cities [4,5] and play a role as significant carbon sinks, mitigating global warming by capturing and storing atmospheric carbon [5,6,7]. Despite their regional importance in providing ecosystem services, páramos face growing threats from climate change, land-use expansion, and fire—one of the most frequent drivers of ecological change [5,8].

The fire of Páramo has been an integral and transformative force in human-páramo interactions for thousands of years (at least 7,000), and this must be acknowledged given its ecological importance [9,10,11,12]. As a result, páramos have burned relatively frequently [13], with fire acting as a natural and anthropogenic disturbance agent that influences the dynamics of vegetation [8,14,15]. Although natural phenomena such as lightning strikes and volcanic eruptions can cause fires, they alone cannot explain the historical fire patterns [11,16]. Instead, human activities have been the dominant cause to date, through deliberate land management or accidental ignition [11,16,17,18]. The ecological effects of páramo fires range from beneficial outcomes, such as the regeneration of fire-resistant species [19], to harmful impacts, including soil erosion, hydrological alteration, and the loss of peat and organic matter essential for carbon storage and water regulation [1,16,20]. These effects of fire and its dynamics are amplified by broader anthropogenic factors that go beyond direct burning.

On a global scale, the occurrence of fires in ecosystems has been extensively investigated and is now recognized as the emerging result of complex interactions among climatic factors, fuel availability and continuity, topography, and human activities [21,22]. These drivers operate and interact on spatial and temporal scales: climate modulates vegetation productivity and flammability, thereby defining broad windows of fire susceptibility; vegetation structure and composition determine fuel continuity; topography shapes local microclimates and influences fire spread; and human activities contribute most ignition sources and alter fuel regimes [23,24].

In contrast to other ecosystems worldwide, knowledge of fire drivers in páramo ecosystems remains scarce and fragmented, largely limited to localized case studies. Existing research highlights the influence of human access, traditional burning practices, and seasonal droughts acting on grass-dominated fuels [7,25,26,27]. However, only a few studies have analyzed fire occurrence patterns or quantified the relative importance of these drivers at broader spatial and temporal scales [7,13,28,29]. This limited understanding leaves key processes unexplored, particularly those explaining the spatial variability of fire occurrence, a gap that becomes increasingly concerning under climate change, which is expected to intensify fire regimes through rising temperatures and expanding human pressures [10,30].

To address this gap, our study analyzes the anthropogenic and biophysical drivers of fire in Colombian páramos. We develop a spatio-temporal analysis of burned area dynamics using the páramo Fire Atlas data [13], for the period 2000–2022, and apply a Random Forest (RF) model to identify the main drivers of the spatial pattern. We hypothesize that the spatial pattern of burned area across Colombian páramos is primarily driven by anthropogenic factors such as land use and accessibility. We further expect that biophysical conditions—including vegetation cover, climate, and altitude—modulate these patterns and contribute significantly to local variation in burned area.

2. Materials and Methods

2.1. Study Area

The study area is located in the Colombian Andes (Figure 1), which concentrates nearly half of the world’s páramo ecosystems [6]. In Colombia, páramos are distributed across the three main mountain ranges and isolated massifs, forming a mosaic of landscapes shaped by strong altitudinal gradients [2]. These ecosystems are characterized by cold temperatures, frequent cloud cover, and intense solar radiation [8,16], conditions that create a distinctive microclimatic variability [31]. Such environmental heterogeneity has driven the evolution of highly specialized flora and fauna, resulting in exceptional endemism and the highest recorded plant diversity among high-mountain ecosystems worldwide [3,32]. Along elevation gradients, vegetation transitions from dense woody cover at lower limits to open grasslands, wetlands, and rocky outcrops at higher elevations [33,34]. These transitions are accompanied by differences in soil depth, moisture, and fuel availability that are critical for understanding fire dynamics [20,25,35].

We use the same páramo area proposed in the páramo fire atlas [13] that considers the boundaries of the WWF Ecoregion map [36], a map of land cover specific to the páramos of the northern Andes [37], and the páramo map created by the Humboldt Institute [6]. This boundary also includes a 500-meter buffer to account for the transition between high Andean forest and páramo and excludes all agricultural areas, which lets us focus on fire impact on nontransformed páramo vegetation. In this study, we use the term "páramo unit" to refer to the official complejo de páramo classification adopted by the Colombian National Environmental System (SINA). Each páramo unit represents a defined high-mountain ecosystem area with a specific name, as recognized in national environmental planning and monitoring frameworks. Similarly, each páramo unit has been geographically separated by the cordillera to facilitate some analysis Table A1.

2.2. Data Sets

2.2.1. Burned Area Data for the Páramos Unit

We examined the burned area using the Páramo Fire Atlas, a burned-area database developed for South America’s páramos providing annual burned area data from 1985 to 2022 at 30 m spatial resolution [13]. Accuracy assessments report low commission errors ( 8%) and relatively low omission errors ( 47%), although some underestimation occurs for small fires. Compared to MODIS MCD64A1, The Páramo Fire Atlas detects a higher proportion of small fires, resulting in approximately four times higher burned-area estimates. Our analysis focused on the 2000–2022 time period because pre-2000 estimates have larger uncertainty due to limited imagery availability (See Figure A1). Figure 2 shows an example of the burned area in the Tota-Bijagual-Mamapacha páramo unit in Colombia.

2.2.2. Explanatory Variables

A set of potential explanatory variables was compiled to represent a broad set of biophysical and anthropic factors hypothesized to influence burned area in the páramo ecosystem. Climatic variables included multi-annual averages of max temperature and accumulation of precipitation, as well as human-related variables relevant to high-Andean environments such as the human footprint, distances to forest, buildings and roads. We also included the density of conflict events, given long-standing internal conflict. Because conflict density showed a highly skewed concentration of near-zero values, it was log-transformed as ${log}_{10}ConflictDens+1e^{-6}$.

On the other hand, protected areas represent a central dimension of conservation in Colombia, where more than half of páramo regions fall under some protection category [38]. We use the database of protected areas in Colombian páramos [39], which are classified into four categories according to their legal restrictions. No protected areas, areas without any degree of protection, high-restriction areas, comprised of national and regional natural parks, wildlife sanctuaries, and flora sanctuaries. Medium-restriction areas included regional and national protected forest reserves, while low-restriction areas encompassed Civil Society Nature Reserves, Regional Integrated Management Districts, and Soil Conservation Districts. To enable a fair comparison between burned areas inside and outside protected areas, burned area values were normalized by the total surface area of each páramo unit. This normalization allowed proportional comparisons of fire occurrence, minimizing biases caused by differences in páramo size or the extent of protected coverage (Figure A2). A detailed summary of each explanatory variable is provided Table 1.

2.2.3. Spatial Analysis

We based the spatial analysis on “The Paramo Fire Atlas” database, using it as the dependent variable. We aggregated the dataset into a 500 × 500 meter grid and calculated the percentage of burned area for each cell. This aggregation reduced the stochastic variability of individual fire events and ensured greater consistency with the resolution of the explanatory variables. We resampled all explanatory variables to the same spatial resolution using bilinear interpolation to maintain consistency.

We evaluated multicollinearity among predictors using Pearson’s correlation coefficient, applying a ±0.65 threshold to exclude one variable from each highly correlated pair, which indicated that elevation and conflict density were correlated with maximum temperature (-0.83) and distance to buildings (-0.69), respectively. Given that exploratory analyses revealed nonlinear relationships between predictors and burned area, we opted not to use traditional parametric approaches such as linear or logistic regression, which rely on assumptions of linearity [47,48,49].

Instead, we developed a predictive model using Random Forest, a non-parametric ensemble learning method well suited to capture nonlinear interactions and complex predictor–response relationships [50]. This approach enabled us to evaluate how biophysical and anthropogenic variables contribute to the spatial patterns of burned area. Variable importance was assessed using the permutation importance metric, which measures the increase in prediction error when a variable’s values are randomly permuted. We report this metric as the percentage increase in mean squared error (%IncMSE). To further explore predictor effects, partial dependence plots were generated to describe the functional relationships between individual predictors and burned area [51]. These plots estimate the average marginal effect of a single variable on the response while holding all others constant, thereby isolating its contribution to the model’s prediction. The RF model was parameterized with 500 trees and three variables considered at each split. All analyses were performed in R, using the terra package for spatial analysis [52] and the randomForest [50] and corrplot packages for statistical analysis [53].

3. Results

3.1. Burned Area

Figure 3 shows the total burned area across each Colombian páramo unit between 2000 and 2022, revealing the persistence of fires over the 22-year period. Burned area distribution varied widely among the three cordilleras: the Eastern cordillera recorded the largest proportion (8.9%), the Central reached intermediate levels (3.5%), and the Western was least affected (0.5%), with near-zero values in páramo units such as Citará and El Duende (Table A1 and Table A2). These results underscore the heterogeneous distribution of burned areas in Colombian páramos, with marked contrasts both among and within cordilleras (Figure A3).

Within this general pattern, the Eastern Cordillera stands out with the most burned páramos units, especially those near major cities—such as Tota-Vijagual (34), Cruz Verde–Sumapaz (9), and Altiplano Cundiboyacense (2)—and in binational systems like Perijá (36), where one-third of its area has burned. In the Central Cordillera, the Chiles–Cumbal páramo (5) shows the highest impact, with one-fifth of its surface burned during the study period.

The largest burned areas were concentrated in a small number of páramos units, with ten areas accounting for most of the burned area (see Figure 4). All of these overlap, at least partially, with the National System of Protected Areas (SINAP). Although a substantial share of the total burned area occurred within protected páramos, comparisons with non-protected zones showed no significant differences in the proportion of burned area (Figure A2 and Table A1).

3.2. Drivers of burned area across the páramo ecosystem

The RF model was used to predict burned area (BA), and the out-of-bag (OOB) statistics showed a mean squared residual of 137.7, corresponding to an RMSE of 11.7 and an explained variance of 54.3%. These results indicate that the model captures patterns in burned area and provides estimates of predictive performance.

Figure 5 presents the ranked importance of each predictor variable (%IncMSE). Precipitation emerged as the most influential driver in burned area within the páramo region (150%), followed by anthropogenic variables such as distance to buildings (138.4%), distance to agricultural areas (134.3%), and the human footprint index (126.3%), and NDVI ranked fifth in importance (122.5%). Variables such as distance to forests, temperature, distance to roads, and protected area category exhibited lower importance in explaining burned area (< 97%).

3.2.1. Biophysical Variables

Annual precipitation was the strongest predictor of burned area variability. Partial dependence plots (Figure 6) show that burned area is 75% higher under drier conditions, particularly when precipitation falls below 1500 mm/year (see Table A3) compared to wetter areas. Burned area is also 80% higher in areas with a higher mean temperature, above 16 °C, relative to cooler páramo zones (<12 °C). Vegetation greenness (NDVI) shows a hump-shaped response: areas with intermediate density (0.13–0.20) nearly double the predicted burned area compared to sparsely vegetated (<0.1) or densely forested areas (>0.25).

Although distance to forests and proximity to roads ranked low in the model results, they still influenced burned area (Figure 6). The burned area increased beyond 2 km from each one, which, in the case of forest, is consistent with the NDVI pattern in which more forested areas experienced fewer fires.

3.2.2. Human Pressure

Agricultural areas and distances to buildings increase the probability of burned areas (Figure 6). Burned area peaks within the first two kilometers (agriculture) and nine kilometers (buildings) of these features and declines with distance, although secondary peaks emerge at greater distances (see Table A3). Both variables are strong predictors of burned areas, and their functional responses indicate average burned areas between 4% and 16%.

In the case of buildings, the burned area increases sharply beyond 8 km. The Human Footprint Index—whose values range from 0 to 100 (where 0 = null impact and 100 = maximum impact)—shows a decreasing trend in predicted burned area between values of 0 and 30 (see Table A3), with the highest burned area concentrated at the lower end (0–20). Beyond index values of 20, variability remains minimal, making the pattern somewhat counterintuitive.

For the protected area category, predicted burned area differs little between protected and unprotected regions. However, within protected areas, zones under medium and high restriction show lower predicted values than those under low restriction, with medium-restriction areas showing the lowest. This variable contributes the smallest effect in the model, with functional responses ranging from 0% to 7% burned area.

4. Discussion

This study shows that fire dynamics in páramo ecosystems result from the combined influence of anthropogenic and biophysical drivers. While human activities have been widely recognized as the main source of ignitions [8,26], our findings demonstrate that biophysical variables also play a critical role in shaping spatial fire patterns (Figure 5). Likewise, our analysis of burned areas across Colombian páramos over the past 22 years reveals strong spatial heterogeneity (Figure 3). Burned area varies significantly among regions, with the Eastern Cordillera—recording the highest frequency and extent—emerging as a critical hotspot (Figure A3).

Contrary to our initial hypothesis, where we expected human variables to explain most of the spatial variability in burned area, the RF model revealed that mean annual precipitation is the most influential predictor. Drier areas exhibit higher burned areas, but the relationship is not linear; instead, precipitation shows a critical threshold, with burned area increasing sharply in regions receiving less than 1500 mm of annual rainfall (Figure 6). This pattern aligns with findings from other Andean ecosystems, where rainfall preceding the fire season and during dry periods regulates fire dynamics [25,26,35]. Overall, these results highlight how moisture availability, interacting with fuel loads and human activities, shapes fire regimes in high-altitude landscapes.

The second most important biophysical predictor was NDVI, which explained a substantial share of burned area variability. The model revealed a unimodal relationship between vegetation greenness (NDVI 0.15–0.25) and burned area (Figure 6), indicating that areas with intermediate vegetation density —such as grasslands dominated by clumped grasses in a continuous layer [54] —are the most susceptible to fire. These grasslands dominate much of the páramo landscape, providing the primary fuel that drives fire behavior [19]. This pattern reflects fuel availability thresholds: low NDVI areas lack sufficient continuous biomass to sustain fire, whereas high NDVI areas usually retain more moisture, which limits ignition [55]. These areas are generally occupied by Andean forests or shrubs [2]. Although the most important variable is biophysical, human pressures have a strong impact, shaping ignition patterns and reinforcing fire pressure across páramos units.

Anthropogenic factors exerted a strong influence on burned area [8]. Proximity to buildings and croplands was closely linked to higher burned area (Figure 6). This pattern reflects greater accessibility for ignition sources, both intentional and accidental. Intentional uses include agricultural expansion, burning to promote regrowth for livestock, deforestation, and the clearing of woody vegetation [19,26]. Accidental ignitions and vandalism have also been reported[26,56]. Historically, páramos have been shaped by the combined influence of fire and grazing as part of their land-use history [8,19,26], as well as domestic fuel use and charcoal production [16]. Cultural practices and local beliefs promote the use of fire, such as burning páramo vegetation during dry periods under the assumption that smoke induces cloud formation and rainfall [16,27]. More recently, demographic pressures, including population growth and the proximity of urban centers, have been linked to heightened fire pressure [20,26].

Despite fire suppression policies—including the prohibition of burning in páramo ecosystems under Colombian law (Law 2111 of 2021; Decree 1076 of 2015; Resolution 358 of 2024), the effectiveness of formal protection measures appeared limited. Our analysis indicates that protection status, as a categorical variable, did not significantly reduce burned area (Figure 5). Although overall burned area did not differ markedly between protected and unprotected areas, zones under medium and high restriction exhibited slightly lower proportions of burned area (Figure 6). Although the National System of Protected Areas (SINAP) provides guidelines for fire prevention and control, its effectiveness is limited by weak institutional coordination, scarce financial resources, and insufficient enforcement capacity. According to the Colombia Country Dossier (2021)[57], only about 14% of the national environmental budget is allocated to SINAP, and its overall management effectiveness index is approximately 65%. These figures reveal persistent gaps between policy design and on-the-ground implementation. Combined with the persistence of local land-use practices that depend on fire, these institutional shortcomings help explain why fire events continue to occur even within protected páramo areas.

Ultimately, our results confirm that burned areas in páramo units do not occur randomly in space or time but rather emerge from predictable dynamics within coupled human–natural systems. Looking ahead, projected reductions in precipitation and rising temperatures are likely to increase flammability, while expanding agriculture, grazing, and urban growth may intensify human pressures [21]. These trends suggest that fire regimes in páramos will be shaped by the combined effects of climate change and evolving socio-economic dynamics. This highlights the urgent need for targeted, adaptive management strategies.

It is crucial to recognize these predictable dynamics, as they enable the design of targeted interventions that address the most influential drivers of fire risk. These insights highlight the urgent need for integrated and spatially explicit fire management strategies in páramo ecosystems. Key actions include implementing early warning systems during dry seasons and promoting sustainable, non-pyric alternatives for pasture renewal where feasible. Monitoring vegetation transitions is equally important to avoid the buildup of continuous, fire-prone fuel matrices. Integrated fire management (IFM) should also incorporate traditional ecological knowledge into evidence-based policies, ensuring strategies are both scientifically sound and socially legitimate. Beyond prevention and control, it is urgent to understand how fire affects biodiversity and ecosystem services, given the limited knowledge available. Addressing this gap is critical to strengthen fire prevention, enhance ecosystem resilience, and align conservation efforts with local practices.

Appendix A.1

Figure A1. The frequency distribution of burned area percentages across all páramo units reveals a strong right-skewed pattern. Most observations correspond to less than 5% of burned area, while only a small proportion exceed 25%. This indicates that, although fire events are widespread, they tend to affect limited portions of the landscape. The long tail of the distribution suggests that extensive burns are rare but occur occasionally.

Figure A1. The frequency distribution of burned area percentages across all páramo units reveals a strong right-skewed pattern. Most observations correspond to less than 5% of burned area, while only a small proportion exceed 25%. This indicates that, although fire events are widespread, they tend to affect limited portions of the landscape. The long tail of the distribution suggests that extensive burns are rare but occur occasionally.

Table A1. Páramo units grouped by Cordillera, with burned and protected proportions.

Páramo	% burned	% protected area	Cordillera
Belmira	2.23	100.00	Western
Cerro Plateado	0.55	18.22
Frontino - Urrao	0.35	84.49
Paramillo	0.04	99.88
Farallones de Cali	0.02	100.00
Tatamá	0.01	100.00
Citará	0.00	56.11
El Duende	0.00	39.77
Chiles - Cumbal	21.34	16.99	Central
Los Nevados	5.49	56.03
Nevado del Huila - Moras	4.06	78.22
Guanacas - Puracé - Coconucos	2.81	24.96
Las Hermosas	2.62	65.08
Chilí - Barragán	1.80	35.73
La Cocha - Patascoy	0.82	34.98
Sotará	0.76	42.05
Sonsón	0.20	43.39
Doña Juana - Chimayoy	0.13	50.85
Perijá	33.11	79.67	Eastern
Tota - Bijagual - Mamapacha	16.59	54.60
Altiplano Cundiboyacense	16.13	11.89
Cruz Verde - Sumapaz	14.44	53.23
Sierra Nevada de Santa Marta	11.00	100.00
Pisba	10.56	29.86
Guantiva - La Rusia	7.34	50.82
Iguaque - Merchán	6.21	54.27
Almorzadero	5.86	0.64
Sierra Nevada del Cocuy	5.33	71.27
Rabanal y río Bogotá	5.22	68.56
Guerrero	3.78	85.34
Jurisdicciones - Santurbán - Berlín	3.31	43.00
Chingaza	1.34	81.53
Tamá	0.12	75.03
Los Picachos	0.00	88.20
Miraflores	0.00	84.57
Yariguíes	0.00	99.36

Table A1. Páramo units grouped by Cordillera, with burned and protected proportions.

Páramo	% burned	% protected area	Cordillera
Belmira	2.23	100.00	Western
Cerro Plateado	0.55	18.22
Frontino - Urrao	0.35	84.49
Paramillo	0.04	99.88
Farallones de Cali	0.02	100.00
Tatamá	0.01	100.00
Citará	0.00	56.11
El Duende	0.00	39.77
Chiles - Cumbal	21.34	16.99	Central
Los Nevados	5.49	56.03
Nevado del Huila - Moras	4.06	78.22
Guanacas - Puracé - Coconucos	2.81	24.96
Las Hermosas	2.62	65.08
Chilí - Barragán	1.80	35.73
La Cocha - Patascoy	0.82	34.98
Sotará	0.76	42.05
Sonsón	0.20	43.39
Doña Juana - Chimayoy	0.13	50.85
Perijá	33.11	79.67	Eastern
Tota - Bijagual - Mamapacha	16.59	54.60
Altiplano Cundiboyacense	16.13	11.89
Cruz Verde - Sumapaz	14.44	53.23
Sierra Nevada de Santa Marta	11.00	100.00
Pisba	10.56	29.86
Guantiva - La Rusia	7.34	50.82
Iguaque - Merchán	6.21	54.27
Almorzadero	5.86	0.64
Sierra Nevada del Cocuy	5.33	71.27
Rabanal y río Bogotá	5.22	68.56
Guerrero	3.78	85.34
Jurisdicciones - Santurbán - Berlín	3.31	43.00
Chingaza	1.34	81.53
Tamá	0.12	75.03
Los Picachos	0.00	88.20
Miraflores	0.00	84.57
Yariguíes	0.00	99.36

Table A2. Distribution of burned area by cordillera.

Cordillera	Total area (ha)	Burned cells	% burned
Central	1,013,731	35,121	3.46
Eastern	1,490,635	133,452	8.95
Western	77,968	362	0.46

Table A2. Distribution of burned area by cordillera.

Cordillera	Total area (ha)	Burned cells	% burned
Central	1,013,731	35,121	3.46
Eastern	1,490,635	133,452	8.95
Western	77,968	362	0.46

Table A3. Thresholds were derived from Random Forest partial dependence plots using the Bai–Perron method (strucchange::breakpoints), with the number of breaks selected by BIC. Values are in the original units of each predictor (excluding categorical variables such as ProtectedAreas).

Variable	Thresholds
DistForest	1.72
Temperature	17.08
Precipitation	1577.41
DistBuildings	9.60
DistRoads	5.60
DistAgriculture	2.14
NDVI	0.14
ProtectAreas	–
Human Footprint	20.50

Table A3. Thresholds were derived from Random Forest partial dependence plots using the Bai–Perron method (strucchange::breakpoints), with the number of breaks selected by BIC. Values are in the original units of each predictor (excluding categorical variables such as ProtectedAreas).

Variable	Thresholds
DistForest	1.72
Temperature	17.08
Precipitation	1577.41
DistBuildings	9.60
DistRoads	5.60
DistAgriculture	2.14
NDVI	0.14
ProtectAreas	–
Human Footprint	20.50

Figure A2. Normalized burned area inside and outside protected areas (PAs; % relative to the available area in each category). Boxplots display the distribution of páramo-level values, paired within complexes to account for variation in total size. Statistical comparison using the Wilcoxon signed-rank test revealed no significant differences between inside and outside fractions (two-sided p = 0.741; one-sided test for inside < outside p = 0.371)

Figure A2. Normalized burned area inside and outside protected areas (PAs; % relative to the available area in each category). Boxplots display the distribution of páramo-level values, paired within complexes to account for variation in total size. Statistical comparison using the Wilcoxon signed-rank test revealed no significant differences between inside and outside fractions (two-sided p = 0.741; one-sided test for inside < outside p = 0.371)

Comparison between Cordilleras

Kruskal–Wallis tests based on normalized values (% burned area relative to total páramo area) indicated significant differences among cordilleras (H = 9.18, p = 0.0101). Post-hoc pairwise Mann–Whitney tests revealed that the Western cordillera differed significantly from both the Central (p = 0.0044) and Eastern (p = 0.0123), whereas no significant difference was detected between Central and Eastern cordilleras (p = 0.2212).

Figure A3. Boxplots of burned area normalized by páramo size (percent of each páramo’s area burned) across the Western, Central, and Eastern Cordilleras. A Kruskal–Wallis test detected differences among cordilleras (H = 9.18, $p=0.010$). Post-hoc pairwise Mann–Whitney U tests showed Western < Central ($p=0.0044$) and Western < Eastern ($p=0.012$), whereas Central vs Eastern was not significant ($p=0.22$). Points represent individual páramos units.

Figure A3. Boxplots of burned area normalized by páramo size (percent of each páramo’s area burned) across the Western, Central, and Eastern Cordilleras. A Kruskal–Wallis test detected differences among cordilleras (H = 9.18, $p=0.010$). Post-hoc pairwise Mann–Whitney U tests showed Western < Central ($p=0.0044$) and Western < Eastern ($p=0.012$), whereas Central vs Eastern was not significant ($p=0.22$). Points represent individual páramos units.