Soil microbial communities vary in composition and functional strategy across soil aggregate size class regardless of tillage

Soils are a complex matrix of microhabitats composed of aggregates, which form due to the interactions of mineral particles, organic matter, and soil biota (Tisdall and Oades, 1982; reviewed by Six et al., 2004). Macroaggregates (>0.25 mm) are held together chiefly by temporary binding agents, such as plant roots and fungal hyphae (Tisdall and Oades, 1982; Wilson et al., 2009). As such, macroaggregates regularly undergo formation, degradation, and reformation and are heavily impacted by agricultural management practices such as tillage (Six et al., 2000; Al-Kaisi et al., 2014). Within macroaggregates, microaggregates (<0.25 mm) form from strong electrostatic interactions between mineral particles and microbially processed organic matter (Oades, 1984; Oades and Waters, 1991). The stability of aggregates and differences in physical and chemical structure impose constraints on residing microbial communities, which may require functional specialization to maximize fitness.

Soil microbial access to resources is controlled by 3 main factors: spatial separation of microbe and resource, protection of substrates via strong electrostatic bonds with the mineral environment, and resource chemistry (Baldock and Skjemstad, 2000; Schmidt et al., 2011; Dungait et al., 2012; Lehmann et al., 2020). Pore space and pore connectivity within aggregates modulate microbial access to resources by controlling water availability and microbial access to nutrients, oxygen, and carbon substrates (Sexstone et al., 1985; Foster, 1988; Rabbi et al., 2016). Macroaggregates have a greater number of large pores and higher pore connectivity, which may provide for a microbial habitat with greater resource availability (Negassa et al., 2015; reviewed by Wilpiszeski et al., 2019). By contrast, microaggregates have smaller, disconnected pores, which can occlude resources from microbial access (McCarthy et al., 2008; Rabbi et al., 2016). In addition to spatial separation of resources, substrates become inaccessible to microbes as a result of organomineral complexes; this is one of the driving mechanisms behind carbon persistence in microaggregates (Chenu and Plante, 2006; Grandy and Neff, 2008; Schmidt et al., 2011). Resource availability is further controlled by carbon chemistry (Lehmann et al., 2020). The majority of carbon contained in macroaggregates is in labile plant-derived compounds, while microaggregates contain more microbially processed carbon (Elliott, 1986; Christensen, 2001; Six et al., 2004; Davinic et al., 2012).

The physical and chemical environment within soil aggregates affects microbial community composition (Trivedi et al., 2015; Bach et al., 2018; Upton et al., 2019) and physiology (Lagomarsino et al., 2012; Bach and Hofmockel, 2014; Trivedi et al., 2017; Pold et al., 2020; Liu et al., 2021). To clarify patterns between community composition and physiology broadly, microbial ecologists have begun to use life history frameworks to explain how microorganisms optimize survival under varying abiotic conditions (Fierer et al., 2007; Wood et al., 2018; Malik et al., 2020). These frameworks reduce the vast phylogenic diversity of microbes into functional groups from which we can make predictions about ecosystem processes (Wallenstein and Hall, 2012). The yield–acquisition–stress (Y–A–S) life history framework proposed by Malik et al. (2020) hypothesizes that trade-offs exist between microbial allocation of carbon resources toward yield, resource acquisition, and stress tolerance. The Y–A–S framework is a microbially focused adaptation of Grime (1977) framework and centers on the emergence of microbial strategies due to resource limitation and stress (Malik et al., 2020). Yield strategists maximize the proportion of resources allocated to biomass and prosper in environments free of resource limitation. Environments where resources are scarce are hypothesized to favor acquisition strategists with metabolic machinery optimized for the uptake of readily available substrates or the breakdown of complex polymers. Habitats with stressful conditions (i.e., low water availability, pH extremes) favor stress tolerators that invest resources into survival mechanisms such as cell membrane integrity or DNA repair at the expense of growth or nutrient acquisition (Schimel et al., 2007; Malik et al., 2020).

Because the Y–A–S life history framework categorizes microbes by their carbon acquisition and allocation, and other functional traits, it is a useful framework to understand microbial contributions to carbon cycling in soils. Thus, applying this theory in the context of agricultural systems may help to predict how carbon is stored, which can provide a path toward improving agricultural sustainability. Further, agricultural soils provide a unique system to test the applicability of the Y–A–S framework at the aggregate scale because field-scale management decisions alter soil structure and redistribute resources. For example, tillage decreases stability in both macro- and microaggregates and accelerates turnover time, which may lead to changes in microbial community dynamics at the aggregate scale (Al-Kaisi et al., 2014). Further, tillage has been shown to reduce carbon relative to untilled soils resulting from a loss of carbon-rich macroaggregates and an increase in the proportion of carbon-depleted microaggregates (Six et al., 2000). Although life history frameworks have been applied previously to soil microbial communities in agricultural systems, these studies were conducted on intact soil (Schmidt et al., 2018; Wood et al., 2018; Malik et al., 2019). Applying the Y–A–S framework—and life history strategies broadly—at the scale of soil aggregates may yield additional insights into microbial community composition and function and their responses to agricultural management, because these microhabitats are unique environments with distinct physical and chemical characteristics.

