Twenty-two sites were selected across three islands (Oahu, Maui, Molokai) within the main Hawaiian archipelago to cover a diversity of soil management, fertility, and taxonomy to maximize the variance in parameters associated with soil health. Preliminary assessments helped categorize sites into 1) current land management/use, 2) land use history, 3) disturbance level, 4) soil order, and 5) predominant mineralogy. Current land management/uses included protected forests (managed to preserve long-term non-native or native forest, greater than 100 years no disturbance from feral ungulates), unmanaged previously intensive agricultural lands (UPIAL, previously monocrop plantation with no current management system, grasses, shrubs, or forest as dominate cover, and less than 100 years no disturbance), pasture (managed with pasture grasses for rearing livestock), organic croplands (no use of chemical pesticides), and conventional croplands (use of chemical pesticides). Land use history indicated simply whether an intensive plantation history (for Hawaiʻi, this is typically sugarcane or pineapple) was present or absent. Disturbance level was defined categorically as low, medium, and high (Table 1). Soil order was according to final GPS coordinates of sample location and NRCS NCSS taxonomic classification (Web Soil Survey). Predominant mineralogy was assigned using taxonomic classification and a diagnostic key (Supplementary Table S1).

The final compilation included an integration of sites across soil types, land use history, and natural versus agricultural landscapes that is representative of Hawaiʻi (Table 2 ). However, as is reflective of reality, some soil types are more represented in some land uses and some current land uses are more likely to be represented in one past land use history or another. Therefore, we purposefully designed this study as a multivariate approach to identifying parameters indicative of soil function, specifically healthy soil function, for the ecosystems of Hawaiʻi (and other similar tropical/subtropical, and volcanic regions). Then, we narrowed down to key soil health indicators that can be linked to inherent drivers to facilitate the next steps of developing an index of soil health and refining the parameters for specific systems with the goal of assisting landowners, managers, and farmers to improve the health and resilience of their lands.

Three replicate samples were collected from each of the 22 sites. Each sample was comprised of five soil cores taken from the 0–15 cm depth of mineral soil using established sampling methods form the Cornell Soil Health Manual (Doran and Jones, 1996; Moebius-Clune, 2016). Briefly, the organic horizon was removed prior to soil coring 0–15 cm of mineral soil at five locations within a 1 m² quadrant. The soil from five cores was homogenized into one sample in a bucket in the field. Three quadrants that each produced one soil sample for a site were located at least 5 m apart within the site. Thus, 66 soil samples were packaged in a cooler with ice and transported to the lab for analysis. Samples were transported to processing and storage facilities at UH Mānoa and subsets were frozen at −20°C and air-dried (<10% moisture). Additionally, a subset for phospholipid fatty acid testing (see below) was kept chilled, not frozen, and shipped immediately under refrigeration to the analysis facility.

Forty-six parameters classified as biological, chemical, or physical and tied to soil function or process were measured for each of the 66 samples ( Supplementary Table S2 ).

The principle components, based on the key soil health indicators covariance matrix, of each sample location was summarized using principal components analysis (PCA in PC-ORD v.7.0) The response matrix included 46 parameters in 66 soil samples. Uneven distribution of data (i.e., a slight horseshoe shape with clear outliers) in the 2D output suggested the need to transform data, which was confirmed by assessing distribution tables of each variable for non-normality (Peck, 2016). Transformations of log, cube root, and square root were tested on each highly skewed variable since all were either positively or negatively skewed with single peaks. Transformed variables were rerun for skewness and the transformation with the lowest skewness value was selected as the best possible transformation. A new PCA using variables transformed for normality showed an improved graphical display regarding spatial distribution of plots and outlier assessment.

Values from the 46 measured parameters served as the main dataset for PCA ordination, and overlays of supplemental environmental data operated as the second matrix including current land management, historical disturbance, soil order, and soil series. Potential variables were individually removed and re-added to determine their influence on the PCA and assess any significant impact on results, while using the second matrix overlays to identify issues with the distributions of data balance or flag potential errors in data manipulation. To avoid potential bias, variables problematic to overall balance were removed. Supplemental environmental data operated as the second matrix including current land management/use, land use history, disturbance level, soil order, predominant mineralogy, and combinations relevant to hypothesis testing.

