Disturbance by soil mixing decreases microbial richness and supports homogenizing community assembly processes.

The spatial heterogeneity of soil's microhabitats warrants the study of ecological patterns and community assembly processes in the context of physical disturbance that disrupts the inherent spatial isolation of soil microhabitats and microbial communities. By mixing soil at various frequencies in a 16-week lab incubation, we explored the effects of physical disturbance on soil bacterial richness, community composition, and community assembly processes. We hypothesized that well-mixed soil would harbor a less rich microbial community, with community assembly marked by homogenizing dispersal and homogeneous selection. Using 16S rRNA gene sequencing, we inferred community assembly processes, estimated richness and differential abundance, and calculated compositional dissimilarity. Findings supported our hypotheses, with > 20% decrease in soil bacterial richness in well-mixed soil. Soil mixing caused communities to diverge from unmixed controls (Bray-Curtis dissimilarity; 0.75 vs. 0.25), while reducing within-group heterogeneity. Our results imply that the vast diversity observed in soil may be supported by spatial heterogeneity and isolation of microbial communities, and also provide insight into the effects of physical disturbance and community coalescence events. By isolating and better understanding the effects of spatial heterogeneity and disconnectivity on soil microbial communities, we can better extrapolate how anthropogenic disturbances may affect broad soil functions.

Introduction

Soil is a staggeringly complex, heterogeneous, and even harsh web of microhabitats that harbors vastly diverse communities of largely uncharacterized microorganisms that drive crucial soil functions such as biogeochemical cycling, organic matter decomposition, and plant productivity (Fierer 2017, Tecon and Or 2017). This diversity is, in part, underpinned by the disconnected nature of soil microhabitats (Treves et al. 2003, Carson et al. 2010) and soil spatial heterogeneity (Tilman 1994, Fierer and Jackson 2006, Rillig et al. 2017). Global changes may increase physical soil disturbances, perhaps through differences in land use (e.g. tillage), weather patterns (e.g. cryoturbation or flooding), or bioturbation (e.g. invasive species). Given reduced heterogeneity, changes to resource availability, and new microbial interactions, how do physical disturbances affect soil microbial diversity? In order to predict changes to soil biodiversity and soil function, it is essential to determine the mechanisms that drive ecological relationships under global change and disturbance.

Ascribing processes and mechanisms of community assembly to observed ecological patterns is a central question within ecology, and a particular challenge within soil communities (Hubbell 2001, Hanson et al. 2012), where ecological patterns are often weakly ascribed to the “black box” of unknown aspects of community ecology (Vellend 2010). An emerging research focus is on determining the relative importance—and thus coexistence—of niche processes (Grinnell 1917, Hutchinson 1957, Chase and Leibold 2003, 2014) and neutral processes (Bell 2001, Hubbell 2001), which both accurately predict ecological assembly patterns at various scales (Adler et al. 2007). The mechanisms by which patterns and community membership unfold are generally categorized into four ecological processes, listed here on a spectrum from niche to neutral: selection, dispersal, diversification, and drift (Vellend 2010, Zhou and Ning 2017). Selection refers to deterministic or niche-based processes dictated by biotic factors, such as intertaxa fitness differences, and abiotic factors, such as environmental filters (Hutchinson 1957). Homogeneous selection describes community assembly under similar conditions or filters, thus decreasing phylogenetic differences between communities (Dini-Andreote et al. 2015). Variable selection occurs when variable conditions produce different selective pressures, thus increasing phylogenetic differences between communities (Stegen et al. 2015). Dispersal describes the movement and establishment of organisms in space, and may occur in soil through processes such as physical disturbance, water percolation, or active dispersal in water films or saturated pores (Zhou and Ning 2017). Homogenizing dispersal increases compositional similarity between communities, whereas dispersal limitation increases compositional differences between communities, perhaps allowing for stochastic demographic changes to community composition, i.e. drift (Stegen et al. 2013). To statistically infer the relative influences of these community assembly processes in soil microbial communities, Stegen et al. (2012, 2013, 2015) have developed a null modeling approach that compares observed phylogenetic distance and dissimilarity metrics between communities to null models of stochastically assembled communities, originally demonstrated with river sediment communities (Stegen et al. 2012).

Despite robustly defined statistical models, characteristics of the soil environment and its inhabitant microorganisms interact in ways that inhibit prediction of their influences on community assembly processes (Evans et al. 2017). The “sparsely populated, frequently dehydrated, maze” of soil offers limited connectivity to its bacterial inhabitants, who typically live in spatially structured biofilms or microbial hotspots (Kuzyakov and Blagodatskaya 2015, Junkins et al. 2022), and may only interact in small communities of perhaps 120 individuals (Raynaud and Nunan 2014). Soil structure, as a source of spatial heterogeneity, can shape microbial communities by furnishing distinct microhabitats (i.e. conditions for variable abiotic selection) or by engendering microbial community isolation (i.e. inducing dispersal limitation; Rillig et al. 2017, Wilpiszeski et al. 2019). Physical disturbance to soil structure may reduce spatial heterogeneity, thus potentially altering microbial interactions and community composition (Mansour et al. 2018). This may give rise to homogenizing community assembly processes, such as homogenizing selection, through more uniform distribution of resources and abiotic conditions; and homogenizing dispersal, through direct movement of organisms.

