Introduction

Agriculture is increasingly expected to contribute to ecosystem services beyond the provisioning of food, feed, fiber, and fuel (Robertson et al. 2014). Climate regulation services from agricultural lands are attractive for their known co-benefits, such as improved soil and water quality (Lal 2010). Reduced soil disturbance, increased use of perennial and woody crops, more efficient use of fertilizer N, and improved animal and manure management have been endorsed as promising strategies for mitigating greenhouse gas (GHG) emissions from agricultural sources (Smith et al. 2008). Applying these GHG mitigation strategies for optimal efficacy requires thoughtful integration of practices in time and space in the context of available resources and anticipated climate change (Wilkes et al. 2013).

Integrated crop-livestock (ICL) systems have been recognized for their capacity to mitigate GHG emissions through the accumulation of soil organic C (SOC) (Franzluebbers et al. 2012; Carvalho et al. 2014). Diverse crop rotations, use of perennial forages, residue retention, and addition of animal manure contribute to increased C inputs relative to outputs, providing the means to increase SOC over time (Silva et al. 2014). Full GHG accounting of ICL systems, however, is largely unknown, as the documentation of soil C stocks, soil-atmosphere methane (CH4) and nitrous oxide (N2O) fluxes, enteric CH4 emissions from livestock, and CO2-related emission from inputs and farming operations require intensive measurement campaigns, often over multiple years. When compiled and expressed on a CO2-equivalent basis, these data provide an estimate of the net global warming potential (GWP) of a production system (Robertson et al. 2000).

Previous GWP evaluations of ICL systems have been conducted in Brazil, where pasture-based production of cattle (Bos taurus) is common (Cardoso et al. 2016). In southern Brazil, Ribeiro et al. (2019) found ICL systems under a black oat (Avena strigosa Schreb.) – soybean (Glycine max L.) rotation and no-till management reduced net GWP under moderate grazing compared to intensive grazing (0.09 vs. 4.92 kg CO2equiv. ha-1 yr-1, respectively). Moreover, their 3.5-year study highlighted the importance of SOC accumulation under ICL management, which offset up to 98% of GHG sources. In an investigation of pasture restoration strategies, Figueiredo et al. (2017) observed conversion of degraded pasture to a crop-livestock-forest integration (CLFI) system resulted in a negative GWP (i.e., net GHG uptake) after accounting for C sinks in soil and woody biomass over a 10-yr period. In aligned research, decreased climate impacts from CLFI systems were accompanied by lower production costs, improved quality of employment, and increased investment in regional economies when compared to conventional systems without integrated enterprises (Costa et al. 2018).

Though climatically distinct from the subtropical conditions of Brazil, the semiarid northern Great Plains of North America shares common attributes of extensive beef cattle production and abundant grazing resources (USDA 2017; Statistics Canada 2016). Non-irrigated cropland is common throughout the region, occupying over 51 Mha (Padbury et al. 2002). Moreover, cropland is often juxtaposed to rangeland, providing agricultural producers the potential to utilize cropland for livestock when rangeland forage conditions are of poor nutritive value or in limited supply. This context favors ICL systems, which by virtue of multiple enterprises, can serve to buffer economic downturns from fluctuating market and weather conditions common throughout the region (Hendrickson et al. 2008).

Prior investigations of ICL systems in the northern Great Plains of North America have evaluated management effects on grain and animal production (Kumar et al. 2019; Tanaka et al. 2005; Karn et al. 2005), soil and water attributes (Singh et al. 2020; Faust et al. 2018; Ryschawy et al. 2017), and GHG flux dynamics (Abagandura et al. 2019; Barsotti et al. 2013). To date, no evaluation has been conducted in the region to quantify ICL system effects on net GWP. Therefore, we sought to determine net GWP for an ICL system in south-central North Dakota USA by coupling long-term assessments of soil C stocks with intensive measurements of GHG emissions from soil, beef cattle, inputs, and field operations over a three-year period. We anticipated GHG emissions from beef cattle inclusion would be offset by enhanced SOC accrual in the ICL system, resulting in a similar GWP to ungrazed cropland. Accordingly, we hypothesized net GWP would not differ between an ICL system with beef cattle and an ungrazed cropping system.

Methods

Site and treatment description

The site was located approximately 5 km south of Mandan, ND USA (46.77, – 100.90; 593 asl) with gentle rolling topography (0–3% slope) and fertile Mollisol soils (Fine-silty, mixed, superactive, frigid Typic and Pachic Haplustolls; USDA 2020a). Soil conditions in the 0–0.3 m depth prior to 2014 were non-saline (0.29 ± 0.05 dS m-1) and slightly acid (pH1:1 = 6.10 ± 0.31), with soil C and N of 20.7 ± 2.3 g C kg-1 and 2.0 ± 0.2 g N kg-1, respectively. Long-term mean annual precipitation and temperature at the site is 456 mm and 4 °C, respectively, with an average 135 d frost-free period. Previous land use at the site was perennial grass pasture [crested wheatgrass; Agropyron desertorum (Fisch. ex. Link) Schult.].