To examine the relationship between microbial function, community composition, and the soil physical environment, we isolated aggregates of 4 different sizes (<0.25, 0.25–1, 1–2, and >2 mm) via a modified optimal moisture sieving method (Bach and Hofmockel, 2014), from soils in a 7-year tillage experiment. We measured 3 metrics of microbial physiology important for soil health and often used to define life history strategy: direct uptake of simple carbon resources using multiple-substrate-induced respiration (MSIR), potential enzyme activity, and carbon use efficiency (CUE). Further, we used predictive metagenomics to infer community stress tolerance within aggregates using traits and genomic markers defined in Malik et al. (2020). We hypothesized that aggregates of different size classes harbor distinct microbial communities that can be classified into life history strategies. Further, we hypothesized that microbial physiology and predicted genes would differ across aggregate size classes in alignment with the Y–A–S framework and that this relationship would be affected by tillage. The hypotheses for each aggregate size class were as follows: (1) aggregates >2 mm will harbor organisms that employ a yield strategy, classified as having the highest CUE, due to high resource availability and a less stressful environment relative to other aggregate size classes; (2) aggregates 0.25–2 mm will harbor acquisition strategists due to associated carbon chemistry that necessitates enzymatic breakdown; and (3) aggregates <0.25 mm will harbor stress tolerators due to low pore connectivity in microaggregates, which may provide stressful conditions (i.e., osmotic/oxidative stress) (Figure 1). Life history frameworks have the potential to improve the predictions of microbial community function and, thus, soil function; this study seeks to examine life history frameworks at a microbially relevant habitat scale to understand how management alters microbial traits and soil function.

Soils were collected from a long-term tillage experiment in Madbury, NH (N 43°10.333′, W 70°56.307′). The soils at the site are Hollis-Charlton and Charlton fine sandy loams and the average annual precipitation is 1,202 mm. The tillage experiment was established in June 2013. Prior to that time, the site had been managed for silage corn and hay for at least 10 years. Six tillage treatments were randomly assigned to plots within each of 4 blocks. The tillage treatments were full-tillage, strip-tillage, and no-tillage, each replicated twice within a block. Full-till (FT) plots were managed by inversion of the soil to 30 cm depth with a moldboard plow followed by the preparation of the seed bed with shallow disking. No-till (NT) plots were not tilled. The experiment was managed as a corn–soybean rotation, with all plots planted to corn in 2013, 2015, and 2017 and soybean in 2014, 2016, and 2018. Herbicides (glyphosate each year except 2015; dicamba in 2015) were applied for weed control prior to planting crops in the NT treatment and several weeks postplanting in both the NT and FT treatments. Liquid dairy manure was applied to each treatment in 2013, several weeks prior to corn planting to maintain soil fertility. In 2016, each of the replicate plots for a given tillage treatment within a block was assigned to an additional treatment: crop seed coated with pesticides or the same crop seeds without the pesticide coating.

Soils were sampled in May 2019, prior to tillage (or herbicide application) and planting to avoid short-term effects of management-related disturbances, as this study was focused on long-term effects of soil management on soil microbial communities. Within each block, we sampled both the FT and NT plots managed with and without pesticide-coated seeds. Soils were sampled using a slide hammer coring device to a depth of 10 cm. Twelve soil cores were collected across each tillage treatment plot and composited in one bag. Samples were stored on ice and transported to 4°C within 6 h. Prior to aggregate isolation, soils were passed through an 8-mm sieve to facilitate the removal of roots and large rocks while still maintaining aggregate structure. The pesticide seed coat treatment did not affect any of the variables measured in this study and these plots have been retained in the analyses.

Aggregates of 4 different size classes (<0.25 mm = microaggregates; 0.25–1 mm = small macroaggregates; 1–2 mm = medium macroaggregates; >2 mm = large macroaggregates) were isolated using an optimal moisture sieving technique with minor modification. This method has been shown to consistently separate aggregates while minimizing disturbance to the soil, making measurements of microbial physiology more reflective of in situ conditions (Bach and Hofmockel, 2014). Soils were dried for 1 week to 0.1-g H₂O g soil⁻¹. Aggregate size classes were isolated by placing repeated batches of 375-g dried soil on stacked sterile sieves with mesh openings of 250 µm, 1 mm, and 2 mm affixed to a vibratory sieve shaker (40 amps for 2.5 min; Retsch AS200) (Tiemann et al., 2015). The vibratory sieve shaker was used to minimize disturbance to the aggregates while still facilitating consistent aggregate separation similar to the procedure of Bach and Hofmockel (2014). A subset of each aggregate size class was dried for the analysis of total carbon, total nitrogen, pH, and electric conductivity (EC). A second subset of 3 g was frozen upon isolation at −80°C for DNA extraction and community classification. A third subset was frozen at −20°C for enzyme analysis. The remainder of each size class was stored at 4°C for physiological analysis via MicroResp™ and CUE. All analyses for this study were conducted within 3 months of initial sampling.

Aggregate pH and EC were determined with a soil to water ratio of 1:5 using Accumet Basic AB15 and Accumet Basic AB30 (Fisher Scientific, Hampton, NH) probes, respectively. Field moisture content was measured by drying 10 g of soil at 105°C for 3 days and recording mass loss. Water holding capacity (WHC) was determined by wetting approximately 50 g of each aggregate to capacity and measuring moisture content (as above) after draining for 48 h. Total carbon and total nitrogen were measured by elemental combustion analysis (Costech EC4010, Valencia, CA).

DNA extractions were performed using the Qiagen DNeasy PowerSoil Kit (Venlo, the Netherlands), following the manufacturer’s protocol with slight modifications to maximize DNA recovery from samples (Geyer et al., 2019). DNA concentrations were quantified using the PicoGreen fluorescence assay (Quant-iT^TM PicoGreen Kit, Thermo Fisher, Germany) on a Biotek Synergy HT microplate reader (Winooski, VT).

