Temporal Dynamics of Bacterial Communities along a Gradient of Disturbance in a U.S. Southern Plains Agroecosystem.

Land conversion for intensive agriculture produces unfavorable changes to soil ecosystems, causing global concern. Soil bacterial communities mediate essential terrestrial ecosystem processes, making it imperative to understand their responses to agricultural perturbations. Here, we used high-throughput sequencing coupled with a functional gene array to study temporal dynamics of soil bacterial communities over 1 year under different disturbance intensities across a U.S. Southern Plains agroecosystem, including tallgrass prairie, Old World bluestem pasture, no-tillage (NT) canola, and conventional tillage (CT) wheat. Land use had the greatest impact on bacterial taxonomic diversity, whereas sampling time and its interaction with land use were central to functional diversity differences. The main drivers of taxonomic diversity were tillage > sampling time > temperature, while all measured factors explained similar amounts of variations in functional diversity. Temporal differences had the strongest correlation with total nitrogen > rainfall > nitrate. Within land uses, community variations for CT wheat were attributed to nitrogen levels, whereas soil organic matter and soil water content explained community variations for NT canola. In comparison, all measured factors contributed almost equally to variations in grassland bacterial communities. Finally, functional diversity had a stronger relationship with taxonomic diversity for CT wheat compared to phylogenetic diversity in the prairie. These findings reinforce that tillage management has the greatest impact on bacterial community diversity, with sampling time also critical. Hence, our study highlights the importance of the interaction between temporal dynamics and land use in influencing soil microbiomes, providing support for reducing agricultural disturbance to conserve soil biodiversity. IMPORTANCE Agricultural sustainability relies on healthy soils and microbial diversity. Agricultural management alters soil conditions and further influences the temporal dynamics of soil microbial communities essential to ecosystem functions, including organic matter dynamics, nutrient cycling, and plant nutrient availability. Yet, the responses to agricultural management are also dependent on soil type and climatic region, emphasizing the importance of assessing sustainability at local scales. To evaluate the impact of agricultural management practices, we examined bacterial communities across a management disturbance gradient over 1 year in a U.S. Southern Plains agroecosystem and determined that intensive management disturbance and sampling time critically impacted bacterial structural diversity, while their interactive effect influenced functional diversity and other soil health indicators. Overall, this study provides insights into how reducing soil disturbance can positively impact microbial community diversity and soil properties in the U.S. Southern Plains.

INTRODUCTION

Rising human populations have resulted in the need for increased land conversion to heavily managed environments for greater food production. Yet, land use change represents one of the largest perturbations to soil ecosystems, significantly impacting both aboveground and belowground communities (1, 2). Whole ecosystem diversity is generally diminished when natural land is converted to agricultural systems, with lasting negative effects on soil health (3). In general, agricultural land use type regulates microbial diversity, plant diversity, and soil physicochemical properties (4–7). The effect of land use on microbial communities has become increasingly important since microbes represent the bulk of biodiversity in terrestrial ecosystems, perform essential ecosystem functions, and are fundamental to ecosystem stability (8–10).

While it has been established that changes in land use shift microbial community structure and diversity, there has been a renewed focus on observing these communities under a gradient of disturbance intensities due to the quickly growing need for sustainable agricultural practices. Different intensities of soil disturbance create unique environments that support microbes with those specific environmental requirements (11). Although terrestrial microbial studies over large spatial scales (6, 7, 12, 13) have demonstrated which soil and environmental factors are important for shaping microbial distribution patterns, they are unable to pinpoint the dynamics required to manage microbial communities at the local level. Agricultural management practices also vary locally, with inputs such as tillage, pesticide and fertilizer use, crop rotation, and residue incorporation directly altering soil microbial biomass (14, 15) and community composition (2, 16). This is critical as there is no ideal community type (17), soil type, or set of soil characteristics (18, 19) when trying to define a functional soil system. By directing attention to gradients of disturbance in a range of land uses commonly found in agroecosystems, local variation can be captured in soil and environmental properties, management type, and plant diversity, which may give insight into the complex dynamics shaping soil communities (11).

Patterns of variability between land use with increasing management disturbance have been studied extensively at single time points, but much less is known about the extent to which land use under a gradient of disturbance intensities interacts with temporal dynamics in altering soil bacterial communities. As seasons transition, variations occur in environmental factors such as solar radiation, temperature, and precipitation, all of which can affect microbial community structure and functions (20–23). Several studies investigating soil microbial community changes in relation to temporal variability have observed community differences in a range of time scales, many of which are associated with shifting environmental conditions (21, 22, 24). These variations in environmental conditions and community structure are often related to land management practices (16, 23, 25) and temporal changes in plant growth and development (26, 27). Specifically, plant growth alters rhizodeposition, promoting microbial activity (28) and modifying community composition by enriching specific microorganisms (29). Expanding on spatiotemporal studies that are specific to the local land use, plant community, and soil conditions are critically needed.

Functional diversity of the soil microbial community is equally important as compositional diversity when examining overall ecosystem diversity. Typically, a high structural and functional microbial diversity is thought to be fundamental to soil health, function, and sustainability by providing functional redundancy critical for ecosystem stability in the presence of stress and disturbance (9, 30–32). The use of functional gene arrays (FGAs) or GeoChip has provided a way to examine relationships between microbial community structure and function by focusing on genes important to microbial processes like biogeochemical cycling and stress responses (33–39). FGAs allow for a thorough analysis of essential ecological questions, especially those concerned with microbial community responses to disturbances (35, 40–44), including soil microbial community responses to land use, land management, and temporal changes (45–47). However, it remains unclear how the functional capabilities of soil microbial communities change under a gradient of disturbance intensities and seasons.

To investigate the effect of land use with increasing management disturbance and season on the temporal dynamics of soil bacterial communities and its underlying mechanisms, we conducted a 12-month field study in agroecosystem land uses with a gradient of disturbance within the U.S. Southern Plains agroecosystem. The agricultural sites included two perennial grasslands and two annual croplands: a native tallgrass prairie (TGP), Old World bluestem (OWB) pasture, no-tillage (NT) canola (Brassica napus L.) field, and conventional tillage (CT) winter wheat (Triticum aestivum L.) field. In this study, we focused on the following questions. (i) Do land use with various types of management disturbance and season shape soil properties? (ii) How do land use and seasonal temporal dynamics interact to influence bacterial community diversity? (iii) What roles do soil and environmental properties play in influencing bacterial community diversity between seasons and under increasing management disturbance? We predicted that soil bacterial diversity would decrease with increasing management disturbance, while seasonal differences would become more discernible with increasing management disturbance. Our results revealed that land use drove differences in taxonomic diversity, while sampling time and its interaction with land use influenced functional gene diversity, and that the biotic and abiotic factors shaping bacterial community diversity also differed spatiotemporally with importance varying with management intensity.