An ICL study was established at the site in 1999 to evaluate long-term outcomes from integrating cattle with annual crops. Since initiation, the study has been conducted in three phases, with each phase reflecting changes to crops and crop sequences as production goals evolved. Crop rotation details associated with Phase I (1999–2005) and Phase II (2006–2014) may be found in Tanaka et al. (2005) and Liebig et al. (2017b), respectively. Beginning in 2015 (Phase III), a three-year adaptable management system was adopted, including intercropped corn (Zea mays L.)/soybean (Glycine max L.) (Year 1), spring wheat (Triticum aestivum L.) interseeded with a mix of perennial and annual cover crops (Year 2), and regrowth of the perennial and any surviving annual cover crops (Year 3) (Fig. S.1). The intercrop was composed of alternating corn and soybean rows with 0.38 m spacing. The interseeded forages included intermediate wheatgrass (Thinopyrum intermedium (Host) Barkw. & D.R. Dewey subsp. intermedium)/timothy (Phleum pretense L.)/alfalfa/hairy vetch/oilseed radish (Raphanus sativus L. var. oleiformis)/red clover/chicory (Cichorium intybus L.). The forages were cross-seeded at a 30º angle into standing spring wheat during early tiller formation. Due to poor production, timothy was included in the cover crop mixture only in 2017. When perennials failed to establish or experienced extreme winter kill in the spring wheat, proso millet (Panicum miliaceum L.) was seeded as a forage the following year. In addition to the grazed rotation outlined above, a three year ‘crop only’ rotation was included using a sequence of corn, spring wheat, and soybean (Fig. S.1). For purposes of this report, the contrasting crop rotations will be referred to as Grazed and Ungrazed treatments.

Crops in each rotation were present every year, and rotations were replicated four times in a randomized complete block design. Plot sizes for Grazed and Ungrazed treatments were 2300 and 500 m2, respectively. No-till planting techniques, coupled with herbicide weed control and application of synthetic fertilizer, were used for the duration of the study following practices used by local producers. Applied fertilizer amounts were the same for Grazed and Ungrazed treatments, and crops were grown under rainfed conditions. A generalized chronology of field activities both treatments may be found in Fig. 1, while additional details concerning production inputs and field operations are included in the Appendix (Table S.1).

Fig. 1
figure 1

Chronology of field activities by crop phase for grazed and ungrazed cropland in an integrated crop-livestock study near Mandan, ND USA, 2016–2019

Full size image

Grazed treatments were seasonally stocked by Angus yearling cattle in the fall of each year at rates of 0.41–0.52 ha steer-1 (0.89–0.93 animal unit equivalent yr-1). Access to crops and crop residue was controlled using electric polywire fencing, which was moved periodically to provide access to forage, either in place or from baled hay harvested during the growing season (as outlined below). Cattle were not back-fenced during grazing periods but permitted to regraze formerly grazed areas. Stocking period duration was 30, 118, and 26 d in 2016, 2017, and 2018, respectively, each in a single grazing cycle. Grazing was extended in 2017 due to favorable weather conditions in the fall coupled with the decision to graze standing corn due to poor grain production. A water source and shelter were located at the end of each replication within Grazed treatments.

Biomass and cattle measurements

Grain yield was measured at crop senescence from a 1.5-m strip in each plot with a Wintersteiger plot combine (Wintersteiger, Inc.; Ried im Innkreis, Upper Austria). Protein content of collected grain was measured on dried samples using a Foss Infratec Nova analyzer (Foss, Inc.; Hilleroed, Demark). Pre-and post-grazing biomass measurements were taken for cover crop, corn stalks and residue, and wheat residues using nine 0.25 m2 frames plot-1 to quantify the biomass available for consumption by beef cattle. Following an adaptive management approach, harvested cover crop bales from the grazed plots were fed alongside cover crops in controlled strips in 2016. In 2017, standing corn was grazed, with supplemental forage from cover crop residues and bales, which were simultaneously fed. In 2018, only grazing of cover crops occurred. Estimates of intake were made using the difference between pre- and post-grazing biomass, corrected for dry matter. These were considered maximum values of intake, since some biomass was lost to trampling, consumption by insects and wildlife, and high winds and not consumed by cattle. In instances when post-grazed sample grids contained greater biomass than pre-grazed samples, dry matter intake was predicted using an equation based on live weight and average daily gain (Minson and McDonald 1987). To obtain average live weight, cattle were weighed two consecutive days following an overnight fast of water and feed at the beginning and end of grazing. The weight of each day was averaged, then the mean liveweight was calculated as the ending weight minus the initial weight, divided by two, then added to the initial weight. Average daily gain was calculated using the difference in final weight minus initial weight, divided by the number of grazing days.

Conversion rate of CH4 as a percent of gross energy intake was estimated for each animal following guidelines outlined by Benchaar et al. (1998) using digestible organic matter, dry basis of the feed, and intake levels as a proportion of body weight using the dry matter intake calculated as previously described. Digestible organic dry matter was obtained for each crop or cover crop by description of percentage of each crop and published nutrient values. Gross energy intake was calculated from the predicted dry matter intake and the gross energy of the crops (data not shown). Values of CH4 production were taken as enteric fermentation for calculation of net GWP (reviewed below).