Bacterial and fungal community composition were analyzed through amplicon sequencing of ribosomal RNA (rRNA). The V4 region of the bacterial 16S gene was targeted using primers 515f/926r and primers ITS1F/ITS2 were used to target fungal communities (White et al., 1990; Parada et al., 2016). The conditions of amplification were as follows: enzyme activation at 95°C for 3 min, 35 cycles of denaturing at 95°C for 30 s, annealing at 55°C for 30 s, and extension at 72°C for 60 s; followed by a final extension at 72°C for 12 min. Amplification of the desired regions was confirmed via gel electrophoresis. 16S and ITS amplicons were sequenced at the UNH Hubbard Center for Genome Studies (University of New Hampshire, Durham, NH) on the Illumina HiSeq 2500 platform (Illumina, San Diego, CA). Average read counts per sample for 16S and ITS amplicons were 44,895 and 89,525, respectively.

Bioinformatic analysis of the bacterial and fungal sequences was conducted in QIIME2 version 2019.4 (Bolyen et al., 2019). For 16S sequences, primers were removed using Cutadapt (Martin, 2011) and sequences were trimmed and denoised using DADA2 (Callahan et al., 2016). Sequences were rarefied to 4,300 reads (Figure S1). Following rarefication, 5 samples were dropped from analysis due to low read counts. Taxonomy was assigned by comparing reads against the GreenGenes database (version 3.8.19). For ITS sequence analysis, primers and conserved regions were removed using ITSxPress (Rivers et al., 2018). The remaining steps follow the 16S methods outlined above through DADA2. ITS taxonomy was assigned using the UNITE database (Nilsson et al., 2019) and sequences were rarefied to 22,000 reads.

Community weighted genome size and rRNA copy number for the bacterial community were predicted via ancestral state reconstruction in R (version 3.6.3) (R Core Team, 2018). Briefly, bacterial amplicon sequences were placed on a reference phylogenetic tree described in Gravuer and Eskelinen (2017) using pplacer (Matsen et al., 2010). Function phyEstimate() in the picante package (version 1.8.2) (Kembel et al., 2010) was used to estimate genome size and rRNA copy number for each 16S sequence based on the phylogenetic distance from a reference sequence with known genome size and rRNA copy number. Amplicon sequence variant (ASV) relative abundance was corrected for rRNA copy number using the script from (Kembel et al., 2012) and community weighted mean (CWM) trait values were calculated using the FD package (Laliberté et al., 2014).

Functional gene prediction for bacterial communities was performed using Phylogenetic Reconstruction of Unobserved Traits 2 (PICRUSt2; Douglas et al., 2020). PICRUSt2 works by predicting the metagenomic composition of 16S amplicons using relative organisms, for which there are fully annotated genomes as a reference. Functional gene abundance, measured as KEGG orthologs (KOs), was predicted via hidden-state prediction (Louca and Doebeli, 2018). The reference database was validated using weighted nearest sequenced taxon index values, which indicate the average phylogenetic distance from query sequences to known reference sequence (Table S1). The KEGG database contains information on experimentally characterized gene function. KEGG gene functions were mapped onto ASVs observed in this experiment using orthologous genes present in both the database and observed ASVs.

MSIR was measured using the MicroResp^TM system (Campbell et al., 2003) with 7 substrates representing various chemical groups: L-arginine and L-sysine (amino acids); D-(+) trehalose, D-(+) glucose, and L-(+) arabinose (carbohydrates); and citric acid and L-malic acid (carboxylic acids). Water (Milli-Q H₂O) was used as a control. Three replicates for each aggregate size class were included. Following a 5-day preincubation, substrates were added (30-mg substrate per g total soil water) to ∼0.362 g of aggregates in deep-well microplates fitted with the MicroResp^TM rubber gasket. A colorimetric CO₂ detection plate was affixed to the rubber gasket and soils were incubated for 6 h at 25°C. Following the incubation, the absorbance of the detection plate was measured, and activity was inferred. CO₂ production was corrected for differences in biomass by subtracting H₂O-induced respiration. Substrate-induced respiration (µg C-CO₂-g dry-soil⁻¹ hr⁻¹) was analyzed on a per substrate basis and as a sum of all substrates.

The potential enzyme activity of 5 hydrolytic enzymes was assessed for each aggregate size class: β-glucosidase (BG), N-acetyl-B-D-glucosaminidase (NAG), Leucine-7-aminopeptidase (LAP), cellobiohydrolase (CBH), and acid phosphate (PHOS), following protocols adapted from Saiya-Cork et al. (2002) and German et al. (2011). Fluorescent standards (MUB and MUC) and substrates (BG, CBH, NAG, LAP, and PHOS) were made within 1 week of conducting analysis and stored frozen at −20°C. Optimal substrate concentration and incubation time were determined using V-max tests conducted on soils from these sites in 2018. For each assay, a sample control, buffer control, and substrate control were measured as well as 16 replicate wells of each aggregate size class by enzyme combination. Fluorescence was measured using a Biotek Synergy HT microplate reader at 360- and 460-nm excitation and emission wavelengths, respectively. Total enzyme activity corrected for differences in microbial biomass was analyzed (µmol enzyme activity µg biomass⁻¹ h⁻¹).