The normally distributed soil health indicator matrix was analyzed using a standardized PCA with a correlation cross-products matrix, which produces correlation coefficients among the variables and further standardizes non-comparable response variables. This method provides a broken stick eigenvalue; the broken stick eigenvalue was less than the actual eigenvalue for the first four axes, therefore these are all presented and interpreted to some degree (Peck, 2016). Rnd-Lambda randomization results agreed with the broken stick method, the last useful axis is four with p = .001 and cumulative variance explained at 71.7%.

Multi-response permutation procedure (MRPP) is a non-parametric multivariate test of differences between groups (Peck, 2016). The A statistic describes effect size with respect to how similar within-group samples are compared to outside the group samples. When A = 1, sample units within each group are identical, when A = 0, groups are no more different than expected by chance. We tested all possible models based on the site attributes—agricultural history, current land use, soil order, mineralogical class, disturbance level, and cropland versus non-cropland. The MRPP was run with Euclidean distance measure on the transformed data. Any contrast with less than two sites (6 sample points) was excluded. In the final model, an adjusted p-value was calculated by dividing the model p-value by the number of pairwise contrasts, the adjusted p-value was used to determine whether a pairwise contrast was significant or not.

The list of 46 potential indicators of soil health was reduced to a short list of key indicators that meet multiple criteria for capturing the breadth of soil health as an ecosystem property, reducing multicollinearity with other variables, and practicality for inclusion in a routine soil test. First, potential parameters were removed if r < .50 or > −.50 with axis 1 in the PCA, which left 26 selected for further assessment. A hierarchical ordering correlogram of the untransformed values for remaining 26 selected parameters was performed in R. Highly correlated parameters were reduced further on the basis of practicality (i.e., combination of cost and difficulty). The final list was cross checked to maintain balanced coverage across biological, chemical, and physical properties. Within the constraints of the original sample design, the key parameters of soil health were compared across a subset of mineral and land use classes to assess their utility as indicators of soil health.

Univariate analyses were conducted to first determine the effect of past land use (PIAL versus none) and then assess the effect of current land use (protected forest, pasture, UPIAL, organic cropland, versus conventional cropland) on each of the 11 soil health indicators. A mixed model ANOVA approach was used to assess the main effect of past or present land use with soil mineralogical classes as the random effect [lmer function in the lme4 package (Bates et al., 2015)]. General linear hypothesis testing [glht function in multcomp package (Lenth, 2020)] was performed to compare group means (Tukey-adjusted). For a subset of the data (including UPIAL, organic cropland, and conventional cropland in LAC and HAC soils), mixed factorial ANOVA was performed to examine the interactive effect of soil mineralogical class and current land use on each soil health indicator with farm/location as the random effect (lmer function in the lme4 package). Group differences were assessed by Tukey multiple comparisons of least square means (lsmeans function in emmeans package). For all post hoc multiple comparisons with the Tukey test, letter groupings were assigned with the cld function (multcomp package (Piepho, 2004)).

Soil was collected from a range of land uses (identified a priori as protected native and non-native forest, managed pasture, unmanaged previously intensive agricultural lands (UPIAL), organic cropland, and conventional cropland) ( Figure 1 ) across a range of soil orders and clay mineralogy (Oxisol, Mollisol, Andisol, Inceptisol, and Vertisol) on three islands. When possible, pairs or triplets of sites were obtained on the same, or related, soil series but different land use. Forty-six parameters across biological, chemical, and physical soil properties were measured and analyzed with a multivariate approach to 1) test for the predominant drivers of inherent (i.e., clay mineralogy, texture) differences in soil on a heterogenous landscape and 2) deduce a key set of indicators that may represent soil health as an ecosystem property. Four significant principle components analysis (PCA) axes cumulatively explained 71.7% of the variance within the soil health indicator dataset. The two dominant axes explained 43.0% and 12.3% of variance, followed by the next two that explained a further 9.0% and 7.4%. Many of the parameters across biological, chemical and physical soil properties, strongly correlated (i.e., r > .5) to the positive or negative side of axis 1 (Supplementary Table S3).