Soil carbon measurements

Soil carbon stocks were determined over a 15-year period (1999–2014) in both treatments as outlined in Liebig et al. (2020b). Briefly, soil cores were collected in the spring of 1999 and 2014 to a depth of 0.91 m using a truck-mounted hydraulic probe. At each sampled location, three cores were composited by depth increment (0–0.08, 0.08–0.15, 0.15–0.31, 0.31–0.61, and 0.61–0.91 m), placed in a double-lined plastic bag, and kept cool until processing. Samples were dried at 35 °C, mechanically ground to pass a 2.0 mm sieve, and split for pH, C and N, and air-dry water content analyses. Soil pH was estimated from a 1:1 soil-water mixture (Watson and Brown 1998), total C and N by dry combustion (Nelson and Sommers 1996), and soil inorganic C (SIC) by volumetric calcimeter for samples with pH ≥ 7.2 (Loeppert and Suarez 1996). Soil organic C was calculated as the difference between total C and SIC. Data were expressed on an oven-dry, volumetric basis using field measured soil bulk density (Blake and Hartge 1986). Soil organic C stocks for both treatments in 1999 and 2014 were calculated on an equivalent mass (EM) basis assuming soil profile mass of 15123 Mg ha-1, approximating the soil mass in the 0–0.91 m depth (Ellert and Bettany 1995).

Gas flux measurements

Soil-atmosphere CH4 and N2O flux was measured from Grazed and Ungrazed treatments from 14 Oct 2016 – 16 Oct 2019 using static chamber methodology (Hutchinson and Mosier 1981). Within each experimental unit, gas samples were collected from a single two-part chamber as described by Liebig et al. (2010). Collar bases for each chamber were occasionally removed to accommodate field activities (e.g., planting, haying, grain harvest). Collars were reinstalled within 24 hr after field activities 1–2 m from original collar placement. Each collar occupied an area of 0.03 m2.

Gas samples from chamber headspace were collected with a 20 ml syringe at 0, 30, and 60 min after installation (10:00 each sampling day). Once collected, gas samples were transferred to 12 ml evacuated Exetainer glass vials sealed with butyl rubber septa (Labco Limited, Buckinghamshire, UK). Gas fluxes were measured every 7–10 d when the soil was not frozen, and every 14–21 d throughout the winter. Concentration of N2O and CH4 inside each vial was measured by gas chromatography using a Varian Model 3800 and Combi-Pal auto-sampler. Description of the gas chromatograph, calibration protocols, and flux rate calculations are described elsewhere (Liebig et al. 2019; Phillips et al. 2009).

Weather measurements

Daily weather data over the course of the study was retrieved from a North Dakota Agricultural Weather Network (NDAWN) station within 0.5 km of the research site (NDAWN 2020). A National Oceanic Atmospheric Administration (NOAA) weather station 5 km north of the study was used to document precipitation received as snow (NOAA 2020).

Net global warming potential and greenhouse gas intensity estimates

Estimates of net GWP and greenhouse gas intensity (GHGI) were inclusive of assessments made at the beginning of fall grazing in 2016 through completion of the 2019 cropping season, but prior to scheduled grazing in 2019. Accordingly, assessments of cattle weight gain and enteric fermentation were associated with years 2016, 2017, and 2018, whereas crop-related assessments were associated with years 2017, 2018, and 2019. This assessment framework was selected to align with measurements of soil-atmosphere gas flux.

Emissions associated with seed, herbicide, and fertilizer production were calculated using conversions developed by West and Marland (2002) and Lal (2004). The same resources were used to calculate field operation emissions and were inclusive of planting, herbicide application, mowing, swathing, raking, baling, and grain harvest. Enteric fermentation was taken as CH4 production (g CH4 d-1) from beef cattle calculated as the product of gross intake energy following guidelines and conversions outlined by Melgar et al. (2019). Differences in SOC stocks between the 1999 and 2014 samplings were taken to represent net soil-atmosphere C exchange for each treatment. Contributions of soil-atmosphere CH4 and N2O flux to GWP were based on annual flux rates using beginning and ending dates in mid-October. It should be noted that the evaluation did not include GHG emissions associated with beef cattle before and after their placement on the experimental plots, nor did it include emissions associated with polywire fence, water pumping, and mineral supplements.

Net GWP was calculated for Grazed and Ungrazed treatments as the sum of emitted CO2 equivalents from eight factors:

$${\text{Net GWP }} = {\text{ SP }} + {\text{ FP }} + {\text{ HP }} + {\text{ FO }} + {\text{EF}} + {\text{ CH}}_{4f} + {\text{ N}}_{2} {\text{O}}_{f} + \, \Delta {\text{SOC}}$$
(1)

where SP, FP and HP addressed seed, fertilizer and herbicide production, respectively; FO addressed field operations; EF was estimated CH4 emission from enteric fermentation (Grazed treatment only); CH4f was soil-atmosphere CH4 flux; N2Of was soil-atmosphere N2O flux; and ∆SOC was SOC change over a 15-yr period.