The mixing and restructuring of microbial communities along with their spatially heterogeneous environments has been termed community coalescence (Rillig et al. 2015, 2016). The coalescence framework, as a potential driver of community composition, is relevant across many microbial environments, such as in aquatic systems where freshwater and brackish communities meet (Rocca et al. 2020), or in soil bioremediation where compost is added to restore contaminated soil (Kästner and Miltner 2016). Soils are vast collections of intermittently connected communities that often move and interact in units, such as in association with a soil particle or aggregate. Coalescence in soil thus occurs during typical soil disturbances, such as bioturbation (Jacquiod et al. 2020), agricultural tillage (Guillou et al. 2019), or cryoturbation (Gittel et al. 2014). Community coalescence events in soil are likely to be spatially fragmented, leaving much of the soil relatively undisturbed. This is one possible mechanism by which soil maintains such high levels of diversity (due to differentially affected subcommunities) and functional resilience (maintained within undisturbed communities) (König et al. 2019).

Here, we explore the effects of community coalescence in physically disturbed soils. What happens to microbial diversity and community assembly processes when soil is mixed to coalesce isolated communities and heterogeneous microhabitats? We hypothesized that well-mixed soil would harbor a less diverse microbial community, and we predicted that community coalescence would decrease richness and result in increasingly homogeneous soil communities dominated by homogenizing dispersal and homogeneous selection. Our goals were to relate differences in richness, compositional (dis)similarity, and relative contributions of community assembly processes—namely dispersal and selection—to physical soil disturbance. To address these goals, we subjected the soil environment to mixing at various frequencies over a 16-week lab incubation, and assessed the outcomes on soil microbial communities and associated community assembly processes using 16S rRNA gene sequencing and statistical models. By isolating and better understanding how spatial heterogeneity affects community assembly processes in soil, we are also better equipped to extrapolate the effects that anthropogenic processes, such as climate change or land use change, may have on broad soil functions. This work also presents relevant considerations for anyone who has homogenized soil for experimental purposes.

Methods

Soil collection

We sought to obtain an unmanaged soil in order to minimize effects of management or prior amendments, and to obtain a soil with moderate or low clay content in order to minimize extracellular DNA (Morrissey et al. 2015) while mitigating soil compaction and stickiness during manipulation. As such, we collected Freeon silt loam soil (a very deep, moderately well-drained, coarse-loamy, mixed, superactive, and frigid Oxyaquic Glossudalf; Luvisols by WRB classification) on 30 August 2018 near Connor’s Lake in Sawyer County, Northern WI, United States (45°44 52.8 N, 90°43 51.6 W, 425 m asl; Figure S1, Supporting Information). Vegetation type was northern mesic forest, early-to-mid seral, dominated by Acer rubrum L. (approximately 80%), Acer saccharum Marsh. (approximately 10%), Betula alleghaniensis Britt. (approximately 5%), and Tilia americana L. (< 5%). A total of two soil cores (1.8 cm diameter) were collected from each of six locations randomly chosen along a 50-m transect. From each of these 12 soil cores, we retained a portion of the A horizon, depth of 15–20 cm, in a large Whirl-Pak bag kept on ice prior to refrigeration. A representative subsample of air-dried soil was submitted for standard analyses and was found to be comprised of 50% silt, 36% sand, and 14% clay (hydrometer method; Bouyoucos 1962); 2.7% organic matter (loss on ignition; Schulte and Hopkins 1996), 2.2% total C, 0.1% total N (C and N by flash combustion), 5.1 pH (1 : 1 water; Richards 1954), 13 mg P kg–1, 22 mg K kg–1 (P and K by Bray-1 method; Bray and Kurtz 1945), 172 mg Ca kg–1, 25 mg Mg kg–1 (Ca and Mg by ammonium acetate method; Thomas 1982), and < 3 mg available N kg–1 (NO3––N + NH4+–N by KCl extraction; Doane and Horwath 2003).

Experimental setup and design

To investigate the effects of mixing and community coalescence on soil microbial community ecological processes and diversity, we incubated small aliquots of soil for 16 weeks, during which sets of these soil aliquots were pooled together, mixed by vortex, and redivided at various frequencies intended to mimic infrequent (e.g. tillage) and frequent (e.g. bioturbation) soil disturbances. To establish the experiment, freshly collected, field-moist soil was gently shaken through an ethanol-sterilized sieve to 2 mm and homogenized, removing any visible roots. We then established eight-tube mixing sets totaling 400 mg soil, with each tube (0.5 ml freestanding tube, catalog number 16466–036, VWR) containing 50 mg soil (± 5 mg; Figure S3, Supporting Information). These eight-tube sets were randomly assigned to mixing treatments, which determined how frequently the soil aliquots in the set would be pooled in one tube (2 ml freestanding tube, catalog number 89004–308, VWR), mixed by vortex, and redistributed for further incubation over the duration of a 16-week incubation period (Fig. 1). The mixing treatments included: two times mixed (2×; soil was manipulated at the beginning of the incubation and again halfway through the incubation), four times mixed (4×; every fourth week), eight times mixed (8×; every other week), 16 times mixed (16×; weekly), and 32 times mixed (32×; twice weekly; Fig. 1). There were four replicate mixing sets for each treatment. To control for the effect of soil disturbance in absence of pooling with other soil, the vortex controls were stand-alone tubes of soil (50 ± 5 mg) that underwent vortex mixing, but were never pooled with any other soil (Fig. 2; n = 8 per mixing treatment). To control for the effects of incubation, there were 32 tubes of soil (50 ± 5 mg), that incubated undisturbed for the duration of the experiment (1×, or control).

Figure 1.

