Carbon balance shift in mountain peatlands along a gradient of grazing disturbance in the tropical Andes (Colombia)

High-elevation cushion peatlands are typical ecosystems of the Andes above 4000 m of altitude, with an important role in hydrology and global carbon sequestration. In Tropical Andean context, grazing livestock is one of the main threats to cushion peatlands, altering the vegetation and the storage carbon function. The aim of this research is to understand how grazing influences cushion peatland functioning by identifying ecological thresholds for carbon balance process. The study was carried out in four Andean peatlands in the northern part of Colombian Andes during 2019–2020. We established 30 plots of 1 m2 where water table level, vegetation cover, and grazing disturbance were monitored. We also measured CO2 fluxes using an infrared gas analyzer connected to a closed static chamber, which registered net ecosystem exchange and respiration data. Considerable variation in the conservation status of Distichia muscoides cushions was found within the sampled peatlands, reflecting an heterogeneous signal of grazing disturbance that is evident at the plot-specific scale. Decreasing water table level was related with changes in dominant vegetation, from compact cushion species to grasses proliferation, exacerbating disturbance effects and carbon emissions. Mixed-effects logistic regression models showed a carbon balance shift, from CO2 sink to net emitter, in plots with high disturbance intensity and low D. muscoides cover. This study provides information for a better understanding of mountain peatlands functioning in the Tropical Andes and underlines the key role of D. muscoides cushions and the water table in carbon balance shift.


Introduction
Peatlands are common ecosystem in tropical mountains, and they have one of the largest carbon density of all the ecosystems, with nearly 1000 ton C ha −1 (Hribljan et al. 2016).Peatland carbon storage requires the accumulation of partially decomposed organic matter over thousands of years, which results from the imbalance between primary production and limited soil decomposition.Human disturbances such as grazing livestock affect the stability of the water table and plant composition, which are key elements that have a direct effect on the ability of the system to accumulate peat (Belyea and Malmer 2004;Moore and Dalva 2006).Our understanding on how grazing disturbance alters carbon accumulation process in peatland ecosystems is still scarce in the tropical Andes.Better understanding of carbon balance shifts once it is pushed past a critical threshold is needed.Management strategies for Andean mountain peatlands rely on the development of indicators that help address the impacts of grazing disturbance (Chimner and Welker 2011).
Carbon storage in soils requires an organic matter input provided by litter or old roots from peat-forming plants, such as Distichia muscoides which is a peat-forming cushion plant with high net primary production rates but it's also dependent of flooding conditions (Cooper et al. 2015).Thus, changes in vegetation composition and water table depth result in the alteration of net carbon assimilation and decomposition rates, (Clymo 1984;Turetsky et al. 2002;Belyea and Malmer 2004).The alteration of the carbon balance happens when stress conditions cause greater respiration loss than the CO 2 uptake, due to peat oxidation and acceleration of decomposition rate.The change in the functioning of the system from a carbon sink to a net emitter of CO 2 corresponds to a cross of an ecological threshold (Turetsky et al. 2002;Moore and Dalva 2006;Larsen and Alp 2015;Toms and Villard 2015).Although studies aimed at understanding ecological thresholds from empirical data are limited and the scope in decision-making is questioned (Dodds et al. 2010;Hillebrand et al. 2020), this approach could be useful to represent critical conditions of carbon balance shift as key indicator oriented toward management strategies to minimize adverse changes in peat accumulation process (Belyea and Malmer 2004;Johnson 2013;Toms and Villard 2015).
Even though mountain peatlands cover small areas, they have large carbon stocks in the soil due to a combination of high mineral content from the geomorphological process of the mountains and the type of vegetation that forms the peat (Hribljan et al. 2017).The steep topography of the Andes, together with a diverse and complex geomorphological setting of high-altitude lakes and glacial valleys, provides adequate conditions for the formation of wetlands and peatlands with unique vegetation (Cooper et al. 2015).At high elevations of Andes mountains (> 4000 masl), peatland vegetation grows in compact cushion of vascular plants such as D. muscoides (Juncaceae) and Plantago rigida (Sklenář and Jørgensen 1999).Since cushion species at high elevations are the main peat-forming plant and play a key role in carbon sequestration, changes in cushions dominance could be catastrophic for the functioning of mountain peatlands.
Andean mountains and in particular their alpine ecosystems are important sources of feed for native and nonnative grazers, including Andean camelids at least for the last 9000 years (Browman 1984), especially in dry Puna region of the central Andes, where peatlands provide the evergreen forage for the livestock in the dry season (Maldonado Fonkén 2014).Grazing plays an important role for biomass production because of the historical relationship between native camelids and cushion peatlands (Garcia et al. 2014).However, in the northern Andes, cushion peatlands are instead affected by introduced grazers such as sheep, horses, and cattle, which also generate a disturbance on plant cover due to their heavier bodies and hard hooves (Maldonado Fonkén 2014; Machaca et al. 2018).Trampling by grazers pug peat soil, removing part of the biomass and the carbon stored in peat due to aerobic decomposition (Salvador et al. 2015).
Here, we hypothesized that human disturbances by high intensity grazing, combined with a lowering water table, could negatively affect the carbon balance of cushion peatlands dominated by D. muscoides, shifting the system from a carbon sink to a carbon source.Therefore, the purpose of this study was to assess how grazing and vegetation influence CO 2 emissions in Tropical Mountain peatlands and the limit at which the system shifts from a carbon sink to a net carbon source in order to provide management tools that protect carbon storage capacity in peatland ecosystems.Our research focused on (1) comparing environmental drivers and vegetation patterns along a grazing disturbance gradient and (2) identifying response variables that affect carbon balance, using the ecological thresholds approach.

