Exploring Linkages between Supporting, Regulating, and Provisioning Ecosystem Services in Rangelands in a Tropical Agro-Forest Frontier

: Rangeland management in former tropical rainforest areas may a ﬀ ect ecosystem services. We hypothesized that management practices like burning and overgrazing reduce supporting (soil quality) and consequently also provisioning (forage productivity and quality) and regulating (nutrient cycling) ecosystem services. We established 31 exclosures in two landscape categories (alluvial soils, low-hills), documented management practices, and assessed 18 soil quality indicators, litter decomposition as a proxy for nutrient cycling, and forage quantity and quality during one year in grasslands of the Lacandon region, southeast Mexico. Path analysis was used to explore direct and indirect e ﬀ ects of livestock management practices on soil-based ecosystem services. Landscape position had direct e ﬀ ects on management practices, and direct and indirect e ﬀ ects on soil properties. Altitude (a proxy for the soil catena, ranging from alluvial soils along the Lacant ú n river to Cambisols and Acrisols in the low-hills) was the variable showing most signiﬁcant negative relations with soil quality and forage production. Decomposition rate was site-speciﬁc and had no relation with landscape position and management. Our study suggests that position on the landscape, which relates to nutrient and water availability, had stronger e ﬀ ects than management practices on forage productivity and quality and drives farmers management practices. authors have read and agreed the version of the manuscript.


Introduction
Agricultural expansion, both for cropland and rangeland, is one of the major drivers of deforestation in the tropics [1,2]. Extensive cattle production occupies more than 27% of rural landscapes across Latin America, and continues to expand [3]. The establishment and maintenance of rangelands has large consequences for local and regional biodiversity. Given that biodiversity plays a critical role in ecosystem functioning, large impacts on ecosystem services can be expected to occur after transformation from forest to grasslands. Ecosystem services have been categorized in four classes, viz., regulating, supporting, cultural, and provisioning services [4][5][6]. Supporting and regulating services are those services crucial for the sustained supply of provisioning services, the products Land 2020, 9,511 2 of 17 that humans obtain from ecosystems (food, feed, building materials, medicine, etc.). Supporting and regulating services are often replaced by external human inputs (mineral fertilizer to replace natural nutrient cycling, pesticides to replace natural pest control) in intensively managed agricultural landscapes, resulting in a disconnection between supporting and regulating services on the one hand and provisioning services on the other. However, for low-income farmers with limited resources, these supporting and regulating services may not be sufficiently replaced and become crucial for provisioning services.
Establishment of extensively managed rangelands in former tropical rainforest areas may cause a decline of supporting and regulating services, and therefore in provisioning services, which subsequently results in a further decline of supporting and regulating services. Deforestation, burning, and rangeland management can lead to a decline of soil quality through erosion, loss of soil organic matter and nutrients, and compaction, thereby reducing the ability of the soil to produce enough forage of sufficient quality, which then feeds back to the levels of soil organic matter and nutrients. Additionally, deforestation in riparian areas reduces water availability for livestock, wildlife, and local human populations [7][8][9]. Here, we focus on a better understanding of the effects of different rangeland management practices, that follow after deforestation, on soil quality and ecosystem services.
In Mexico, livestock ranching is currently the most common agricultural activity, occupying 56% of national territory [10,11]. In Marqués de Comillas, located in the Lacandon region, Chiapas, southeast Mexico, rangelands are the major anthropogenic disturbance, representing around 50% of the area [12]. Currently, many of the pasturelands in Marqués de Comillas are not managed efficiently. Local farmers have faced important deficiencies in technical assistance, training, and access to technology. Often livestock, as a private tradeable commodity, functions as a substitute for money savings, allowing for sales under immediate need for cash [13]. Apparently, current management practices fail in their intent to contribute to rural employment and to provide good levels of meat and milk production. Although the economic value of bovine production has increased since 2006, local communities have not improved their living conditions and income. Some studies reported that local farmers of Marqués de Comillas perceive that their lands are 'poor', and face degradation due to livestock ranching [14,15]. Additionally, the establishment of mono-dominant patches of bracken fern (Pteridum aquilinum) was identified as an important indicator of land degradation and a reason for field abandonment. Landowners perceived that invasion by bracken fern is favored by the recurrent use of fire to control the establishment of secondary vegetation, and by declining soil fertility [16]. In an attempt to prevent further decline in ecosystem services, and make cattle ranching more 'profitable' or 'sustainable', farmers rotate cattle over large areas, leaving parts of the rangelands in fallow periods to recover their productivity [15]. However, these practices drive further rangeland expansion at the expense of native forests at larger spatial scales. Sustainable use of rangelands could prevent further deforestation in agro-forest frontier areas, and thereby also contribute to maintaining biodiversity in the remaining forest fragments. Conservation of biodiversity and ecosystem services will ultimately depend on our ability to sustainably manage human-modified ecosystems.
There is an urgent need to understand how these rangelands respond to management practices, how they contribute to human well-being, and how people could use them in a sustainable way. A variety of effects of conversion of forest to rangeland has been documented. A study in the Atlantic Zone of Costa Rica found significant changes in soil properties associated with conversion of forest to pasture; namely an increase in pH, bulk density and concentrations of NH 4+ , and a decrease in concentrations of NO 3− and rates of N-mineralization [17]. Other studies have shown changes in the soil biota. In the Brazilian Amazon, land-use change resulted in large reductions in fungal richness and an increased prevalence of generalist taxa [18]. In the Mexican rainforest, pasture establishment was associated with increased vulnerability of the soil microbial communities to seasonal changes [19]. However, few studies have so far investigated how differences in rangeland management affect soil physical, chemical, and biological quality, and related ecosystem services underpinning rangeland productivity in the humid tropics. The objectives of our study were to (1) characterize current livestock management practices in the geomorphologically diverse landscape of the Marqués de Comillas study area, and (2) understand how these practices affect soil quality and related supporting, regulating and provisioning ecosystem services. More specifically, we explored the direct and indirect effects of livestock management practices on soil-based supporting and regulating ecosystem services and their linkages to two provisioning ecosystem services: forage productivity (quantity of palatable plants) and forage quality (forage protein concentration), supplied by these grasslands (Figure 1). We established 31 plots in different landscape categories (alluvial areas and low-hills) in the Marqués de Comillas study area, documented management practices, and measured soil quality indicators, litter decomposition potential, and forage productivity and quality.
rangeland productivity in the humid tropics.
The objectives of our study were to (1) characterize current livestock management practices in the geomorphologically diverse landscape of the Marqués de Comillas study area, and (2) understand how these practices affect soil quality and related supporting, regulating and provisioning ecosystem services. More specifically, we explored the direct and indirect effects of livestock management practices on soil-based supporting and regulating ecosystem services and their linkages to two provisioning ecosystem services: forage productivity (quantity of palatable plants) and forage quality (forage protein concentration), supplied by these grasslands (Figure 1). We established 31 plots in different landscape categories (alluvial areas and low-hills) in the Marqués de Comillas study area, documented management practices, and measured soil quality indicators, litter decomposition potential, and forage productivity and quality.
We hypothesized that (1) the rangelands located away from the river in low-hills show lower soil quality and ecosystem services supply than the rangelands in alluvial areas; (2) the more intensively managed rangelands show a decrease in soil quality, supporting services (decomposition rate as a proxy for nutrient cycling), and provisioning services (forage productivity and quality); and (3) the rangelands with higher soil quality have higher decomposition rates (direct effect) and higher forage productivity and quality (indirect effect) than rangelands with lower soil quality.

