1 Introduction

Specialization of commercial agroecosystems and the corresponding economies of scale have led to impressive gains in productivity and labor-use efficiency in recent decades. However, the long-term sustainability of highly specialized systems and concentrated agricultural landscapes is increasingly under scrutiny. Dependence on a small number of agricultural commodities (Puma et al. 2015), unsustainable mining of water and soil resources (Foley et al. 2011), loss of agrobiodiversity (Tilman et al. 2006), and volatility in commodity markets (Russelle et al. 2007) make these systems more vulnerable to climate unpredictability and changing resource availability. Concern that these vulnerabilities could endanger critical ecosystem services to and from agriculture entails the need for production systems that reliably maintain or boost productivity without further environmental externalities.

Integration of crops and livestock in commercial production systems intensifies production by leveraging synergies between multiple agroecosystem components rather than relying on increased input levels. However, integrating commercial crop and animal production involves tradeoffs at several scales. At the farm and regional scales, socio-economic motivators such as risk reduction and land-use efficiency on the one hand are balanced by high upfront capital costs and managerial complexity on the other (Garrett et al. 2017). While these tradeoffs are relatively well understood, much less is known about the biophysical tradeoffs at the field scale, i.e., between ecosystem functions driving animal production and those driving crop production. From a management perspective, characterizing these tradeoffs is important for developing best management practices and informing adaptive management strategies under disparate environmental conditions. These potential tradeoffs also have implications for the long-term sustainability and resilience of integrated systems, i.e., their ability to continue providing desired ecosystem services under adverse conditions (Folke 2006).

Rotating livestock into productive cropland can have effects on agroecosystem properties, including soil physical, chemical, and biological attributes, that carry over into the subsequent cropping season (Carvalho et al. 2018). For example, inter-seasonal interactions between livestock grazing, grass physiology, and soil properties can affect soil water availability and consequently the growth and development of crops in rotation. Reductions in stored soil water in grazed areas due to increased evapotranspiration rates from grazed winter forages could negatively impact crop resilience if water conditions are limiting in the following season. On the other hand, improved soil structure and thus capacity for root expansion and exploration in grazed areas could have the opposite effect, bolstering subsequent crop resilience to water stress (Cecagno et al. 2017).

This study pairs a detailed within-season analysis of soil-plant-water dynamics with over 16 years of results from a long-term, no-till integrated beef/soybean system in Rio Grande do Sul, Brazil (Fig. 1). The dynamics of plant-soil interactions remain under-researched in this diversified system, but evolution in soil physical, chemical, and biological properties due to crop-livestock integration at the site has been well documented. With regard to soil physical properties, moderate grazing of winter forage cover crops has been shown to improve soil aggregate structure while leaving bulk density and porosity unchanged compared with ungrazed cover forages (Conte et al. 2011). With regard to soil chemical properties, soil total organic carbon and nitrogen and particulate organic carbon and nitrogen all have been shown to increase in the moderately grazed vs. ungrazed system (Assmann et al. 2014). Furthermore, acidity declines more after liming in moderately grazed systems, while basic cation use efficiency increases (Martins et al. 2014b). Evidence of soil biological responses to grazing in this system has been seen in increased microbial functional diversity and decreased microbial activity in grazed compared with ungrazed plots (Chávez et al. 2011), although seasonal dynamics in microbial biomass C, N, and P are similar between the two systems (Souza et al. 2010a). The mechanisms of these changes range from altered nutrient recycling from dung and urine patches, to root growth dynamics during the winter pasture season, to residue quantity and quality, among others. Despite these differences in soil properties and substantial year-to-year variation in soybean yields, differences in soybean yields between the grazed and ungrazed systems have yet to be detected except under higher-than-recommended grazing intensities (Carvalho et al. 2018).

Fig. 1
figure 1

a Beef steer in annual black oat (A. strigosa) and Italian ryegrass (L. multiflorum) pasture under moderate grazing intensity (20-cm sward height). b Pre-harvest soybean at the end summer growing season (March–April), showing volunteer regrowth Italian ryegrass in the understory. Photos: P. Nunes