Visual representation of mixing treatments over the course of the 16-week soil incubation. The mixing treatment represents the number of times each mixing set was pooled together, mixed by vortex, and redivided. The 1× controls (n = 32) were incubated undisturbed.

Figure 2.

Experimental setup depicting one pooled mixing set in the 4× mixing treatment. Note: the incubation interval varied in length depending on the mixing frequency. At the conclusion of the 16-week incubation, DNA was extracted from all soil in each incubation tube for 16S rRNA gene amplicon sequencing. See Figure S2 (Supporting Information) for an expanded version of this figure, detailing all treatments.

At respective times of mixing, soil from the eight tubes within a given mixing set was combined in one larger tube. Gravimetric moisture was restored to approximately 25% using autoclaved Milli-Q water, and the pooled soil was agitated using a vortex mixer (catalog number 02215365, Fisher Scientific) fitted with a horizontal tube holder at speed seven for 5 seconds. Following vortex mixing, the soil was evenly divided back into the eight incubation tubes (Fig. 2). The 1× and 2× treatments underwent monthly moisture correction, on an individual tube-basis (without pooling or mixing), to mitigate excessive soil drying.

The cap of each incubation tube had one 1/32″ hole for air exchange, drilled at a 45° angle (for the vortex controls, an intact cap was used during vortex mixing). All tubes and caps were autoclaved prior to use. Tubes were incubated in two identical dark incubation boxes at room temperature and > 95% relative humidity to reduce soil drying (Figure S3, Supporting Information). Temperature and relative humidity were continuously monitored in each incubation box (data not shown). The incubation boxes were frequently opened for treatment manipulation, and thus kept aerated.

In order to characterize the microbial community at the time of soil sampling, we also retained 32 × 50 mg (± 5 mg) soil samples, which were frozen at −80 °C without incubation (“Initial”). At the conclusion of the experiment, all tubes were frozen at −80 °C prior to DNA extraction. An electrode deionizer (catalog number 05.8091.100, Haug North America, Mississauga, ON, Canada) and antistatic nitrile gloves were used to minimize static attraction and repulsion.

DNA extraction and 16S rRNA gene sequencing

Using the DNeasy PowerLyzer PowerSoil Kit (Catalog No. 12855, Qiagen, Germantown, MD), total genomic DNA was extracted from all soil within each incubation tube, which ranged from 30 to 55 mg soil per tube at the end of the incubation. Care was taken to transfer all soil residue and DNA through a series of washes with PowerBead Solution in conjunction with vortex agitation. Complete library preparation details can be found in the Supporting Information. Briefly, the 16S rRNA genes of extracted DNA were amplified in triplicate using PCR. Variable region V4 of the 16S rRNA gene was targeted using forward primer 515f and reverse primer 806r (Walters et al. 2016). Amplified DNA was normalized and purified, prior to paired-end 250 base pair sequencing on an Illumina MiSeq sequencer at the UW-Madison Biotech Center. To obtain high coverage, the same library was sequenced twice under identical conditions, and total reads were pooled for each sample after processing as described next. Sequencing data was processed using a QIIME2 (Bolyen et al. 2019) pipeline, with DADA2 (Callahan et al. 2016) as the operational taxonomic unit (OTU, or amplicon sequence variant)-picking algorithm, and taxonomy assignment using the SILVA 132 reference database (Quast et al. 2013, Yilmaz et al. 2013). This yielded 16 687 633 demultiplexed sequences, which was reduced to 12 961 153 after denoising, with a mean length of 238 base pairs (SD = 5.5). Excluding blanks, a total of 9264 OTUs were identified. Amplicon sequences are available in the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) under accession PRJNA820861. Our primers targeted both bacteria and archaea, but because our communities were dominated by bacteria (> 99.2% of total reads), for simplicity, we will simply refer to bacteria when discussing communities in this manuscript. Over 96% of archaeal reads represented the phylum Crenarchaeota.

Community assembly process assignments

