The high plain subregion of Colombia is between 100 and 300 m above sea level. Precipitation ranges between 1,500 and 2,500 mm/year, distributed between April and November, with an average temperature of 27°C and maximum and minimum values of 36 °C and 21°C, respectively (Amézquita et al., 2013; Rivera et al., 2013).

The soils in this area are classified by the US Soil Taxonomy as oxisols and are relatively well-drained, with very low fertility and high aluminum saturation. The soil pH ranges from 4.0 to 4.5, with very low organic matter concentration (1.5%–2.0%), phosphorus (range of 1–2 mg kg^-1 Bray II), and exchangeable bases (Ca, Mg, and K), the latter being usually less than 0.1 cmolc kg^-1 and exchangeable aluminum saturations between 60% and 90% (Pérez-Lopez et al., 2019). From the physical point of view, the soils are naturally dense, with the bulk density ranging from 1.4 to 1.65 g cm^-3. The texture ranges from clay loam to sandy textures; they are fragile to mechanization processes, with high susceptibility to compaction, erosion, and surface sealing, and have low infiltration and water storage capacity (Pérez-Lopez et al., 2019).

The study was conducted between 2012 and 2014 at AGROSAVIA’s Carimagua Research Center of the Colombian Agricultural Research Corporation in the municipality of Puerto Gaitán (4° 19′01″N; 72° 04′59″W) in the Meta (Colombia). A 70-ha plot with Urochloa humidicola was selected with indicators of moderate degradation, such as low forage yield of less than 700 kg DM ha^-1year^-1, soil compaction, and the presence of termite mounds. A general review of the plot showed the need for pasture renovation and infrastructure adaptation to ensure enhanced management for breeding cattle production. The management strategies evaluated are given as follows:

A closed-chamber technique (CCT) was used to determine the net GHG fluxes: methane (CH₄), nitrous oxide (N₂O), and carbon dioxide (CO₂). Sampling was performed monthly and always in the morning hours in each system during 2013 and 2014 for a total of 18 samplings. Three rings were inserted into the soil in each plot and remained throughout the study period. At the time of sampling, the chamber lid was placed, and 20 mL air sample was taken at 0, 20, and 40 min after placing the lid. A 15-mL sample was dispensed in vacuum tubes and sent to the greenhouse gas laboratory at CIAT. Concentrations of CO₂, N₂O, and CH₄ taken in 2013 were determined in a Shimadzu GC-14A gas chromatograph, equipped with a flame ionization detector for methane analysis and an electron capture detector for nitrous oxide analysis. For the determination of CO₂ concentrations, a Qubit Systems S151 gas analyzer with infrared detection (IRGA) was used. In 2014, a Shimadzu GC-2014 was used that read all three gases simultaneously.

The fluxes were calculated following the GRACEnet methodology (Parkin and Venterea, 2010). Concentration in ppm (v/v) was transformed into mass units according to the ideal gas equation. Once the mass was calculated, the flux was estimated using the HMR function (Pedersen et al., 2010; Pedersen A. R., 2015) by the R 3.2.1 software (R Core Team, 2015). If the concentrations vs. time relationship generated a nonlinear trend, the Hutchinson and Mosier (1981) equation was used. When the trend was linear, the linear regression equation was used. Flows outside the range reported in the literature were discarded from the analysis. The area under the curve was interpolated to estimate the cumulative emission per period. Cumulative emissions of the three gases were converted to CO₂-equivalent values (t ha −¹ period −¹) using global warming potential values of 25 and 298 for methane and nitrous oxide, respectively.

SOC was determined from soil chemical and physical sampling, i.e., C concentration (%) and bulk density. Sampling was performed at two depths: 0–10 and 10–20 cm. The C content was determined by applying the following equation: T O C   k g / m 2 = ρ * f C * P m , (2)

where ρ = soil bulk density (kg m^-3), fC = carbon concentration (kg C kg soil−¹) by Walkley & Black, (1934), and Pm = sampling depth (m).

During the dry season (December–March) and rainy season (April–November), the following parameters were evaluated in the pastures.

• Forage yield was quantified by placing a frame (0.5 m × 0.5 m) in 10–20 sites according to the pasture management condition, and the grass was cut with a sickle or scissors. Subsequently, the green weight of the forage was determined using an electronic scale. Forage moisture was quantified after drying at 70°C for 3 days in a forage drying oven and used to calculate forage dry matter content.

