Abstract
Integrated crop-livestock systems can have subtle effects on soil quality over time, particularly in semiarid regions where soil responses to management occur slowly. We tested if analyzing temporal trajectories of soils could detect trends in soil quality data which were not detected using traditional statistical and index approaches. Principal component and cluster analyses were used to assess the evolution in ten soil properties at three sampling times within two production systems (annually cropped, perennial grass). Principal component 1 explained 33% of the total variance of the complete dataset and corresponded to gradients in extractable N, available P, and C : N ratio. Principal component 2 explained 25.4% of the variability and corresponded to gradients of soil pH, soil organic C, and total N. While previous analyses found no differences in Soil Quality Index (SQI) scores between production systems, annually cropped treatments and perennial grasslands were clearly distinguished by cluster analysis. Cluster analysis also identified greater dispersion between plots over time, suggesting an evolution in soil condition in response to management. Accordingly, multivariate statistical techniques serve as a valuable tool for analyzing data where responses to management are subtle or anticipated to occur slowly.
1. Introduction
Integrated crop-livestock systems (ICLS) are recognized globally for their contributions to improve agricultural sustainability [1–4]. Recent emphasis on conservation agriculture, climate-smart agriculture, and sustainable intensification has underscored the potential role of ICLS to create more productive and resilient agricultural systems [5, 6]. An inherent emphasis on multiple enterprises makes ICLS well-suited for future growing conditions, where production synergies between enterprises can serve to enhance adaptability to increasingly variable weather and market conditions, while concurrently minimizing input costs [7, 8].
Integrated crop-livestock systems can be more management and labor intensive than single-enterprise production systems [4]. Accordingly, producers need to know if investments in implementing ICLS translate into improvements in key response metrics. Relevant information would utilize metrics characterizing agroecosystem sustainability and, ideally, do so quickly following adoption of management practices to assess trajectory over time [9]. Such information could serve to justify investments in ICLS or, if necessary, provide quantifiable guidance for adjusting management to more effectively meet production, economic, and/or environmental goals.
Soil quality serves as an important response metric to agroecosystem management given its foundational role in affecting agricultural and environmental outcomes through impacts on ecosystem services [10]. Formally defined, soil quality refers to the capacity of soil to function [11]. Accordingly, soil quality contributes to numerous ecosystem services within agricultural landscapes (e.g., water retention/filtration, climate regulation, and biodiversity conservation) [12].
The central role of soil to deliver key ecosystem services has served as a focal point for the development of numerous assessment methods. Over the past 30 years, soil quality assessments have evolved from the evaluation of a minimum dataset of soil properties [13] to elaborate indices involving scoring functions based on management goals across multiple scales [14–16]. As agricultural practices are intensified to meet anticipated production and energy needs [17], soil quality assessments will be essential to ensure deployed practices are sustainable.
Previous soil quality evaluations of ICLS have focused on treatment effects on single soil physical, chemical, or biological properties or indices of multiple properties [4, 18–21]. Outcomes from these evaluations have helped understand soil responses to ICLS, as results are typically easy to interpret and convey, thereby improving understanding by clientele. Despite this benefit, there are unique attributes associated with ICLS that may make traditional assessment of soil property dynamics, whether by single properties or an index, problematic. Soil properties can change slowly in ICLS, particularly in semiarid regions where soil responses to management occur slowly [22]. Accordingly, long-term research may be necessary to assess ICLS effects on soil quality, yet short-term funding cycles may make it difficult to conduct such research. Secondly, traditional Soil Quality Index methods “merge” outcomes from multiple soil properties to provide an aggregated numerical rating [11, 16]. In doing so, traditional indexing methods may mask nuances inherent to the data and, thus, miss important spatiotemporal trends.
Given this context, we sought to analyze near-surface soil quality indicator data collected over a six-year period from an ICLS in a semiarid region using principal component analysis (PCA). These data were previously analyzed using traditional statistical and index approaches to evaluate soil property responses to residue management, frequency of hoof traffic, season, and production system [18, 23]. We hypothesized a reanalysis of the same data with multivariate statistical techniques would identify trends in soil property dynamics not previously observed.
2. Materials and Methods
2.1. Site and Treatment Description
The ICLS experiment was located at the USDA-ARS Northern Great Plains Research Laboratory near Mandan, ND, USA (46°46′ W, 100°54 N). The experimental site is within the temperate steppe ecoregion of the USA, with a semiarid climate characterized by long, cold winters and short, hot summers. Mean annual precipitation and temperature at the site are 414 mm and 4°C, respectively, and the average frost-free period is 131 days. Gently rolling uplands (0 to 3% slope) characterize the topography of the experimental site, and predominant soil types are Temvik-Wilton silt loams (fine-silty, mixed, superactive, frigid Typic, and Pachic Haplustolls) [24].
A thorough treatment description can be found elsewhere [18]. Briefly, two 6.0 ha crested wheatgrass [Agropyron desertorum (Fisch. ex. Link) Schult.] pastures were sprayed with glyphosate [N-(phosphonomethyl) glycine; 0.7 kg a.i. ha−1] twice in mid-May and converted to an annual cropping sequence of oat/pea (Avena sativa L./Pisum sativum L.), triticale/sweet clover (Triticum aestivum x Secale cereale/Melilotus officinalis L.), and corn (Zea mays L.). Beginning in 2007, the crop sequence was changed to oat/alfalfa (Medicago spp.)/hairy vetch (Vicia villosa Roth)/red clover (Trifolium pratense L.), brown midrib sorghum-sudangrass (Sorghum bicolor L. Moench)/sweet clover/red clover, and corn also using no-till planting techniques. Each phase of the three-year cropping sequence was present in both pastures, which were used as replicates. Management decisions related to seeding, fertilizer, and weed control by herbicides followed recommended practices by area producers.
The oat and triticale/sorghum crop mixtures were harvested for grain from mid-August to early-September with the straw spreader removed from the combine, which created a swath of crop residue for winter grazing. The corn was swathed for forage in mid- to late-September. Each crop strip was split into three residue management treatments that included no residue removal (Control; 0.05 ha), residue removal with a baler (Hayed; 0.05 ha), and residue removal by grazing with livestock (Grazed; 1.69 ha). The Control and Hayed treatments were randomly assigned within each crop strip. For the Grazed treatment, the swathed crop residues from the cropping sequence represented winter forage for ten 4–6-year-old nonlactating bred Hereford cows, due to calve in late March. Plot areas within a crop phase were 27 m × 23 m (0.06 ha) for the Control and Hayed treatments and 54 m × 312 m (1.68 ha) for the Grazed treatment. Stocking rate within a crop phase for the Grazed treatment was 0.2 ha cow−1 (5.9 cows ha−1), corresponding to winter grazing practices used by local producers.
Grazing commenced in mid-November and ended in mid-February, with the oat crop mixture grazed first, triticale/sorghum crop mixture second, and corn last. Access to crop swaths was controlled using electric fences oriented at right angles to the swaths. Fences were moved daily to provide access to fresh forage. A shelter and “frost-free” water fountain were located at the end of each pasture within the Grazed treatment.
Two perennial grass pastures were used for comparison to the annually cropped treatments. The pastures, each 6.0 ha, were composed of a mixture of native and introduced cool-season perennial grasses. Similar to Grazed treatments, swathed grass from the perennial pastures was used as winter forage from 1999 through 2002 for the same number of cows managed in the annually cropped treatments. Within each perennial pasture a nongrazed strip was split into Hayed and Control treatments as outlined above. Drought-induced limitations in forage in 2002 and 2003 resulted in the perennial grass treatments being hayed but not grazed in 2003. From 2004 through 2008, the perennial grass treatment was not swathed but lightly grazed with ten Hereford or Angus cows from mid-October to mid-January.
2.2. Soil Sample Collection, Processing, and Analyses
Soil samples were collected in April of 2002, 2005, and 2008 after the swathed crops had been grazed but not replanted. Within the annually cropped pastures, samples were collected from nine subplots (three per crop phase) in the Control and Hayed treatments, oriented randomly in each treatment but between crop rows. Samples from the Grazed treatment were collected from two transects differing in frequency of hoof traffic, also between crop rows. Nine subplots (three per crop phase) were established in each transect perpendicular to crop swaths approximately 100 m (representing high-traffic; HT) and 200 m (representing low-traffic; LT) from the shelter and water source. Sampling protocol for the western wheatgrass pastures followed that in the annually cropped pastures, with the exception of fewer subplots in each treatment (three in the Control, Hayed, and Grazed hoof traffic transects). Within each subplot, six soil cores were collected from the 0 to 7.5 cm depth using a 35 mm (i.d.) step-down probe and composited. Each sample was saved in a double-lined plastic bag, placed in cold storage at 5°C, and analyzed within 6 wk of collection.
Soil samples were dried at 35°C for 4 d and then ground by hand to pass a number 10 (2.0 mm) sieve. Identifiable plant material (>2.0 mm diameter, >10 mm length) was removed during sieving. Electrical conductivity and pH were estimated from a 1 : 1 soil-water mixture [25, 26]. Soil NO3-N and NH4-N were determined from 1 : 10 soil-KCl (2 M) extracts using cadmium reduction followed by a modified Griess-Ilosvay method and indophenol blue reaction [27]. Plant-available soil P was estimated by bicarbonate extraction [28]. Potentially mineralizable N was estimated from the NH4-N accumulated after a 7 d anaerobic incubation at 40°C [29]. Total soil C and N were determined by dry combustion [30], and as soil pH was <7.2 within the surface 7.5 cm depth, total soil C was considered equivalent to SOC. All data were expressed on an oven-dry basis prior to data analyses.
2.3. Data Organization
Prior to statistical analyses, data were organized to facilitate characterization of all sampled locations as well as sampled locations specific to the annually cropped treatments. Data group 1 (S1) was composed of 10 soil properties from 40 plot locations in both annually cropped and perennial grass treatments (Table 1). Within S1, there were 24 locations in the annually cropped treatments and 16 locations within the perennial grass treatments. Data group 2 (S2) focused on the 24 locations within the annually cropped treatments, with an equivalent number of grazed and ungrazed locations (12 each). For both S1 and S2 data groups, data were included from each sampling time: 2002, 2005, and 2008.
2.4. Statistical Analyses
Individual plot dynamics were quantified by combining a series of multivariate analyses using a method developed by Dolédec and Chessel [31] and adapted by García-Martínez et al. [32] and Ryschawy et al. [33]. The method distinguished impacts of external (the system environment) and internal (system structure and functioning) factors, while considering the evolution of individual plots over time. We applied this statistical methodology to S1 to analyse plot trajectories of near-surface soil quality indicators (Table 1).
To analyse plot trajectories, we used a particular type of PCA (Principal Component Analysis) called Within-Class Analysis (WCA) to consider evolution of individuals (plots) while not considering the general date effect on different plots. We thus considered the deviation of each soil property in a plot from the average across plots, considering time period effect. The first step of the WCA consisted in standardizing the data. For each soil property (column), data were centered by subtracting the mean of the column to the values and scaled by dividing the values of each column by their standard deviation (with a column mean of 0 and a standard deviation of 1). Within-Class Analysis used an R statistical software package dedicated to Analysis of Environmental Data (Exploratory and Euclidean Method; ade4 package) [34]. Within-Class Analysis used orthogonal PCAs where factors were used as covariables to partition rows (time was used as the main factor to partition rows considering each date separately for plots). The analysis characterized differences between plot locations once a general common trend over time was eliminated. We considered the three sampling times (2002, 2005, and 2008) as the factor partitioning rows.
Main factors within data group S1 were summarized by PCA. We then performed a WCA resulting in a data frame W1 containing transformed data for each of the 40 plot locations at the three different dates considered, as the deviation between individual and date-class mean coordinates on within-analysis axes. The results of the WCA therefore represented interplot location trajectories.
To summarize soil trajectories for each individual plot location for all dates, we performed a final PCA on the WCA individuals considered at each date. Per the Kaiser criterion, we thus considered the main factors with eigenvalues >1 as summarizing observed changes. A Hierarchical Cluster Analysis (HCA) (with squared Euclidean distance and Ward’s aggregation method) was carried out on the four main factors of the PCA to establish a typology of plot locations according to their change over time.
To analyse changes in plot locations elucidated by HCA, we calculated means and standard deviation of variables for each measured variable. All statistical analyses were performed using R 3.1.0 software [34]. We repeated this statistical methodology on data group S2 to evaluate soil property dynamics based on grazed (Grazed) and ungrazed (Control, Hayed) areas within the 24 plot locations assigned to annually cropped treatments.
3. Results
3.1. Analysis of All Sampled Locations
The first four principal components accounted for 82.3% of the variability based on WCA of the whole sample (data group S1). Associations between soil properties were evident following the WCA, as three distinct groupings were observed graphically (Figure 1). Soil nitrate (NO3N), ammonium (NH4N), and available phosphorus (P) appeared linked, as did soil carbon to nitrogen ratio (CNRATIO), soil organic C (SOC), total nitrogen (TN), and electrical conductivity (EC). Soil pH (PH) and potentially mineralizable N were also linked. Only soil bulk density (SBD) was not grouped with another soil property.
Contributions of measured soil properties to the four principal components are provided in Table 2. Principal component 1 explained 33% of the variability and corresponded to gradients of available N and P (NO3N, P) and CN ratio levels. Principal component 2 explained 25.4% of the variability and corresponded to gradients of PH, SOC, and TN. Principal component 3 explained 13.6% of the variability and corresponded gradients in EC and PMN. Principal component 4 explained 10.3% of the variability and was primarily influenced by EC and NH4N.
Dispersion between plots showed a trend toward greater dissimilarity in soil properties over time (Figure 2). Based on standard deviation of mean values, plot differences were closer in composition in 2002 and 2005 than in 2008, with the latter showing substantial dispersion at the end of the 6 yr study period (Table 3).
3.2. Analysis of Production Systems
Cluster analysis found distinct groupings based on production system, as annually cropped treatments and perennial grasslands were clearly distinguished graphically (Figure 3). Such trends were evident by comparing soil properties in each treatment between 2002 and 2008 (Table 4). Plots under annual cropping had a greater accumulation of mineral nutrients (NO3N, NH4N, and P) compared to perennial grassland, whereas perennial grasslands had greater accumulation of organic matter (SOC, TN), increased mineralizable N (PMN), and a lower dispersion between plots compared to annually cropped treatments.
3.3. Analysis of Grazing Effects under Annual Cropping
In the analysis of grazing effects (data group S2), 83.4% of the variability was explained based on WCA. Principal component 1 explained 29.8% of the variability and corresponded to gradients in total carbon and nitrogen accrual rates (SOC, TN, and CNRATIO), whereas principal component 2 explained 24.7% of the variability and corresponded to gradients of mineral nutrient accrual rates (NO3N, NH4N, and P) (data not shown).
Hierarchical cluster analysis revealed four clusters that were differentiated by grazing and treatment replication (Figure 4). Treatment replications were distinguished graphically along the horizontal plane, with one replication (Field A) left of the vertical plane and the other replication (Field B) right of the vertical plane. Differences in soil properties between treatments were limited to Field B, where soil pH was lower and available P was higher under grazing (Table 5).
4. Discussion
Traditional statistical and index-based approaches for assessing soil quality may provide a limited characterization of soil property dynamics [35]. Multivariate statistical analyses such as PCA provide a useful method for screening a diverse collection of functionally relevant soil properties to identify potential data trends [36]. These techniques provide an objective approach to extract and weight information in complex datasets [37], thereby providing a potentially valuable tool for researchers seeking to link the status of soil properties, soil function, and agroecosystem management.
Analysis of whole sample dynamics reflected the predominant influence of soil chemical properties on total variance in the dataset. Soil properties most influenced by additions of N and P through fertilization (e.g., NO3N and P) accounted for substantial variance between plots. Annual change in nutrient-related properties would be expected to vary based on weather-driven factors affecting plant growth and nutrient loss [38]. Electrical conductivity, a relative measure of the total quantity of ions in soil solution, also accounted for an appreciable amount of variance in the dataset.
Previous PCA evaluations characterizing soil quality have found organic matter-related properties to account for a large amount of variance in data [35, 36]. A similar outcome was observed in this evaluation, as SOC, TN, and CNRATIO accounted for a moderate amount of variation among soil properties across principal components. The limited measurement timespan (6 yr), coupled with the inherent difficulty in detecting change in properties with a large background signal due to previous management as a perennial grassland [39], may have compromised the usefulness of organic matter-related properties in this evaluation. In general, soil properties associated with organic matter are slow to change in semiarid agroecosystems, sometimes requiring a decade or longer for management effects to be detected [22].
Despite limited change in some soil properties in this evaluation, the aggregate visual assessment of data dispersion proved useful in detecting a general evolution of the dataset. Such dispersion was confirmed upon inspection of standard deviation values for soil properties within each year, as mean values of all soil properties, except CNRATIO, had the greatest standard deviation in the last year of the evaluation. Coupled graphical and tabular data characterizations seem particularly useful for datasets like the one analyzed in this evaluation, as visual assessments infer possible data trends generally, while tabular assessments provide confirmation specifically.
Outcomes from the cluster analysis of the complete dataset showed data groupings based on production system, with perennial grass and annually cropped treatments separated visually. Such separation was evident with cluster analysis despite limited observed differences in mean values of soil properties between production systems (Table 4) [18]. Moreover, multivariate analysis and graphical representation of these previously analyzed data clearly showed limitations associated with the aggregated Soil Quality Index (SQI) approach to assessing soil quality, as differences in the overall SQI between production systems were not observed in any year [18]. From a soil management perspective, outcomes from the cluster analysis underscored the value of perennial grass systems to enhance soil quality through increased accrual of soil carbon and nitrogen, while minimizing levels of nutrients susceptible to loss. Such attributes are important features of sustainable agricultural systems [40] and are central to agroecosystem resilience to climate-induced stressors and improvements in environmental quality from agricultural production [41, 42].
In contrast to whole production system analyses, outcomes were less distinct when data were restricted to annually cropped fields with and without grazing. Cluster analysis revealed an overarching effect of replication on data groupings, of which grazed and ungrazed plots lacked homogeneity. It is possible that the limited number of data points, coupled with a need for greater time-in-treatment to resolve management effects on soil properties, resulted in indeterminate responses to PCA. It is important to note that no differences in soil properties were found between grazed and ungrazed systems in any year when the same data were analyzed using traditional statistical analyses [18]. Future evaluations of grazed and ungrazed treatments in both annually cropped and perennial grassland systems are warranted given the important role of livestock to affect biomass productivity through influences on soil condition [43].
Collectively, statistical methods used in this evaluation contribute to an ongoing evolution of soil quality assessment. As the interpretation of individual soil properties can be complicated by multiple response functions [16], use of aggregate approaches, such as PCA and HCA, potentially offer novel insights into how soil responds to management [14]. Use of statistical approaches capable of resolving changes in soil condition in the near-term can be particularly useful in regions where alterations in individual soil properties require a decade or more to be detected [22]. Similarly, these approaches may be valuable in detecting subtle changes in soil induced by systemic drivers (i.e., climate). The capacity to detect shifts in soil condition quickly under the context of anticipated climate change will be increasingly valuable, particularly for dryland agriculture [44]. Timely management interventions will be essential in the future to mitigate soil degradation and maintain soil function.
5. Conclusions
Use of multivariate statistical analyses revealed trends in data not previously detected using traditional statistical and index approaches for assessment of soil quality. Principal component and hierarchical cluster analyses provided a helpful means to visually discriminate between annually cropped and perennial grass treatments, while concurrently identifying a clear trend toward greater dispersion in the data over time. Based on the typology used in this evaluation, apparent multivariate statistical techniques represent a valuable tool for analysis of data collected at different time periods where changes in response variables are anticipated to occur slowly.
Application of these techniques to datasets in other semiarid regions may contribute toward identifying trajectories in soil responses to applied treatments in the near-term. Such information could justify research expenditures, and therefore continuation of experimental treatments until soil responses can be verified with traditional statistical approaches and/or indices. Conversely, these techniques could be used to inform possible changes in applied treatments if visual discrimination is not detected.
Abbreviations
| CNRATIO: | Soil carbon to nitrogen ratio |
| EC: | Electrical conductivity |
| ICLS: | Integrated crop-livestock systems |
| NH4N: | Soil ammonium |
| NO3N: | Soil nitrate |
| P: | Available P |
| PCA: | Principal component analysis |
| PH: | Soil pH |
| PMN: | Potentially mineralizable N |
| SBD: | Soil bulk density |
| SOC: | Soil organic carbon |
| TN: | Total nitrogen |
| WCA: | Within-Class Analysis. |
Disclosure
The US Department of Agriculture (USDA) prohibits discrimination in all its programs and activities on the basis of race, color, national origin, age, disability, and, where applicable, sex, marital status, family status, parental status, religion, sexual orientation, genetic information, political beliefs, and reprisal, or because all or part of an individual’s income is derived from any public assistance program (not all prohibited bases apply to all programs.). USDA is an equal opportunity provider and employer. Mention of commercial products and organizations in this manuscript is solely to provide specific information. It does not constitute endorsement by USDA-ARS over other products and organizations not mentioned.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors gratefully acknowledge the contributions of Justin Hartel, Marvin Hatzenbuhler, and Robert Kolberg for planting and swathing crops, Clay Erickson, Gordy Jensen, and Curt Klein for managing cattle, and John Bullinger, Jason Gross, and Holly Johnson for assisting with soil sampling, processing, and analyses.