Factors contributing to GWP were converted to CO2 by multiplying by 3.67 (44 g CO2 mol-1 ÷ 12 g C mol-1). Carbon dioxide equivalents for soil-atmosphere CH4 and N2O flux were calculated assuming direct GWP of 1 kg CH4 ha-1 = 28 kg CO2 ha-1 and 1 kg N2O ha-1 = 265 kg CO2 ha-1 (100-yr time horizon) (Myhre et al. 2013). Results for net GWP were expressed as kg CO2equiv. ha-1 yr-1. Values were related to crop and animal productivity by dividing net GWP by total yield, protein, and energy. Total yield was taken as the sum of corn, spring wheat, and soybean grain yield for both treatments, along with animal weight gain for the Grazed treatment. Grain protein was calculated using measured values (outlined above), while meat protein was assumed 19.8% for trimmed retail cuts of raw beef (USDA 2020b). Energy values for corn, spring wheat, soybean, and beef were assumed 365, 339, 446, and 224 kcal 100 g-1, respectively (USDA 2020b). Calculated quotients for total yield, protein, and energy reflected the GHGI for each treatment.

Data analyses

Values for grain yield, soil-atmosphere CH4 and N2O flux, and SOC stocks under Grazed and Ungrazed treatments were analyzed within and across years using appropriate models in SAS (SAS Institute 2014). Factors contributing to GWP and GHGI were analyzed using PROC GLIMMIX in SAS with treatments and replicates as fixed and random effects, respectively (Schabenberger 2005). The slice option in the LSMEANS statement was used to generate P-values for comparisons, and differences among means were documented at P ≤ 0.05. Variation of means was expressed using standard error (Steel and Torrie 1980).

Results

Weather conditions

Monthly precipitation and temperature over the course of the study was typically wetter and cooler than long-term means (Fig. S.2). Near-continuous snowfall in the winter months was coupled with lower than normal temperatures Jan-March, contributing to wet and cool soil conditions following spring snowmelt (Liebig et al. 2020a). Monthly precipitation during the growing season (April-Sep) was typically less than the long-term (107 yr) mean in April-May, only to rebound with greater than normal precipitation in later summer months. Precipitation was deficient only in 2017 during the first four months of the growing season, with 122 mm less precipitation April-July compared to long-term monthly means.

Production metrics

Differences in grain yield within and across crops highlighted production tradeoffs associated with the ICL system used in this study (Table 1). The Grazed treatment had one less grain harvest compared to the Ungrazed treatment over the three-year rotation, as the third year of the Grazed treatment was allocated to perennial forage. Additionally, demand for available soil water by crops was likely higher in the Grazed treatment compared to the Ungrazed treatment due to the inclusion of intercropped forages in the former. Moreover, corn was grazed and not harvested for grain in 2017 due to negligible production caused by deficient moisture conditions early in the growing season (April-July), contributing to less grain yield overall for the Grazed treatment. Accordingly, grain yield was 1.5–3.0 times greater across crops in the Ungrazed treatment compared to Grazed treatment (Table 1). Within crops, corn yield was 1.4–3.2 times greater in the Ungrazed than Grazed treatment, with significant differences between treatments observed for each year and averaged across years. Spring wheat grain yield was greater under the Ungrazed treatment compared to the Grazed treatment in 2018, though differences were not observed in other years nor averaged across years.

Table 1 Grain yields for grazed and ungrazed cropland in an integrated crop-livestock study near Mandan, ND USA, 2017-2019
Full size table

Like grain yield, sum of protein and energy in grain across crops was significantly greater under the Ungrazed treatment compared to the Grazed treatment each year and averaged across years (Table 1). Within crops, protein and energy in corn grain was significantly greater in the Ungrazed treatment compared to the Grazed treatment in 2018 and 2019. Differences in spring wheat grain protein between treatments was limited to 2019, where protein was greater under Grazed cropland than Ungrazed cropland (455 vs. 379 kg ha-1; P = 0.0279). In 2018, energy in spring wheat grain was greater under the Ungrazed treatment compared to the Grazed treatment (13393 vs. 11190 Mcal ha-1; P = 0.0324). Mean protein and energy values across years for corn and spring wheat were not different between treatments.

Daily cattle weight gain in the Grazed treatment differed across years, with values greatest in 2017 (1.1 ± 0.1 kg animal d-1), intermediate in 2018 (0.8 ± 0.1 kg animal d-1), and least in 2016 (0.2 ± 0.1 kg animal d-1) (P < 0.0001) (Table 2). The extended time spent grazing standing corn in 2017 under favorable weather conditions contributed to greater total weight gain compared to 2016 and 2018. When expressed on an areal basis, total animal weight gain was 10.4, 170.4, and 31.9 kg ha-1 in 2016, 2017, and 2018, respectively (P < 0.0001). Total protein and calories in meat followed similar trends as animal weight gain across years (Table 2).

Table 2 Cattle allotments, grazing schedule, weight gain, and accumulated protein and energy in meat from beef cattle for an integrated crop-liveastock study near Mandan, ND USA, 2016–2018
Full size table

Factors contributing to global warming potential

Input production and field operations