CUE was measured using the ¹⁸O water tracing method from Geyer et al. (2019). Aggregates were preincubated at 40% WHC and 25°C for 5 days. For each sample, paired incubations were conducted with either a 20 at% ¹⁸O-water enrichment or unlabeled deionized water. Aggregates were incubated at 25°C for 24 h. Then, 3 mL of vial headspace was sampled via syringe and injected into a LI-COR 6252 benchtop CO₂ analyzer (Lincoln, NE). Samples were analyzed for δO¹⁸ quantification via temperature conversion elemental analysis (PyroCube Elementar Analysis System, Hanau, Germany) at the UC Davis Stable Isotope Facility. Due to machine failure during O¹⁸ quantification, CUE analysis includes 33 datapoints spanning all 4 aggregate size classes and both tillage treatments (<0.25, n = 8; 0.25–1, n = 7; 1–2, n = 8; >2 mm, n = 9). For comparison, all other physiological analyses include 64 datapoints (n = 16 per aggregate class).

Microbial biomass was calculated via the substrate-induced respiration method of Anderson and Domsch (1978) with modification. Microbial biomass can be calculated from the following conversion: 1-µL CO₂ h⁻¹ glucose-induced respiration corresponds to 40-µg microbial biomass carbon over short incubations (1–3 h) at 22°C (Anderson and Domsch, 1978). Glucose-induced respiration from MicroResp™ measurements (µL C-CO₂ g⁻¹ dry soil h⁻¹) was converted to microbial biomass C (µg C g⁻¹ dry soil h⁻¹) using the equation from Anderson and Domsch (1978):

Soils were incubated at 25°C for 6 h rather than 22°C for 1–3 h. As such, incubation conditions were modifications from the original Anderson and Domsch method and may represent an overestimation of microbial biomass. However, we expect that these modifications would affect samples uniformly, and thus, the relationship between treatments would be maintained. To validate measurements of microbial biomass, linear regression was performed to test the relationship between total carbon and microbial biomass as these two measurements are known to be highly correlated (Figure S2) (Anderson and Domsch, 1989). To further validate biomass calculations, linear regression was performed to test the relationship between soil DNA (ng g⁻¹ dry soil) and microbial biomass (Figure S3). Calculations for CUE follow equations from Geyer et al. (2019) adapted from Spohn et al. (2016).

Statistical analyses were conducted in R (version 3.6.3; R Core Team, 2018). Normality and homoscedasticity of all data were assessed using Shapiro–Wilks and Levene’s test, respectively. pH did not meet the assumptions of normality and was log transformed prior to analysis. Differences in abiotic characteristics among aggregates and between tillage treatments were determined using two-way analysis of variance (ANOVA). Post hoc Tukey HSD tests were performed to determine which aggregate size classes were driving statistically valid differences. Statistical differences in MSIR and enzyme activity between tillage treatments and aggregate size classes were determined using permutational multivariate ANOVA (PERMANOVA) in the “vegan” package (Oksanen et al., 2018). Experimental block significantly impacted MSIR and enzyme activity, so permutations were constrained by block to minimize this effect. As tillage significantly impacted both MSIR and enzyme activity, individual substrates and enzymes across aggregate size classes were analyzed separately for NT and FT treatments via ANOVA following log transformation to meet the assumptions of normality.

Differences in microbial communities and predicted gene composition across aggregate size classes and tillage treatments were assessed with PERMANOVA and visualized using nonmetric multidimensional scaling with Bray–Curtis dissimilarity. Indicator ASVs in bacterial communities were identified using the “indicspecies” package in R. Bacterial and fungal microbial communities were distinct between experimental blocks, so PERMANOVA randomizations were constrained by block to account for this effect. KOs were analyzed using STAMP v2.1.3 (Parks et al., 2014) for differential abundance of predicted genes across aggregate size classes using White’s nonparametric t test with Benjamini–Hochberg false discovery rate (FDR) correction. FDR was set at 0.05.

Soil aggregate distribution differed between the FT and NT treatments. Macroaggregates (<2 mm) were significantly reduced within FT plots relative to NT plots. In both FT and NT plots, the 0.025–1 and <2 mm aggregates made up the majority of the aggregate fractions (Figure S4). Total carbon did not vary across aggregate size classes within either tillage treatment (Table S2; F = 1.917; P > 0.05). However, under NT, total nitrogen was greater in the <0.25 mm aggregate size class compared to in the 0.25–1 and 1–2 mm size classes (Table S2; F = 4.74; P < 0.01).

Microbial community composition differed between tillage treatments and aggregate size classes; however, there was no interaction effect. Bacterial community composition was modified by tillage treatment (PERMANOVA; F = 6.81; P < 0.001) and differed by aggregate size class (PERMANOVA; F = 1.45; P = 0.01) (Figure 2a). Differences in community composition between aggregate size classes were driven largely by the >2 mm and <0.25 mm size classes (P = 0.036). Diversity metrics of both evenness and Shannon’s diversity were unaffected. WHC, pH, and EC were strongly correlated with bacterial community composition (pH, P = 0.005; EC, P = 0.005; and WHC, P = 0.001). Indicator species analysis identified two ASVs in the order Actinomycetales associated with the <0.25 mm aggregate size classes. Two ASVs in the phylum Verrucomicrobia and Actinobacteria were associated with the 0.25–1 mm aggregate size classes. Additionally, there were 5 indicator species in the >2 mm aggregate bacterial communities spread across 5 phyla (Table S3). Like bacterial communities, fungal communities differed at both the field (PERMANOVA; F = 5.27; P < 0.001) and aggregate scale (PERMANOVA; F = 1.45; P < 0.01; Figure 2b). However, fungal communities in the <0.25 mm aggregate size class differed significantly from those of both the 1–2 and >2 mm aggregate sizes. All environmental variables tested (pH, EC, WHC, total N, total C, and C:N) were correlated with fungal community composition in the NT plots.