Soil chemical and physical variables were analyzed considering a completely randomized design. The fixed factors were the systems (management strategies), time, and depth, and the random factors were associated with the replicates. For the analysis of the data, an analysis of variance model with repeated measures over time was used when the temporal autocorrelation was significant; in this case, different covariance matrices were evaluated using Proc GLIMMIX of SAS version 9.4 (SAS, 2016). Normality and homoscedasticity assumptions were evaluated based on Q–Q plots and model residual vs. predicted plots.

For the construction of indicators associated with the soil’s physical (water regulation), chemical (soil fertility), and biological (biodiversity) variables, a principal component analysis (PCA) was performed using the ADE4 1.4.5 library (Dray and Dufour, 2007) from R software (R Core Team, 2019). Calculations were carried out following the methodology proposed by Velásquez, Lavelle, and Andrade (2007).

An initial PCA allows for identification of the variables that best discriminate among the different management strategies. Variables with significant contribution (>50% of the maximum value) to either of the first two principal component axes were selected, and their contribution to PCA axes 1 and 2 was multiplied by the overall variability explained by each PCA axis to generate a weight factor for each variable. Values for each variable were then multiplied by their corresponding weight factor and summed to generate a raw sub-indicator value using the following formula: I i 1 = F 1 × α I a + β I b + γ I c ·   ·   · + F 2 × α I I a + β I I b + γ I I c ·   ·   · . I_i1 = value of the ecosystem service indicator at plot 1.F₁ = % of the value of the variance explained by axis 1 of the PCA of the indicator data group I _i .αI, βI, and γI = respective contributions of variables a, b, and c to factor I.a, b, and c = values of the variables measured at plot 1.

Finally, the indicators were rescaled between 0.1 and 1 to compare ES such as soil fertility, water functions, and biodiversity (soil macrofauna) associated with management strategies. For the climate regulation indicator, the CO₂-equivalent (t ha −¹ period −¹) was estimated from GHG fluxes (Rodríguez et al., 2013; Lavelle et al., 2014) and rescaled between 0.1 and 1 for comparison, and stock SOC values were rescaled between 0.1 and 1.

The analysis of variance of the indicators was performed considering a complete randomized design where the fixed factors were the management strategies and time (years) and the interaction strategies x time. The random factor was associated to the replicate.

Soil Ca, Mg, and P concentrations were affected by management, time (year), and sampling depth. The highest Ca, Mg, and P concentrations were found at 0–10 cm depth, where systems received corrective applications between 1 and 1.5 t ha−¹ of dolomite lime and fertilizers. The concentrations of Ca, Mg, and P ranged from 1.02 to 1.97 cmolc kg−¹, 0.24 to 0.48 cmolc kg−¹, and 6.1 to 10.9 mg P kg−¹ in the pasture–crop system during the time of study. The concentration of Ca, Mg, and P ranged from 1.05 to 2.56 cmolc kg−¹, 0.23 to 0.59 cmolc kg−¹, and 1.8 to 7 mg P kg−¹, respectively, in the tree management system. The concentration of Ca, Mg, and P ranged from 0.46 to 1.8 cmolc kg−¹, 0.14 to 0.37 cmolc kg−¹, and 2.1 to 2.6 mg P kg−¹, respectively, in the till fertilizer system, and the values of Ca, Mg, and P ranged from 0.68 to 0.88 cmolc kg−¹, 0.39 to 0.59 cmolc kg−¹, and 1.1 to 1.6 mg P kg−¹, respectively, in the forest patch management system. According to the PCA, the behavior of soil chemical properties is explained by a total (components 1 and 2) of 56.2%, 50.4%, and 46% for the years 2012, 2013, and 2014, respectively. Trees and pasture–crop rotation management had the highest values of nutrient concentrations.

In general, during the 3 years of evaluation, there was a separation of pH and Al in the two depths sampled and a positive relationship of some elements such as P, Ca, and Mg in the systems with the addition of chemical fertilizers (Figure 1).

For the soil fertility indicator, there were significant differences regarding the time (year) x system effect. In 2012, the pasture–crop rotation system that included fertilizer additions and amendments had the highest indicator values compared to the other management systems (Figure 2). In 2013, the pasture–crop system still had the highest indicator values. The areas with the lowest values corresponded to plots that had lower rates of fertilization such as pasture areas managed only with tillage and fertilization (till–fert). In 2014, the trend was the same as in the previous year; the plots with the highest values corresponded to those that received amendments and fertilizers (pasture–crop and trees).

