2.1 Soil sampling

Topsoil samples for this study were collected on mineral soils (i.e., <300 mg⋅g⁻¹ SOM) through two distinct programs: the Ontario Topsoil Sampling Program (OTSP) and the Ontario Soil Survey (OSS). The OTSP is a collaborative project between the Ontario Ministry of Agriculture, Food and Rural Affairs (OMAFRA) and the School of Environmental Sciences, University of Guelph, with the goal to characterize soil physical, chemical, and biological properties in the surface layer of targeted agricultural soils across southern Ontario, Canada (Fig. 1), to better understand soil health. The OSS, on the other hand, is an OMAFRA program primarily designed for the production of digital soil maps, with ongoing projects in Peterborough County and Ottawa (Fig. 1). In both programs, the Ap (horizon disturbed by human activity, also called plough layer; Canadian System of Soil Classification (Agriculture and Agri-Food Canada 1998)) horizons were sampled, but in the OTSP, sampling depth was terminated at 30 cm. Sampling design for both programs are described briefly below.

Fig. 1.

Map of the Ontario Topsoil Sampling Program, the Ottawa soil survey, and the Peterborough soil survey study areas with sampling locations. Topographic base maps courtesy of the Ontario GeoHub, Ontario Ministry of Natural Resources and Forestry. Map projection: MNR Lambert Conformal Conic, GCS North America 1983.

