Ecological insights into soil health according to the genomic traits and environment-wide associations of bacteria in agricultural soils

Abstract

Soil microbiomes are sensitive to current and previous soil conditions, and bacterial 'bioindicators' of biological, physical, and chemical soil properties have considerable potential for soil health assessment. However, the lack of ecological or physiological information for most soil microorganisms limits our ability to interpret the associations of bioindicators and, thus, their utility for guiding management. We identified bioindicators of tillage intensity and twelve soil properties used to rate soil health using a 16S rRNA gene-based survey of farmland across North America. We then inferred the genomic traits of bioindicators and evaluated their environment-wide associations (EWAS) with respect to agricultural management practice, disturbance, and plant associations with 89 studies from agroecosystems. Most bioindicators were either positively correlated with biological properties (e.g., organic matter) or negatively correlated with physical and chemical properties. Higher soil health ratings corresponded with smaller genome size and higher coding density, while lower ratings corresponded with larger genomes and higher rrn copy number. Community-weighted genome size explained most variation in health ratings. EWAS linked prominent bioindicators with the impacts of environmental disturbances. Our findings provide ecological insights into bioindicators of soil properties relevant to soil health management, illustrating the tight coupling of microbiome and soil function.

Fig. 1. An overview of the general relationship between bioindicators of biological, physical, and chemical classes of soil properties.

Among all correlated bioindicators of soil health properties shown in these co-occurrence network diagrams, the majority of positive correlations where with biological ratings (A and B), while the majority of negative correlations where with physical and chemical ratings (C and D). Networks were divided based on whether indicators exhibited positive (A and B) or negative correlations (C and D) with health ratings. In (A) and (C), nodes are colored according to health rating class and, in (B) and (D), according to whether a taxon is represented by a described species. The relationship among bioindicators was visualized in a network to highlight the prevalence of key taxa (see labeling of nodes in B and D); differences in the relationships of bioindicators with soil health classes (in A and C); and the high number of uncultured/unclassified bioindicators (in B and D). Nodes represent bioindicator OTUs aggregated to their lowest resolved taxonomic rank and scaled by the total number of OTUs. Edges represent co-occurrence of indicator OTUs for one or more of the same health rating. Edge weights are scaled by the number of co-occurring OTUs common between nodes. In (A) and (C), nodes were colored based on majority rules according to the number of OTUs representing a given health class. Classes were hyphenated when no majority was achieved.

Fig. 2. Correlations among genomic traits and also total health score.

In (A), correlations were based on average trait scores weighted by the relative abundance of taxa-specific traits values (i.e., community-weighted data) in the soil health data. In (B), the same calculation was made from the genomic database used to assign trait values to taxa. This side-by-side comparison illustrates that the relationships among traits in community-weighted data partially reflected the existing relationships observed in the genomic data (Mantel statistic r = 0.64; p = 0.01). The strength of each Pearson’s correlation corresponds with color intensity as indicated by the scale provided. Circle area corresponds to the inverse of p value with non-significant values indicated by a small, colorless circle.

Fig. 3. A summary of correlations between soil health ratings and community-weighted traits.

Community-weighted genome size explains the most variance in overall health ratings. All Pearson’s r > |0.3| are shaded blue and all significant correlations are shown in bold.

Fig. 4. The relationship between tillage, soil health, and community weighted genomic traits.

Tillage was significantly different according to variation in overall soil health rating, community-weighted rrn copy number, and community-weighted genome size. Significant differences are denoted with asterisk based on t-tests (** p < 0.01; *** p < 0.001).

Fig. 5. Community-weighted genome size explains significant variation in active carbon content.

Points were colored on the basis of active carbon rating in soil health assessment. Differences in the proportion of unclassified taxa assigned traits was not correlated with active carbon rating (r = 0.06, p = 0.1; see Fig. S6).

Fig. 6. Variation in the relative abundance of six key bacterial taxa indicates that the effects of active carbon rating and tillage differ with respect to genome size.

A set of six taxa were selected to represent extremes of genome size from three of the largest (left size, indicated by upward arrow), to three of the smallest (right size, indicated by down arrow). Collectively, these six taxa comprised 14% of all reads and their relative abundance has a strong impact on relationships indicated in Figs. 4 and 5. Active carbon ratings were divided into categories that range from very low (0–20) to very high (80–100). Pairwise statistically significant differences (p < 0.05), according to Tukey HSD, are denoted by lettering.

Fig. 7. Environment-wide associations of the taxa indicated in Fig. 6.

Each plot shows indicator OTUs for (A) plant association (i.e., bulk vs. rhizosphere soil), (B) soil disturbance and (C) management practices. In (B) and (C), indicator values were assigned as positive or negative based on whether a factor was a reference or treatment (e.g., no-till vs. till) with details of designations provided in Table S3. Only OTUs shared among the AgroEcoDB and soil health data were included. Pairwise statistically significant differences (p < 0.05), according to Tukey HSD, are denoted by lettering.