Study Site
The study was conducted in Marqués de Comillas, a part of the Lacandona rainforest region, Chiapas, southeast Mexico (16°01′N, 90°55′W) ( Figure 2). The climate is warm-humid with an average annual temperature of 24 °C. Annual rainfall is around 3000 mm, with a dry period from February through April, with less than 100 mm rainfall per month [20]. A total of four geopedologic land units were identified in the approximately 50 km 2 area, based on variation in geomorphological and We hypothesized that (1) the rangelands located away from the river in low-hills show lower soil quality and ecosystem services supply than the rangelands in alluvial areas; (2) the more intensively managed rangelands show a decrease in soil quality, supporting services (decomposition rate as a proxy for nutrient cycling), and provisioning services (forage productivity and quality); and (3) the rangelands with higher soil quality have higher decomposition rates (direct effect) and higher forage productivity and quality (indirect effect) than rangelands with lower soil quality.

Study Site
The study was conducted in Marqués de Comillas, a part of the Lacandona rainforest region, Chiapas, southeast Mexico (16 • 01 N, 90 • 55 W) ( Figure 2). The climate is warm-humid with an average annual temperature of 24 • C. Annual rainfall is around 3000 mm, with a dry period from February through April, with less than 100 mm rainfall per month [20]. A total of four geopedologic land units were identified in the approximately 50 km 2 area, based on variation in geomorphological Land 2020, 9, 511 4 of 17 and pedological conditions within the larger landscape context of the study area: (1) fluvial terraces, (2) flood plains, (3) low-hills, and (4) karst mountains [21,22]. We combined fluvial terraces and flood plains into one category of alluvial areas, as both are characterized by deposition of fine sediments due to regular flooding of the Lacantún river. Alluvial areas represent 16% of the land and the main soil types there are Fluvisols, Cambisols and Luvisols. The low-hills represent 70% of the area, and consist of claystone-sandstone, siltstone-sandstone and sandstone-conglomerates. The main soil types found in the low-hills are Cambisols associated with Acrisols [21]. The karst mountains are mostly steep and soils are shallow. As these latter soils are unsuitable for pastures, they were excluded from our study. The two landscape categories (alluvial soils and low-hills) are easily recognized by farmers and are part of the discourse of local people when talking about soil management, and were therefore included in this study. pedological conditions within the larger landscape context of the study area: (1) fluvial terraces, (2) flood plains, (3) low-hills, and (4) karst mountains [21,22]. We combined fluvial terraces and flood plains into one category of alluvial areas, as both are characterized by deposition of fine sediments due to regular flooding of the Lacantún river. Alluvial areas represent 16% of the land and the main soil types there are Fluvisols, Cambisols and Luvisols. The low-hills represent 70% of the area, and consist of claystone-sandstone, siltstone-sandstone and sandstone-conglomerates. The main soil types found in the low-hills are Cambisols associated with Acrisols [21]. The karst mountains are mostly steep and soils are shallow. As these latter soils are unsuitable for pastures, they were excluded from our study. The two landscape categories (alluvial soils and low-hills) are easily recognized by farmers and are part of the discourse of local people when talking about soil management, and were therefore included in this study. Several vegetation types characterize the area, with rainforest still being the dominant one. Extensive areas of former old-growth forest have been converted to agricultural fields (croplands) and later to rangelands. The landscape matrix is highly dynamic and under high anthropogenic pressure. In 2011, 35% of the landscape in this region was covered by rangelands, while the rest of the land was covered by mature forest (52%), secondary forest (5%), croplands (4%), forestry plantations (<3%), cocoa shade farming (<1%), and oil palm plantations (<1%) [7].
Mayans and other human groups, the original inhabitants of the area, abandoned it more than 500 years ago. Human recolonization started in the early 1970s, when governmental programs encouraged immigration. The Marqués de Comillas municipality has 28 ejidos, a unique land tenure system of common property that emerged from the Mexican Revolution to distribute lands to poor landless peasants [23,24]. A typical productive unit of an individual farmer in this region has an average size of 30-50 ha, from which 20-40 ha are usually devoted to cattle ranching [25]. Several vegetation types characterize the area, with rainforest still being the dominant one. Extensive areas of former old-growth forest have been converted to agricultural fields (croplands) and later to rangelands. The landscape matrix is highly dynamic and under high anthropogenic pressure. In 2011, 35% of the landscape in this region was covered by rangelands, while the rest of the land was covered by mature forest (52%), secondary forest (5%), croplands (4%), forestry plantations (<3%), cocoa shade farming (<1%), and oil palm plantations (<1%) [7].
Mayans and other human groups, the original inhabitants of the area, abandoned it more than 500 years ago. Human recolonization started in the early 1970s, when governmental programs encouraged immigration. The Marqués de Comillas municipality has 28 ejidos, a unique land tenure system of common property that emerged from the Mexican Revolution to distribute lands to poor landless peasants [23,24]. A typical productive unit of an individual farmer in this region has an average size of 30-50 ha, from which 20-40 ha are usually devoted to cattle ranching [25].