RESULTS

Changes in soil and environmental variables across land use and management gradient.

Over the 1-year sampling period, all the measured soil properties were significantly affected (P < 0.001) by land use (Fig. 1a). Only properties influenced by climate (soil water content [SWC]) and management practices (topsoil nitrate [TopN]) significantly differed by season (P < 0.05). Notably, soil organic matter (OM) and soil total nitrogen (TN) decreased as management disturbance increased. The OM differed by season only for CT wheat. SWC significantly decreased (P < 0.05) in CT wheat, with the lowest SWC of all sites being observed in CT wheat during January, whereas all other land uses had lows in SWC during summer months (Fig. 1b). Minimum daily air temperatures occurred in December 2016 and January 2017, and maximum air temperatures occurred in August 2016 and July 2017. The greatest monthly rainfall was recorded in April 2017 (227 mm) and the lowest monthly rainfall in November 2016 (15 mm). Both croplands had higher averages of TopN that significantly differed from those of the grasslands less management disturbed (P < 0.05). Elevated levels of TopN were present during summer and fall in both croplands. Ammonium (NH4) was significantly (P < 0.05) lower in NT canola than other land uses.

FIG 1

Soil chemistry within each land use type across a 1-year sampling period. (a) Soil chemistry averages for each land use across 1 year. (b) Factors that significantly differed by season across the whole land use gradient shown by season within land uses. Error bars represent the standard deviations. Letters represent significant differences of P < 0.05 between pairwise land use comparisons or seasons within land use. The same letter indicates no significant difference.

Impact of land use and seasonality on soil bacterial communities.

To determine the effect of land use and season over the sampling period, α-diversity was calculated for the bacterial communities (Table 1). For all land uses, seasonal variability had a greater impact on α-diversity than land use, with all indices significantly different between seasons (P < 0.001). Overall, bacterial richness was lower in summer and fall than winter and spring, but no single land use had a more diverse or rich community throughout the sampling period. Shannon diversity was the lowest in the fall across all land use types, and fall was significantly different from other seasons (P < 0.05). The two fields that differed the most as far as management disturbance, TGP and CT wheat, were compared separately to see if land use differences were observed when focusing on the most different land use types and levels of management disturbance (Table 1); interestingly, significant differences were still driven by season.

TABLE 1

Bacterial community structural and functional differences in α-diversity based on land use and sampling time

Alpha diversity	Effect	16S				GeoChip
		All fields		TGP and CT		TGP and CT
		F	P	F	P	F	P
Chao 1	Field	0.517	0.672	0.896	0.350
	Season	12.05	<0.001	5.196	0.005
Observed OTUs	Field	0.573	0.635	1.401	0.244
	Season	11.16	<0.001	3.663	0.022
Pielou	Field	1.891	0.138	1.738	0.195	0.410	0.529
	Season	3.494	0.020	0.280	0.839	1.790	0.192
Shannon	Field	0.426	0.735	0.415	0.523	0.939	0.343
	Season	9.129	<0.001	1.970	0.137	1.772	0.195

Correlations based on analysis of variance. Season for GeoChip represents sampling time since only three times were used across the data set. TGP, native tallgrass prairie; CT, conventionally tilled winter wheat.

The effect of land use and season on the β-diversity of soil bacterial communities was examined using principal-coordinate analysis (PCoA) based on the Bray-Curtis distance metric. The bacterial community structure of soils was visibly separated by land use, with CT wheat generally isolated from the other land uses (Fig. 2). NT canola and OWB pasture were the most similar in community structure. Each field had observable temporal differences in community structure, and the visible temporal differences increased with increasing management disturbance. Permutational multivariate analysis (PERMANOVA) supported the PCoA plot (Table 2), indicating that the structure of the bacterial community was significantly shaped by land use (P = 0.001, R2 = 0.2949) and season (P = 0.001, R2 = 0.1067), but the effect of season was not as strong as that of land use. When comparing just the TGP and CT wheat, the same significant differences were observed, but the effect of land use (P < 0.05, R2 = 0.3614) on the bacterial community structure was even greater.

FIG 2

Principal-coordinate analysis of Bray-Curtis dissimilarity for soil bacterial communities (16S rRNA gene) showing the differences in four fields along a management disturbance gradient over a 1-year sampling period. (a) Differences in community structure between land uses. Land uses include conventionally tilled (CT) wheat, no-till (NT) canola, Old World bluestem (OWB) pasture, and native tallgrass prairie (TGP). (b) Differences in community structure separated by season.

TABLE 2

Effect of land use and season on bacterial community structures

Distance metric	Effect	16S				GeoChip
		All fields		TGP and CT		TGP and CT
		R 2	P	R 2	P	R 2	P
Bray-Curtis	Field	0.2949	0.001	0.3614	0.001	0.0842	0.005
	Season	0.1067	0.001	0.1140	0.014	0.1363	0.034
	Field × Season	0.1102	0.104	0.0783	0.112	0.1049	0.011
Jaccard	Field	0.2057	0.001	0.2394	0.001	0.0798	0.004
	Season	0.0876	0.001	0.1125	0.018	0.1299	0.042
	Field × Season	0.1303	0.032	0.0928	0.083	0.1054	0.014

The 16S permutational multivariate analysis of variance (adonis) model was set up as dissimilarity ∼ field × season. 16S analysis was done by including 4 fields: TGP, OWB pasture, NT, and CT. It was also performed using only prairie and CT since GeoChip included only those two fields. The GeoChip permutational multivariate analysis of variance (adonis) model was set up as dissimilarity ∼ field × season + block with permutation constrained by block to deal with the effect of data on multiple arrays.

Responses of soil bacterial community composition across the land use gradient.