Carbon dioxide emissions associated with seed and fertilizer production were similar across years for each crop phase in Grazed and Ungrazed treatments (Table 3). The only substantial deviation in CO2-associated emissions from seed and fertilizer production across years happened in the cover crop phase in 2019, where the planting of proso millet was necessary due to the failure of the perennial cover crops to establish the previous year. Herbicide-associated CO2 emissions were variable across years within crop phases due to the adaptive nature of weed control under no-till management. Carbon dioxide emissions from herbicide production in the Grazed treatment ranged from 70 to 110, 40 to 62, and 0 to 70 kg CO2equiv. ha-1 yr-1 for corn-soybean, spring wheat + cover crop, and cover crop phases, respectively. For the Ungrazed treatment, CO2 emissions from herbicide production ranged from 70 to 74, 33 to 46, and 62 to 75 kg CO2equiv. ha-1 yr-1 for corn, spring wheat, and soybean phases.

Table 3 Carbon dioxide emissions associated with production inputs and field operations for grazed and ungrazed cropland in an integrated crop-livestock study near Mandan, ND USA, 2017–2019
Full size table

Field operations were generally more frequent and variable across years in the Grazed treatment compared to the Ungrazed treatment due to additional planting and multiple haying operations in the former. Carbon dioxide emissions from field operations in the Grazed treatment ranged from 69 to 105, 119 to 169, and 44 to 87 kg CO2equiv. ha-1 yr-1 for corn-soybean, spring wheat + cover crop, and cover crop phases, respectively. Field operation emissions in the Ungrazed treatment were mostly consistent across years, with means or ranges of 80, 94–104, and 71 kg CO2equiv. ha-1 yr-1 for corn, spring wheat, and soybean phases.

When averaged across years, CO2 emissions associated with seed, fertilizer, and herbicide production in the Grazed treatment were 107, 139, and 53 kg CO2equiv. ha-1 yr-1, respectively (Table 5). In the Ungrazed treatment, the same were 67, 139, and 61 kg CO2equiv. ha-1 yr-1. Mean CO2 emissions associated with field operations in the Grazed and Ungrazed treatments were 98 and 83 kg CO2equiv. ha-1 yr-1, respectively.

Methane emission from enteric fermentation

Estimated CH4 emission by enteric fermentation from beef cattle in the Grazed treatment was greater in 2017 compared to 2016 and 2018 (Table 4), corresponding to the extended grazing period in 2017. When expressed on an areal basis, CH4 emission by enteric fermentation was 2.7 ± 0.1, 40.9 ± 3.3, and 4.7 ± 0.1 kg CH4-C ha-1 in 2016, 2017, and 2018, respectively. When averaged across years and expressed on a CO2–equivalent basis, enteric fermentation accounted for 602 kg CO2equiv. ha-1 yr-1 in the Grazed treatment (Table 5).

Table 4 Estimated CH4-C production from beef cattle for an integrated crop-livestock study near Mandan, ND USA, 2016-2018
Full size table
Table 5 Effects of grazed and ungrazed cropland on greenhouse gas sources and sinks, net global warming potential (GWP), and greenhouse gas intensity (GHGI) per unit product, protein, and energy for an integrated crop-livestock study near Mandan, ND USA
Full size table

CH4 and N2O flux

Both treatments exhibited negative soil-atmosphere CH4 flux over the study (Fig. 2A), with mean annual values of – 1.10 ± 0.08 and – 0.87 ± 0.07 kg CH4-C ha-1 yr-1 for Grazed and Ungrazed treatments, respectively. Repeated measures analysis of hourly soil-atmosphere CH4 flux found greater CH4 uptake rate during the first two years of the study compared to the third (– 15.4 vs. – 7.8 µg CH4-C m-2 h-1; P ≤ 0.0001), and greater CH4 uptake during summer (Jun-Aug) compared to other seasons (– 18.8 vs. – 9.3 µg CH4-C m-2 h-1; P ≤ 0.0001) (Table S.2). No differences in hourly soil-atmosphere CH4 flux were found between Grazed and Ungrazed treatments (– 14.4 vs. – 11.9 µg CH4-C m-2 h-1; P = 0.1837), nor among crop phases (– 14.0, – 12.1, and – 12.7 µg CH4-C m-2 h-1 for corn-soybean/corn, spring wheat + cover crop/spring wheat, and cover crop/soybean phases, respectively; P = 0.1819). When expressed on a CO2-equivalent basis, annual uptake of CH4 by soil was – 41 and – 32 kg CO2equiv. ha-1 yr-1 in the Grazed and Ungrazed treatments, respectively (Table 5).

Fig. 2
figure 2

Annual soil-atmosphere CH4 and N2O flux for grazed and ungrazed cropland in an integrated crop-livestock study near Mandan, ND USA, 2016–2019. Error bars reflect standard error of the mean. Annual values are reflective of dates beginning and ending in mid-October

Full size image