Study Design, Plot Establishment and Plot Characteristics
In total, 31 exclosures (3 m × 3 m) were established in individual paddocks across four ejidos, Boca de Chajul, Playón de la Gloria, Galacia, and El Pirú, to measure biomass productivity and forage quality. The paddock selection process aimed to capture diversity in terms of landscape position, and past and actual management practices. Most rangelands in Marqués de Comillas are nevertheless located in the low-hills. Therefore, 26 of the exclosures were established in this landscape category. All exclosures were established on actively managed rangelands. Researcher and landowner decided together on the exact location of the exclosure after a joint visit to the farm. Given the heterogeneity of plant cover in the paddocks, care was taken that the exclosures were located such that the plant community structure inside and outside the exclosure was similar at the beginning of the cattle exclusion period in May 2018. For each exclosure we collected data on the presence of trees, biomass productivity, decomposition of standard material, and biological, chemical, and physical soil quality indicators.

Characterization of Livestock Management Practices
Semi-structured interviews were conducted with individual farmers to obtain information about the historical and current management practices (interview questions are provided in Appendix A). During a field visit with the owner we collected additional information about management decisions. Given that sometimes the official owner was not the one who made decisions or managed the land, additional in-depth interviews were conducted when necessary (for example, to other members of the family or to employees in charge of the rangeland). During the interview, the interviewee was asked to draw a sketch of the rangeland and to locate fences, water bodies, fragments of tree vegetation, and roads. The majority of farmers (28 of 31) practiced an incipient rotational system, which consists of reallocation of livestock from one paddock to another to allow the pastures to recover after a period of grazing (this was assessed in the path model through the variable grazing time (GT), were GT = 1 corresponds to continuous grazing). Implementation of rotational grazing systems is a relatively recent phenomenon, with the general objective of increasing productivity by ensuring that palatable plants receive sufficient resources and by enabling livestock to harvest forage more efficiently [26]. The rotational system (average grazing time per paddock and direction of the rotation between areas) was represented in the sketch, as well as the history of establishment (different parts of the rangeland were established in different years, as the rangeland area generally expanded over time). We also asked detailed questions about duration of the productive activity, fire use, agrochemical use, machinery use, stocking rate, duration of grazing period, grazing intensity (expressed as stocking rate and grazing periods), and type of production (calves, or adults for meat). When possible, information was corroborated by cross-checking. For example, years since the rangeland was established were cross-checked by Landsat images modelling, and size of farms was measured using ArcGis software.
We also calculated the Ecological Disturbance Index (EDI) according to Zermeño et al. (2015) as a general metric for the intensity of disturbance that considers various management practices: [27]. The EDI is based on size (SZ) (field surface in hectares under cattle ranching), duration (D) (years since the rangeland was established) and severity of agricultural management (SV). SV is based on five criteria, viz., fire frequency (FF), frequency of herbicide use (HF), use of tractor frequency (TF), grazing severity (GS), and tree cover (TC). FF refers to the number of burning events in the history of field use; HF is quantified as the number of agrochemical applications per annual productive cycle; TF is the number of times a tractor was used per year; GS refers to the stocking density (SD) (number of livestock heads per hectare) adjusted by grazing time (GT) (proportion of time under pasture use in an annual cycle); TC is a parameter that adopts values from 0 to 2, depending on the relative cover of trees within a given field: 2, tree cover 0-25%; 1, 25-50%; and 0, >50%.