The soil bacterial community was dominated on average by several bacterial taxa across all land use types commonly found in soils, including Actinobacteria (20.16% to 24.67%), Proteobacteria (21.20% to 24.36%), Firmicutes (9.68% to 18.03%), Acidobacteria (9.20% to 13.53%), and Chloroflexi (4.23 to 9.90%) accounting for over 75% of the relative abundance of each system (see Fig. S2 and Table S1 in the supplemental material). The lower-relative-abundance phyla were comprised of Bacteroidetes, Gemmatimonadetes, Planctomycetes, and Verrucomicrobia, all of which made up at least 1% of the bacterial community in all land uses over the sampling period (Fig. S2 and Table S1).

Relative abundance of dominant bacterial phyla in soils separated by land use and season. Values are based on proportional frequencies of DNA sequences classified at the phylum level. Fields include conventionally tilled (CT) wheat, no-till (NT) canola, Old World bluestem (OWB) pasture, and native tallgrass prairie (TGP). Download FIG S2, TIF file, 2.0 MB.

Average relative abundance (%) of abundant phyla in soil from conventionally tilled (CT) winter wheat, no-till (NT) canola, Old World bluestem (OWB) pasture, and native tallgrass prairie (TGP) systems from August 2016 to July 2017. Download Table S1, DOCX file, 0.02 MB.

Significant differences in the relative abundance of the main bacterial taxa between types of land use were evident at the phylum level (P = 0.001). The greatest number of significant differences of relatively abundant phyla were between CT wheat and TGP communities, while NT canola and OWB pasture had the smallest number of significant community differences (detailed in Table S2 in the supplemental material). Within individual land uses, the phylum relative abundance of the bacterial communities differed significantly by season (P < 0.01), except for OWB pasture (P = 0.066). Markedly, significant changes in many phyla across seasons were unique to specific land use types (detailed in Table S3 in the supplemental material). At a lower taxonomic level, roughly 20% of the operational taxonomic units (OTUs) in each land use were not present in any of the other land uses, with the smallest percentage of shared OTUs between CT wheat and the TGP (see Fig. S3a in the supplemental material). Within all land uses, the greatest percentage of unique OTUs was observed during the spring (Fig. S3b to e), corresponding to warming air temperatures, rainfall peaks, and resuming plant growth. CT wheat had the most shared OTUs during the fall and spring, which were the wheat growing seasons. TGP also had the most shared OTUs during the peak growing season for warm grasses during summer and fall. In general, as the amount of management disturbance increased between land uses, the bacterial communities became increasingly different at the phylum and OTU levels.

Unique OTUs across different land uses (a) and seasons (b to e). The percentage of unique OTUs found in each season is presented using a Venn diagram. Download FIG S3, TIF file, 2.2 MB.

Pairwise comparison P values for relative abundance of phyla significantly different between land use types from August 2016 to July 2017. P values were false-dicovery rate (FDR) corrected. Land use types include conventionally tilled (CT) winter wheat, no-till (NT) canola, Old World bluestem (OWB) pasture, and native tallgrass prairie (TGP). All phyla used in pairwise comparisons were previously found to be significant by field using PERMANOVA with 999 permutations and the Bray-Curtis dissimilarity metric. Download Table S2, DOCX file, 0.02 MB.

F values of phylum relative abundance significantly different between seasons. The last timepoint in July was left out of the analysis because it stood alone for the summer 2017 season. All phyla were individually compared across seasons using PERMANOVA with 999 permutations and the Bray-Curtis dissimilarity metric. Asterisks indicate significant differences: *, P < 0.05; **, P < 0.01; *** P < 0.001. Italicized values indicate a P value of 0.05 < P < 0.1. Download Table S3, DOCX file, 0.02 MB.

Effect of soil and environmental factors on soil bacterial diversity.

While α-diversity indices were not significantly different across the land management disturbance gradient, they were affected by local soil properties that differed between land uses. The richness and diversity of the CT wheat bacterial community significantly (P < 0.05) decreased when there were high levels of TopN present. Similarly, the diversity and evenness significantly (P < 0.05) decreased in the presence of elevated TopN in NT canola. Only the richness in the TGP was significantly influenced by SWC. There were no detectable relationships between soil variables and α-diversity for OWB pasture.

The influence of management, soil, and environmental factors on β-diversity was determined using the Mantel test and multiple regression on a distance matrix (MRM), with correlations to individual taxa within land uses shown in Table S4 in the supplemental material. Overall, bacterial community differences were significantly driven by tillage > SWC > sampling time > minimum air temperature according to the Mantel test (Fig. 3a). Between fields (spatially), soil factors shaped by land use were significantly important, including OM > TN > SWC. In comparison, sampling time had the greatest correlation to temporal community differences as well as OM and TN. Since potentially significant correlations among factors could exist, MRM was further used to determine the contributions of different environmental factors to shaping bacterial community structure. In general, differences based on MRM were similar to those of the Mantel test, with tillage > sampling time > minimum air temperature being significant (Fig. 3b). Notably, spatially related factors had a stronger impact (R2 = 0.62) on the bacterial community than temporally related factors (R2 = 0.46). The same soil factors were key to spatial difference based on the MRM and Mantel test (Fig. 3c), with average air temperature also having moderate importance. Sampling time again had the strongest relationship with temporal differences (Fig. 3d). The factors average rainfall > TN > TopN were also key to temporal community dynamics.

FIG 3

Effects of soil and environmental factors on soil bacterial community structure. (a) Correlations on overall, spatial, and temporal differences based on Mantel test. (b) Multiple regression on distance matrix (MRM) on overall community structure. (c) MRM for spatial differences in community structure. (d) MRM for temporal differences in community structure. Bars with diagonal lines represent negative regression coefficients. Gdist, geographical distance between sampling sites. Significance is expressed as follows: ***, P < 0.001; **, P < 0.01; and *, P < 0.05.

Spearman correlations between relatively abundant phyla and soil properties within land use types. Correlation coefficients with P < 0.05 are indicated in boldface, coefficients with 0.05 < P < 0.1 are indicated in italics, and data were omitted where the P value is >0.1. Download Table S4, DOCX file, 0.02 MB.