Treatment effects on soil-atmosphere N2O flux were previously reviewed (Liebig et al., 2020a). Briefly, mean values of annual N2O flux were not different between treatments over the evaluation period (P = 0.9039), with fluxes of 2.6 and 2.5 kg N2–N ha-1 yr-1 for Grazed and Ungrazed treatments, respectively (Fig. 2B). Among crop phases of the Grazed treatment, annual N2O flux was greatest under corn-soybean (3.5 kg N2O-N ha-1 yr-1), intermediate in spring wheat + cover crop (2.7 kg N2O-N ha-1 yr-1), and least under cover crop (1.5 kg N2O–N ha-1 yr-1). For the Ungrazed treatment, corn and wheat phases exhibited elevated annual N2O flux (both 2.8 kg N2O–N ha-1 yr-1) relative to soybean (1.9 kg N2O–N ha-1 yr-1). Major flux events during spring thaw (March-April) had a disproportionate effect on N2O flux in both treatments, accounting for 34–54% of cumulative flux across crop phases over the course of the study. When averaged across years and expressed on a CO2–equivalent basis, soil-atmosphere N2O flux equated to 1062 and 1046 kg CO2equiv. ha-1 yr-1 for Grazed and Ungrazed treatments, respectively (Table 5).

Soil organic carbon

Profile SOC stocks did not differ between treatments in 1999 and 2014 (Fig. 3). Prior to initiating the study in 1999, SOC stocks were 179.8 ± 3.2 and 183.1 ± 8.3 Mg C ha-1 in the Grazed and Ungrazed treatments, respectively. In 2014, SOC stocks increased to 184.1 ± 2.7 Mg C ha-1 in the Grazed treatment and 187.9 ± 7.3 Mg C ha-1 in the Ungrazed treatment, with changes most pronounced in the 0–0.08 m depth (data not shown). At this surface depth, increases in SOC concentration were complemented by large increases in soil bulk density over time, an outcome presumed to be caused by reduced root biomass and increased vehicular traffic between 1999 and 2014 (Liebig et al., 2020b). The outcome resulted in an accentuated SOC stock increase at 0-0.08 m, but when considered over the entire soil profile, change in SOC stocks over time was modest for both treatments (4.4 ± 3.4 Mg C ha-1 for Grazed, 4.8 ± 6.9 Mg C ha-1 for Ungrazed; P = 0.9742). Associated annual changes in SOC stocks were 0.29 and 0.32 Mg C ha-1 yr-1 for Grazed and Ungrazed treatments, respectively. Accordingly, Grazed and Ungrazed treatments were net CO2 sinks with uptake capacity of – 1074 and – 1162 kg CO2equiv. ha-1 yr-1 (Table 5).

Fig. 3
figure 3

Soil organic C stocks for grazed and ungrazed cropland in 1999 and 2014 for an integrated crop-livestock study near Mandan, ND USA. Error bars reflect standard error of the mean. Soil carbon stock change based on a soil mass of 15123 Mg ha-1, approximating the 0–0.91 m depth

Full size image

Net global warming potential and greenhouse gas intensity

For the Grazed treatment, the relative impact of summed factors contributing to net GWP decreased in order of (1) SOC stock change, (2) soil-atmosphere N2O flux, (3) enteric fermentation, (4) fertilizer production, (5) seed production, (6) field operations, (7) herbicide production, and (8) soil-atmosphere CH4 flux (Table 5). Summed factors contributing to net GWP for the Ungrazed treatment decreased in relative impact of (1) SOC stock change, (2) soil-atmosphere N2O flux, (3) fertilizer production, (4) field operations, (5) seed production, (6) herbicide production, and (7) soil-atmosphere CH4 flux. Factors differing in their CO2-equivalent contribution to net GWP between treatments included seed production (Grazed > Ungrazed), field operations (Grazed > Ungrazed), and soil-atmosphere CH4 flux (Ungrazed > Grazed). Additionally, enteric fermentation in the Grazed treatment contributed 602 kg CO2equiv. ha-1 yr-1 absent in the Ungrazed treatment. Soil-atmosphere N2O flux negated nearly all the CO2equiv. sink capacity accrued from SOC stock change in the Grazed treatment (1062 vs. – 1074 kg CO2equiv. ha-1 yr-1). A similar outcome with N2O flux and SOC stock change was observed in the Ungrazed treatment (1046 vs. – 1162 kg CO2equiv. ha-1 yr-1).

Summing across factors and years, net GWP was positive for Grazed and Ungrazed cropland (946 and 200 kg CO2equiv. ha-1 yr-1, respectively; P = 0.0331), implying net CO2equiv. emission to the atmosphere by both treatments (Table 5). Measures of mean GHGI were significantly greater under the Grazed treatment compared to the Ungrazed treatment for grain and beef (142 vs. 0 g CO2equiv. kg product-1; P = 0.0141), protein (1113 vs. 21 g CO2equiv. kg protein-1; P = 0.0093), and energy (26 vs. 0 g CO2equiv. Mcal-1; P = 0.0164).

Discussion

Treatments evaluated in this study allowed for direct assessment of GHG tradeoffs associated with including livestock on cropland in a semiarid region. As typical of ICL systems, the specific management practices evaluated in this study were tailored to site conditions and available resources (Liebig et al. 2017a). However, practices were generally representative of ICL systems in the northern Great Plains, in that cropland was used for mixed production of forage and grain, beef cattle were used as livestock, and the system design was inherently flexible so as to allow for adaptive responses to variation in weather and markets (Hendrickson et al. 2008; Hanson et al. 2007).