Because winter grazing treatments were expected to differentially impact starting soil water conditions, as well as soil physical and chemical parameters influencing plant available water and uptake in the soybean season, this study monitored soil water content across the growing season using several complementary techniques. We paired soil water measurements with indicators of plant water status and physiology, hypothesizing that, even in favorable climate years, the effects of livestock integration on soil water storage and nutrient cycling carry over into the subsequent cropping season and potentially affect crop yield–determining factors. Our overall objective was to describe agroecosystem processes that may ultimately influence resilience to variable and unpredictable climates and resource availabilities. Specifically, we intended to characterize these soil-plant dynamics during a normal-to-above-average precipitation year, both to establish a baseline of system behavior against which responses in less favorable years could be compared and to fill a gap in the understanding of this system. Most studies on crop-livestock integration worldwide only draw on a few years of data, but longer-term approaches are needed to describe system trajectories and potential tradeoffs from livestock integration. This study is unique in its ability to provide insight into the potential mechanisms driving agroecosystem functions in this widely implemented integrated crop-livestock system in southern Brazil after 16 years of winter grazing. Our results are important both for the purposes of adaptive management planning and for long-term sustainability assessments in diversified production systems.

2 Materials and methods

2.1 Experimental site

The experimental site is located at Espinilho Farm in São Miguel das Missões district, Rio Grande do Sul, Brazil (28°56′ 14″ S 54° 20′ 52″ W, 465-m elevation). The Planalto physiographic region has a warm, humid subtropical climate with an average annual temperature of 19 °C and average annual precipitation of 1850 mm, distributed unevenly and primarily in the summer months. The soil is an Oxisol (Rhodic Hapludox): deep, well drained, and clayey in texture, with 540, 270, and 190 g of clay, silt, and sand, respectively, per kilogram of soil. The site’s soils are characterized by low pH levels (3.7–4.7), moderate aluminum and base saturation, and moderate to high potassium and phosphorus availability.

Since 2001, the site has been managed as an integrated annual beef-soybean rotation under no-till. The experimental design is a randomized complete block with 3 replicates for the grazing treatments and 2 replicates for the ungrazed control. For this study, one control plot was divided in two to complete the replicate design, despite some degree of non-randomization introduced. Permanent plots are subjected to 4 grazing intensity treatments in the winter. We sampled only the moderate grazing intensity at 20–30 cm sward height in addition to the controls. This moderate grazing treatment represents the best management practice for this site as indicated by soil health and animal performance (Kunrath et al. 2014). Plot areas ranged from 0.05 ha for ungrazed control plots to 2.1 ha for grazed plots. The varying areas facilitate sward height management using the put-and-take method (Mott and Lucas 1952), in which animals are added or removed from plots to maintain target sward heights. However, the smaller ungrazed plots also result in an increased likelihood of oversampling or partial randomization, and consequently non-independence of measurements. These limitations are related to a continuous, long-term experimental protocol originally designed 16 years ago as a grazing trial; making modifications to the design at this stage would have introduced unwanted inconsistencies into the extended dataset.

Each winter (May–October), experimental plots are seeded with a mixture of black oat (Avena strigosa Schreb.) and Italian ryegrass (Lolium multiflorum Lam.) over the naturally reseeded ryegrass from the previous season. The grazed pastures are set to beef steers under continuous, put-and-take stocking for an average of 120 days year−1 (minimum 99 days, maximum 150 days) after the forage has accumulated 1500 kg ha−1 of dry matter (usually between June and July). Put-and-take stocking results in variable stocking rates across seasons, typically ranging from 750 to 1100 kg live weight ha−1. A total of 250 kg ha−1 urea is applied in two doses during the winter pasture season, the first when the pasture reaches stage V3-V4 and the second just before animals enter the plots, approximately 1 month later. KCl (160 kg ha−1) and triple superphosphate (130 kg ha−1) are applied at pasture sowing. In the summer between November–December, the pastures are desiccated with 1.5 kg ha−1 Roundup WG and direct-seeded with soybean. Soybean is typically sown between November 15 and December 15 and harvested after approximately 120 days. In 2016, inoculated seed (Brasmax Tornado 6863 RR) was sown on November 13–14, 2016, and harvested 139 days later, on April 1–2, 2017. Fungicides and insecticides were applied as needed throughout the season, along with 1 L ha−1 Super K400 foliar fertilizer applied together with the first fungicide dose.