To determine if a given mixing treatment increased the influence of any community assembly process as compared to unmixed control condition, we adapted a null-modeling method (Stegen et al. 2012, 2013, 2015; R code at https://github.com/stegen/Stegen_etal_ISME_2013) that estimates the influence of selection or dispersal using the abundance-weighted beta-mean nearest taxon distance (βMNTD; the mean phylogenetic distance between each OTU in one community and its closest relative in another community; Fine and Kembel 2011), and Bray–Curtis dissimilarity (BC-Dis), respectively (Bray and Curtis 1957). Unlike applications of this model to different field-based communities, we used the 1× control soil condition as the baseline scenario, creating null distributions for both βMNTD and BC-Dis based on every pairwise comparison amongst the 1× controls (496 comparisons total). Observed βMNTD and BC-Dis values from mixing treatments (112 within-pooled mixing set comparisons for each mixing treatment) were then compared to the null distributions to determine the relative effects of selection and dispersal. The 1× control samples, which were incubated undisturbed after the initial soil homogenization and mixing, represented stochastic community assembly in absence of mixing-induced selection or dispersal pressure. Thus, our inferred community assembly processes are always relative to the incubated 1× control. Detailed methods for community assembly processes assignment follow.

To identify a dominant influence of selection, the null distribution values of βMNTD (βMNTDNull) were arranged in ascending order and the 95% confidence interval (CI) was nonparametrically identified by finding the 0.025 and 0.975 quantiles. We then took the observed βMNTD (βMNTDObs) values for every possible pair of communities within a mixing set and compared that to βMNTDNull. Homogeneous selection was identified in comparisons for which βMNTDObs was below the 95% CI of βMNTDNull, indicating lower mean phylogenetic distance between community members than observed in the null. Variable selection was identified in comparisons for which βMNTDObs was above the 95% CI, indicating higher mean phylogenetic distance. Comparisons that fell within the 95% CI of βMNTDNull values were considered to lack a dominant influence of selection, and were subsequently tested for the influence of dispersal.

To identify a dominant influence of dispersal, the null distribution values of BC-Dis (BC-DisNull) were arranged in ascending order and the 95% CI was nonparametrically identified by finding the 0.025 and 0.975 quantiles. We then took the observed BC-Dis (BC-DisObs) values for every possible pair of communities within a mixing set for which selection was not identified, and compared these values to the BC-DisNull. Homogenizing dispersalwas identified in comparison for which BC-DisObs was below the 95% CI of BC-DisNull, indicating a higher level of similarity between community compositions than was observed in the null condition; and dispersal limitation was identified in comparisons for which BC-DisObs was above the 95% CI, indicating lower similarity. Comparisons that fell within the 95% CI for both metrics were considered undominated by any particular community assembly process.

Data analysis

Data analysis was performed in R (R Core Team 2018), version 4.0.3, using ggplot2 (Wickham 2016) for data visualization. Unless otherwise noted, the experimental unit is a tube, and may be referred to as “sample” or “community.” To describe richness, we used the weighted linear regression model of OTU richness estimates, which weights observations based on variance, to calculate 95% CIs for treatment means using the betta() function in breakaway for R (breakaway::betta; Willis and Bunge 2014), interpreting only treatments with nonoverlapping 95% CIs. Beta diversity was visualized for Bray–Curtis dissimilarities (Bray and Curtis 1957) of relative abundance data using principal coordinates analysis (PCoA) created with phyloseq::ordinate (McMurdie and Holmes 2013). To test for a significant effect of mixing treatment on community composition, we used permutational multivariate analysis of variance (PERMANOVA) to partition Bray–Curtis dissimilarity matrices among sources of variation using vegan::adonis (Anderson 2001). A significant result (P < 0.05) was subjected to post hoc pairwise comparisons, adjusting P-values using the Benjamini–Hochberg method (Benjamini and Hochberg 1995) to identify significant differences between the 1× control and mixing treatments. To test if treatments differed in their dispersion relative to the 1× control, we used homogeneity of multivariate dispersions (PERMDISP; vegan::betadisper; Anderson 2006). To quantify the degree to which tubes differed across pooled sets but within the same mixing frequency, we calculated Bray–Curtis dissimilarities (vegan::vegdist; Oksanen et al. 2019) for all pairs of tubes from different pooled sets, but the same mixing frequency, using analysis of variance (ANOVA) and a post hoc Dunnett’s Test (Dunnett 1955) to test for significant treatment differences relative to the 1× control. When assigning community assembly processes as described above, selection was inferred by βMNTD (picante::comdistnt; Kembel et al. 2010). Dispersal was inferred as described above by calculating Bray–Curtis dissimilarities (phyloseq::distance) on OTU relative abundances.

After evaluating our key questions, we assessed differential abundance to identify significant treatment-driven shifts in relative abundances of taxa. For this analysis, we compared each treatment to the 1× control (excluding taxa with mean relative abundance < 0.00001) and subjected these datasets to a beta-binomial regression model and “Wald” hypothesis test in corncob::differentialTest (Martin et al. 2021), which controls for the effect of the mixing treatment on dispersion. We report the µ value, which is the coefficient used to estimate relative abundance in the corncob model and is proportional to the fold-change in relative abundance between the treatment and control. To further understand changes in community composition, we sought to test the relationship between mixing treatment and mean predicted rRNA gene copy number, which may correlate with potential growth rate, by calculating the weighted mean predicted 16S rRNA gene copy number for each sample (Nemergut et al. 2016) and compared treatments using ANOVA and post hoc testing, as described above. Predicted rRNA gene copy numbers were obtained using the ribosomal RNA operon database (rrnDB; Stoddard et al. 2015). The R code used to perform these analyses and to create the following figures is available at https://github.com/jaimiewest/Soil-Mixing.

Results

Soil mixing decreased bacterial richness

Increased mixing frequency decreased bacterial richness in pooled mixing set samples (Fig. 3), with the most frequently pooled and mixed soil treatments (16× and 32×) demonstrating lower richness than the 1× controls, 2×, 4×, and 8× treatments, as well as the initial soil community. However, the stand-alone unpooled vortex controls maintained a consistent level of richness not statistically different from that of the 1× controls, regardless of mixing frequency.

Figure 3.

Community-level OTU richness, by pooled mixing treatment (closed points), or frequency of vortex mixing (open points). “Initial” represents richness of freshly collected soil that did not undergo incubation. Error bars represent 95% CI (± 1.96 * SE).

Pooled and mixed soil communities became more similar to each other while diverging from unmixed controls

The treatment-driven clustering pattern apparent in the PCoA ordination illustrates the importance of mixing frequency on soil microbial community composition data (Fig. 4, P = 0.001, R 2 = 0.73, PERMANOVA). Though mixed soil communities were significantly different from the unmixed 1× controls (Fig. 5A, P < 0.001 for all treatments), communities within a given pooled mixing set became more similar to each other with mixing (i.e. decreased dispersion, Fig. 5B, P < 0.0001, PERMDISP; and P < 0.0001 for all treatments, Tukey’s HSD). To this end, we can visually identify subclustering of pooled mixing sets within the 16× and 32× treatments (Fig. 4, e.g. the two eight point clusters in the upper left corner of the plot). The initial, unincubated samples were included in the PCoA in order to gauge the overall effect of the lab incubation on soil communities, which is much smaller than the effects of mixing.

Figure 4.

PCoA of Bray–Curtis dissimilarities of community relative abundances, colored by mixing treatment. Each point represents one tube community. Vortex controls (open points) were mixed but never pooled with other soil. Initial communities (gray points) represent the community present in freshly collected, unincubated soil.

Figure 5.

Bray–Curtis dissimilarities of bacterial community composition. (A) Boxplots summarize dissimilarity compared to 1× controls, quantified for each possible pairwise comparison between a treatment tube and a 1× control. The 1× boxplot represents pairwise comparisons amongst all 32 1× controls. (B) Boxplots represent the dissimilarity amongst tubes within each pooled mixing set, or amongst the vortex controls for a given mixing frequency. Note that there are four pooled mixing sets per treatment, and comparisons are made only within mixing set. (C) Boxplots represent the dissimilarity amongst all tubes within the same mixing treatment, excluding pairs of tubes from the same pooled mixing set (i.e. excluding the comparisons summarized in panel B). Asterisks represent statistically significant treatment differences from 1× based on (A) PERMANOVA, (B) PERMDISP, and (C) ANOVA: *** = P < 0.001, ** = P < 0.01, and * = P < 0.05. Statistical significance for vortex control treatments are reported in the text.

Vortex control communities were also significantly different in composition from the 1× controls, though to a lesser magnitude than the pooled and mixed soil treatments (Fig. 4, open points and Fig. 5A in black; P < 0.001, R 2 = 0.33, PERMANOVA; and P < 0.01 for each treatment compared to 1×). However, vortex controls did not become more similar to each other within a mixing treatment (compared to 1× controls; Fig. 5B in black; P = 0.15, PERMDISP).

Because we found that the pooled mixing set communities became more similar to each other with mixing, yet the unpooled vortex controls did not become more similar to each other given the same mixing treatment, we wanted to determine if there was an overall effect of the mixing treatment on community (dis)similarity amongst communities subjected to the same mixing regime, but not mixed together—i.e. comparisons of tubes from the same treatment, excluding tube pairs from the same pooled mixing set (Fig. 5C). Compared to dissimilarity amongst 1× tube communities, we found a significant treatment effect (P < 0.0001, ANOVA), with significant decreases in pairwise Bray–Curtis dissimilarities at 2× and 4× (P < 0.01, Dunnett’s), and a significant increase at 8× (P < 0.0001, Dunnett’s). Dissimilarity amongst communities at 16× and 32× was not significantly different than dissimilarity amongst 1× tube communities (Fig. 5C). Note that we tested for significant differences in Bray–Curtis dissimilarities using ANOVA because the exclusion of comparisons yielded an incomplete distance matrix, rendering PERMDISP inapplicable.

Community assembly

Soil mixing altered the relative dominance of ecological community assembly processes (Fig. 6). Pairwise comparisons within the pooled mixing sets (Fig. 6A) demonstrated that homogeneous selection dominated community assembly at the highest mixing frequencies. With less frequent soil mixing, there was a greater proportion of undominated comparisons, with 73% and 31% undominated at 2× and 4×, respectively. Homogenizing dispersal was most dominant at 4×, with 44% of comparisons, yet this growing proportion was overtaken by homogeneous selection as mixing frequency increased.

Figure 6.

The dominant community assembly processes as driven by mixing frequency (A) amongst communities within pooled mixing sets (one mixing set is illustrated in the inset), or (B) between each possible pooled sample + vortex control combination (of the same mixing frequency; one mixing set + one vortex control is illustrated in the inset). The community assembly processes are assigned using a null modeling approach, detailed in the text. The influence of selection is determined using the β-mean nearest taxon distance, and the influence of dispersal is determined using Bray–Curtis dissimilarity. The null models to which each metric is compared were created using the 1× controls; thus community assembly is inferred only relative to the 1× control condition.

To identify differences in community assembly attributable to whether soil was pooled or not, we tested comparisons between each possible pooled sample and vortex control pair within mixing frequency. A primary effect of soil agitation at a given frequency might result in homogeneous selection between pooled samples and vortex controls, however, we found little evidence for this mechanism. There was a dominance of dispersal limitation at 8×, and an increasing dominance of variable selection as mixing frequency increased from 8× to 16× and 32× (Fig. 6B). Comparisons were largely undominated at 2× and 4×.

Taxonomic composition shifted with increased soil mixing

In order to better understand community coalescence and to identify key taxa associated with community assembly processes in soil, we explored shifts in community composition related to the soil mixing treatments. The Supporting Information contains an expanded version of this section. The 1× controls had the highest phylum-level relative abundances of Proteobacteria, Acidobacteria, Chloroflexi, Verrucomicrobia, and Actinobacteria, which comprised over 80% of mean relative abundance, and also reflected the phylum-level composition of the initial soil communities (Figure S4, Supporting Information). After frequent soil mixing, over 80% of mean relative abundance were taxa from the phyla Actinobacteria and Proteobacteria, with one genus, Nocardioides (Propionibacteriales), comprising almost 30% of the mean relative abundance at 32×, and with one particular Nocardioides OTU emerging as the most abundant OTU in each of the 32× communities. The dominance of several OTUs at high soil mixing frequencies is apparent in the stark differences in cumulative mean relative abundance curves across mixing treatments (Figure S5, Supporting Information). The four most abundant OTUs at 32× (detailed in the Supporting Information), and the 10 most abundant OTUs at 16× comprised over 50% of cumulative mean relative abundance for respective treatments; this same proportion of relative abundance was comprised of almost 100 OTUs in the infrequently or unmixed treatments. Unlike their pooled soil counterparts, the vortex controls generally resembled the 1× controls in their phylum-level mean relative abundances across vortex mixing frequencies (Figure S4, Supporting Information).

When assessing taxonomic differential abundance relative to 1× controls, we found 392 taxa with positive differential abundance (i.e. enriched in mixing treatments). To make our assessment more tractable, we focused on taxa with the biggest responses (µ > 1.0), and only considered enriched taxa with mean relative abundances greater than 0.002 (0.2%) following enrichment. Fewer than 10% of the total enriched taxa met these criteria (Table S1 and Figure S6, Supporting Information). The most relatively enriched OTU (highest differential abundance estimate) at 32× relative to 1× was from the family Nocardioidaceae (Actinobacteria), which had no 100% ID matches in the NCBI nucleotide database. There were four additional OTUs from the genus Nocardioides that were considered relatively enriched, including the most abundant OTU found in every 32× community, referenced above. In this case, it is possible that these similar OTUs are, in fact, different copies of the rRNA gene that exist within a singular organism’s genome, and this is an instance of splitting a single genome into separate OTUs (Schloss 2021). Though some OTUs increased in relative abundance monotonically with increasing soil mixing frequency, other taxa peaked in relative abundance at moderate frequencies of soil mixing (Figures S7 and S8, Supporting Information). The non-monotonic responses to mixing may be attributable to bacteriophage predation once bacterial populations hit certain levels, following a “killing the winner” hypothesis, as has been observed in Winogradsky columns (Esteban et al. 2015).

There were 2152 total taxa with negative differential abundance (i.e. depleted) in the pooled and mixed treatments compared to 1×, which greatly exceeded the number of enriched taxa. Similarly to our approach for enriched taxa, we focused on the taxa with the strongest negative responses (µ < −1.0), and only retained depleted taxa that were not extremely rare to begin with (mean relative abundances greater than 0.002 in the 1× treatment). Fewer than 3% of the total depleted taxa met these criteria (Table S2 and Figures S9 and S10, Supporting Information).

In the vortex control communities, we found a similar number of enriched taxa with positive differential abundance, and some overlap with the OTUs found to be enriched in the pooled and mixed soil treatments (Figure S11, Supporting Information), but notably fewer depleted taxa (Figure S12, Supporting Information), with just four OTUs depleted across treatments, after filtering out the very rare or weakly responding taxa as described above. These four OTUs were also depleted in the mixed soil treatments.

Predicted weighted mean 16S rRNA gene copy number increased with mixing

The weighted mean predicted 16S rRNA gene copy number was statistically different across mixing treatments (ANOVA, P < 0.001), and increased with mixing frequency from a mean value of 2.09 for 1× to a value of 2.51 for 32× (Figure S13a, Supporting Information). Each treatment, except 2×, was significantly higher than 1× (P < 0.05 for 4×; P < 0.001 for 8×, 16×, and 32×; Dunnett’s). Notably, the proportion of OTUs for which a predicted 16S rRNA gene copy number was available also increased with rate of mixing (Figure S14, Supporting Information); about 30% of 1× OTUs as compared to 64% of 32× OTUs had matching genera in the rrnDB. Nocardioides sp., which comprised over 30% of relative abundance in 32×, largely accounted for this difference in OTU copy number availability by treatment (Table S3, Supporting Information). Further, with a predicted mean copy number of 2.62, Nocardioides sp. heavily weighs on this analysis, given the high proportion of OTUs for which we do not have a predicted gene copy number. To test the influence of Nocardioides on this analysis, we removed it from the calculation and found that the trends remained significant (ANOVA, P < 0.001; Figure S13b and Table S3, Supporting Information).

Discussion

The aim of this study was to determine the effects of community coalescence via physical disturbance on community composition and ecological community assembly processes in the patchy, disconnected, and heterogeneous soil environment. Consistent with our hypotheses, we found that more frequently pooled and mixed soil harbored less rich bacterial communities, with community assembly marked by homogeneous selection. The findings from this study impact our understanding of how physical disturbance affects soil communities and contribute to our growing understanding of the vast bacterial diversity observed in soil.

Bacterial richness and community coalescence

Using soil subsamples that would generally be considered homogeneous and highly similar in community composition (see “Initial” samples in Fig. 4), we demonstrated the effects of community coalescence by pooling together and mixing soil at various frequencies. The 20% decrease in bacterial richness at 32× (Fig. 3) may be attributable to competitive advantage under changing resource availability, as suggested by taxa enrichment (Figures S6 and S7, Supporting Information), or attributable to selection for stress-resistant organisms when abiotic conditions shifted beyond organismal tolerance (Rillig et al. 2015, Castledine et al. 2020). These results closely mirror those of a meta-analysis (Rocca et al. 2019), which found that alpha diversity in soil decreased by a mean of 20% across a variety of environmental disturbances encompassing a range of stressors. In our study, some combination of stress and competition could reasonably decrease richness, and, being a closed system, we would not anticipate sources of increased richness (e.g. from speciation) over the relatively short incubation interval.

The maintenance of richness in vortex controls (Fig. 3), with few depleted taxa (Figure S12, Supporting Information), supports the possibility that diversity, and thus dissimilarity (Fig. 5B), amongst these closed communities, may have been maintained by heterogeneous resource availability. For instance, an idiosyncratic fragment of organic matter in one tube could contribute to the rapid growth and selection for a particular community. As a stand-alone tube, this remains an isolated community. Conversely, if this tube belonged to a pooled mixing set, this community would be subsequently dispersed throughout the soil in all tubes of the mixing set, thus decreasing within-set dissimilarity (as seen in Fig. 5B for the pooled mixing sets) and likely contributing to both homogeneous selection and homogenizing dispersal (Fig. 6A). Due to the small individual mass of each sample in this study, we had insufficient soil to analyze postmixing edaphic characteristics, but future work could attempt to correlate resource availability with community composition.

The vortex controls may be analogous to soil aggregates, which can host isolated communities under variable selection due to patchy resource availability (Rillig et al. 2017, Wilpiszeski et al. 2019). While we would expect that some microbes may continue to remain isolated in protected soil pore spaces, or manage to persist due to priority effects during coalescence events (Castledine et al. 2020), our results suggest that, under frequent mixing conditions, the swift ascendency of a few taxa generally outweighs other mechanisms that might maintain diversity, with a parallel outcome of decreased richness.

Soil mixing selects for fast growth

Community coalescence is an often-overlooked form of disturbance (Mansour et al. 2018; e.g. Rocca et al. 2019), but the increased interconnectedness, forced chance encounters, and potential for degradation of refuges that all characterize coalescent communities in soil may help to describe the frequently observed phenomena of emergence and enrichment of previously rare taxa under rapidly changing biotic and abiotic conditions (Allison and Martiny 2008). For example, new coalescent communities were distinct and dissimilar from the 1× controls (Figs 4 and 5), and initially rare OTUs became abundant after frequent soil mixing (Figure S7, Supporting Information). Further, another study found that emergent rare taxa comprised over half of the observed OTUs in mixed brackish water coalescent communities, with many of these rare taxa becoming highly abundant at times (Rocca et al. 2020). Rare taxa impart phylogenetic plasticity to the microbiome, which can enable functional resilience during periods of transition (Jousset et al. 2017, Jia et al. 2018).

As a potential mechanism for enrichment of rare taxa, our results indicate a mixing-driven increase in mean predicted gene copy number (Figure S13, Supporting Information), suggesting that this trait imparts a selective advantage to soil organisms under frequent coalescence. These apparently mixing-loving, or at least mixing-tolerant, microbes are likely generalists that can translate available resources into fast growth, as was demonstrated in a lagoon coalescence experiment, where a diverse bacterial community of oligotrophic specialists was overcome by copiotrophic generalists (Beier 2021). Another study comparing wastewater communities in static vs. shaken conditions found that fast-growing organisms were enriched in the unstructured, shaken environments, whereas structured, unshaken environments favored organisms that invest in metabolite-mediated life strategies, presumably by maintaining proximity between expensive enzymes and their producers (Junkins et al. 2022). Thus, less structured environments, such as frequently mixed soil, seem to put organisms that rely more on extracellular metabolism at a disadvantage, instead favoring fast growers.

The prevalence of several OTUs of the genus Nocardioides (Figures S6 and S7, Supporting Information), which has a relatively higher predicted 16S rRNA gene copy number than the community predicted mean copy numbers for either the initial (unincubated) soil or the 1× control in this study, supports previous work that found this genus to be relatively enriched in coalescent soils of a bioremediation study (Wu et al. 2019), and in the high-disturbance earthworm drilosphere soil compared to the undisturbed bulk soil in a no-tillage wheat experiment (Schlatter et al. 2019). Nocardioides has also been associated with straw mineralization (Bernard et al. 2012), extracellular DNA degradation (Morrissey et al. 2015), and rapid atrazine mineralization (Topp et al. 2000). These examples suggest that Nocardioides may be a generalist that thrives in coalescent communities by translating available resources into fast growth. Other enriched or abundant organisms at high mixing frequencies also carry higher predicted 16S rRNA gene copy numbers (Figure S13b, Supporting Information), indicating that fast growth is a generally important trait under frequent mixing.

While shaped by coalescence, soil communities remain distinct

The dominance of homogeneous selection under frequent mixing and the absence of a dominant community assembly process under infrequent mixing (Fig. 6A) suggests that large populations (here, the OTUs that become relatively abundant with frequent mixing) tend to be governed by deterministic forces, and small populations (here, rare taxa that persist in the less-frequently mixed sets) are more subject to stochasticity and drift (Hanson et al. 2012). However, the specific community composition in a given tube was not driven by mixing frequency alone—in fact, pairwise comparisons between each pooled sample and vortex control combination (Fig. 6B) demonstrate evidence for variable selection. Critically, this result first highlights that the homogeneous selection identified between tubes from pooled sets is not simply due to selection for communities adapted to the soil being physically agitated—rather, there are outcomes that are specifically the result of pooling previously isolated communities during mixing. This indicates that different drivers of community composition govern the pooled samples vs. the vortex controls, despite the same mixing frequency. This may be due to variations in resource availability and biotic interactions between the vortex controls and the pooled coalescent communities, given the relatively smaller volume of soil in each one-tube vortex control as compared to an eight-tube pooled mixing set. For example, scale of mixed soil may engender differences in proximity of extracellular enzymes and metabolites to their producers (Junkins et al. 2022) and the potential associated differences in predicted 16S rRNA gene copy number, discussed above.

While distinct in composition from unmixed controls (Fig. 5A), the exact changes in community composition due to mixing varied from one pooled set to the next, as illustrated by the high level of dissimilarity across pooled samples from different mixing sets of the same mixing treatment (Fig. 5C). With moderate mixing, more stochastic community assembly processes observed at 2× and 4× (Fig. 6A) produced a mixing set-agnostic response by which we see increasingly similar community composition regardless of whether comparisons are made within mixing sets (Fig. 5B) or across mixing sets (Fig. 5C). However, as mixing frequency increased, dissimilarity decreased within mixing sets (Fig. 5B), yet dissimilarity across mixing sets remained high (Fig. 5C). Together, this emphasizes that, while community coalescence likely selects for mixing-adapted taxa (e.g.Nocardioides, amongst other strong responders; Figure S7, Supporting Information), the specific outcomes of community composition will differ, likely depending upon small differences within starting communities, or resource variability at microsites (Wilpiszeski et al. 2019) that are accentuated by repeated coalescence.

Another notable observation lies in the comparison at 2× between Fig. 6A (pairwise comparisons within mixing sets) and Fig. 6B (pairwise comparisons of each pooled sample + vortex control combination): these comparisons only differed in their treatment and handling (pooled vs. unpooled mixing) at a singular mixing event, halfway through the incubation. However, we see a sizeable difference in the outcome, with almost 30% of pairwise observations within the 2× pooled mixing sets demonstrating homogeneous selection or homogenizing dispersal (Fig. 6A), whereas comparison vs. the 2× vortex controls were largely undominated (Fig. 6B). This highlights how one soil mixing event may produce a change in the dominant community assembly process, suggesting that even subtle or infrequent soil coalescence events, such as annual tillage, could substantially shift community composition and its driving processes on a small scale.

Disturbance disrupts mechanisms that maintain soil bacterial diversity

Generally speaking, soil is largely undisturbed. That said—at relatively small scales, soil fauna burrow and consume soil. Root growth displaces soil, ultimately creating pore space as roots senesce. Soil microbes themselves contribute to aggregate formation, organo-mineral associations, and other miniature soil “disturbances.” Natural and anthropogenic disturbances, such as cryoturbation and tillage, can be significant across a landscape. These disturbance events are largely fragmented, point disturbances, and occur perhaps only occasionally in any given location. In this experiment, we demonstrated that even infrequent soil coalescence can have an impact on community composition and community assembly processes, while frequent community coalescence events resulted in significant losses of bacterial richness and the introduction of deterministic selective processes. We expect the selective processes at work are likely biotic, as we see sharp increases in the relative abundances of likely copiotrophic bacteria such as Nocardioides, in the absence of typical environmental selection filters (e.g. pH, temperature, and moisture). At a field scale, however, these results imply that high levels of diversity would likely be maintained despite mixing events. For example, a heavy rainfall that facilitates pore connectivity or a tillage event that mixes soil over a short distance may select for fast-responding taxa while decreasing richness on a small scale, but high diversity will likely persist across the landscape, as exemplified by the different mixing sets in our study, and the stand-alone vortex controls. As such, our findings generally support the hypothesis that both soil heterogeneity and spatial disconnectivity underpin the high diversity of the inhabitant microbial communities in soil (Fierer and Jackson 2006, Portell et al. 2018).

Methodological considerations

There are several important methodological considerations to this experimental study, which we detail in the Supporting Information. Briefly, these include consideration of sequencing depth, DNA from dead or dormant taxa, and the specific role of dormancy in this study.

Future directions

One future direction might be to assess microbial community function in soil undergoing natural coalescence events to elucidate the complicated relationships between microbial community diversity and function (Raynaud and Nunan 2014, Young and Bengough 2018). We could predict that frequent mixing decreases potential functional breadth due to decreased richness and the dominance of several opportunistic OTUs. However, due to high functional redundancy in soil microbial communities (Louca et al. 2018), whether there would be meaningful impacts from such a reduction may be questionable. Further, rare taxa, which we found to be characteristic of frequent community coalescence [see also Allison and Martiny (2008)] can impart functional resilience to the microbiome (Jousset et al. 2017, Jia et al. 2018), and therefore, we might also predict that soil function is maintained despite decreased richness. Another direction could be to test large-scale diversity and functional resilience of spatially fragmented or isolated coalescence events, as our results indicate that diminished diversity may only play out on smaller scales. Finally, another extension of this work could be to study the effects of soil mixing on fungi, which play an important role in soil structure and function (Crawford et al. 2012); there are likely particular implications of disturbance by mixing for filamentous fungi that connect habitats (Cairney 2005).

Conclusions

Community coalescence in what may be considered a homogeneous soil demonstrates that the bacterial community can change considerably with mixing to support potentially fast-growing bacteria, as richness otherwise declines. Homogeneous selection and homogenizing dispersal were the predominant community assembly processes in frequently pooled and mixed soil, whereas less disturbed soil was undominated by any particular community assembly processes. Despite strong mixing effects, initial differences in community composition and resource distribution appear to be important for final, mixed community compositions, as demonstrated by differences between physically disturbed (vortex controls) vs. pooled and disturbed samples. Our results generally suggest that soil heterogeneity, preserved in relatively unmixed soil, underpins the vast microbial diversity characteristic of soil environments.

Click here for additional data file.