Each land use was then examined separately to determine if soil and environmental variables contributed equally to the variation in bacterial community structure. Based on canonical correspondence analysis (CCA) (see Fig. S1 in the supplemental material), variance partitioning was used to determine if climate variables, nitrogen measurements, or other soil properties explained the most variations in community structure (Fig. 4). Nitrogen measurements had the largest impact on the CT wheat bacterial community and interacted with the other soil properties and climate variables (Fig. 4a). The majority of the variation of bacterial communities was explained by the soil properties for the NT canola site (Fig. 4b), with nitrogen measurements and climate variables having a minor interaction. For OWB pasture, the variations explained by all groups were similar, with soil properties and nitrogen measurements explaining almost the same amount of variation (Fig. 4c). The distribution of the variation explained for TGP was similar to what was observed at the OWB site (Fig. 4d), with all sets of variables having a relatively equal impact on the bacterial community structure. The greatest interaction of all variables was observed for the TGP. Therefore, management and sampling time had the greatest impact on diversity, and management greatly impacted the importance of different soil and climatic factors in shaping bacterial communities within fields.

FIG 4

Variation partition analysis (VPA) of bacterial community structure explained by soil properties, nitrogen measurements, and climate variables for each land use. Variable groupings include soil variables (SWC and OM), nitrogen measurements (TopN, NH4+, and TN), and climate (rainfall and temperature) variables. Total nitrogen (TN) was included only in the CT wheat nitrogen measurements based on the CCA results. Bacterial community data are based on 16S rRNA gene sequencing data.

Canonical correspondence analysis (CCA) for each land use type. Redundant variables were removed if present (VIF > 15). (a) Conventionally tilled (CT) wheat, (b) no-till (NT) canola, (c) Old World bluestem (OWB) pasture, (d) native tallgrass prairie (TGP). All CCA formula values are significant (P < 0.05). Download FIG S1, TIF file, 2.5 MB.

Functional community differences between tilled cropland and native prairie.

Functional diversity at three distinct times across the 1-year data set was investigated using functional gene array for the CT wheat and TGP. For the α-diversity of the functional community, no significant differences were found based on land use or sampling time when looking at evenness and diversity (Table 1). Significant differences were observed for the functional β-diversity of the two fields. Land use, sampling time, and the interaction of the two factors were all significant in shaping functional diversity (Table 2). Interestingly, sampling time and the interaction of sampling time and land use had a stronger effect on the functional structure than land use alone.

To investigate the significant differences in functional community structure, response ratios were used to compare relative gene abundances between land use and sampling time. When compared by land use, all genes that were significantly different were greater in the TGP. Comparisons of sampling times between land use had significant differences in genes involved in carbon cycling, organic remediation, nitrogen, and metal homeostasis. Surprisingly, very few functional gene differences were observed when comparing the two fields during August (see Fig. S4 in the supplemental material). The greatest significant differences between the functional community structures of CT wheat and the TGP occurred during January (Fig. 5a), when the relative abundances of almost all genes that were significantly different were greater in the TGP. Fewer significant differences in function were observed in May (Fig. S4) than in January.

FIG 5

Functional differences between conventional tilled (CT) wheat cropland and native tallgrass prairie (TGP). (a) Differences in relative gene abundance during January between TGP and CT wheat functional community based on response ratio. All genes greater than 0.0 were greater in the TGP community. All genes present are indicated as significantly different by 90% confidence interval, 95% confidence interval (*), and 99% confidence interval (**). (b) Canonical correspondence analysis (CCA) examining the relationships between soil and environmental factors on function community structure using GeoChip data. (c) VPA of functional community structure explained by soil properties (SWC and OM), nitrogen measurements (TopN and NH4+), and climate variables (rainfall and temperature).

Differences in relative gene abundance during August and May between the tallgrass prairie (TGP) and CT wheat function community based on response ratio. All genes greater than 0.0 were greater in the TGP community. All genes present are significantly different by the 90% confidence interval, 95% confidence interval (*), and 99% confidence interval (**). Download FIG S4, TIF file, 2.6 MB.

Soil and environmental factors also impacted the functional potential of CT wheat and TGP. The local diversity of the CT wheat field was impacted by OM, SWC, and temperature measurements. Evenness and diversity significantly decreased (P < 0.05) as OM increased, while both indices significantly increased (P < 0.05) with increasing air temperature. An increase (P < 0.05) in local diversity was also associated with increasing SWC. No significant relationships were identified in the TGP land use. CCA was also used to explore the impact of soil and environmental factors on the functional community composition (Fig. 5b). Nitrogen measurements (TopN and NH4+) appeared to have a greater impact on CT wheat, and soil properties (SWC and OM) appeared to be more important for the functional structure in the TGP. Average air temperature (AvgTemp) had a more important relationship with the function of the TGP community than the CT functional community. Overall, OM grouped closer to the first axis, where functional differences were observed based on land use, while AvgTemp and NH4+ were closer to the second axis, where functional differences were separated by sampling time. SWC and TopN appeared to influence functional structure the most based on the interaction between land use and sampling time. Similarly, variance partitioning analysis (VPA) showed a large amount of the variation in the functional structure to be unexplained (Fig. 5c), and soil properties, nitrogen measurements, and climatic factors all explained a comparable amount of variation. Soil properties greatly interacted with climatic factors and nitrogen measurements. While all soil factors and climatic variables were important to variations in functional diversity, more links to these factors and local functional diversity were observable for the highly managed CT wheat field.

Links between structural and functional community types in tilled cropland and native prairie.

The Mantel test was used to examine the relationship between taxonomic, phylogenetic, and functional bacterial community structure, focusing on key functional groups (e.g., carbon cycling, nitrogen cycling, and virulence). When considering overall community interactions for the TGP, almost all functional groups had a negative relationship with taxonomic diversity (see Fig. S5a in the supplemental material). In comparison, most of the relationships between functional groups and phylogenetic diversity were positive, with a moderately significant association with methane cycling genes (r = 0.229, P = 0.07). The strength of the relationship between functional groups and phylogenetic diversity was significantly greater (P < 0.001) on average than that with taxonomic diversity for the TGP. For individual sampling times, a similar pattern was observed where the strength of the relationships between functional groups and phylogenetic diversity was greater than taxonomic diversity. For the CT wheat field, the majority of relationships between functional groups, taxonomic diversity, and phylogenetic diversity were negative (Fig. S5b). The most significant relationships were found between functional groups and taxonomic diversity for CT wheat, including carbon cycling (r = 0.26, P = 0.05), methane cycling (r = 0.31, P = 0.02), and organic remediation (r = 0.24, P = 0.05). The strength of the relationships was significantly greater (P < 0.001) between taxonomic diversity and functional groups on average. The correlation strength for each sampling time for CT wheat was greater between taxonomic diversity and functional diversity for individual sampling months as well. Primarily, functional structure was more strongly associated with phylogenetic structure in the TGP compared to CT wheat, where functional structure was more closely related to taxonomic structure.