Measurements were collected during three periods of the summer 2016–2017 cropping season (November to April), representing early (V2-R1), late-mid (R1-R5), and late (R5-R7) soybean growth stages. Each grazed plot was paired with an ungrazed control plot for a total of three experimental blocks. Only a limited number of sensor structures could be transported and assembled at the site. All measurements were taken from the same experimental block where the sensor structures were positioned for 5 consecutive days to ensure adequately representative sampling and to capture a range of weather conditions for each block while avoiding confounding interactions between block and weather conditions. Sensors would then be moved to the next block for a further 5 days, for a total of 15 days of sampling per period. Soil samples were collected twice during each 5-day window. Eight and four samples were collected in grazed and ungrazed plots, respectively, for 2 soil layers (0–10 cm and 10–20 cm), according to a crossed transect sampling scheme. The unbalanced sampling design was adopted to account for differences in area between grazed and ungrazed plots.

The site received a total of 672 mm of precipitation during the summer growing season with a relatively uneven distribution; 38% of precipitation occurred during the month of January alone (257 mm). Total monthly precipitation was about 69% greater in January than the 30-year average for the region. Mean daily temperature was 23.8 °C for the season, about 14% higher than the 30-year average.

2.2 Soil moisture and hydraulic properties

Soil samples were dried to constant weight at 70 °C and weighed to determine gravimetric water content. Volumetric water content was calculated by correcting for bulk densities previously determined for grazed and ungrazed plots in the 5–10 and 10–20 cm layers (1.30, 1.28, 1.33, and 1.36 g cm−3 for grazed 5–10 and 10–20 cm and ungrazed 5–10 and 10–20 cm, respectively; Cecagno et al. 2017).

Granular matrix electrical resistance sensors (WATERMARK model 200SS, Irrometer Company, Riverside, CA) for measuring soil matric potential were installed at the beginning of the soybean season and remained in place until harvest. Each of 2 randomly selected plots per treatment contained 2 sensors at 20-cm depth and 2 sensors at 40-cm depth (a total of 4 sensors per treatment per depth) that recorded soil water content at 10-min intervals.

Water retention curves for each treatment were determined using the simplified evaporation method (Peters and Durner 2008) with a HYPROP® precision mini-tensiometer setup (UMS GmbH, Munich, Germany). Undisturbed soil samples were collected in January 2017 from the 0–10 cm soil layer using 250-mL PVC rings, stored in plastic wrap, and processed within approximately 10 weeks of collection. Samples were saturated and then subjected to a week-long dry-down period at room temperature (28 °C) with a fan due to the high clay content of the soil. The HYPROP® measured water potential and sample weight continuously during the dry-down period, and pressure head in the tensiometer base was calculated using HYPROP-FIT® v.3.5.1.13951. Water retention curves were fit using the modified van Genuchten model for soil hydraulic properties given in Eq. 1 below:

$$ \theta ={\theta}_{\mathrm{r}}+\frac{\theta_{\mathrm{s}}-{\theta}_{\mathrm{r}}}{{\left[1+{\left|\alpha h\right|}^n\right]}^m} $$
(1)

where θ is the volumetric water content (mm3 mm−3), θr is the residual volumetric water content, θs is the saturated volumetric water content, h is the soil water potential (kPa), α is a scaling parameter inversely proportional to the mean pore diameter (cm−1), and n and m are shaping parameters for soil water characteristics: m = 1–1/n, 0 < m < 1.

2.3 Plant measurements

Canopy development and light-use dynamics were captured at 10-min intervals using normalized difference vegetation index (NDVI) and photochemical reflectance index (PRI) sensors (SRS-Nr and SRS-Pr, METER Group, Inc., Pullman, WA) with a half-angle field of view of 18°. Sensors were mounted 10 ft above ground level and 45° from surface normal resulting in an elliptical visible canopy area of approximately 10 m2. NDVI and PRI sensors were coupled with upward-facing hemispherical sensors to record downwelling radiation at 630 nm and 800 nm for NDVI sensors and 570 nm and 531 nm for PRI sensors. Indices were calculated as follows:

$$ NDVI=\frac{\alpha {R}_{800}-{R}_{630}}{\alpha {R}_{800}+{R}_{630}} $$
(2)
$$ PRI=\frac{\alpha {R}_{531}-{R}_{570}}{\alpha {R}_{531}+{R}_{570}} $$
(3)

where α is the ratio of incident red (R800) to incident near-infrared light (R630) for NDVI and the ratio of incident 570 nm (R570) to incident 531 nm light (R531) for PRI.