Land use, specifically the legacy of intensive cultivation, predominated over soil type to influence the indicators of soil health. Visualization of axis 1 and 2 of the PCA showed that regardless of current land use and soil type, sites with a history of long-term intensive cultivation clustered independently from other land uses within forest and pasture classifications (Figure 2). Sites without a plantation agricultural history (i.e., 80+ years of sugarcane or pineapple) separated out from those sites with intensive land use history along axis one. The top five strongest (i.e., r > .95) parameters driving the sites with no intensive agricultural history toward the negative side of axis 1 were high Gram negative bacteria, total phospholipid fatty acids (PLFA), organic carbon (OC) concentration, total N concentration, and actinomycetes. The negative side of axis 1 was also driven by high values of many additional biological parameters not listed as well (Supplementary Table S3). In contrast, those sites with a plantation history had low concentration of those parameters negatively related to axis one, and high bulk density (BD), dissolved OC (DOC) to dissolved organic N (DON) ratio, actinomycetes to bacteria ratio, clay concentration, and crystalline Fe oxides. Axis 2 did not provide a clear separation among the past or current land use classifications.

To a certain degree, some current land uses correspond to areas with an intensive agricultural past and others to areas without due to land availability and suitability. For example, the sampling sites classified as conventional cropland and UPIAL (by its own definition) reside exclusively within the cluster defined by an intensive agricultural history. Likewise, all the sampled protected forests occurred in areas without the agricultural past. However, organic croplands and pastures resided within both types of areas with and without the intensive agricultural past. Visualization of axes two and three showed no separation among land use history or current status (Supplementary Figure S1).

Soil type, as defined by soil order and more broadly by mineralogical class (i.e., high activity clays, low activity clays, poorly and non-crystalline minerals, and histic) was not strongly associated with axis 1 (Figure 3). Both high and low activity clays aligned with the positive side of axis 1 associated with high BD, DOC:DON, clay and crystalline minerals. However, high activity clays included Vertisols and Mollisols, and associated with the positive side of axis two, driven by high extractable Ca²⁺, K⁺, P, and Na²⁺, and sand concentration. Low activity clays included Ultisols and Oxisols, and associated more with the negative side of axis 2, driven by high ratio of pyrophosphate to hydroxylamine extractable Al, crystalline Fe oxides, concentration of mega size class water stable aggregates, soil hardness at the surface layers, and BD. Histic soils fell on the negative side of axis 1 which correlated highly with total OC, total N, Actino:Bacteria, and AM Fungi. Visualization of axis 2 and 3 helped further separate out the histic and poorly and non-crystalline mineral (PNCM) groups (Supplementary Figure S1). Particularly with respect to the Andisols and andic Inceptisols having high concentration of poorly and non-crystalline minerals (Al_H+0.5Fe_H), PNCM separation from other soils was driven by silt and sand concentration, fungi to bacteria ratio, and concentration of “mega” size class water stable aggregates.

Multiple multi-response permutation procedure (MRPP), a non-parametric multivariate test of differences between groups (Peck, 2016), models were run among the categorical classifications and selected combinations for hypothesis testing to identify indicators of soil health (Supplementary Table S4). The combination of agricultural history and disturbance level was the second most significant contrast (adjusted p = .008, A value = .2636). The significant pairwise comparisons indicate that those sites classified as PIAL-medium are different than those PIAL-high and none-high; none-low is different than PIAL-medium and none-high. These results indicate first that sites with a plantation history and medium disturbance classification were more like one another than to those with a high disturbance level (i.e., currently in intensive cultivation), regardless of whether there was plantation past land use or not. Second, undisturbed sites with no plantation history were more like one another than to those with a plantation history and medium current disturbance classification or those with no plantation history but currently under intensive practices. The combination of agricultural history and current land use was also among the significant contrasts tested and yielded similar results to the interaction with plantation history and disturbance level. In the pairwise comparisons, those sites classified as PIAL and currently unmanaged (UPIAL) are different from PIAL and currently in conventional cropland, and from sites with no plantation history and currently in organic cropland, pasture, or protected forest. The simpler current land use model had the third highest A value (.2583) and 10 pairwise contrasts that separated from one another along the dominant PCA axis.