Soil Quality and Litter Decomposition
We took four samples with a fixed volume from 0 to 25 cm depth in each plot and bulked them into one composite sample. Samples for the determination of pH, and total C, N and P were taken in October 2018, samples for the determination of available ammonium, nitrate and orthophosphate in March 2019, at the height of the dry season. Composite samples were stored in self-sealed plastic bags and kept refrigerated until being processed in the lab. A subsample was taken and frozen for fatty-acid analysis.
Chemical soil quality was determined through pH (in 1:10 soil: water solution); soil organic matter percentage by wet combustion [28]; total C and N with a AEO Flash 2000 elemental analyzer (Thermo Scientific) after HCl to remove carbonates; total P after macro-Kjeldahl digestion [29] and colorimetry with a Bran-Luebbe Auto Analyzer III [30]; and inorganic N (ammonium, nitrate) after 2M KCl extraction [31,32] and available P (orthophosphate) after Mehlich 3 extraction and colorimetry [33]. Physical soil quality was assessed through particle size analysis [34], sand-free water-stable macro-and microaggregates, bulk density, and water infiltration rate measurements. Micro-and macro-aggregates, were measured from an undisturbed soil sample taken in each plot using a metal cylinder (10 cm diameter, 0-25 cm depth) in October 2018. Samples were sealed and stored in jars to protect their structure until being processed in the lab. Sand-free water-stable aggregates were quantified in three fractions: macro-aggregates (0.250-1 mm and 1-4 mm) and micro-aggregates and particles smaller than 0.250 mm [35]. Bulk density was also determined from two similar 0-25 cm soil cores taken with a 10-cm diameter PVC cylinder. The core was sealed and stored until being processed in the lab where the cylinder was emptied, the soil weighed, oven-dried, and weighed again. Soil infiltration rate was measured near the exclosure plot with a single-ring infiltrometer (20 of cm of diameter, 14 cm height aboveground, and 2 cm belowground) according to Verbist et al. (2010) [36]. We made three measurements in each plot and calculated the average. Soil microbial functional group abundances were quantified through fatty-acid analysis from 3 g finely ground soil from the frozen composite sample of each plot, according to Frostegård et al. (2011) [37]. The fatty acids in phospholipid and neutral lipid fractions were identified with MIDI software and quantified from their retention times in relation to the internal standard (fatty acid methyl ester 19:0). Seventeen biomarkers were selected. In total, three biomarkers were used as indicators of fungal biomass, 18:2 ω6,9 (PLFA fraction) [38], 18:1 ω7 (NLFA fraction) and 16:1 ω5c (NLFA fraction), the latter specifically for arbuscular mycorrhizal fungi (AMF) [39]. A total of 14 PLFAs represented bacterial biomass, 12:0 3OH, i14:0, 16:0 2OH, 16:0 3OH [40], i15:0, a15:0, i16:0 [41], i17:0, a17:0, cy17:0, cy19:0 ω8c [39], i17:0 3OH [42], 10 Me17:0, and 18:0 2OH [40]. We also calculated the fungi: bacteria ratio as well as the Gram-positive: Gram-negative bacteria ratio by adding the abundances of the relevant biomarkers.
Decomposition rate of a standard material (to allow comparisons between plots) was measured as a proxy for the nutrient cycling potential of each plot (referred in the text as 'litter decomposition rate'). Decomposition rates were determined using the litterbag method [43]. Fresh leaf material was collected from the 31 study sites, dried and mixed to create a homogeneous dry mixture. Litterbags (mesh size of 2 mm in the upper side and 0.25 mm mesh in the bottom side) containing an initial dry weight of six grams of the standard material mixture were placed within each exclosure from June 2018 (month 0) to June 2019 (month 12). In total, four bags were collected from each plot in July 2018 (month 1), October 2018 (month 4), December 2018 (month 6), and July 2019 (month 12). After litterbag collection we carefully brushed off the soil residues, removed manually in-grown roots, and then oven-dried and weighed the remaining material. After weighing, the material was incinerated in a muffle oven to correct for soil contamination in the bags. Decomposition rate was calculated as the rate constant k (per day) based on the simple exponential model: dN dt = −kN, with N being the mass remaining at collection time, and k the decomposition rate [44].