Leaf water potential was determined using the Scholander pressure chamber method at pre-dawn (4:30 am–6:00 am UTC− 3:00) and at solar noon (11:00 am–1:00 pm UTC− 3:00). For both pre-dawn and afternoon measurements, the newest, fully expanded leaf was harvested on 16 plants in a grazed plot and 8 plants in an ungrazed plot by cutting at the base of the petiole, immediately placing in a humidified plastic container, and measuring within 20–60 s. Stomatal conductance was recorded simultaneously with the midday water potential measurements on similarly positioned leaves of an additional 16 (grazed) and 8 (ungrazed) plants using a handheld porometer (SC-1, METER Group, Inc., Pullman, WA).

Potential evapotranspiration (ETo) was calculated by the Penman-Monteith FAO-56 method using data from an on-site meteorological station (VantagePro2™, Davis Instruments, Vernon Hills, IL). Missing data points for temperature were interpolated from linear regression with the nearest automatic meteorological station in Cruz Alta, RS, Brazil, approximately 78 km from the experimental site (28° 36′ 12.38″ S, 53° 40′ 24.95″ W, 427 m ASL).

Soybean yield was determined following standard protocol at full maturity, 139 days after sowing. All plants along five randomly selected 2-m transects (0.9 m2 area per transect, 4.5 m2 per plot) were harvested by hand, threshed, and cleaned. After weighing, grain yields were adjusted for 13% moisture content.

2.4 Statistical analysis

All analyses were performed in R statistical computing software v.3.4.1. A linear multi-level model approach was used to analyze volumetric water content, stomatal conductance, leaf water potential, mean midday maxima for NDVI and PRI, and crop yield. The volumetric water content model included sampling date as a repeated measure, grazing treatment, sampling depth, and their interaction as fixed effects, and the three-way interaction among treatment, depth, and sampling date as a random effect. The stomatal conductance model included soybean growth stage as a repeated measure; ETo as a linear covariable; growth stage, grazing treatment, and their interaction as fixed effects; and block-treatment interaction and growth stage–block interaction as random effects. Leaf water potential was represented by the same model with the addition of time of day (morning or afternoon) as a fixed effect. The crop yield model included treatment as a fixed effect and block-treatment interaction as a random effect.

For NDVI, mean midday maxima (measurements from between the hours of 11:00 am UTC− 3:00 and 1:00 pm UTC− 3:00) were fit with a mixed linear model with grazing treatment, growth stage, and their interaction as fixed effects and block-treatment and growth stage–block interactions as random effects. Mean midday minima for PRI were fit with the same model plus a covariable of net solar radiation, independently testing each of the three growth stages to correct for non-homogeneity of variances among stages. Significance levels were corrected using the Sidak correction for three independent tests. Model selection for covariables of potential evapotranspiration, date, and net solar radiation was performed using a generalized additive model with smoothing splines of linear and non-linear functions of the predictor variables. All response variables were transformed where necessary to meet assumptions of residual normality, homogeneity of variance, independence, and homogeneity of slopes. Type III analysis of variance was performed for each of the multi-level models with Kenward-Roger approximation for denominator degrees of freedom. In the case of significant interaction terms, a Tukey-adjusted post hoc test was performed on pairwise comparisons at 95% confidence levels.

For temporally dynamic soil water retention and soil matric potential variables, models were fit for both pooled (grazed and ungrazed treatments combined as fixed effect) and unpooled (each treatment as fixed effect) data. For water retention curves, shaping parameters of a dynamic process model were optimized by PORT routines in the “stats” package of R to minimize root-mean-square error. A lag element was included in the model to account for residual autocorrelation and non-independence across time. The response variable (volumetric water content) was normalized with respect to the overall mean to account for different initial saturated soil water contents. To test whether the van Genuchten model parameters were significantly different between treatments, pooled vs. unpooled linear models for water retention with the optimized van Genuchten model and the lag element as explanatory variables were compared using the Akaike information criterion (AIC) for model fit. For soil matric potential, linear models representing the slope of dry-down segments between rain events with a fixed effect of growth stage (no lag element) were compared using the same pooling method.

3 Results and discussion

3.1 Soil moisture content