The microbial acquisition strategy was gauged by 2 metrics: activity in response to carbon substrates and investment into enzyme production. Microbial uptake and acquisition of simple carbon molecules was measured via MSIR of 7 carbon substrates across 3 compound classes: amino acids, carboxylic acids, and carbohydrates. We used the microbial respiration induced by these simple substrates as a proxy for activity and analyzed the microbial response to all 7 substrates in total (i.e., as the sum of respiration across substrates) and individually. Aggregate size (PERMANOVA; F = 1.77; P < 0.05) and tillage (PERMANOVA; F = 11.99; P < 0.001) impacted community response to carbon substrates with consistently higher activity in NT plots. Although there were differences in aggregate microbial community composition, MSIR was largely consistent across aggregate size classes with few exceptions. Glucose-induced respiration decreased with increasing aggregate size. Arabinose, also a carbohydrate, elicited an opposite response albeit with a much lower magnitude. Microbial response to citric acid was affected by aggregate size (F = 4.301, P < 0.01) with the highest respiration in 1–2 mm aggregates and lowest in aggregates <0.25 mm (Table S2).

Extracellular enzyme activity was measured to explore microbial community investment into nutrient and carbon acquisition; specifically, we measured 2 carbon-acquiring (BG and CBH), 2 nitrogen-acquiring (LAP and NAG), and 1 phosphorus-acquiring (PHOS) enzymes. As a result of differing microbial biomass between tillage treatments (F = 27.92; P < 0.001) and aggregate size classes (F = 10.96; P < 0.001), enzyme activity was relativized per µg biomass (Steinweg et al., 2013; Malik et al., 2019; Figure 3). We analyzed enzyme data jointly (i.e., the treatment effect on the responses of all 5 enzymes together as multivariate data via PERMANOVA) and individually (i.e., each of the enzymes via ANOVA). Aggregate size class and tillage treatment did not impact enzyme activity when analyzed jointly via PERMANOVA. However, when individual enzymes were analyzed, patterns began to emerge. Biomass relativized CBH activity was significantly lower in tilled plots (0.21 ± 0.019-nmol activity µg biomass⁻¹) relative to no-tilled plots (0.32 ± 0.022-µmol activity µg biomass⁻¹) (F = 15.662; P < 0.001). As aggregate size class increased both CBH and BG activity tended to increase however, not significantly (P = 0.08 and P = 0.09 respectively). NAG activity was 1.5 times higher in aggregates >2 mm relative to aggregates <0.25 mm (F = 5.25, P < 0.01). Nitrogen acquiring enzyme activity (i.e., sum of LAP and NAG) was significantly higher in aggregates >2 mm and significantly greater than the 0.25–1 mm and <0.25 mm size classes however was unaffected by tillage.

CUE was analyzed to understand microbial investment into the production of new biomass (i.e., yield in the Malik et al. [2019] framework). Tillage had no significant effect on CUE. However, CUE was affected by aggregate size (F = 3.49; P = 0.029), which was highest in aggregates <0.25 mm (0.27 ± 0.032) and lowest in aggregates >2 mm (0.17 ± 0.016) (Figure 4).

CWM genome size and rRNA copy number were analyzed as a measurement of community stress tolerance and growth, respectively. Predicted CWM genome size was nearly impacted by tillage in the two-way ANOVA model (P = 0.0802). In the FT system, predicted genome size was significantly larger in the >2 mm size class and ranged from 4.46 ± 0.039 to 4.72 ± 0.095 Mbp in <0.25 mm and >2 mm aggregates, respectively (F = 5.029; P = 0.018). However, under NT, CWM genome size was greatest in the 0.25–1 mm aggregate size class (F = 4.21; P = 0.015; Figure 5). rRNA gene copy number was unaffected by either aggregate size or tillage treatment.

Traits conferring stress tolerance were inferred from predicted functional pathways obtained from KOs predicted using PICRUSt2 (Douglas et al., 2020), a tool which allows for inference into the potential genes, and therefore the potential functions, of the bacterial community. KOs will hereafter be referred to as “predicted genes.” Statistical analyses were performed on log-transformed predicted genes. Predicted genes varied by tillage and aggregate size class (PERMANOVA; till treatment: F = 5.41; P < 0.001; aggregate size class: F = 2.18; P < 0.001). There was also a significant interaction of tillage and aggregate size class (F = 1.68; P = 0.003). In the <0.25, 1–2, and >2 mm aggregate size classes, KO abundance differed between till treatments; however in the 0.25–1 mm aggregate size class, tillage had no effect (P = 0.116).

Post hoc pairwise comparison revealed that the predicted gene composition of the <0.25 mm aggregate size class differed from that of the 0.25–1 mm (P < 0.01), 1–2 mm (P < 0.01), and >2 mm (P < 0.05) aggregate size classes. The top 1% of predicted genes with the highest variance between aggregates were subsampled from an initial pool of 7,258 genes and differential abundance of the 73 high variance genes retained was calculated between the <0.25 mm aggregate size class and all other aggregate sizes (Figure 6). Differences in gene composition between samples were driven largely by the abundance of genes encoding fatty acid synthesis/metabolism. Further, stress tolerance genes that encode for cold shock proteins and sigma 70 factors were also enriched in microaggregates.