Biomass Production and Forage Quality
Biomass production of palatable species and forage protein concentration were used to estimate the provisioning service of forage productivity and quality. Cattle-grazing exclosures were established by fencing 3 m × 3 m in plots and divided into four 1-m 2 subplots. We left out the outermost 50 cm of each exclosure to avoid border effects. A total of two subplots were harvested monthly (treatment 1) and the other two were harvested only at the beginning and at the end of the experiment (14 months) (treatment 2). Treatment 1 referred to a simulation of continuous grazing, where all plant biomass present was removed. Palatable and nonpalatable biomass were separated after careful examination and classification. Treatment 2 represented the cumulative biomass production in the absence of grazing effects and allowed us to evaluate if grazing affected plant biomass production. In total, four out of a total number of 434 monthly biomass samples, were discarded as cows had entered the exclosure after damaging the barbed wire (in different plots, at different times). Productivity estimates might therefore slightly underestimate total productivity. All vegetation (grasses, forbs and tree seedlings) was clipped at approximately 5 cm height with garden scissors or a small machete. Vegetation was separated into two groups: palatable and non-palatable species. Palatability of plants was established on the basis of local knowledge; when needed, information was corroborated with the literature [45]. Collected material was dried in a handmade oven until it reached constant weight and stored for further analysis. Biomass productivity with and without monthly clipping was expressed as cumulative dry weight (g) per unit area over 12 months. Forage quality was assessed as protein concentration, which is a factor that determines cattle weight gain. It was measured at the beginning of the experiment in June 2018, using a representative subsample of palatable vegetation in each plot. The subsample was ground into fine powder to determine total nitrogen after Kjeldahl method and calculate protein concentration, as its N concentration multiplied by 6.25 [46].

Data Analyses
We developed a conceptual or path model ( Figure 1) to test (1) whether landscape and history of use affected current livestock management practices and soil quality; (2) if livestock management practices had a direct effect on supporting, regulating and provisioning services; (3) the effect of biological, chemical and physical soil quality on decomposition and forage productivity and quality. Path analysis is a form of multiple regression analysis that evaluates hypothesized causal models by examining the relationships between a dependent variable and a set of independent variables. Given the sample size (31 plots) relative to the number of variables (39) collected in this study, we grouped the variables into seven categories: landscape position, age of the plot (as a proxy for history of management), current management practices, chemical, physical and biological soil quality, decomposition rate (as a proxy for nutrient cycling), and provisioning ecosystem services (forage productivity and quality). Landscape position included three separate variables: landscape category, distance from the river, and altitude. Distance from the river represents the beeline distance between the exclosure plots and the Lacantún river; altitude (as a proxy for landscape position) was calculated as the vertical elevation of the plot above sea level (masl).
For each variable category we selected 1-2 variables to be included in the path analysis. To reduce the number of variables we performed two Principal Component Analyses (PCA), one for management and one for soil quality ( Figure S1; Table S1). To avoid collinearity, a variance inflation factor (VIF) was calculated. When two or more variables within categories showed high VIF values, factors were eliminated and the PCA was run again. We kept the variables that resulted in higher percentage of variance explained. The first two components were used in subsequent analyses to determine the direct and indirect effects of management and soil properties on ecosystem services. Variables for final analysis were selected based on the following criteria: (1) After PCA, axes that explained at least 60% of variation were picked; (2) if two or more individual variables showed a strong negative correlation, one of them was eliminated; (3) variables that showed a strong correlation with two axes were eliminated. All data were normalized after transformation (mean = 0, standard deviation = 1) so that all responses Land 2020, 9, 511 8 of 17 had equal weight. Once individual variables were selected, standard β coefficients were calculated using multivariate regressions to evaluate the relationships between predictors and response variables.
We used simple regression models to (a) illustrate the statistically significant (p ≤ 0.05) correlations according to multivariate analyses (grazing time vs altitude, pH vs altitude, stocking rate vs altitude, pH vs grazing time, AMF vs altitude, pH vs grazing time, Gram+:Gram-ratio vs altitude, forage protein content vs AMF); (b) to examine the relation between the Ecological Disturbance Index and the measured provisioning ecosystem services (forage productivity and forage quality); (c) to examine the relationship between soil properties and management variables and the three proxies used to estimate ecosystem services: forage productivity, forage protein concentration, and decomposition. We used the scores of the first components of the two PCAs described above, which best explained the variation in soil properties (PO 3− , NH 4+ , NO 3− , OM, pH, sand content) and management (grazing time, herbicides, stocking rate, meat production, fire) to integrate the variation in soil quality and management.
To compare the effect of landscape category on forage productivity, we conducted one-way ANOVA and post-hoc multiple comparisons (after Bonferroni correction) by using the average values of dry biomass of palatable species collected in a year of sampling.
All statistical calculations (linear regressions, multiple regressions (method enter), PCA and one-way ANOVA) were completed with IBM SPSS statistical software version 25.