The most distinct group was the UPIAL land use class, which was significantly different from all other land uses (Figure 4). Pasture sample units spanned axis 1 yet were significantly different than UPIAL and organic croplands in multivariate space indicating a latent interaction with mineralogy unable to be explored further within the constraints of our dataset. Pastures were not different from protected forests and the two groups showed a lot of overlap in the 2-dimensional visualization of axis 1 and 2, also signaling a potential influence of similar mineralogy. In the case of conventional agriculture, the sample units were so dispersed within the group, that no differences emerged between it and the other current land uses, except for UPIAL. However, organic croplands were similar enough to one another to be significantly different than UPIAL, pastures, and protected forests.

The combination of cropland and mineral classifications also was among the significant contrasts tested (Supplementary Table S4) and provided additional insight into the nature of the interaction of minerals with land use by sub-setting the dataset to reduce the complexity of the five current land use classes to simply cropland (organic and conventional) and not cropland (UPIAL, pasture and protected forest). The non-cropland PNCM sites were different from non-cropland HAC and LAC. For both HAC and LAC, those sites in croplands were different from those not in cropland. Within the constraints of the dataset, which has greater coverage of HAC and LAC across the cropland/not cropland designations, being in cropland affected soil health indicators for both HAC and LAC. Additionally, HAC in cropland was different from LAC not in cropland, and vice versa.

Eleven dynamic soil parameters emerged from a multi-step dimension reduction process as indicators of soil health across diverse land uses, histories, and soil types. First, at its foundation PCA is a dimension reduction approach, and 26 parameters correlated strongly ( ≥ .5 or ≤ −.5) with axis one. A correlation matrix of the untransformed values showed covariance among many of those 26 parameters (Supplementary Figure S2). From this covarying block, consideration of the practicality of the parameter’s inclusion in a rapid, accessible soil health index (i.e., cost and difficulty) and coverage of biological, chemical, and physical parameters further reduced the list to CO₂ burst, HWEC, PMN, total OC %, and WHC. Among the parameters not included in the block: ß-glucosidase, ß-glucosaminidase, and water stable mega aggregates (mega WSA), the DOC to DON ratio, actinomycetes to bacteria ratio, and BD remained. Because actinomycetes to bacteria ratio is not feasible in a rapid soil health test, it was removed from the final list.

With the criteria of strong relationship to PCA axis 1, non-covariance, practicality, and inclusion of biological, chemical, and physical parameters, 11 parameters emerged as potential indicators of a soil health gradient across the soils and ecosystems in Hawaiʻi (Table 3). Summary values show the range, mean, and median of each indicator across the dataset (Table 4). For contextualization, all indicators except pH and BD are greater in soils without an intensive plantation history than in those with. For several of these indicators (total OC %, CO₂ burst, ß-glucosidase, ß-glucosaminidase, PMN, and HWEC), protected forests are greater than conventional croplands, while the DOC to DON ratio was lower. Similar differences were also present for pasture compared to conventional cropping except for a few (PMN and HWEC). Among those indicators with significant contrasts, total OC %, CO₂ burst, PMN, and HWEC are the same for UPIAL versus conventional cropland, but ß-glucosidase and ß-glucosaminidase are greater while the DOC to DON ratio was lower for UPIAL than conventional cropland. There were fewer differences between organic and conventional cropping, but organic cropping had greater total OC%, and ß-glucosaminidase than conventional.