Soil water content was consistently higher in the ungrazed treatments than in the grazed treatments across measurement methods (Fig. 2). Grazing treatment had a significant effect on volumetric water content (F = 135.9, p < 0.001), as did sampling depth (F = 4.9, p = 0.027; Fig. 2a). Water content was approximately 9% higher in the ungrazed plots compared with the grazed plots (average of depths), and approximately 2% higher in the 0–10 cm layer than in the 10–20 cm layer (average of treatments). Sampling date did not significantly impact volumetric water content according to a conservative degree of freedom test to account for non-independence of consecutive measurements (F = 10.4, p = 0.166).

Fig. 2
figure 2

a Soil volumetric water content (mm3 mm−3) for the 0–10 cm and 10–20 cm layers in grazed (brown, solid lines) and ungrazed (green, dashed lines) treatments. Error bars represent ± 1 standard error from the mean. b Soil water matric potential (kPa) recorded at 60 min intervals with a granular matrix soil water potential sensor at depths of 20 cm and 40 cm. Shaded regions represent ± 1 standard error from the mean. Soybean growth stages indicated below x-axis (early, V2-R1; late-mid, R1-R5; late, R5-R7)

Likewise, ungrazed treatments maintained matric potentials an average of 8 kPa (25%) higher than grazed treatments at both 20- and 40-cm depths across the season. Ungrazed treatments also dried down more slowly than grazed treatments as indicated by a 7% less negative slope of matric potential decline between rain events (F = 4.59, p = 0.03; Fig. 2b).

Soil water retention at 10-cm depth was also greater in ungrazed treatments than in grazed treatments, meaning that volumetric water content was higher in soils from the ungrazed treatments at the same level of suction from evaporative water loss. When water retention curves were fit with the optimized van Genuchten retention function, the ungrazed treatment curve maintained an average of 3% greater volumetric water content beyond 2000-kPa suction (n = 4 ungrazed, n = 3 grazed; Fig. 3), representing the more saturated end of the pressure gradient. Model fit showed significant improvement (smaller AIC) with the inclusion of a grouped treatment effect (AIC = − 40,640 for grouped vs. AIC = − 40,494 for pooled treatment effects; α = 0.05).

Fig. 3
figure 3

Soil water retention curves for grazed (brown, solid line) and ungrazed (green, dashed line) measured by the simplified evaporation method (faded points) and fit with a modified van Genuchten model for soil hydraulic properties optimized for minimum root mean square error. Volumetric water contents were normalized by overall mean to account for different initial water contents among samples

These results suggest that different soil environments in the two treatments significantly impacted the dynamics and availability of soil water throughout the growing season. Direct effects of grazing on soil moisture dynamics include shifts in residue accumulation and soil physical properties. Grazing decreases the amount of residual biomass available at the start of the growing season, which ranges from 6000 to 8000 kg dry matter ha−1 in ungrazed plots to 2500–4000 kg dry matter ha−1 in moderately grazed plots (Carvalho et al. 2018). The surplus of residual biomass in ungrazed treatments keeps the soil surface cooler, reducing evaporation and maintaining water availability in superficial soil layers (Martins et al. 2016). While the hoof action of grazing animals can reduce infiltration rates by compacting surface soil layers, large amounts of residual material even in grazed treatments in this experiment would likely buffer this effect. In any case, compaction from grazing animals is usually limited to the most superficial soil layers and can be ameliorated at depth by root action during the cropping season (Bell et al. 2011; Cecagno et al. 2017). The effects of increased root action in grazed plots on soil physical properties is also evident from past studies at the site, which report significantly larger mean weight diameters and lower bulk density in grazed plots at depths > 10 cm (Souza et al. 2010b; Conte et al. 2011; Cecagno et al. 2017).

Soil organic carbon is often credited with improving soil hydraulic properties due to its favorable impact on soil structural physical attributes such as aggregation and bulk density (Saxton and Rawls 2006). Percent soil organic matter in the 0–10 cm layer was higher in the grazed plots in 2016–2017 (4.0% and 3.2% for grazed and ungrazed plots, respectively), consistent with previous years’ results at the site (Assmann et al. 2014). Higher organic matter proportions in the grazed plots despite retaining only half as much surface residue as ungrazed plots suggests that manure recycling, plant aerial biomass recycling, and plant root deposition during the ryegrass phase may contribute more to the soil organic pool than residue retention in this system (Wilson et al. 2018). Nevertheless, our results indicate that residue retention is more important for soil water content and availability in this system than soil organic carbon.