Rangeland Characteristics and Management
The 31 study plots (paddocks) represented different ages since establishment of livestock management, ranging from 4 to 40 years of livestock management. Altitude of the plots varied from 143 masl, at 0.17 km from the Lacantún river to 186 masl at 13 km from the river. The size of the farms varied between 5 and 50 ha (Table 1). A total of five plots were located on alluvial soils and 26 plots were located on low hills. Average stocking rate was 2.24 heads ha −1 , varying from 0.1 to 9 heads ha −1 . The majority of farmers rotated cattle among different paddocks (28 of 31). All farmers burned the land when the pastures were established, but only two farmers reported to have used fire in the last 5 years. Several herbicides were reported, belphosphate, glyphosate, picloram, and paraquat being most frequently used. Tractor use is a recent practice in the area and was reported in only three plots. There were two types of production: calf production, where farmers sell the calves when they stop drinking milk, and meat production, where farmers fatten the cattle until it reaches a certain weight (Table 1).

Soil Quality and Biomass Productivity
Soil properties showed high variation among plots. We observed in general low levels of available nutrients, low to medium levels of soil organic matter, and large variations in pH, bulk density, infiltration rate, stable macroaggregates, and bacteria: fungi ratios (Table 2). Litter decomposition rates showed a similar variation, close to one order of magnitude (0.003-0.016 day −1 ). The annual cumulative productivity of palatable species was also highly variable (114-1418 kg ha −1 ) in the non-grazing treatment, and in the monthly cut treatment, which was overall lower than the non-grazing treatment (18-852 kg ha −1 ). The protein concentration was relatively low (4.46-15.37%) ( Table 2).

Effects of Management on Soil Quality and Ecosystem Services
The results of the path analysis are depicted as the simplified final model with standardized β coefficients ( Figure 3; Table S2). There was a significant effect of altitude (a proxy for the landscape context in relation to the river) on pH, with rangeland plots in the alluvial soils having a higher soil pH than plots in the low hills ( Figure S2a). Altitude was not significantly related with soil organic matter content. Altitude was related negatively with the ratio of Gram+: Gram-bacteria ( Figure S2b) and positively with the relative abundance of AMF ( Figure S2c). Pastures located at low altitude also had higher stocking rates ( Figure S2d) and longer grazing times ( Figure S2e) than pastures in the low hills. Stocking rate was directly and positively related with pH ( Figure S2f) while the relationship between grazing time and pH was not significant. None of the soil biological, chemical, or physical factors were significantly correlated with decomposition rate. The abundance of AMF was negatively related with forage protein content. Altitude showed a tendency to be associated with higher productivity of palatable species, but the relation was not significant. No relation was detected between the age of the plot and management variables.
between grazing time and pH was not significant. None of the soil biological, chemical, or physical factors were significantly correlated with decomposition rate. The abundance of AMF was negatively related with forage protein content. Altitude showed a tendency to be associated with higher productivity of palatable species, but the relation was not significant. No relation was detected between the age of the plot and management variables.
After applying linear regression models, our data did not show any significant relation between the composite index Ecological Disturbance Index (EDI) and forage productivity and forage quantity. Additionally, no relation was found between the scores of selected PCA axes for management and soil properties and forage productivity, forage quality, and organic material decomposition.

Landscape Position Drives Farmers Management Practices, Soil Quality and Forage Quantity but not Quality
We identified seven management practices that are common to all farmers. Nevertheless, the intensity and incidence differ among the study plots, showing that each farmer has a unique system of livestock management partly driven by the location of his/her land. As hypothesized, the position of rangeland plots in the landscape had a strong effect on several soil characteristics, forage quantity and quality, and thereby on management decisions, i.e., grazing intensity ( Figure 3). Sites close to the Lacantún river had a higher pH and were more fertile, due to annual flooding, in terms of soil available P, OM and particle size distribution, than pastures on the low-hills as observed also by Suazo-Ortuño et al. (2015) and Navarrete (2018) [16,21]. As a consequence, farmers with pastures close to the river maintained higher stocking densities (number of cattle heads per ha) and longer grazing times than farmers on the low-hills because their plots tended to have a higher forage quantity ( Figure 3). Nevertheless, the ANOVA showed no significant differences among landscape categories, probably due to the high variability between individual paddocks (F = 0.90, p = 0.767).
Although higher fertility in alluvial soils could theoretically provide better options for intensification of cattle ranching practices close to the river, we did not observe a significant relationship between landscape position and forage productivity (only a marginally significant trend). This weak relation is probably due to the fact that livestock ranching in Marqués de Comillas is low input and even the more fertile lands cannot sustain a high production in the absence of After applying linear regression models, our data did not show any significant relation between the composite index Ecological Disturbance Index (EDI) and forage productivity and forage quantity. Additionally, no relation was found between the scores of selected PCA axes for management and soil properties and forage productivity, forage quality, and organic material decomposition.