For the subset of soils with adequate representation across mineralogy (HAC and LAC) and current land use (UPIAL, organic cropland, and conventional cropland) there were significant interactive effects among several soil health indicators, including total OC %, CO₂ burst, ß-glucosidase and ß-glucosaminidase, PMN, HWEC, and mega-WSA (Figure 5). For LAC soils, many indicators were consistently lower for conventional than organic croplands, including OC %, CO₂ burst, ß-glucosidase, ß-glucosaminidase, PMN, and mega-WSA. In general, UPIAL tended to have lower values than organic management, but greater than conventional cropland, which was especially apparent for total OC %. In contrast to LAC soils, very few significant effects of land use on HAC soils were detected, and only PMN and HWEC was less in conventional than in organic croplands while mega-WSA were greater. Of those indicators without significant interactions, mineralogy was significant for the DOC to DON ratio and WHC, and in both cases HAC was greater than LAC (Figure 6). Current land use was significant for the DOC to DON ratio (conventional > UPIAL), WHC (organic > UPIAL and conventional croplands), and BD (UPIAL > organic croplands) (Figure 7). Mineralogy and land use had no detectable effect on soil pH.

Soil health conceptualized as a measurable ecosystem property based upon key indicators provides a powerful tool for monitoring progress in restoration projects or implementation of best management practices to promote sustainable agroecosystems. A new paradigm of soil organic matter dynamics, which is central to soil health, is driving the development of new compartmental models tied to measurable soil parameters (Blankinship et al., 2018). Soil organisms, particularly microbial communities that are proximately responsible for the flow of nutrients, C, and energy in the soil ecosystem, rely on accessible organic matter for metabolic substrate. Many of the emergent process-based ecosystem models are microbial models [e.g., (Wieder et al., 2015; Robertson et al., 2018)] and hold promise to improve both decision support tools and earth system projections.

In policy and programs intended to incentivize maintaining and aggrading soil health across multifunctional landscapes and diverse stakeholders, expectations must be gauged accordingly. Long-term, intensive monocrop agriculture, which in Hawaiʻi was predominantly pineapple and sugarcane plantations established post-Western contact, leaves a detrimental legacy on indicators of soil health. The adverse effects on soil biological properties and microbial communities persists following both abandonment and land use/management change to practices consistent with soil health management principles (e.g., perennial grasses or crops, organic matter inputs, and no or reduced tillage). Especially because the legacy of intensive cultivation history may carry-over into the success of land-based initiatives now and into the future, it is important to understand the resultant differences in baseline conditions as well as the limitations to improvements in soil health indicators when building decision support tools and programs.

Soil health metrics interwoven with process-based ecosystem models that underpin decision support tools used by policy makers may also assist in accessing aid and improved economic outcomes that are critical to success in overcoming adoption barriers. Ecosystem functions such as greenhouse gas emission, C storage, nutrient transformation, biomass production, and regulation of hazards and extreme events link directly to key services contributing to human wellbeing (Adhikari and Hartemink, 2016; Saco et al., 2021). Conventional soil organic matter models (e.g., RothC, DNDC, EPIC, and CENTURY) embedded within established decision support tools that assist land managers and policy makers alike contain some of these ecosystem functions. However, their ability to simulate (sub) tropical or volcanic soils which have very different properties from the temperate or continental soils for which they were developed currently limits their usefulness. Alternatively, microbial models such as (MEMS v1.0) developed from the Microbial Efficiency-Matrix Stabilization framework incorporate measurable parameters that constrain C pools sizes and modulate fluxes (Robertson et al., 2018). Measurable pools and rate modifiers in MEMS v1.0 include microbial biomass and turnover, dissolved organic matter, sorption/desorption dynamics, and exoenzyme activity. The overlap between these parameters and the indicators of soil health suggest that compartmental models designed to simulate soil organic SOM matter dynamics and nutrient and GHG fluxes may benefit from the integration of soil health measurements into their initialization and projections for (sub) tropical soils.