Sample locations for the OTSP were selected using a conditioned Latin hypercube sampling (cLHS) approach (Minasny and McBratney 2006), which seeks to generate a maximally stratified sampling plan across a set of given environmental covariates for the area of interest. Given the goal of understanding soil health, a combination of covariates to ensure stratification of the sampling locations across topographic (e.g., digital elevation model (DEM), curvature, slope gradient, topographic position index, topographic wetness index, deviation from maximum elevation), environmental (e.g., soil series, gamma-radiometric (potassium, thorium, uranium), biodiversity (i.e, ecodistricts), and agricultural (e.g., annual crop inventory (Agriculture and Agri-Food Canada 2018)) characteristics of southern Ontario was selected. The source DEM was the provincial DEM for Ontario at 30 m spatial resolution (Ontario Ministry of Natural Resources and Forestry 2021). All covariates not generated from the DEM were resampled to match the grid of the DEM using either bilinear (continuous) or nearest neighbour (categorical) interpolation. Sample size specified to the cLHS algorithm was 500, and these locations were named “seed points”. For each seed point, three samples were collected to represent upper, mid, and lower slope positions in landscapes with sloping terrain. In level landscapes, three samples were collected with a minimum 50 m separation between samples. Samples were collected by excavating a small soil pit (30 cm by 30 cm) to the base of the topsoil and sampling from all sides of the soil pit with a soil knife, which was cleaned between soil pits. Sampling was restricted to either the depth of the topsoil, or a maximum depth of 30 cm for sites where the topsoil exceeded 30 cm (see   Fig. S1 (cjss-2021-0145suppla.docx) for distribution of topsoil sampling depth). Soil was homogenized by hand and a subsample put into a cooler during transportation and immediately refrigerated at 4 °C until processed in the laboratory for evolved CO₂ and NH₃ analyses (see the detailed methodology in Section 2.2). The OTSP is ongoing; soil sampling occurred in 2019–2022, with about 50% of samples collected during summer 2019. Only samples collected in 2019 and 2020 were included in our analysis at this time (largely due to uncertainty of soil sampling at the time of statistical analysis and lack of additional funding in 2021). For this study, soil analyses of interest were particle size, evolved CO₂ and NH₃, PMN, and SOM (see the detailed methodology in Section 2.2).

Sample locations for the OSS were also selected using the cLHS algorithm from covariates representing topographic, environmental, and agricultural characteristics of the study areas. In the case of the OSS, the study areas are much smaller than the OTSP; therefore, covariates were all generated at 10 m spatial resolution. Since the soil survey activities were not limited to topsoil sampling of agricultural fields, only Ap horizons were extracted from the OSS datasets for this analysis. Topsoil samples for the OSS were collected in the months of May to November over a 4-year period (2016–2019). Unlike the OTSP, sampling of topsoil horizons in the OSS was not limited to the upper 30 cm when such horizons were thicker than 30 cm (see  Fig. S2 (cjss-2021-0145suppla.docx) showing distribution of Ap horizon depth). Samples were collected from all sides of the soil pit and homogenized. The OSS project is not a soil health project; particle size analysis and SOM are routine analyses for soil survey activities and were utilized in this study on soil health indicators (see details on methodology in Section 2.2).

2.2 Soil analyses

Solvita CO₂-burst and Solvita labile amino N (SLAN), commercially available soil health tests, are indicators of soil microbial activity. Evolved CO₂, as measured by Solvita CO₂-burst, was evaluated following the updated protocol released in 2019 by the Woods End^® Laboratories Inc., Mt. Vernon, ME. Briefly, a 30 cc scoop (about 25–35 g) of 2 mm sieved and oven-dried (40 °C) soil subsample was used (Brinton 2019). Each soil subsample was transferred to a 50 mL plastic vial. Depending on the weight of the soil subsample, 9–10 mL of distilled water was added to each vial with a transfer pipette (Brinton 2019). Solvita CO₂-burst paddles (Woods End^® Laboratories Inc., Mt. Vernon, ME) were inserted into each soil subsample. The vials (consisting of soil and paddles) were placed in 475 mL glass jars; jars were sealed and left undisturbed for 24 h at room temperature (22 °C). Using a digital colorimeter reader, evolved CO₂ was recorded and expressed as mg CO₂–C kg⁻¹ soil.

For evolved NH₃, as measured by SLAN, a 4 g subsample of 2 mm sieved and oven-dried soil was used; the soil subsample was added to plastic vials containing 10 mL of 2N NaOH. The vials, along with soil subsample and SLAN paddles (Woods End^® Laboratories Inc., Mt. Vernon, ME), were transferred to 250 mL glass jars. Like evolved CO₂, glass jars were sealed and left undisturbed for 24 h at room temperature (22 °C). Evolved NH₃ was recorded using the digital colorimeter reader and expressed as mg NH₃–N kg⁻¹ soil.

Soil organic carbon was determined as the difference between total carbon and inorganic carbon. Prior to quantification of inorganic carbon, soil samples were combusted at 475 °C for 3 h. Using the dry combustion method (LECO 828 Series Carbon/Nitrogen Analyzer), total carbon and inorganic carbon were determined. As described by Pribyl (2010), SOM was determined by multiplying soil organic carbon values using the conversion factor of 1.724, commonly known as the “van Bemmelen factor”, which assumes that organic matter is 58% carbon. Particle size analysis was done using the pipette method with hydrogen peroxide pretreatment to remove SOM as described in Sheldrick and Wang (1993).

PMN was determined following the procedures outlined in the Cornell University CASH Laboratory Standard Operating Procedures by Schindelbeck et al. (2016). Briefly, a fresh soil sample (8 g) sieved through 2 mm is anaerobically incubated for 7 days (Schindelbeck et al. 2016). Nitrogen mineralization potential is estimated by taking the difference between the preincubation and postincubation ammonium–N concentration (Schindelbeck et al. 2016; Rogers et al. 2018). This 7-day anaerobic incubation was added to this study after its initial conception. As a result, fresh soil samples were not available for all sampling locations; however, air-dried soil samples were available. To allow an adjustment to the PMN results for the air-dried samples, an experiment was conducted on 30 samples to calculate a correction factor for converting PMN-dry to PMN-fresh basis. PMN was determined for the 30 soil samples in both fresh and air-dry conditions, and a linear regression was developed with y intercept set to zero, and all samples adjusted according to PMN-fresh = PMN-dry × 0.352 (R² = 0.89, P < 0.0001; Fig. 2). Hereafter, only the fresh-basis PMN data are reported and discussed.

Fig. 2.

Linear regression of potentially mineralizable N (PMN) determined on soil samples in both air-dry and fresh conditions used to establish a correction factor from dry- to fresh-basis PMN.

2.3 Scoring functions

For all the tested soil health indicators, a “more is better” approach was used to provide a scoring value, where high score values indicate greater soil health (Fine et al. 2017; Wu and Congreves 2022). Soil health scores were categorized into ratings of very low (0–20), low (20–40), medium (40–60), high (60–80), and very high (80–100) classes, and to facilitate interpretation, a five-color scale consisting of red, orange, yellow, light green, and dark green was applied to each rating (Fine et al. 2017). Scoring functions were determined for each soil health indicator using the following procedures. First, data outliers were removed using the interquartile range (IQR) technique. For each soil property, the IQR of the data was calculated as the difference between the 75th and 25th percentiles. The upper limit of the data to be retained for further analysis were calculated as 75th percentile + 1.5 × IQR and the lower limit as 25th percentile − 1.5 × IQR, also known as Tukey's inner fences (Hoaglin 2003). Once the outliers were removed, the individual observations were grouped based on two different systems (Table 1): (i) into one of three soil textural groups (coarse, medium, or fine) as outlined in CASH (Moebius-Clune et al. 2016) and (ii) into one of five soil textural groups (clays, clay loams, loams, sandy loams, and sands) based on the OASHCS recommendations (Ontario Ministry of Agriculture, Food, and Rural Affairs 2018). The nonparametric Kruskall–Wallis (Kruskal and Wallis 1952) rank sum test was used to detect an effect of texture group, as it is commonly used as the alternative to a one-way analysis of variance for nonparametric data (Brunner et al. 2020), followed by the post hoc Dunn's test (Dunn 1961) with Bonferroni P-value adjustment to determine which texture groups differed significantly. Brunner et al. (2020) showed that unequal sample sizes can result in erroneous results from the Kruskal–Wallis rank sum test when comparing more than two groups, and that an alternative is to use pseudo-ranks. The pseudo-ranks technique was therefore applied using the “pseudorank” package (Happ et al. 2020) in R version 4.1.0 (R Core Team 2020). Where the Kruskal–Wallis and Dunn's tests confirmed significant differences between texture groups, the scoring functions were developed by texture groups, which were significantly different. If no significant differences existed between texture groups, the data were combined. Scoring functions were determined by calculating the mean and standard deviation of the data and generating a cumulative normal distribution, which becomes the scoring function (Moebius-Clune et al. 2016) and assigns a score ranging from 0 to 100.

Table 1.

Textural groups of the soil groups investigated in the Ontario's Agricultural and Soil Health Conservation Strategy (OASHCS) and Comprehensive Assessment of Soil Health (CASH).

2.4 Data visualization and correlation analysis

Regardless of soil textural groups, a principal component analysis (PCA) was conducted in SAS (SAS Institute version 9.4, Cary, NC, USA) on the measured values (without outliers) of the soil health indicators. Based on the scree plots, the first two PCs were selected; the eigenvalues of the soil health indicators within the first two PCs were visualized in two-dimensional plots to understand the interdependence among soil health indicators and better summarize the trends of variation in the dataset. Moreover, Pearson correlation analysis (at P < 0.05) was conducted to further examine the relationship between the soil health indicators.

3.1 Effect of soil texture on soil health indicators

Due to the known effects of texture on soil parameters, each dataset was subdivided into soil textural groups based on three (i.e., the same textural groups as CASH) or five groups (as recommended by OASHCS). Regardless of the number of textural groupings (three or five), our results demonstrated a significant effect of soil texture on all the tested soil health indicators (P < 0.05; Tables 3, 4,  S1 (cjss-2021-0145suppla.docx), and  S2 (cjss-2021-0145suppla.docx)), which was consistent with other jurisdictions (Stott et al. 2013; Fine et al. 2017; Nunes et al. 2021).

Table 3.

Sample size, mean, and standard deviation (SD) for the chemical and biological indicators of soil health sampled from topsoil within three soil textural groups in Ontario.

Table 4.

Sample size, mean, and standard deviation (SD) for the soil chemical and biological indicators of soil health sampled from topsoil within five soil textural groups in Ontario.

For evolved NH₃, evolved CO₂, and SOM concentration, all three textural groupings (i.e., CASH approach) were significantly different from each other (Table 3), whereas for PMN concentration, significant differences were observed between coarse and medium, and coarse and fine textural groups but no differences were detected between medium and fine textural groups (Table 3). In fact, this aligns with CASH, where the PMN soil health indicator was evaluated based on a single scoring function curve because no significant differences were observed in PMN between the three texture groups (Moebius-Clune et al. 2016).

The multiple comparison of means showed that the trends across five texture groups differed, and the trends were inconsistent between the soil health indicators (Table 4). For example, evolved NH₃ mean values for clay and clay loam groups were not significantly different from each other but were significantly different from the loam and sandy loam groups; however, the sandy group was not significantly different from any of the other four soil textural groups. A somewhat similar trend was observed for the PMN indicator, with clay, clay loam, and loam not being significantly different from one another but being significantly different from the sand group (Table 4). In addition, for PMN, the sandy loam group only differed significantly from the loam and sand groups. The difficulty in separating the sands from the other soil texture groups is likely a result of the relatively small sample size for the sandy group (ranging from 17 to 24) for the evolved NH₃, evolved CO₂, and PMN indicators where the sample size is an order of magnitude smaller than the other texture groups. For the SOM indicator, clay and clay loam groups were not significantly different, but loam, sandy loam, and sand groups were all significantly different from one another (Table 4). In the present study, the dataset for SOM was larger (n = 2781) than the other indicators (n = 752–777). Regardless, the SOM concentration was not different between clay (440 samples) and clay loam (593 samples; Tables 4 and  S2 (cjss-2021-0145suppla.docx)), which suggests that there was no need for five texture groups. Lastly, when looking across all four indicators, the five texture groupings recommended by OASHCS generally indicated that the clay and clay loam groups were similar, as were the sandy loam and sand groups. When considering how the textures are grouped (Table 1), mean separation of the five texture grouping is remarkably similar to the three texture groups; hence, scoring functions were developed using three texture groupings (i.e., CASH) and forms the basis of further results and discussion.

All the soil health indicators, except PMN, followed the pattern of a decrease in mean values with an increase in the coarseness of soil, which is consistent with Amsili et al. (2021), Nunes et al. (2018), Fine et al. (2017), and van Es and Karlen (2019). As expected, greater concentration of SOM in fine soils compared to coarse soils reflects the increased potential of soil carbon retention in fine-textured soils (von Lutzow et al. 2006; Samson et al. 2020; Amsili et al. 2021). Likewise, high concentration of evolved CO₂ and NH₃ from fine-textured soils indicates a greater microbial activity with the fineness of soil texture, which is perhaps related to greater soil carbon in fine-textured soils (Amsili et al. 2021). Unlike all other indicators, concentration of PMN was 5.2% greater in medium- than fine-textured soils (Table 3). It is unclear why the PMN trend contrasts with the other indicators but Bechtold and Naiman (2006) observed a higher fraction of mineralizable and labile nitrogen in medium- than fine-textured soils.

3.2 Scoring functions

We developed scoring functions for each soil health indicator using three soil textural groups (Fig. 4). Like CASH, we used a cumulative normal distribution function, based on the mean and standard deviation of the soil health indicator, to provide a soil health score and rating for each indicator evaluated in this study. Our results showed very little overlap in scoring function curves among soil textural groups (Fig. 4), which confirmed that separate scoring functions by each soil textural group should be developed for each soil health indicator. This contrasted with previous CASH research where soil respiration and PMN were not different among soil textural groups; hence, the same scoring functions for all the soil textural groups were used (Fine et al. 2017; Amsili et al. 2021). This discrepancy might be due to methodological differences, such as depth of soil sampling (surface 15 cm vs. topsoil ∼25 cm ( Fig. S1 (cjss-2021-0145suppla.docx))) and incubation time to measure evolved CO₂ (i.e., 4 days with CASH compared to 1 day in our study).

Fig. 4.

Cumulative normal distribution scoring functions developed using the mean and standard deviation of evolved CO₂ (a), evolved NH₃ (b), potentially mineralizable N (c), and soil organic matter (d) sampled from topsoil in Ontario. Points represent the empirical cumulative distribution function while lines represent the cumulative normal distributions.

As expected, the soil health scoring and ratings varied with the soil textural group (Table 5). For instance, evolved NH₃ rating was very low when concentrations ranged from 154 (coarse texture) to 211 mg NH₃–N kg⁻¹ (fine texture; Table 5). For SOM, measured values between 42 (coarse texture) and 51 mg⋅g⁻¹ (fine texture) represented high ratings (Table 5). Very low to low soil health ratings are expected to indicate a constraint in soil functioning and identify the need to make appropriate soil and crop management decisions to improve soil health (Fine et al. 2017); thus, it is important to have soil health scoring and ratings based on soil texture and regionally (i.e., Ontario) specific.

Table 5.

Soil health scores, ratings, and the corresponding chemical and biological soil health indicator values sampled from topsoil within three soil textural groups in Ontario.

3.3 Relationship among indicators

Based on the scree plots, two principal components (PCs) were selected, which explained 68.9% of the total variability in the dataset. About 50% of the variance was accounted for by PC1, where the eigenvectors of the soil health indicators followed the order of evolved NH₃ > PMN > evolved CO₂ > SOM (Table 6). On PC2, which explained 18.6% of the variance, SOM had the greatest eigenvector whereas PMN had the least (Table 6). Typically, each PC is represented by soil physical, chemical, or biological indicators (Fine et al. 2017; Wu and Congreves 2022). In our study, the biological soil health indicator, evolved NH₃, had the greatest influence on soil health assessment in PC1, whereas the chemical soil health indicator SOM was the most important on PC2 (Table 6). Unlike our study, previous findings by Fine et al. (2017) and Zebarth et al. (2019) indicated a high influence of soil respiration (biological) and SOM (chemical) on PC1. A soil health study by Chahal and Van Eerd (2019) found a high influence of soil respiration on PC1 and SOM on PC2. Furthermore, no clear groups among the soil textures were observed between PC1 and PC2 (Fig. 5a). Among the three soil textural groups, the two-dimensional plots revealed a high degree of overlap (Fig. 5a).

Fig. 5.

Principal component (PC) analysis demonstrating the variability in the dataset between soil textural groups (a) and factor groupings (b) of chemical and biological soil health indicators, sampled from topsoil in Ontario, on PC1 and PC2. PMN, potentially mineralizable N; SOM, soil organic matter.

Table 6.

Eigen analysis of chemical and biological soil health indicators sampled from topsoil in Ontario.

Across all soils, biplots of the eigenvectors of soil health indicators on PC1 and PC2 indicated two distinct groups: (i) evolved CO₂, PMN, and evolved NH₃ and (ii) SOM (Fig. 5b). Although SOM and evolved CO₂ fell into the same quadrant in the biplot, evolved CO₂ was more closely associated with evolved NH₃ and PMN than SOM. Similar results were observed by the correlation analysis, where the correlation coefficient of evolved CO₂ with evolved NH₃ (r = 0.35, P < 0.0001) and PMN (r = 0.33, P < 0.0001) was greater than with SOM (r = 0.29, P < 0.0001). Contrary to Marshall et al. (2021), our study results suggest that the amount of readily available soil nitrogen has a greater influence on soil respiration (soil microbial activity) than the amount of substrate (represented by SOM), which is consistent with Zhu et al. (2016). As expected, a positive relationship of evolved NH₃ with PMN (r = 0.44, P < 0.0001) was observed, further confirming PMN as an indicator of soil microbial activity and soil health (Sharifi et al. 2007). From the correlation and PCA, SOM had a positive but weak relationship with evolved NH₃ (r = 0.30, P < 0.0001), evolved CO₂ (r = 0.29, P < 0.0001), and PMN (r = 0.28, P < 0.0001). Although weak but positive association of SOM with PMN, evolved NH₃ and CO₂ confirm the multifaceted role of SOM in influencing mineralization and soil microbial activity (Maestrini et al. 2014; Marshall et al. 2021).

Although simplified soil health assessments are promoted as cost- and time-efficient, comprehensive assessments consisting of a suite of soil health indicators might be informative. For instance, in our study, despite the correlations among the soil health indicators, the degree of orthogonality revealed by PCA indicates an independent effect of the indicators other than evolved NH₃ and PMN. We did not assess any soil physical indicator in this study, which is a limitation. It is expected that adding other indicators to the soil health database might provide additional information related to soil functioning and contribute to improving soil health interpretation. Most soil health assessments sample surface soil to 15 cm depth but we were interested in characterizing topsoil (median depth was approximately 25 cm). One might expect greater concentrations in soil sampled from shallower depth soils (i.e., 15 cm), especially in no-till systems, but this is not always the case (Wu and Congreves 2022). This study is not promoting (nor excluding) any commercially available soil health test. Our goal is to enhance understanding of soil health assessment and to facilitate the process and development of a provincial soil health test.