Grazing may also have an indirect impact on soil moisture dynamics through shifts in plant regulation of soil water storage across rotation periods. Grazing animals continually defoliate forage biomass, which triggers root mortality and turnover in the short term but stimulates compensatory shoot production and continued root growth throughout the season (Fulkerson and Donaghy 2001). As a result, moderately grazed areas accumulate greater overall herbage mass in the course of the season despite having less standing biomass at any given time (Carvalho et al. 2018). In addition to supplementing the soil organic matter pool through root deposition, this process also results in a longer period with actively transpiring vegetation during winter. In contrast, ungrazed plots experience only one winter growth cycle to about 50-cm standing biomass, followed by lodging and senescence. The ungrazed plots may therefore have a lower cumulative plant water requirement over the length of the winter, and consequently greater stored soil water at the beginning of the soybean growing season.

Furthermore, decreased water extraction by plant roots in the ungrazed plots during the cropping season may contribute to slower rates of soil dry down observed in these plots relative to grazed plots. Soils in the Planalto region are often acidic and prone to aluminum toxicity. Following liming for acidity correction in 2010, ungrazed plots remained more acidic and more saturated by aluminum than grazed plots, especially at depths below 10 cm (Martins et al. 2014a). These results were attributed to less root action and therefore decreased nutrient cycling and bioporosity in the ungrazed treatments. A pH of 3.7 (compared with 4.7 in grazed plots), combined with aluminum saturation of 32.5% (compared with 3.1% in grazed plots) in the 2016–2017 season, would be sufficient to severely limit root expansion and water extraction in soybeans (Silva et al. 2001). Soil chemical properties limiting water extraction by soybean roots may thus interact with the above differences in water loss through evaporation to explain consistently wetter soil conditions in ungrazed plots despite periodic rewetting events.

In summary, animal grazing of winter forage in rotation with summer soybean decreased soil water availability throughout the cropping season. This difference is likely due to a combination of decreased residue retention and more extensive exploitation of the soil profile by actively transpiring grazed forages in the winter. Conversely, soils in ungrazed plots had consistently higher water content and were slower to dry down, possibly due to (1) reduced root extraction capability of aluminum-impaired roots or (2) reduced need for root exploration due to sufficient water content in surface layers to meet crop demand. These differences have been inconsequential for soybean productivity in favorable rainfall years and even in two extreme drought years previously recorded at the site (Martins et al. 2016). In one of these years, the authors suggested that redistribution of water storage to deeper soil layers in grazed plots points to a factor that may ultimately favor soybean resilience when water becomes limiting, as water stored in deeper soil layers is less susceptible to losses (Cecagno et al. 2017).

These results indicate that tradeoffs among multiple underlying processes may be at work in determining the productivity of soybean crops in integrated systems. In particular, the improvements in soil physical properties observed over the years in the grazed plots may offset their lower soil water retention and the higher water use by winter pastures to allow grazed plots to achieve similar yield outcomes to ungrazed plots. Although water is more abundant in the ungrazed plots, it may be less available to crops, setting the stage for opposing water cycling functionalities in the two systems. However, different outcomes might be expected under more prolonged drought scenarios or in regions with less reliable climates, suggesting the need for continued monitoring of soil water dynamics in an assortment of integrated systems across climatic zones. For management purposes, these results also suggest that soil water content measurements may be less informative than plant water status or overall plant available water when determining water budgets and yield predictions for integrated systems.

3.2 Plant water status

We hypothesized that altered soil physico-chemical properties stemming from grazing activity could carry over into plant physiological processes. Our results do not support this hypothesis with regard to crop water status, though they do for crop time to maturation and radiation use efficiency (detailed in Sect. 3.3). Various methods for determining evaporative demand indicated that atmospheric conditions during our study period remained between 0.2–2 kPa for the length of the growing season, peaking in early-mid November and declining as the season progressed. This range is well below the 3 kPa vapor pressure deficit thresholds for physiological stress in soybeans. An analysis of covariance between potential evapotranspiration and leaf water potential indicated a significant difference in afternoon leaf water status between grazed and ungrazed treatments when averaged across the season (p = 0.03; Fig. 4). However, this difference was driven by more negative leaf water potentials in the ungrazed treatment during early growth stages (V2-R1), a difference that largely disappeared in later growth stages. Stomatal conductance was not significantly different between treatments regardless of growth stage or potential evapotranspiration (p = 0.19), although there was a non-significant trend of higher stomatal conductance in the grazed compared with ungrazed treatment.