In the process of aggrading soil health, landscapes regain resilience through improved soil functions. For Hawaiʻi and other (sub)tropical and volcanic regions, land-based management relating to conservation, biocultural restoration, climate action for mitigation/adaptation, and increased local food production strive for their specific goals. But, these practices also contribute more broadly to sustainability when viewed through the lens of improved soil health and the associated expansive network of co-benefits and regulation services. Understanding, representing, and projecting outcomes associated with soil health is critical to incorporating their full value into complex watershed-based management and interdisciplinary social-ecological forecasting that link directly to building resilient, climate ready landscapes and communities.

Key drivers for the inherent differences in the soil health indicators, including land use history, current land use practices, and mineralogy, must be understood and integrated into the development of a soil health index. Any index should set minimum and maximum benchmarks and weigh parameters according to equitable standards. Therefore, the state of each driver (e.g., timeline of intensive use history, time since implementation of current land use and management, and predominant mineralogy) must be ascertained and recorded in databases designed for syntheses of soil health into the future.

The legacy of a plantation history is a strong driver of differences in soil health indicators, but greater complexity associated with current land use, management, and soil type also is present and important to understand while developing a robust soil health index for the (sub)tropics. For example, the level of disturbance in current management practices and, outside of croplands, mineralogy both also affected soil health indicators. Results suggest that soil health may differ inherently for high versus low activity clays and whether a system is cultivated intensively for food production (cropland) under conventional versus organic management affects soil health regardless of mineralogy.

Importantly, the soil health indicators of pastures were not different from protected forests, suggesting that well-managed grazing lands may be similar in its soil health status as protected forests. However, the most distinct group was the unmanaged, previously intensive agriculture lands (UPIAL land use class), which was significantly different from all other land uses. Unmanaged abandoned agricultural lands are more similar to each other than to sites that remain in intensive cultivation. But, they are also more similar to each other than to sites without plantation history and currently are in organic croplands, pasture, or protected forest. Upon further inspection, the univariate analysis suggests that, while the relationship between organic versus conventional cropland was largely consistent across soil type for most soil health indicators, abandoned cropland was more variable. For UPIAL, in some cases, soil health indicators fell in between organic and conventional croplands, while sometimes aligning more closely to conventional or organic for other indicators. This finding further highlights the imprint that intensive agriculture may have on the health of a soil and demonstrates the constraints to rebuilding soil health upon the cessation of soil disturbance without proactive management strategies.

A suite of readily measured indicators emerged out of a first-principle approach to determining a holistic indication of soil health across a range of soils and systems in Hawaiʻi encompassing much of the diversity in the tropics and subtropics. These indicators integrate biological, chemical, and physical properties with key functions associated with soil C and nutrient cycling, water relations, and generally, the provisions of a soil environment conducive to a diverse soil organismal community. These indicators are consistent with current measures and assessments of soil health (Rinot et al., 2019), but developed with a more organic and equitable process, without carry over of ingrained bias (Lehmann et al., 2020). In the development of a soil health index, indicators may be weighted differentially for systems. Further, cropping systems should be paired with additional fertility testing and nutrient management for optimal environmental and yield outcome.

Every land use may improve its soil health status within a reasonable range of expectations when considering land use history, current land use, and mineralogy. Land use history had a predominate effect on all but two of the 11 soil health indicators significant. We also detected greater complexity in the interactions between land use and mineralogy, suggesting that we must capture such nuances when we interpret changes in soil health indicators. Therefore, our results provide robust evidence for the contextual inclusion of land use history, current land use, and mineralogy as we develop a soil health index.

The measurement of soil health as a dynamic ecosystem property is only possible by properly identifying the right suite of indicators specific to a region and metering that measurement to appropriate benchmarks for a system defined by past land use, current land use, and mineralogy. Moving forward, providing a soil health index of (sub)tropical and volcanic soils will help to assist currently underserved producers and land managers improve the health and productivity of their lands and simultaneously reap the co-benefits of a healthier environment and society. Within this framework, fair and equitable programs can be established to improve economic outcomes as well as C neutrality goals.