Relationship between functional diversity, taxonomic diversity, and phylogenetic diversity for tallgrass prairie (TGP) and CT wheat based on the Mantel test. The Mantel test used Pearson correlations. Taxonomic diversity was based on the Bray-Curtis distance metric, and the UniFrac weighted distance metric was used for phylogenetic diversity. GeoChip functional data were grouped by probes into the main gene categories. Download FIG S5, TIF file, 2.8 MB.

To look closer at relationships between structure and function, functional groups were compared to OTUs. Several OTUs had a significant relationship (P ≤ 0.05) with functional groups that were present in both fields (see Table S5 in the supplemental material) with the direction of the relationship often varying between fields. For example, all the OTUs linked to carbon cycling, carbon degradation, and methane cycling in both fields had a negative correlation for the TGP and a positive correlation for CT wheat, likely related to their importance to the processes directly or indirectly within each field. OTUs significantly correlated with N fixation present in both fields were all Bacillus spp. and were present in a much greater abundance in CT wheat, generally having highs in August and May compared to highs in just May for the prairie. One OTU (Massilia sp.) significantly correlated with organic remediation in both fields and had a higher abundance overall in the TGP, with peaks in August, whereas abundance was greatest in May for CT wheat. Significant OTUs were then compared within fields to determine if OTUs were significantly correlated with more than one function. When examining the functional groups of carbon cycling, nitrogen cycling, and organic remediation, the TGP had almost 60 OTUs that significantly correlated with all three processes, compared to only 3 in CT wheat, potentially reflective of how specialized the community is in the CT wheat field. The three OTUs in CT wheat all belonged to a different phylum, while the significant OTUs in the TGP included many species in the orders Rhizobiales, Bacillales, and Solirubrobacterales.

OTUs significantly correlated (P ≤ 0.05) to functional gene groups present in both the tallgrass prairie (TGP) and conventionally tilled (CT) wheat field. Correlations are based on Spearman correlations, with the direction of relationship presented as positive (+) and negative (−). Download Table S5, DOCX file, 0.02 MB.

DISCUSSION

Impact of land use and seasonality on soil properties.

Land use change, management intensification, and season have different effects on soil properties and thus impact microbial communities in different ways. In this study, we examined how the soil ecosystem was affected by an increasing amount of management disturbance across four land uses commonly found in the U.S. Southern Plains. Land use change and management intensification modify the soil environment and generally reduce soil quality (48), as illustrated by the decrease in OM, TN, and SWC under tillage management compared to other land uses (49). Reducing management disturbance resulted in several soil properties being indistinguishable between land uses, further signifying that removing intensive management improves vital soil properties (1). Meanwhile, it had been previously observed that land uses under comparable amounts of management resulted in similar edaphic properties when cropland and non-cropland soil properties were compared (3, 50), which may explain parallels between properties in NT canola and OWB pasture, which received similar yearly management. Only soil properties related to climate and management significantly differed by season. Lower SWC was evident in times of low monthly rainfall or increased daily temperatures. For the croplands, soils exhibited highs of TopN during summer and fall due to fertilizer application, which is expected as management practices in agricultural fields are largely seasonally dependent (23, 51). Overall, even though land use was the greatest determinant of soil properties, sampling time was also key for explaining differences in soil properties, especially as management disturbance increased.

Impacts of land use and seasonality on soil bacterial diversity.

Determining how soil microbial community diversity is impacted across time and space is crucial for preserving soil health against continued environmental changes. The α-diversity and β-diversity of bacterial communities were distinctively altered by land use and season. As has been observed in a similar study comparing land use types and temporal dynamics (21), season had the most significant impact on α-diversity, with different land uses having greater diversity at various times of the year. However, the interactive effect of season and land use on belowground diversity remains unclear as most studies emphasize spatial variability over temporal variability (52). On large spatial scales, variation in α-diversity is not significantly explained by land use but rather soil properties (7), with moisture and nutrient availability generally being the most notable factors (20, 53, 54). Increases in TopN in the croplands decreased α-diversity (55, 56), while SWC influenced α-diversity in the TGP. Both properties generally vary over shorter periods of time, making them potentially better predictors of seasonal microbial community changes (6). Even with the documented impact of season on α-diversity, it is thought that the importance of temporal dynamics is underestimated due to the presence of relic DNA (57), the response of different taxa to environmental changes (52), and the lack of focus on living/active cells (58, 59). Given that no land use had the greatest α-diversity throughout the whole 1-year period and α-diversities were influenced by soil properties that vary seasonally, it is important to assess temporal dynamics when trying to determine differences in the microbial community.

The effects of land use were far more critical for regulating the β-diversity of the bacterial communities across the management gradient. The β-diversity of all land use types differed from that of the TGP (Fig. 2), with tillage management having the most significant impact (Fig. 3). Between land uses, several soil factors and air temperature were critical for differences in bacterial diversity (Fig. 3), while the distance between sites had no significant effect on the smaller scale of our study. Within fields, as management decreased, less variation in β-diversity was explained by the measured properties, of which the relative importance became more evenly distributed (Fig. 4). No group of variables was the most important to variations in TGP bacterial diversity. In comparison, slightly greater importance of nitrogen and other soil properties was found to be associated with variations in OWB pasture, possibly due to changes in microbial community composition and diversity from management disturbances in grasslands (60, 61). More variation in bacterial diversity was explained in the croplands. While both croplands were fertilized, soil N content was far more important to bacterial diversity in CT wheat, presumably because fertilizer was applied with no residue cover and directly incorporated through tillage. Soil properties that increased under NT management like SWC and OM explained more variations in the NT canola, supporting that reduced management increases carbon storage and moisture availability (62, 63). Sampling time was also a significant driver of diversity differences (Fig. 3), with rainfall and soil nutrients again having considerable influence. This is consistent with previous studies where climate variables, soil moisture, and nutrient availability dictated temporal changes (6). Although several factors were exclusive to shaping bacterial diversity based on time or space, SWC, OM, and TN continually appeared to be notable factors impacting the bacterial communities (64–67), with land use type being critical to explain differences in diversity, especially compared to the native system. It should be noted, that this was the first time canola was planted on the NT cropland, which had previously been a long-term winter wheat system. While plant species can influence the microbial communities, many other factors in croplands likely outweigh the introduction of a new crop. In agricultural systems, crops are cultivated in various soils being impacted by the soil type, soil properties, and land management, often reducing the importance of the rhizosphere microbial community for plant growth compared to native ecosystems (28). Soil properties have also been shown to override the influences of crop type on soil bacterial communities (68), with land use and management strongly shaping soil properties (1, 21). Additionally, a mesocosm experiment using soil collected from long-term monoculture cropping systems determined that the cropping history of the soil was the main factor determining the microbial community composition when a new crop was introduced (69). Together, these points help emphasize that the plant type during this single growing season was likely not responsible for the overall observed differences.