Fig. 4
figure 4

Relationship between potential evapotranspiration and afternoon leaf water potential in soybeans in previously grazed (green, dashed lines) or ungrazed (brown, solid lines) plots at three growth stages (early, V2-R1; late-mid, R1-R5; late, R5-R7) and averaged across the season

Although soil water content and retention were consistently higher in the ungrazed treatments, we observed more negative leaf water potential early in the season in this treatment compared with the plants grown after grazed pasture. This result may be explained by the impacts of greater residue cover on emergence time in the ungrazed treatment; early-season plants in the ungrazed plot may have been slightly younger with a less extensive root system during this drier-than-average period. Alternatively, increased residue loads in the ungrazed plots are also correlated with decreased initial soybean populations compared with the grazed plots (Caetano 2017), resulting in lower density stands and potentially less intraspecific competition for water resources later in the season. While measurements of evaporative demand indicated relatively greater potential water stress during the early vegetative stages of soybean growth in the ungrazed plots, they were well within the limits of normal physiological functioning during the critical grain-filling stages (January–February). High stomatal conductance in plants from both grazing treatments is consistent with lack of water stress.

Past results at this experiment have shown that, even under water-limiting conditions, the difference in soil water content between grazing treatments does not affect relative crop physiological responses. During a severe drought year (less than 40% of average annual rainfall), moderately grazed plots had significantly lower soil water content than ungrazed plots, but both maintained water within the plant-available range and presented no noticeable difference in afternoon leaf water potential or ultimate soybean productivity (Martins et al. 2016). The authors suggested that the ability of crops in the moderately grazed treatments to maintain physiological functions under low soil moisture is a product of the higher-level, systemic interactions at work in integrated crop-livestock systems. For instance, despite lower soil water content early in the season, water availability to the plant may be maintained in grazed plots by improved soil structure (Souza et al. 2010b) and chemical properties (Martins et al. 2014a), resulting in improved water storage at depth and ability of roots to explore the profile. Ultimately, these properties would allow crop plants in previously grazed plots to extract available water more effectively, maintaining system function despite lower absolute soil water content. Thus, apparent negative impacts on a factor at the field level may not translate to noticeable effects at the system level, especially when environmental conditions are optimal.

3.3 Crop performance and canopy dynamics

Soybean yields in 2017 were similar between treatments at 4208 kg ha−1 and 4157 kg ha−1 in the grazed and ungrazed plots, respectively (p = 0.90), about 18% higher than the regional average for soybean yields in Rio Grande do Sul state (CONAB 2017).

While NDVI dynamics for both grazed and ungrazed treatments trended together for most of the growing season (December–February), a significant drop in NDVI corresponding with senescence occurred in the ungrazed treatment approximately 2 weeks earlier than in the grazed treatment (p = 0.01; Fig. 5a). Coupled with this earlier senescence were significantly higher PRI reflectance signatures from the ungrazed treatment (p = 0.03; Fig. 5b) beginning as early as January (mid-season). Likewise, diurnal PRI trends in ungrazed plots were higher than in grazed plots across days, especially in the early- and mid-season growth stages (p < 0.001; Fig. 5c).

Fig. 5
figure 5

Seasonal midday maxima and minima for a normalized difference vegetation index (NDVI) and b photochemical reflectance index (PRI) of soybean canopy in previously grazed (brown, solid lines) or ungrazed (green, dashed lines) plots. c Diurnal trends for PRI of soybean canopy in early- (V2-R1), mid- (R1-R5), and late-season (R5-R7) growth stages of previously grazed or ungrazed plots