Study site
The study was conducted at the high elevations of El Cocuy National Park, in the northern Andes of Colombia (Fig. 1).A large portion of the protected area overlaps with an Indigenous territory, as well as private lands managed by farming families who frequently graze their livestock within the protected area.The area is characterized by an annual mean temperature of 6.4 °C and a bimodal precipitation regime, with two wet periods in May and October with an average rainfall of 868.4 mm (López-Moreno et al. 2022).We randomly selected four peatlands dominated by cushions of D. muscoides, combined with graminoids such as Agrostis breviculmis and Calamagrostis effusa and the moss Breutelia chrysea (Fig. 2).All sites had more than five meters of peat, were located at a similar elevation, and had a long history of livestock use, with evidence of sheep grazing disturbance (Table 1).The peatlands were visited on five occasions during the 2019-2020 period to include the annual climatic variation.

Vegetation sampling and environmental factors
We located a total of 30 1 m 2 plots randomly positioned and stratified across the main plant and ground cover gradients observed on the selected peatlands.Plant composition was sampled in all the plots by visual estimation using a 1 m 2 quadrat with 10 × 10 cm sub-quadrants (Mueller-Dombois and Ellenberg 1974), where plant and ground covers were estimated to the nearest 1% for each species.From plant composition information, ground cover was classified in cushions of D. muscoides, hollows, dead cushions, or exposed peat-, grass-, and moss-dominated sites.Plant species were identified following standard herbarium identification protocols, and species names were assigned following the international plant name index (IPNI 2023.https:// www.ipni.org/ accessed on February 28, 2023).
Due to the differential sheep grazing intensity within individual peatlands, disturbance was not homogenous throughout each site.Thus, disturbance was characterized at plot level and nearby area using the percentage of exposed peat cover and the evidence of sheep alteration such as excrements, hoof marks, trails, and/or herbivory.The type of specific disturbance was classified according 1 3 to a semi-quantitative intensity scale from 1 to 6, where the presence of excrement represented the lowest intensity while pugged and peat implied the most severe (Online Resource: Fig. S1).The disturbance category was developed from combining exposed peat cover and the sheep alteration scale.For water table level, we installed 30 monitoring wells that were located near the 1m 2 plots.Water level was measured manually during each visit to the study site.When the water table was above the peatland surface, a positive value was recorded and, if it was below the surface, a negative value was recorded.The water table level was monitored for each plot at each visit to the study site.Also, in each visit, a water sample was collected to measure pH and electrical conductivity (EC).Average monthly air temperature and precipitation were provided by the Instituto de Hidrología, Meteorología y Estudios Ambientales (IDEAM) from a meteorological station closest to the study sites.

CO 2 measurements
CO 2 flux was measured using an infrared gas analyzer (EGM4; PP Systems, Amesbury, USA) connected to a closed transparent static chamber of 60 cm of diameter and a height of 50 cm (Hutchinson and Mosier 1981).For each measurement, soil and air temperature, relative humidity, atmospheric pressure, and photosynthetically active radiation using a PAR sensor attached to the inside of the chamber were measured (Apogee MQ-200 Quantum meter).Using a total of 105 flux measurements with a duration of 124 s, each was performed in two consecutive phases: (1) a phase where the chamber receives natural solar radiation, registering net ecosystem exchange (NEE), and a second phase where the chamber does not receive radiation, registering autotrophic and heterotrophic respiration (RE).Gross primary production (GPP) was calculated as the difference between net ecosystem exchange and ecosystem respiration for each pair of light and dark, flux measurements.

Statistical analysis
Parametric and non-parametric analyses of variance (one way-ANOVA and Kruskal-Wallis) were performed to identify the differences between the environmental and flux variables among the sampled peatlands.Plant composition was analyzed using a cluster analysis with UPGMA method as well as analysis of similarity between sites and vegetation plots based on the Bray-Curtis dissimilarity matrix among samples.Additional MDS analysis was used to analyze the homogeneity differences in plant composition and the main environmental gradients associated with the changes in species composition.The main environmental patterns were analyzed using principal components analysis (PCA).A repeated-measures analysis with mixed effects was applied using lme4 package to quantify the relationships between explanatory variables (vegetation cover and water table) and CO 2 fluxes, where plot and period were incorporated in the models as random effects, and maximum likelihood estimation was performed.CO 2 fluxes was the dependent variables with GPP, NEE, and ER analyzed separately.Five NEE outliers were omitted due to equipment problems.
The relationships between vegetation and water table and NEE were analyzed using mixed-effects binomial logistic regression models to identify the ecological thresholds (McCullagh and Nelder 1989;Toms and Villard 2015).NEE was coded as a binary variable, converting negative values to values of 0 (net ecosystem uptake) and positive values to values of 1 (net ecosystem emission).The threshold was determined when the probability of the system going from sink to emitter crossed 0.5 probability (Toms and Villard 2015).The selection of the logistic models that better explained the changes in NEE was made based on the significance of the model and AIC value.Residuals were checked for normality and that homogeneity of variance assumptions was met, and significance was determined at an alpha value of 0.05.Also, average and standard error of NEE values were calculated to evaluate the carbon emissions behavior above and below the determined threshold for all the sampled plots.All

Vegetation patterns and environmental factors
A total of 43 species were registered, including 13 bryophytes and 30 vascular plants.Plot vegetation was dominated by D. muscoides with a mean 54% cover in all study sites (F4,130 = 35.2,p = 0.001; Table 2).Exposed peat was more common in Curial alto and Corral Grande (15 and 16% cover, respectively) and commonly associated with sedges (Carex sp.), grasses (Agrostis breviculmis), rosettes such as Werneria pygmaea, herbs that were common on disturbed non-peatland areas such as Lachemilla nivalis and Oreomyrrhis andicola, and mosses such as Bryum argenteum and Breutelia spp.Five cover types were classified in the plots: 1) D. muscoides, 2) mosses, 3) grasses, 4) exposed peat, and 5) other vascular plants (Table 2).
Water table monitoring showed that site and time of the year had a significant effect on water table depth, being Curial site the wettest, followed by Corral Grande, Chapetona, and Curial Alto the driest.Water table levels ranged from 10 to 33 cm below the soil surface and were very variable throughout the year, with a mean depth of − 19.6 ± 12.2 cm during the wet season, and -30.6 ± 10.3 cm during the dry season (January and February), when some piezometers were completely dry (Table 2).Plots with less exposed peat showed more stable water table levels throughout the year (F4,130 = 35.2,p = 0.001).
Water table showed a significant positive relationship with D. muscoides cover from the mixed model (Pseudo-R2 = 0.5, p = 0.005) and a negative relationship with grasses (Pseudo-R2 = 0.2, p = 0.06).Plots with water table above 20 cm depth had up to 50% D. muscoides cover, whereas plots with lower water table (< 20 cm from the surface) exhibited an average increase of 15% cover of grasses (Fig. 3).A negative relationship between water table depth and the percentage of exposed peat was also observed (Pseudo-R2 = 0.38, p = 0.03), except in Chapetona site, where almost no exposed peat was recorded.
The ordination diagram shows that there was a plot-scale heterogeneity regardless of the peatland.Plant composition changed across the water table and disturbance gradients.Plots with higher exposed peat cover had higher dominance of species associated with drier conditions such as Agrostis breviculmis, Lachemilla nivalis, Bryum argenteum, and Oreomyrrhis andicola.In contrast, plots with higher water table level and low disturbance intensity were dominated by D. muscoides and Plantago rigida and had less exposed peat cover (Fig. 4).Exposed peat cover was strongly correlated with disturbance categories (r = 0.89) and was the environmental factor that better explained variance between plots.