Enzyme activity and CUE were negatively correlated. Linear models determined CUE was negatively correlated with total carbon acquiring, total nitrogen acquiring, and total phosphorus acquiring activity, suggesting that there is a trade-off in microbial investment to growth versus carbon and nutrient focused resource acquisition (Figure 7).

To highlight microbial investment into yield versus stress tolerance strategies, linear regression was used to compare CUE and CWM genome size. Previous studies have shown that genome size is reduced in environments with high acidity, aridity, heat, and salinity, which suggests that smaller genomes may be a potential stress tolerance strategy (Cortez et al., 2022; Simonsen, 2022). CUE was negatively correlated with genome size (Figure 8), suggesting that having a larger genome may come at the expense of CUE. However, this does not support the hypothesized trade-off between yield and stress tolerance.

Most studies observe microbial function and composition using bulk soils after sieving and homogenizing soil structure. However, this may obscure the complex interplay between microbial communities and the soil physical environments in which they reside. The structural complexity of soils includes innumerable microhabitats including inter- and intra-aggregate spaces (Foster, 1988; reviewed by Wilpiszeski et al., 2019). The results presented here suggest that observing microbial communities at the finer spatial resolution of aggregates can yield insights into the variability of microbial function, which can be understood using microbial life history strategies frameworks. Our findings show that soil aggregates harbor unique microbial communities with different preferences for specific low molecular weight carbon compounds and that microbial CUE differs with aggregate size. Additionally, genes conferring stress tolerance vary by aggregate size class. Taken together, these results suggest that trade-offs within the Y–A–S life history strategy framework apply to soil aggregate microbial communities and that these relationships persist regardless of tillage.

At the field scale, tillage changed aggregate size class distribution across this study. Tilled plots had 5.5% more microaggregates <0.25 mm and the proportion of aggregates >2 mm was significantly reduced. While carbon distribution in agricultural soil generally varies depending on both the scale of study (aggregate vs. whole field) and agricultural management (Grandy and Robertson, 2007; Trivedi et al., 2015; Bach et al., 2018), bulk soils from this tillage experiment do not show differences in total carbon with tillage possibly due to the relatively recent adoption of the NT treatment. Further, we found there was no difference in total carbon between aggregate size classes, which may have been a result of the aggregate isolation method chosen (Xu et al., 2017).

Both bacterial and fungal community structure were impacted by tillage and aggregate size class. Previous studies have found that tillage drives changes in community structure largely due to changes in the chemical and physical composition of soils (Carbonetto et al., 2014). Community structure in aggregates is driven by both the quality of resources (Davinic et al., 2012) and physical structure of the aggregate environment including pore size and connectivity within aggregates, which can alter microbial access to carbon (Sexstone et al., 1985; Mummey et al., 2006; Ebrahimi and Or, 2016). Differences in substrate quality and physical structure have been shown to alter microbial community composition in aggregates previously (Trivedi et al., 2017). We found no interaction between aggregate size and tillage, suggesting that communities within aggregates are affected by tillage in a consistent manner regardless of aggregate size. This highlights the importance for studying aggregate-scale interactions in agricultural soils.

MSIR differed between aggregate size classes; however, there were no consistent patterns across substrate classes. Glucose-induced respiration accounted for a quarter of total MSIR in <0.25 mm aggregates. In agreement with Lagomarsino et al. (2012), glucose respiration was inversely related to aggregate size in both FT and NT treatments. However, arabinose (also a carbohydrate) elicited an opposite response. Higher arabinose utilization in larger aggregate size classes may be the result of microbial communities that are accustomed to processing coarse plant biomass as L-arabinose is an important component of plant cell walls (Seiboth and Metz, 2011). For example, fungi are responsible for degrading plant cell wall components, such as cellulose, hemicellulose, and pectin (Aro et al., 2005). Therefore, higher arabinose utilization in larger aggregates could potentially be attributed to a larger fungal community. Microbial activity in response to citric acid additions accounted for 25% of total MSIR in aggregates greater than 1 mm and increased with increasing aggregate size contrasting the findings of Lagomarsino et al. (2012). Given the role of plant roots as binding agents in macroaggregate formation (Tisdall and Oades, 1982), microbial communities in larger aggregates may be more adept at utilizing organic acids, like citric acid, given that they constitute a large proportion of plant root exudates (Tawaraya et al., 2014).

In addition to direct uptake of simple compounds, microbes also maintain resource demands by producing enzymes to break down complex polymers. When analyzing the additive effects of nitrogen and carbon acquiring enzymes, investment into total nitrogen acquiring enzymes was higher in >2 mm aggregates relative to in 0.25–1 mm aggregates, suggesting that macroaggregate communities may harbor communities that favor resource acquisition. One possible explanation for this trend is that high enzyme production in macroaggregates could be a result of favorable pore structure and diffusion gradients found in >2 mm aggregates. Enzymes released exogenously diffuse into the environment to reach their target substrate; thus, resource acquisition via enzyme production should be favored in environments with high pore connectivity, such as macroaggregates (Rabbi et al., 2016). Despite often having more complex carbon in need of enzymatic depolymerization, microaggregates may be unfavorable environments for diffusion due to low porosity making return on investment into enzymes low (Rabbi et al., 2016). Individual carbon acquiring enzymes (BG and CBH) were unaffected by aggregate size. CBH, however, was significantly lower in FT plots, which could be a result of changes in fungal community composition due to soil disturbance (Helgason et al., 2010; Kyaschenko et al., 2017).