While much is still unknown about the relationship between taxonomic/phylogenetic and functional diversity, it is widely believed that increased diversity, including functional diversity, sustains soil functions and creates greater resilience to disturbance and stress (70, 71). Taxonomic/phylogenetic and functional diversity can also be differentially affected by soil and environmental properties. Based on results from the FGA analysis, land use and sampling time were both central in shaping the functional diversity of the CT wheat and TGP field, although land use alone had less of an effect than sampling time or the interaction of sampling time with land use. The reduced effect of land use on functional diversity is likely due to shared taxa between communities leading to more similar functional traits (72, 73) and the redundancy of many biogeochemical gene families across microbial groups (74). TGP functional diversity was associated with greater SWC, OM, and air temperature, and CT wheat functional diversity was associated with higher N content. Available N has been shown to significantly impact the active bacterial community and increase the number of taxonomic and phylogenetic groups that specialize in using N compounds (58). We also attempted to uncover the correlations between taxonomic/phylogenetic diversity and functional diversity, although deciphering such correlations is not straightforward. Functional diversity had stronger correlations to taxonomic diversity than to phylogenetic diversity in the CT wheat field, whereas in the TGP, functional diversity had stronger relationships with phylogenetic diversity. It is possible that the CT wheat community remains more phylogenetically similar over time, while the taxonomic community changes more rapidly. These types of patterns have been previously observed and suggested as warning signs of biodiversity loss due to environmental changes (75, 76) resulting from intensive management practices in agroecosystems.

Impacts of land use and seasonality on soil bacterial community composition.

Throughout our study, the greatest management disturbance resulted in the greatest impact on the bacterial community, as shown by the results of tillage treatment at both the phylum and OTU levels. The impact of land management, especially tillage, on bacterial community composition has been extensively documented (1, 77, 78), and although less studied, season has considerable influence on composition as well (21, 23, 79). For all land uses, the most unique OTUs were present during the spring season. During spring, air temperatures begin to rise and rain increases. Temperature and moisture not only impact the physiological activity of bacterial communities, but also regulate plant activity, including rapid growth and increasing root exudates (80, 81). Such large seasonal changes are likely responsible for differences in community composition observed between land uses as well as the increase in bacterial richness during the spring season. Monitoring changes in microbial composition over time and in response to management is one of the best ways to determine sustainable agricultural practices as it can indicate early potential changes in soil functionality, although it is necessary to remember there is not one optimal microbial community composition.

To examine the functional gene community composition, relative gene abundances of the whole communities were compared between CT wheat and the TGP. Between the two land uses, the abundances of all genes that significantly differed were always greater in the TGP. Such differences are believed to be reflective of microbial functional gene abundance and diversity (46), although gene presence does not necessarily mean the gene is being actively transcribed. More distinct differences in gene abundances between land uses were apparent when comparing specific sampling times (Fig. 5; Fig. S4). In general, seasonal microbial community differences are usually more evident in agricultural soils compared to native soils (21) due to seasonal management practices and plant activity. The greatest differences occurred during January, when plants in both fields were generally not active, air temperatures reached yearly lows, and CT wheat had the lowest SWC. The importance of soil water content in regulating microbial activities is well known, with soil water content being a key abiotic factor linked to functional diversity (82). Furthermore, the greater ground cover (i.e., residues) during the winter in the TGP may help alleviate the stress of the colder temperature on the microbial community, with greater plant litter amounts also increasing water infiltration and reducing soil evaporation (83). Therefore, the effects of reduced SWC and reduced ground cover could lead to decreased microbial diversity and activity under CT wheat. The smallest number of differences in the functional gene community was observed in August. The tallgrass prairie mainly consists of warm-season grasses; therefore, the plant community is in peak growth during this time, likely releasing nutrients to support microbial activity. In comparison, the CT wheat field is tilled during the summer fallow season to incorporate residues for decomposition providing organic carbon and nitrogen, again likely resulting in increased microbial activity (84, 85). Even though there were clear differences in the functional diversities of the microbial communities in relation to land use and sampling time, it is equally necessary to survey changes in functional gene abundance as shifts in diversity alone do not always result in differences in the biogeochemical functional ability of the soil microbial community (86, 87).

Conclusions.

Environments in agroecosystems are continually modified due to land use and management practices that can, directly and indirectly, influence soil bacterial communities. Soil communities are exposed to variability in space and time, making no single biotic or abiotic factor the sole reason for shifts in bacterial community composition, raising the need for continued research on a range of agricultural systems. In this study, we investigated the effects of land use and sampling time on the structural and functional diversity of bacterial communities as well as the interactions with soil and environmental factors in four land uses in the U.S. Southern Plains. First, our results indicated that land use, especially with intensive management, had the greatest impact on taxonomic diversity, while sampling time and time within a specific land use were more important for differences observed in functional diversity. Next, soil nutrients, particularly nitrogen, and soil water content were determined to be critical for variations in community taxonomic and functional diversity across land management and sampling time. Last, functional diversity was also reduced under intensive management, with species likely being more specialized in function due to fertilizer usage and more strongly linked to taxonomic diversity than phylogenetic diversity. Although the impacts on functional and structural diversity may have different relationships with land use and sampling time, it is clear that both types of diversity are important for structuring the interactions of edaphic properties, climatic factors, and bacterial communities. The results contribute to the idea that preserving microbial diversity should be one of the main focuses of sustainable agriculture. While these observations may be regionally specific, we recommend sampling around management practices (e.g., August) as sampling in relation to a specific management practice or environmental change likely provides the most insight when trying to determine the impact on soil health. This is one reason why microbes show great promise as a soil health indicator as they can respond to disturbance before plant communities and soil properties. Additionally, we further recommend the use of no tillage as it increased the total nitrogen, organic matter, and water content in the soil, in comparison to CT management, which increased the reliance on nitrogen inputs, generating a less diverse and likely more specialized bacterial community. Moving forward, continued monitoring of changes in bacterial communities within local land uses’ corresponding natural and anthropogenic disturbances will likely be most useful when trying to make informed decisions about managing soil health and ecosystem services.