Carbon fluxes
The studied peatlands had an average carbon assimilation rate of −0.04 ± 2.1 g CO 2 m −2 h −1 .Carbon fluxes were positively related to soil temperature, solar radiation, and water table level.GPP ranged from − 0.7 to − 0.22 g CO 2 m −2 h −1 (Table 1).GPP was influenced by PAR, with a significant response function during the hours of highest radiation (around 1900 µmol m-2 s -1 ; Fig. 5).Carbon assimilation values decrease at low water table levels and increase with increasing water table level (Pseudo-R 2 = 0.29, p = 0.04; Fig. 5).In contrast, ER responded oppositely to increasing water tables (Pseudo-R 2 = 0.45, p = 0.5).Also, ER had a strong relationship with soil temperature, since ER values increase with soil warming (Pseudo-R 2 = 0.35, p = 0.12; Fig. 5).D. muscoides exposed peat and grass cover had a significant relationship with NEE.Carbon assimilation decrease in exposed peat-dominated plots, where positive or low NEE values were recorded (Pseudo-R 2 = 0.13, p = 0.01).In contrast, plots with high D. muscoides cover had negative NEE values, which means high-carbon assimilation rate (Pseudo-R 2 = 0.15, p = 0.0001).Also, high cover of grass was related with positive NEE values.The logistic regression models showed a variation on carbon balance values related to the loss of vegetation cover of the cushions at the plot scale.In both cases when D.muscoides cover decreased and when the exposed peat area increased, a threshold point was found for which there was a greater probability that respiration exceeded assimilation.Plots with less than 33% of D. muscoides cover had a high probability of a net carbon emission, reaching a 0.78 probability of being net carbon emitters when D. muscoides cover was completely lost (Fig. 6).In the same way, plots with more than 25% cover of Fig. 5 Linear relationship between net ecosystem exchange (NEE; top) and water table level (cm); gross primary production (GPP; center) versus photosynthetically active radiation (PAR; µmol m −2 s −1 ); and ecosystem respiration (ER; bottom) and soil temperature (°C) for all sampled peatlands of the Cocuy National Park exposed peat had also a high probability of net carbon emission (Fig. 7).Also, average and standard error of NEE values were related with changes in D. muscoides and exposed peat cover according to the NEE patterns above and below the determined thresholds.