System Effects on Net GWP and GHGI

The positive net GWP for grazed cropland corresponded to a similar outcome for ICL systems in subtropical conditions (86–4919 kg CO2equiv. ha-1 yr-1; Ribeiro et al. 2019). Outcomes from aligned ICL research in a wetter portion of the northern Great Plains (eastern South Dakota) suggest the likelihood of positive net GWP from grazed cropland given appreciable positive fluxes of both N2O and CH4 (Singh et al. 2020). The positive, though significantly lower net GWP for ungrazed cropland corresponded with other studies of semiarid cropping systems (116–1948 kg CO2equiv. ha-1 yr-1; Liebig et al. 2019; Tellez-Rio et al. 2017), though a study on a similar soil type in northeast Montana found rainfed cropping systems to have negative net GWP (– 292 to – 2107 kg CO2equiv. ha-1 yr-1; Sainju et al. 2014).

In this study, the difference in net GWP between grazed and ungrazed cropland was largely driven by emissions from enteric fermentation by beef cattle. This single factor accounted for 80% of the difference in net GWP between treatments and served to significantly increase GHGI attributes under grazed cropland. Broadly, adaptive management attributes of the Grazed treatment highlighted how net GWP could vary widely by year depending on forage availability and quality, as well as weather conditions. For instance, the 2016 grazing period was characterized by low forage production, low forage quality, and a short grazing period, resulting in a net GWP of – 30 kg CO2equiv. ha-1. In 2018, elevated forage production and quality contributed to improved cattle growth, with a resulting positive net GWP (770 kg CO2equiv. ha-1). In 2017, the decision to graze standing corn coupled with mild weather conditions – thereby extending the grazing period – intensified grazing effects on net GWP (2097 kg CO2equiv. ha-1). In contrast, annual variation in net GWP was less in the Ungrazed treatment (e.g., – 595, 582, and 614 kg CO2equiv. ha-1 yr-1 in 2017, 2018, and 2019, respectively). The large fluctuation in net GWP under the Grazed treatment underscores the sensitivity of GHG outcomes to grazing management in ICL systems. More broadly, it is important to note that net GWP accounting of beef cattle on cropland in this study does not reflect GHG outcomes from other agroecosystems used by cattle at other times of the year. Previous estimates of net GWP on nearby native rangeland were found to be negative at moderate and heavy stocking rates (– 783 and – 618 kg CO2equiv. ha-1 yr-1, respectively; Liebig et al. 2010). Accordingly, typical ICL management in the region, as reflected by grazing on native rangeland during the growing season coupled with grazing on cropland later in the year, would serve to reduce overall net GWP.

Methane Emission from Cattle

Daily CH4 production estimates for enteric fermentation across years were comparable to measured assessments in a three-year study of ICL systems in a subtropical environment (151 vs. 122 g CH4-C animal d-1, respectively) (Pontes et al. 2018) and literature-derived values for ‘good pasture’ in the U.S. (134 g CH4-C animal d-1) (Westberg et al. 2001), underscoring the validity of the approach used in this study. Resulting estimates elucidated the strong influence of grazing duration and diet quality on enteric fermentation emissions from beef cattle in ICL systems. The relatively short grazing periods in 2016 (30 d) and 2018 (26 d) yielded modest enteric fermentation emissions, with total CH4 production for both periods  < 5 kg CH4-C ha-1. In contrast, enteric fermentation emissions in 2017 were over eight-fold greater than other years, an outcome caused by a lengthy grazing period (118 d) due to limited snowfall and moderate temperatures coupled with access to standing corn. Daily CH4 production estimates in 2017 were over 2 and 4 times greater than 2018 and 2016, respectively, highlighting the tradeoff of including standing corn as a supplement to forage when weather conditions are favorable. However, when daily CH4 production estimates were normalized to mean daily weight gain, values in 2017 (241 g CH4-C kg-1) were intermediate to 2016 and 2018 (321 and 153 g CH4-C kg-1, respectively) due to significantly higher weight gain in 2017, underscoring the importance of framing results in the context of production. Though there are notable caveats to comparing studies with differences in stocking rate and grazing duration, enteric fermentation emissions from this study were less than those observed in an ICL study evaluating grazing intensities of winter black oat (Avena strigosa Schreb.) (952–2382 kg CO2equiv. ha-1 yr-1; Ribeiro et al. 2019), but greater than nearby historically native pasture and seeded perennial forage (176–563 kg CO2equiv. ha-1 yr-1; Liebig et al. 2010).

Inputs and field operations