MATERIALS AND METHODS

Site description and field sampling.

The study was conducted at the United States Department of Agriculture, Agricultural Research Service, Grazinglands Research Laboratory at El Reno, OK (35°34.1′N, 98°03.6′W; 414 m above sea level), from August 2016 to July 2017. Soil was collected from four sites: native tallgrass prairie (TGP), introduced (OWB) pasture, NT canola field, and CT winter wheat field. The grasslands and croplands were approximately 2.7 km apart. All four sites are included in the Southern Plains site of the Long-Term Agroecosystem Research (LATR) Network (88, 89). The soil type for the research area was Bethany silt loam (a fine, mixed, superactive, thermic Pachic Paleustolls). The study area has a temperate continental climate, with summer months being characteristically hot and dry with a 30-year (1980 to 2010) average daily maximum and minimum air temperatures of 22.5°C and 8.8°C, respectively, and rainfall mostly occurring in May to June and September to October, with an average annual rainfall of 860 mm (90–93).

Native tallgrass prairie consists of mainly warm-season native mixed grasses. This mixture includes big bluestem (Andropogon gerardii Vitman.), little bluestem [Schizachyrium scoparium (Michx.) Nash], Indiangrass [Sorghastrum nutans (L) Nash], and switchgrass (Panicum vergatum L.). Old World bluestem (Bothriochloa spp.) is a monoculture pasture site that was established well over 20 years before this study was conducted. Both pasture sites had deep soils (>1-m depth) and high water-holding capacity (90, 94). During the sampling period, the TGP was grazed by beef cattle for approximately five of the sampling months at a stocking density of 0.13 head/ha for 30 days and 0.83 to 1.06 head/ha for the remaining months. The OWB pasture was grazed for roughly eight of the sampling months at a stocking density ranging from 0.65 to 0.94 head/ha. Prescribed burns of the pasture sites are on a 4-year rotation, with the most recent burning occurring in February 2014. The OWB pasture is managed using annual fertilizer and herbicide treatments (91), while the native TGP is not fertilized treated. In these pastures, vegetation generally greens up in April and enters the senescence phase toward the end of October, with peak growing season occurring during the May to June period (95).

As a cool-season crop, winter wheat is the dominant cultivated ecosystem in the U.S. Southern Plains, generally converted from native tallgrass prairies. Winter wheat has been planted in the study sites under CT management since the late 1990s. In Oklahoma, winter wheat fields are managed for multiple purposes, such as grain-only, graze-grain, and graze-out. The CT wheat field was grain-only during the 2015 to 2016 growing season and graze-out (no grain production; cattle grazing from November through May) during the 2016 to 2017 growing season. Each year, the seedbed was prepared for planting using a chisel plow treatment to a depth of 31 cm, which resulted in complete disturbance of soil and residue mixing (84). The NT canola field was converted from a CT wheat field in 2015. It was grain-only wheat during the 2015 to 2016 growing season and on canola rotation during the 2016 to 2017 growing season as a part of the 4-year crop rotation (canola, grain-only wheat, graze-grain wheat, and graze-out wheat). It was the first year canola had been grown in the NT plot. Both croplands were left fallow from June to September, being fertilized and planted between late August and mid-October. Crops had fall and spring growing seasons and were dormant during the winter months. The CT wheat site was harvested in early June, and the NT canola site was harvested in late June. Detailed management data have been previously published (96).

Soil sampling began in August 2016 for all sites. Soil sampling was conducted every 2 weeks during the fall and spring months and once a month during summer and winter months, resulting in 20 sampling times per field between August 2016 and July 2017 for a total of 80 soil samples. During each soil sampling time point, eight cores roughly 20 m apart were taken in a random-walking pattern throughout each field at a depth of 0 to 15 cm using a 2.5-cm-diameter soil probe. Soil cores were pooled and homogenized at each sampling time for a representative sample of each plot. Soils were sieved to 2 mm to remove debris and stored at −80°C until analysis.

Soil properties and climate data.

Weather data, including average monthly rainfall, maximum air temperate, average air temperature, and minimum air temperature, were gathered from an Oklahoma Mesonet weather station (http://www.mesonet.org/index.php/weather/local/elre) in El Reno (ELRE), Oklahoma. The Mesonet tower is located on the native TGP site used in this study (35°32.9′N, 98°02.2′W). Soil chemical analyses were performed at the Oklahoma State University Soil, Water and Forage Analytical Laboratory (https://agriculture.okstate.edu/departments-programs/plant-soil/soil-testing/). Tests included topsoil nitrate (TopN), soil organic matter (OM), soil total nitrogen (TN), and ammonium (NH4+). Gravimetric soil water content (SWC) was determined by oven drying for ≥24 h at 65°C or until the weight no longer changed (84).

Soil DNA extraction, PCR amplification, and sequencing.

Microbial genomic DNA was extracted from 0.25 g of soil using the Quick-DNA fecal/soil microbe miniprep kit (Zymo Research, Irvine, CA, USA) according to the manufacturer’s instructions, and DNA was eluted with sterile water. For each soil sample, four technical replicates were extracted. DNA was quantified with the Qubit double-stranded DNA (dsDNA) BR assay kit (Thermo Fisher Scientific, Waltham, MA, USA), as described in the manufacturer’s instructions. DNA dilutions of 2 ng/μL were prepared for use in PCR. PCR was performed using a two-step barcoding protocol (97). For the first DNA amplification, primer pair M13-tagged 341F (5′-GTAAAACGACGGCCATACGGGNGGCWGCAG-3′) and 785R (5′-GACTACHVGGGTATCTAATCC-3′) were used (98). The second PCR step used a barcoded version of the forward primer and the 785R primer stated above. The PCR for amplification used a 50-μL mixture containing 0.1 μL of each primer (100 μM), 2 μL of template DNA, 25 μL of Phusion high-fidelity PCR master mix with HF buffer, and 23 μL of water. The PCR conditions were a preliminary denaturation phase at 95°C for 5 min then 30 cycles at 95°C for 30 s, 55°C for 30 s, 72°C for 60 s, and final extension at 72°C for 10 min. PCR products were checked on a 1% agarose gel for amplification and purified using a QIAquick PCR purification kit (Qiagen, USA) before the barcoding reaction. The PCR barcoding step used a 30-μL reaction mixture containing 0.15 μL of reverse primer, 1 μL of the barcode forward primer, 5 μL of purified PCR product, 15 μL Phusion high-fidelity PCR master mix with HF buffer, and 8.85 μL of water. The PCR conditions for the barcoding reaction were 95°C for 5 min, followed by 6 cycles at 95°C for 30 s, 55°C for 60 s, and 72°C for 90 s, and a final extension at 72°C for 10 min. PCR product was then pooled and further purified before sequencing using an Illumina MiSeq platform (Illumina, USA) at the Oklahoma Medical Research Foundation (Oklahoma City, OK).

Sequence analyses.

Raw FASTQ files were checked for quality with FASTQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc) then demultiplexed and processed using QIIME (version 1.9.0) (99). Low-quality sequences and chimeras were removed. Operational taxonomic units (OTUs) were clustered at 97% sequence similarity using the UCLUST function in QIIME. Technical replicates for each sample were combined to increase sequence number per sample and get one representative sample per time point. The OTU representative sequences were taxonomically identified using the SILVA 16S database. Sequences were rarefied to 12,000 sequences per sample based on the lowest sequence number sample to use for α-diversity, calculated using the vegan package (100) in R version 4.0.3 (101). β-Diversity was calculated using the vegan package using the unrarefied data.