Substrate-independent CUE methods capture the measures of in situ microbial efficiency that are more representative of the constraints imposed by the environment than can be obtained by the addition of substrates (Geyer et al., 2019), making them a good method for assessing microbial strategies (e.g., yield) in response to the aggregate environment. The values reported in this study fall within the same range of other studies using substrate-independent CUE (Geyer et al., 2019; Malik et al., 2019). Contrary to our hypothesis that aggregates >2 mm would harbor yield strategists, CUE was lowest in macroaggregate communities, and this relationship was not impacted by agricultural management. It is possible that high yield is advantageous to microorganisms in microaggregates, where carbon resources are scarce due to physical separation (Dungait et al., 2012). Despite emerging patterns, it is hard to draw broad conclusions as microhabitat variation in CUE is generally understudied (Anthony et al., 2020).

Agricultural management and resultant changes in abiotic soil properties impact the genomic potential of soil microorganisms (Trivedi et al., 2015). Predicted metagenomic features differed across bacterial communities between aggregate size classes; however, this relationship was unaltered by agricultural management. Specifically, CWM genome size increased significantly with increasing aggregate size. This relationship may be the result of differences in the stability of aggregates as microbial habitats. Microaggregates are more stable environments with slower turnover times and longer persistence (Al-Kaisi et al., 2014). Studies have shown that stable environments select for smaller, streamlined genomes while variable environments tend to select for organisms with larger genomes and a greater breadth of metabolic capabilities to adapt to changing conditions (Bentkowski et al., 2015). Thus, differences in stability between macroaggregates and microaggregates may impose selection pressures that drive genome size. Further supporting this idea, the relationship between genome size and aggregate size was made more apparent by agricultural management. In aggregates >2 mm, genome size was greater for FT than NT (4.72 ± 0.095 and 4.54 ± 0.023 mega base pairs, respectively). Tillage decreases macroaggregate stability (Al-Kaisi et al., 2014). As such, larger genome sizes in FT macroaggregate communities could be a stress tolerance response to decreased macroaggregate stability relative to aggregates under NT. Aggregates <0.25 mm had similar genome sizes in NT and FT soils (4.48 ± 0.012 and 4.46 ± 0.039 mega base pairs, respectively), which is consistent with the idea that microaggregates are more stable regardless of tillage (Al-Kaisi et al., 2014). rRNA gene copy number is a genomic feature included in most studies of bacterial life history strategy (Fierer et al., 2007; Krause et al., 2014; Schmidt et al., 2018) as it is positively correlated with bacterial growth rate (Klappenbach et al., 2000; Roller et al., 2016). However, unlike other life history frameworks, within the Y–A–S framework, growth rate is an emergent property of the microbial community and not a strategy in and of its self (Malik et al., 2020). Consistent with this, we did not see any trends in rRNA copy number between aggregate size classes or tillage and there were no apparent trade-offs with any of the physiological metrics studied.

The functional genes of a microbial community determine functional capabilities and the emergent community life history strategy (Wood et al., 2018). Gene composition of aggregates <0.25 mm differed from that of other aggregate size classes, supporting the notion that microaggregates may select for microorganisms employing a stress tolerance life history strategy. Stress tolerance is largely conferred by cell membrane chemistry, given that this is the only barrier of defense between a cell and the surrounding environment (Russell et al., 1995). Of the 8 predicted genes significantly enriched in the <0.25 mm aggregate size class, 4 were related to fatty acid synthesis or metabolism: long-chain acyl-CoA synthetase, acyl-ACP dehydrogenase, 3-oxoacyl-reductase, and acetyl-coa c-acetyltransferase. Fatty acid synthesis and metabolism is important for stress tolerance as bacteria have been shown to alter the composition and distribution of fatty acids along cell membranes in response to stressors such as pH (reviewed by Guan and Liu, 2020). Genes encoding for cold shock proteins (Csp) were also enriched in aggregates <0.25 mm. Despite their name, many Csp gene families are also important in microbial response to osmotic, oxidative, and starvation stress (reviewed by Keto-Timonen et al., 2016). Thus, their enrichment in aggregates <0.25 mm further supports the notion that microaggregates harbor a stress tolerant bacterial community. Sigma factors are RNA polymerase subunits responsible for transcription initiation (Paget, 2015). Sigma 70 contains 4 subgroups of sigma factors, 3 of which are involved in response to nutrient limitation, oxidative, and osmotic stress (Paget, 2015). Predicted genes coding for RNA polymerase sigma 70 factors were significantly enriched in microaggregates, which further bolsters the evidence of stress tolerance. While microorganisms exist in complex communities in nature, most previous work to link genes and the functions they encode to stress tolerance have been explored in single organism, culture-based studies. Therefore, our application of these inferences should be further tested as accessibility to community-level traits increases due to advances in metagenomics.

The ultimate goal of the Y–A–S framework is to predict microbially driven ecosystem function and response to environmental change by grouping organisms based on shared function in response to environmental conditions. Microbial function often cannot be predicted by genes alone (Pold et al., 2020) because of in situ trade-offs that are better assessed with physiological measurements. CUE was negatively correlated with enzyme activity for both carbon and nitrogen acquiring enzymes. This trade-off between yield and investment into resource acquisition is consistent with the Y–A–S framework (Malik et al., 2019; Malik et al., 2020). However, the trade-off between CUE and phosphorus acquiring enzymes was less apparent. There were no discernable patterns of acquisition related to simple carbon molecules within aggregate sizes. This lack of relationship between CUE and MSIR may suggest that the genomic machinery for simple carbon uptake does not come at an overall cost to yield (Geyer et al., 2020). Enzyme activity and CUE measurements in this study represent one time point, and we cannot infer how these relationships change with fluctuating environmental conditions.