Differences in individual factors contributing to GHG emissions highlighted how inclusion of forages and livestock under adaptable management can lead to greater emissions in ICL systems. Emissions associated with seed production were 49% greater under Grazed cropland compared to Ungrazed cropland. Here, inclusion of intercropped soybean in the corn-soybean phase and the cover crop mixture in the spring wheat + cover crop phase of the Grazed treatment added additional seed-associated emissions each year. Also, the reseeding of forage during the third year of the study of the cover crop phase (proso millet) created an emission burden absent in the Ungrazed treatment. While additional emissions from seed production seem unavoidable in ICL systems due to the need for nutritious forage throughout the crop rotation, it is possible to reduce the contribution of this factor by limiting selected forages to those that optimize production and quality for GHG mitigation purposes, or by extending the period of the perennial phase during the crop rotation cycle (Hammond et al. 2014).

Differences in field operation emissions between treatments largely reflected additional activities in the Grazed treatment associated with soybean interseeding and forage production. The number of tractor and combine passes across all crop phases in the Grazed treatment ranged from 13-18 passes yr-1, whereas 11–12 passes yr-1 were needed in the Ungrazed treatment. The seeding of soybean and cover crops in the Grazed treatment added an additional tractor pass in the corn-soybean and spring wheat + cover crop phases, respectively. Moreover, three to six tractor passes were required the subsequent year in the cover crop phase, as forage harvest was inclusive of mowing, swathing, and baling. As such, opportunities exist to reduce field operation emissions in ICL systems through direct grazing of standing forage. This was done in 2017 for corn that was grazed as feed, thereby eliminating one combine pass typically conducted for grain harvest. While grazing cover crops was a possibility to reduce field operations in this study, using the available forage as hay was employed instead as forage resources for cattle on nearby grazing land were abundant and of good quality during summer. This decision allowed cover crop forage to be used later in the year when the quality of forage on grazing land typically declines.

Managing GHG Sources

Strategic management of inputs and livestock offer opportunities to reduce GHG sources from grazed cropland. As reviewed above, direct grazing of a limited diversity of high-quality forages can mitigate GHG emissions through fewer field operations and lower seed inputs. Additionally, more efficient utilization of available N can have multiple mitigation benefits, inclusive of fewer field operations, lower fertilizer inputs, and lower soil-atmosphere N2O flux, though the efficacy of the latter requires efficient N cycling and maintenance of aerobic soil conditions (Eckard et al. 2003). Moreover, use of nitrification inhibitors to reduce N2O flux from ammonia-based fertilizers and urine patches shows promise for improving N use efficiency (Di and Cameron 2002) as does certain forage species with biological nitrification inhibition capacity (Simon et al. 2020), but evaluations in semiarid cropping systems are lacking. Finally, numerous opportunities exist to reduce enteric fermentation emissions from ruminant livestock in grazed cropping systems, though some options have proven more effective and economically viable than others (e.g., culling unproductive animals and improving diet quality vs. supplementing forage with plant secondary compound or manipulating rumen microbial populations) (Eckard et al. 2010). Substituting ruminant for non-ruminant livestock can also result in lower enteric fermentation emissions, though tradeoff emissions elsewhere in the production cycle must be considered (Grossi et al. 2019).

Enhancing GHG sinks

Enhancing the GHG sink capacity of both grazed and ungrazed cropping systems requires management to minimize C losses with concurrent increases in soil C inputs. To this end, management practices that limit soil disturbance are critical to mitigating C losses to mineralization and soil erosion (Lal 2010). As treatments in this study were managed using no-till practices, soil C losses through these pathways were minimized. Opportunities to increase SOC stocks in semiarid cropping systems include incorporation of perennial grasses and woody plant species (Ledo et al. 2020). Abundant evidence supports the potential for large increases in SOC stocks under highly productive forage and/or biofeedstock grass species (Franzluebbers et al. 2012; Schmer et al. 2011), while there is complementary evidence for SOC stock increases when transitioning from cropland or pasture to silvopastural systems (De Stefano and Jacobson 2017), with potential for soil carbon storage in deeper soil depths (> 0.75 m) (Howlett et al. 2011). Moreover, incorporation of woody plant species can generate additional marketable products (e.g., nuts, fruit, timber), and create microenvironments that increase crop production while mitigating animal stress during weather extremes (Krishnamurthy et al. 2019; Vallejo et al. 2012). As such, expanding the number and diversity of enterprises in semiarid cropping systems may provide significant contributions to economic and environmental sustainability over the long-term, and could serve as an effective strategy to sustainably intensify agricultural production. However, integrating GHG sources and sinks across enterprises would be incumbent to a comprehensive assessment of net GWP and GHGI.

Conclusions

Livestock integration on cropland can be a viable strategy to sustainably intensify agriculture. Outcomes reported here, however, suggest otherwise. Within the scope of GHG-related metrics, inclusion of beef cattle on semiarid cropland resulted in a significantly greater GHG burden compared to ungrazed cropland. Future research should focus on the design of ICL systems with net CO2-equivalent sink capacity. Designs that improve nutrient-use efficiency and increase soil C accrual are needed. As ICL systems are uniquely characterized by the dynamic allocation of resources in space and time, effective designs should accommodate an adaptive management framework.

In addition to novel designs, future evaluations of ICL production systems must be inclusive of an assessment portfolio spanning production, environment, and human dimension domains. Only through inclusive evaluations can land managers adequately assess the many tradeoffs within and between domains when considering ICL systems as a sustainable intensification strategy.