Landscape Position Drives Farmers Management Practices, Soil Quality and Forage Quantity but Not Quality
We identified seven management practices that are common to all farmers. Nevertheless, the intensity and incidence differ among the study plots, showing that each farmer has a unique system of livestock management partly driven by the location of his/her land. As hypothesized, the position of rangeland plots in the landscape had a strong effect on several soil characteristics, forage quantity and quality, and thereby on management decisions, i.e., grazing intensity ( Figure 3). Sites close to the Lacantún river had a higher pH and were more fertile, due to annual flooding, in terms of soil available P, OM and particle size distribution, than pastures on the low-hills as observed also by Suazo-Ortuño et al. (2015) and Navarrete (2018) [16,21]. As a consequence, farmers with pastures close to the river maintained higher stocking densities (number of cattle heads per ha) and longer grazing times than farmers on the low-hills because their plots tended to have a higher forage quantity ( Figure 3). Nevertheless, the ANOVA showed no significant differences among landscape categories, probably due to the high variability between individual paddocks (F = 0.90, p = 0.767).
Although higher fertility in alluvial soils could theoretically provide better options for intensification of cattle ranching practices close to the river, we did not observe a significant relationship between landscape position and forage productivity (only a marginally significant trend). This weak relation is probably due to the fact that livestock ranching in Marqués de Comillas is low input and even the more fertile lands cannot sustain a high production in the absence of external inputs. It is also a very low-profit (or even non-profitable) activity, with a benefit: cost ratio of 0.24-0.99 for small farms with less than 60 hectares [7,47]. The alluvial soils cover 16% of the area but only 1% of the rangelands is found in the floodplains [22]. Instead, alluvial soils are preferred for agriculture (mainly maize, and beans) and, more recently in some ejidos, for oil palm plantations. Most agricultural activity in the floodplains is geared towards the production of subsistence and marketable food crops. Recently, African oil palm (Elaeis guineensis) is increasing in the floodplains to the extent that it might further reduce rangelands in this area. Since low-hills are perceived as poor soils, and farmers cannot afford fertilizers or chemicals for weed and pest control because of the low revenue, they rarely decide to use them for crop production. Instead, they use them for extensive cattle ranching requiring minimum investment and representing a low value but stable asset to be used when required. For example, El Pirú, located far away from the Lacantún river, only possesses land in the low-hills. The famers of this community rent land in Flor del Marqués, the neighbor community close to the river, to produce subsistence crops and their land is mainly dedicated to cattle production, ecotourism and, more recently, oil palm production.

Soil Properties Differed Widely with Landscape Position and among Plots
Altitude, a proxy for distance from the river and thus also for the influence of flooding and the groundwater table, affected mainly soil pH and microbial functional groups (Figure 3). Neutral to basic pH values in alluvial soils are likely explained by basic-cation accumulation in surface layers. We found a negative relation between altitude and the ratio Gram+: Gram-bacteria ( Figure 3), with the highest values of that ratio observed in the floodplains. It has been proposed that Gram-bacteria, which are consider as copiotrophs, use more labile, plant-derived carbon sources such as exudates, while Gram+ bacteria, considered as more oligotrophs, use more recalcitrant carbon sources [48]. However, our data did not support this pattern. Ho et al. (2017) noted that the simple classification of Gram+ and Gram-bacteria may be too crude to capture changes along multiple environmental gradients [49]. Furthermore, Gram+ bacteria are less sensitive to stress conditions as many of them form resistant walls [50] which may help survival under conditions of flooding. Observations by Unger et al. (2009) that the Gram+: Gram-ratio increased under flooding would support this hypothesis [51] Altitude had a positive effect on the abundance of AMF (Figure 3) likely due to a higher P availability in alluvial pastures. Negative relationships between phosphate availability and AMF biomass and functioning have been commonly reported [52]. Furthermore, intermittent anaerobic conditions in alluvial areas could have negatively affected the obligate aerobic AMF [53]; consistent with that explanation we observed gley layers and water accumulation in plots closest to the river.
Stocking rate was the livestock management practice with the highest explanatory values for soil quality, followed by grazing time and the only soil variable affected was pH ( Figure 3). These results were unexpected. Grazing intensity combines number of heads per hectare (stocking rate) and grazing time, and thus reflects cattle influence on biomass removal, biomass recovery, dung and urine deposition and trampling that in turn can have cascading effects on soil physical, chemical and biological properties beyond pH. However, grazing intensity did not have an effect on organic matter ( Figure 3). Soil OM in these rangelands was low, with the highest soil OM values measured being approximately half the values measured in successional and old-growth forest sites of the same region [54]. These results suggest that this land use has led to soil impoverishment. Conversion of forest to pasture often results in large losses of soil organic matter and nitrogen [21,55]. However, there is less literature on how duration of land under pasture further impacts on soil organic matter levels. A meta-analysis by Zhou et al. (2017) showed a significantly negative correlation between pasture age and soil organic carbon and nitrogen [56]. These arguments may also partially explain why landscape position differences to some extent translated into differences in soil quality but ultimately had no impact on litter decomposition rates and on the productivity of palatable plants.