Although there were apparent trade-offs between yield and acquisition strategies, there was conflicting support for the trade-offs between stress tolerance and yield strategies. There was a negative relationship between genome size and CUE (R² = 0.28; P = 0.003), which suggests that larger genomes may come at a cost to yield. However, microaggregate communities had high CUE and an abundance of genes conferring stress tolerance, which may suggest that stress tolerator communities do not have inherently lower yield or that the Y–S trade-off may exist only when presented with unfavorable abiotic conditions. CUE is a dynamic trait dependent on both genomic constraints on physiology and the conditions of the surrounding environment (Domeígnoz-Horta et al., 2020). Six of the 8 predicted genes enriched in microaggregate communities were potentially related to stress tolerance. The remaining 2 genes were related to ABC membrane transporters responsible for uptake of simple molecules, which could put place microaggregates along the stressed and resource limited axis of Y–A–S framework (Malik et al., 2020). Further work needs to be done to verify the S–A and S–Y trade-offs in microaggregates under varying levels of abiotic stress.

By analyzing both physiological and metagenomic characteristics of aggregate microbial communities, we were able to assign life history strategies within the Y–A–S framework. In support of our hypothesis, we were able to classify microaggregates communities as stress tolerant based on their predicted metagenomic profile. However, the microaggregate community also had the highest CUE, which may be the result of a greater number of oligotrophic microorganisms (Trivedi et al., 2017). Large macroaggregate communities >2 mm were classified as acquisition strategists based on increased enzyme activity relative to the other aggregate size classes. Physiological responses of small and medium microaggregate (0.25–2 mm) communities did not show a strong tendency toward any particular life history strategy, suggesting that either these microbial communities do not face the same selection pressures or perhaps that communities in mid-sized aggregates are more generalist in support of a multistrategy approach. In addition, the work presented here provides support for trade-offs in resource allocation hypothesized in Malik et al. (2020). By validating this framework using soil aggregate communities, we may be able to incorporate these functional assignments into models that account for microbial strategies when predicting soil C dynamics (Wieder et al., 2015). For example, tillage shifts the distribution of aggregates, thus changing the proportions of microbial communities classified by Y–A–S tolerance, which has implications for maintaining healthy, carbon-rich agricultural systems. Using microbial life history strategies to understand the functional outcomes associated with changes in microbial communities provides insight into the impact of management activities on soil function.

Sequencing data were deposited in the NCBI Sequence Read Archive under accession number PRJNA803717. Datasets and R scripts for analyses conducted in this study can be found in the Ernakovich lab repository at https://github.com/ErnakovichLab/aggregate_life_history_strategies.

The supplemental files for this article can be found as follows:

Figure S1. Rarefaction curve showing amplicon sequence variants (ASVs) observed at various sequencing depths.

Figure S2. Linear regression of microbial biomass versus total carbon.

Figure S3. Linear regression of microbial biomass against soil DNA measured.

Figure S4. Aggregate size class distribution in no-till (left) and full-till (right) treatments.

Table S1. Community weighted NSTI values for samples included in PICRUSt2 analysis.

Table S2. Soil and microbial community properties for no-till and full-till aggregate size classes.

Table S3. Bacterial indicator species in aggregate size classes.

The authors would like to thank the New Hampshire Agricultural Experiment Station for funding this research. LTB would like to thank all the members of the Ernakovich, Frey, and Grandy labs for their input and support, in particular, Stacey Doherty, Nathan Blais, Sean Schaefer, and Mel Knorr. Thanks also to Evan Ford and John Palmer for managing the field experiment.

Funding was acquired by JGE and RGS and provided by the New Hampshire Agricultural Experiment Station. This is Scientific Contribution Number 2927. This work was also supported by the USDA National Institute of Food and Agriculture Hatch Project 1016134 to JGE and a USDA NIFA AFRI Competitive Grant (No. 2017-67013-26594) to RGS.

The authors declare no competing interests with the work presented in this study.

Contributed to conception and design: LTB, RGS, ASG, JEM, JGE.

Contributed to acquisition of data: LTB, NGW, JEM.

Contributed to analysis and interpretation of data: LTB, JEM, JGE, KMG.

Drafted and/or revised this article: LTB, RGS, ASG, JEM, NDW, KMG, JGE.

Approved the submitted version for publication: LTB, RGS, ASG, JEM, NDW, KMG, JGE.

How to cite this article: Bernhardt, LT, Smith, RG, Grandy, AS, Mackay, JE, Warren, ND, Geyer, KM, Ernakovich, JG. 2022. Soil microbial communities vary in composition and functional strategy across soil aggregate size class regardless of tillage. Elementa: Science of the Anthropocene 10(1). DOI: https://doi.org/10.1525/elementa.2022.00023

Domain Editor-in-Chief: Steven Allison, University of California Irvine, Irvine, CA, USA

Guest Editor: Elsa Abs, University of California Irvine, Irvine, CA, USA

Knowledge Domain: Ecology and Earth Systems

Part of an Elementa Special Feature: The World of Underground Ecology in a Changing Environment