Functional gene array.

The GeoChip 5.0S array containing ∼60,000 probes per array (35) was used for functional gene analysis. Microarray analysis was conducted following previously described protocols (34, 35). In short, three time points were chosen from the 1-year sampling period to represent various sampling seasons from TGP and CT wheat. These two fields were chosen for further examination because CT wheat croplands are the most common type of land conversion of native prairie systems. The DNA extracted for sequencing was also used for this part of the study. Each DNA sample was purified using Agencourt AMPure XP (Beckman Coulter, CA, USA) bead purification following the manufacturer’s protocol. The quality of the DNA was determined using a Nanodrop ND-1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE) based on the A260/A280 and A260/A230 ratios, and DNA concentrations were quantified again using Qubit dsDNA BR assay kit (Thermo Fisher Scientific, Waltham, MA, USA). DNA was labeled using random priming and cyanine dye, purified using Qiagen QIAquick purification kit per manufacturer’s instructions, and dried. Resuspended labeled DNA was mixed with hybridization solution and pipetted into the center of the well of the gasket slide (Agilent), and the array was assembled and sealed and placed into a rotisserie hybridization oven to hybridize in the presence of 10% formamide at 67°C for 24 h. Once hybridization was complete, slides were washed and imaged using a NimbleGen MS200 microarray scanner (Roche NimbleGen, Madison, WI, USA).

GeoChip data were normalized and quality filtered with methods modified from previous ones (34, 35). Probes flagged as outliers (bad spots) were removed from all samples. Then, for each array, the sum of the signal intensity was calculated, and the maximum sum value was applied to normalize the signal intensity in each array, producing a normalized value for each spot in each array. Normalized data were further denoised as follows. A probe signal is counted as low rank in a sample if a raw signal is <500, a signal/noise ratio (SNR) is <2, a signal/background ratio (SBR) is <1.3, a coefficient of variation (CV) is >0.8, a signal is <99% of detected probes, or a signal is <50% of designed probes. Only the probes showing low-rank signals across all samples were removed as noise. Then, probe signals with an SBR of <1.1 were filled with zeros, considered undetected.

Statistical analyses.

Differences in measured soil properties were compared across land use and season using the Kruskal-Wallis rank sum test. The effects of land use type and seasonal sampling time and their interactions with α-diversity and β-diversity indices were analyzed using R. The main effects and interactions of land use and season on α-diversity indices were tested using avop in the “lmPerm” package (102) in R. Tukey’s post hoc test was used when significant values (P < 0.05) were returned for one of the main effects or interactions. Analysis of variance (ANOVA) was used to compare the effect of soil properties to α-diversity within land use types. Principal-coordinate analysis (PCoA) was performed using Bray-Curtis distance metrics. The statistical significance of effects of season and land use on β-diversity was assessed by permutational multivariate analysis of variance (PERMANOVA) using adonis in the vegan package (100). The same analysis was used for α-diversity and β-diversity of the functional community. Bacterial community composition differences were compared by field using adonis and pairwise field comparisons. The same method was used for within-field seasonal differences of bacterial community composition. To link β-diversity to measured soil and environmental factors, the Mantel test and multiple regression on distance matrix (MRM) were modified as previously described (103). In the modified Mantel test and MRM, β-diversities of spatial, temporal, or all pairwise comparisons were, respectively, subjected to a linear mixed model with random effect of intercepts in different seasons. The significance test was based on constrained permutation of samples considering the repeated measurement. The factors in the MRM models were forward selected based on adjusted R2. Canonical correspondence analysis (CCA) was conducted to determine the effect of soil properties and environmental factors on the bacterial community within land use and for the overall functional data set. CCA models were considered significant when P was <0.05 and redundant variables had been removed (variance inflation factor [VIF] of >15). Each variable was additionally checked for significance within each model. Variance partitioning analysis was conducted for bacterial community composition of each field based on CCA results. Variables were separated into three groups representing nitrogen measurements (TopN, NH4+, and TN), soil variables (OM and SWC), and climate factors (rainfall and air temperature). To examine differences in relative gene abundance based on GeoChip data, genes present in at least 50% of the samples across treatment were used. Response ratios were determined using an online available MicroArray functional gene microarray analysis system (http://ieg.ou.edu/microarray/) based on 90%, 95%, and 99% confidence intervals for land use and sampling time (104, 105). The Mantel test was also used to examine correlations between functional diversity, taxonomic diversity, and phylogenetic diversity. Correlation coefficients were compared using a two-sided t test. Spearman’s correlation was used to look at relationships between functional groups and OTUs using cor.test in R.

Data availability.

16S rRNA gene sequences were deposited in the Sequence Read Archive (SRA) database under BioProject accession no. PRJNA816491.