NDVI signatures are associated with leaf-area index and biomass production due to the reflectiveness of leaves to near-infrared light (Gamon et al. 2015). PRI, on the other hand, captures changes in reflectance in a much narrower wavelength region: the green and yellow areas of the spectrum associated with carotenoid absorption and shifts in the xanthophyll cycle. Seasonal PRI dynamics are thus influenced by changing chlorophyll to carotenoid pigment concentration ratios, while diurnal PRI trajectories track photochemical and non-photochemical quenching mechanisms, making the latter an indicator of light-use efficiency in photosynthesizing organisms (Porcar-Castell et al. 2012; Gamon et al. 2015). Diurnal trends for PRI typically reach a minimum towards solar noon, when photosynthetic photon flux density is at its peak (Gamon et al. 2015). Therefore, seasonal PRI signals indicate a higher relative chlorophyll concentration in soybeans grown in ungrazed plots, while diurnal PRI signals suggest higher light-use efficiency in ungrazed plots, particularly during the critical grain-filling period (R5-R7). Taken together, these indices point towards faster physiological maturation and earlier senescence in soybeans grown in ungrazed relative to grazed plots.

Differences in soil water availability are not enough to explain altered light-use efficiency and time to maturation, as water was not limiting in either treatment. Furthermore, accelerated maturation in ungrazed plots appears to contradict observations of delayed emergence due to heavier residue retention compared to the grazed plot (Caetano 2017). Earlier yellowing of plants in ungrazed plots could indicate nitrogen limitation due to waterlogging, decreased rhizobium activity, or nitrogen immobilization by excess carbon inputs from retained residues. However, these proposed mechanisms are not supported by the higher chlorophyll to carotenoid ratios and increased light-use efficiency indicated by PRI. Instead, higher chlorophyll concentrations in the ungrazed treatment despite near-saturated soil water conditions suggest earlier remobilization of resources to grains and thus reduced time to maturation in the ungrazed plots.

The ability of moderately grazed plots to maintain levels of crop productivity equivalent to ungrazed plots is notable given these alterations in plant physiological processes. These results appear to hold regardless of agrometeorological conditions and/or negative differences in variables that presumably affect individual plant performance (e.g., soil water availability, maturation time). Nevertheless, the lengthier dry-down period that may be required for soybean crops in the grazed plots has the potential to affect management operations, especially if different plots ripen at different times, and should be taken into consideration in management plans for integrated systems.

In summary, soils in previously grazed plots had lower overall soil water content and retention, but this difference did not affect plant water status. Soybean plants in previously grazed plots stayed green longer and used radiation less efficiently but produced equivalent yields to soybeans in ungrazed plots. Our results suggest that tradeoffs in soil water depletion by grazed pastures on the one hand are countervailed by improvements in soil quality and plant available water capacity on the other. These processes tend to stabilize the performance of crops in rotation, leading to equivalent production despite altered physiological dynamics and soil environments. These insights inform management of integrated systems under varying environmental conditions by bridging our understanding of grazing animal effects on soil conditions with the carry-over effects on crop productivity over the long term. In particular, they highlight the importance of good grazing management. Avoiding overgrazing in these systems not only ensures adequate residue coverage to prevent surface layer degradation from trampling but also stimulates continuous pasture root action and the associated co-benefits for soil physical and chemical properties. Our results serve as a springboard for further investigation into the interactions at play in integrated systems, whether for adaptive management purposes or for assessing resilient responses to environmental stresses in the field, laboratory, or model simulation setting.

4 Conclusions

We showed that 16 years of beef cattle integration into soybean systems in southern Brazil has shifted the dynamics of key biophysical parameters, including soil water content and retention and crop growth and phenology, compared with ungrazed systems. Despite these differences, soybean yields were equivalent between treatments, suggesting inter-seasonal tradeoffs between pasture water use, soil physico-chemical properties, and crop water use in grazed compared with ungrazed plots. These results raise questions as to other possible mechanisms (e.g., microbial processes) governing integrated crop-livestock system productivity and potential biophysical buffering capacity, especially under less optimal environmental conditions. Future work should focus on determining the soil, plant, and microbial attributes that most contribute to the ability of grazed plots to maintain crop productivity even under altered soil water dynamics, and the frequency with which we might expect yield costs or benefits in these systems over the long term.

Beyond the processes explored here, there are multiple other interacting factors and trade-offs that govern overall outcomes for crop productivity and system function in integrated crop-livestock systems, including crop species in rotation, soil type, weather conditions, tillage regime, grazing management practices, and microeconomic and socio-cultural factors. This study’s aim was to add an intensive examination of within-season soil water and canopy dynamics during an above-average precipitation year to the growing database of grazing animal impacts on agroecosystem functions in integrated crop-livestock systems. Such efforts will contribute to the optimal management, sustained productivity, and resilience of these increasingly important intensive production systems under a variety of enviro-climatic conditions.