Discussion
The variation in species composition along the disturbance gradient suggests that the signal for grazing disturbance is heterogeneous within the peatlands, and it is evident at plot scale rather than on the whole peatland as a unit.The variation in vegetation cover along the disturbance gradient and water table fluctuation could lead different trajectories of plants composition.The relationship between low water table and high grass cover showed how graminoids proliferate at dryer areas, suppressing cushions species and favoring the proliferation of species more adapted to grazing use and drier condition such as Agrostis breviculmis, Lachemilla sps, and Oreorrhymmys andicola (Benavides et al. 2013ab;Loza Herrera et al. 2015;Greenwood et al. 2019).Water table drawdown could be also favored by trampling and grazing, by reducing cushions surface humidity (Machaca et al. 2018;Zhong et al. 2020), which could also accelerate graminoids rise and peat decomposition.Also, oral communications with local people and field observations showed that sheep did not forage Distichia sp.itself, but rather the graminoids growing over the cushions, evidencing that sheep grazing pressure only occurred in drier areas with palatable plant cover.This interaction could lead to positive feedback between graminoids proliferation favored by water table drawdown and increasing disturbance, opening new questions such as whether the presence of grasses is explained, to a greater extent, by disturbance or by the low water table.
Our results for net ecosystem uptake were in a range of − 0.001 and − 0.67 g CO 2 m −2 h −1 , much lower than those reported for cushion-dominated grazed peatlands in Ecuador (− 1.25 g CO 2 m −2 h −1 ) (Sánchez Garces et al. 2017).Besides both studies being carried out at similar altitudes (> 4200 m) and hydrological conditions, the Ecuadorian peatland barely had exposed peat cover as well as noncushion plant abundance, which could explain the carbon balance differences (Urbina and Benavides 2015;Sánchez Garces et al. 2017).Water table also showed a spatial relationship with carbon balance values, with high ER values and low positive NEE values (lower C uptake) at low water table level (Zhong et al. 2020).Water table drawdown increases the oxygen diffusion within the unsaturated peat and promotes peat oxidation by an increase in soil microorganism ER activity, producing carbon loss (Moore and Dalva 2006;Jaatinen et al. 2008;Hribljan et al. 2014).
This research allowed us to understand the influence of water table fluctuations on the instantaneous CO 2 fluxes, highlighting the importance of short-term variation understanding to infer the effects of water table drawdown if it is more prolonged in time.Besides the explored methods were limited to explain the variability of peat accumulation process at large temporal scales, drier conditions could increase decomposition rate of peat and grasses proliferation, leading to a rapid loss of stored carbon and driving to an ecosystem shift from carbon sink to a net carbon emitter.Although cloud cover was not considered in the analysis, it has a significant effect on carbon balance shift, having a positive impact on CO 2 assimilation rate (Planas et al. 2020).However, the one-year monitoring of CO 2 fluxes expects to lead us to include annual climatic variation and rule out the effect of cloudiness over the NEE results.
With one-year observations, the logistic regression models showed a carbon unbalance related to the disturbance signal at a plot scale, from which respiration exceeded assimilation.Carbon balance shift from negative to positive NEE was observed both for the decrease in D. muscoides cover, and the increase of exposed peat cover, pointing out an ecological threshold from C sink to C source (Pelletier et al. 2015).Areas with less than 33% of D. muscoides cover and more than 25% of exposed peat cover are likely to have crossed the ecological threshold and they are losing the stored carbon, compromising the persistence of the ecosystem (Sánchez Garces et al. 2017).
For ecological thresholds models, single drivers were analyzed following the best statistical models, confounding variation for the identification of threshold responses (Simmonds et al. 2017).However, other abiotic variables could interfere, highlighting the importance of considering complementary drivers that also affect or interact with NEE such as water table level, microtopography, water chemistry, and climate variability (Belyea and Malmer 2004).Furthermore, grazing pressure combined with environmental changes could bring forward the carbon balance shift at lower disturbance intensity (Machaca et al. 2018), which could lead us to devise management strategies that delay threshold points such as develop cushion restoration and restrict livestock access at drier season.
Little is known about how ecological processes response to threshold cross or the capacity of vegetation and humidity to recover after certain disturbance intensity (Benavides et al. 2013a,b).Evidence of disturbance effect will be dependent on the spatial scale measure, where the response variable and the environmental drivers selected will determine the identification of a threshold for carbon balance shift (Belyea and Malmer 2004).Despite better understanding of the resilience of the systems is required, small spatial and temporal scale processes such as NEE could be an earlywarning signal for ecosystem decline (Eppinga et al. 2009).