Land Use History Was Not Related to Actual Management Practices and Soil Quality Indicators
Surprisingly, age of the plots did not show a relationship with any management practices, and with ecosystem services, suggesting that the ecosystem services measured in our study have still not declined with time under pasture use. Our results contrast with a similar study conducted in plots with similar age range in the dry forest regions of Mexico where it was found that grazing time up to 40 years promoted forage production but decreased forage protein concentration [57]. The authors suggested that farmers apply regular burning to increase productivity in the short term (possibly because production is P-rather than N-limited) but that that practice reduces forage quality, because of N losses through volatilization in hot fires. In our study area, fire was only very infrequently used after the rangelands had been created. All plots showed low availability of essential macronutrients N and P, hence explaining the low forage productivity and low forage protein concentration. Nutrient availability thus limits forage and/or milk production despite adequate water availability in the region. Our results support previous observations that cattle raising is hardly productive in this study region.
In our study we were interested in the diversity of ways in which farmers managed their pastures and how landscape, and past and current management relate to soil quality and indicators of ecosystem services. This additional layer of farmer diversity may have obscured significant relations between grazing history, soil quality and ecosystem services that would have appeared with a plot selection that was exclusively focused on the temporal dimension. The lack of significant relationships with history fits, however, with farmer perceptions on the efficacy of their practices, as they have various ways to manage their rangelands over the long term, none of which are clearly superior over other sets of management practices.
Our data did not show any significant relation between EDI and productivity and forage quality. A major difference between our study and that by Zermeño et al. (2015) is that the latter authors focused on forest regeneration potential, whereas we focused on the ecosystem services delivery by the grasslands themselves [27]. A further reason why in our study EDI was not useful is that we focused on pastures only, whereas Zermeño et al. (2015, 2016) compared widely different land use types [27,58]. We, therefore, suggest that a separate analysis of the different components of management, rather than deriving a composite indicator as EDI, is a better approximation to understand the long-term delivery of provisioning services within these rangelands.

Implications and Recommendations for Sustainable Management
We hypothesized that time under rangeland management would have a negative effect on soil quality and provisioning ecosystem services. Nevertheless, the results did not support such relationship, suggesting that farmers are so far able to maintain the (marginal) value of their pastures. However, the decline in SOM is a warning sign as it may likely compromise ecosystem services in the future.
In our study, the quality of forage was low even in the youngest, more fertile, alluvial soils. Also, we found that the relative abundance of AMF had a negative effect on the protein concentration of the pasture. While AMF may be stimulated by low P availability, situations in the low-hills in which both P and N availability are very low may result in N immobilization in the AMF mycelium [52,53,59], thereby exacerbating N-limitation by plants and thereby reducing forage quality and protein content ( Figure 3). Low forage quality does not allow productivity of cattle to be sustainable. Farmers do not apply fertilizers of any kind to these pastures so natural low fertility is expected to worsen with time. One adult cow needs 5 tons of forage per year so one hectare of these pastures could not maintain a single cow. The best plots in our study could sustain 0.3 and the worse 0.03 adult cows ha −1 with the forage biomass measured. This low forage productivity may explain why the region has focused on calf production, that requires less forage that maintaining adult cows. Current trends are to send calves to feedlots as soon as they reach the desired weight. According to Sagarpa (2014), Mexico produced an average of one head of cattle per 12.69 ha −1 of pasture per year, the state mean for Chiapas requires on average 1.80 ha head −1 year −1 (but with a wide interval of 0.80-18.90 ha head −1 year −1 ) [60]. Thus, the plots evaluated in this study have a productivity, both in forage and cattle, within the range reported for this region.
The majority of farmers have limited economic resources and keep the activity mostly as a self-subsistence strategy (savings) and, as stated in the interviews, a way of expressing their identity. Given the low rewards of extensive cattle ranching, oil palm plantations promoted and subsidized by state governmental programs are growing in these low-fertility soils. In oil palm plantations farmers compensate for the low soil fertility by heavily using fertilizers. The use of fertilizers, herbicides and pesticides is widespread due mainly to the credits and subsidies that currently exist for this profitable crop in contrast to subsistence crops. From this point of view, extensive cattle farming would require sustainable intensification practices to increase productivity [61][62][63] and supporting farmers with governmental training programs to increase profitability and farmers resilience.

Conclusions
Our results showed that management practices related to cattle ranching in Marqués de Comillas have so far not strongly degraded soil-related supporting and regulating services. Nevertheless, the delivery of provisioning services has been very poor. In human-modified landscapes, engaging local actors and sustainable practices could prevent further poverty, environmental degradation, and deforestation, and allow the transition towards more sustainable cattle production.