Total porosity was more strongly correlated with mesoporosity and inversely related to microporosity (Figure 3) in all years; however, the strength was not the same across years. Although texture is related to porosity, the contribution to principal components is minimal (Figure 3). Mesoporosity had variations originated by management systems and soil use change in the years of evaluation; its values ranged from 4.98% to 15.6%; the highest values were found for the first year, mainly due to the action of tillage in soil management involving pastures and annual crops (pasture–crop).

The water function indicator was possibly (alpha<0.1) affected by the management system. The highest values were observed in the pasture–crop and forest patch management systems, which had values above 0.5 for the years 2012 and 2014. In 2013, this indicator decreased in all systems; however, the forest patch and pasture–crop systems still had the highest values (Figure 4).

In 2012, the Blattodea order, which is composed primarily by termites, had the highest density (individuals m^-2) at a depth of 0–10 cm, representing 35% (4,613 individuals m^-2) of the population, followed by earthworms with 28% (3,632 individuals m^-2). At the depth of 10–20 cm, this trend was similar, where termites were also the dominant organisms; however, the population decreased with depth, corresponding to 15% with 2032 individuals m^-2, followed by ants (Hymenoptera) and beetles (Coleoptera).

According to the PCA, the first two components of soil macrofauna variables explained 55.1, 56.0, and 54.6% of the variability for the years 2012, 2013, and 2014, respectively (Figure 5). In 2013, a high correlation between termites and total macrofauna abundance (individuals m^-2) was observed with Blattodea being the order with the highest density, with 17,984 individuals m², representing 72% of the total abundance. For 2014, even though there was a decrease in density, termites were the taxonomic unit with the highest correlation with total abundance and 9,139 individuals m^-2 representing 69% of the population.

In 2012, the till–fert system tended to differentiate from the other management systems with an average value corresponding to 0.32 (Figure 6). In 2013, the forest patch system had a higher value than others, with an indicator value of 0.6. For 2014, no significant differences were found among systems.

The goal of the systems was to improve zootechnical indicators of the cattle breeding system with financial sustainability and animal welfare. Degraded U. humidicola pastures, without intervention during the rainy season (April–November), had a forage yield of 435–839 kg DM ha−¹ and yields ranging from 389 to 671 kg DM ha−¹ during the dry season. The renovation of U. humidicola by means of the crop–grassland strategy increased forage yield from 1,324 to 1929 kg DM ha−¹ in the rainy season and from 869 to 1,152 kg DM ha−¹ in the dry season. When pasture recovery was carried out conventionally with soil tillage and fertilization, yields ranged from 1,232 to 1893 kg DM ha−¹ in the rainy season and from 749 to 904 kg DM ha−¹ in the dry season (Table 3). The forage quality of the U. humidicola pastures showed crude protein (CP) between 3.1% and 6.6% and forage degradability between 51.6% and 63.6% in the dry season. In the rainy season, CP ranged between 5.3% and 11%, and degradability between 60.7% and 77.5%. Regarding corn forage quality, it was found that the CP contents varied between 7.0% and 8.5%, with degradability between 58.5% and 62.1%.

Forage production from annual crops also showed outstanding results. Maize reached green forage yields of 35.1 and 25 t ha−¹, with the differences due to the genotypes used. Forage sorghum Corpoica JJT18 had yields that varied between 39.6 and 50 t ha−¹. Sorghum also has the possibility of generating a second crop by regrowth, which can represent 36%–45% of additional forage to conserve. Regarding soybean forage production, the Corpoica Taluma 5 variety reached yields between 18.5 and 26.5 t ha−¹, the Panorama 29 variety, between 17.1 and 26.4 t ha−¹, and the Corpoica Sabana 7 variety, 27.7 t ha−¹. Corpoica JJT18 sorghum provided CP values of 6.0–7.6% and degradability between 58.1% and 69.7%. For soybean varieties, CP values ranged from 12% to 24% and degradability from 56.7% to 71.1% (Table 4).

The management and environmental conditions prevailing during the study favored the growth and development of forest species. Diameter at breast height (DBH) recorded in 2015 showed that the highest values were reached by Acacia mangium with 58 cm, followed by Piptadenia opacifolia with 45 cm, and finally, Cassia grandis with 25 cm. With the differential development of the species, ground shade projection related to the growth and architecture of the canopy that generates the microclimate for animal comfort was similar for the first two species and less for the last one. Species growth also plays an important role in C fixation in the wood biomass accumulated, soil and water conservation, biodiversity, and the beauty of the rural landscape.