Conclusion
This study provides information for a better understanding of the functioning of mountain peatlands of the Tropical Andes, crucial in carbon storage and water regulation (Salvador et al. 2015), and underlines the key role of D. muscoides cushions and the water table in the peat accumulation process.Considerable variation in the conservation status of D. muscoides cushions was found within the sampled peatlands, showing an heterogenous signal of grazing disturbance and an effect on the carbon balance at microscale level.Since carbon balance determines a key function of peatlands, high-resolution variables like NEE could be a useful response variable to anticipate responses in largerscale processes (Eppinga et al. 2009).This hypothesis needs further exploration through additional data for a better characterization of spatial patterns on carbon balance behavior.
Rather than providing specific number of vegetation cover at which the system crosses non-return state, these results could guide decision-makers about consistent monitoring of specific response variables, and appropriate spatial and temporal scale at which thresholds may arise (Huggett 2005).

Fig. 1
Fig. 1 Location of study sites in El Cocuy National Park, in the northern Andes of Colombia

Fig. 3
Fig.3 Linear relationship between water table level (cm) and vegetation cover type (%) for all sampled peatlands of the Cocuy National Park.(Top): D. muscoides cover; (center): grass cover and (bottom): exposed peat cover

Fig. 6
Fig. 6 Probability of net carbon emission according to D. muscoides cover.Red line represents the ecological threshold where there is a 50% probability of a positive carbon balance for all sampled peatlands of the Cocuy National Park

Fig. 7
Fig. 7 Probability of net carbon emission according to exposed peat cover.Red line represents the ecological threshold where there is a 50% probability of a positive carbon balance for all sampled peatlands of the Cocuy National Park

Table 1
Descriptions of the studied cushion-dominated peatlands sampled at El Cocuy National Park Geographical coordinates, altitude and number monitored plots