Introduction

Understanding the mechanisms determining variation in community composition (beta diversity) across space and time is a major goal of community ecology [1,2,3,4,5,6,7]. Growing evidence has indicated that various processes can generate variation in community composition among localities. For example, community composition might vary solely due to dispersal limitation [4, 8] or stochastic processes, such as stochastic dispersal [9] and demographic stochasticity [10]. In addition, abiotic factors, as important determinants of community composition, could favor the persistence of different species along an environmental gradient [11,12,13]. However, the extent to which variation in community composition is determined by environmental filtering also depends on local biotic interactions (e.g., competition, predation, mutualism, and parasitism) [14,15,16]. Theory predicts that environmental selection interacting with local biotic interactions can result in either a single stable state or alternative stable states (hereafter referred to as ASS, also known as multiple stable states, multiple stable points, multiple stable equilibria, etc.), or mosaic cycles, where species replace one another in a cyclic or more complicated sequence through time [17,18,19]. The concept of ASS suggests that there can be more than one stable state that communities might approach during assembly or after disturbance, even under the same environment and with the same species pool [20,21,22,23,24]. A single stable state implies a fixed species composition given certain environmental conditions, whereas the other two scenarios can give rise to substantial variation in community composition.

Soil bacterial communities are among the most diverse systems in nature [13, 25]. Great advances have been made in understanding compositional variation in such communities in recent years [26,27,28,29]. Numerous studies have confirmed that environmental selection is the primary force determining the species composition of soil bacterial communities [12, 13, 25, 30]. Nonetheless, stochastic processes are also suggested to be main drivers of community divergence, which are considered as opposites to deterministic processes responsible for compositional variation unexplained by environmental factors [12, 28, 29]. However, it is worth noting that an implicit assumption in inferring stochastic processes from biogeographic patterns of community composition is that there is only a single stable state for communities under identical environments and species pools. The drawback of this assumption is that it totally neglects local biotic interactions that may promote alternative states of communities [18, 19], thus leading to overestimation of the effects of stochastic processes evaluated from the variation partitioning approach.

Coexisting soil bacterial species are engaged in a variety of interactions that may serve as local processes shaping community composition, including cross-feeding, quorum sensing, interference (toxin secretion) and exploitation competition [31,32,33,34,35]. A change in the abundance of a single species is expected to affect the fitness and population dynamics of species that directly interact with it, and such small deviations have the potential to be amplified indirectly through networks of species interactions. In addition, the effect of a change in species abundance on community composition via biotic interactions might be species specific. For example, closely related species are more likely to share similar ecological characteristics [36,37,38], leading to intense competition for resources [39] and antagonistic interactions [40]. Thus, an increase in the abundance of different species in separate local communities might specifically lead to a decrease in the abundance of their closely related clades, resulting in species-specific community divergence.

Whether or not community divergence mentioned above may persist is determined by how species abundance is further regulated in local communities. While frequency dependence and the underlying stabilizing/destabilizing processes have usually been discussed within the context of coexistence [41, 42], they are clearly linked to community divergence as well [18, 19, 43]. If stabilizing processes (also known as negative frequency dependence) were universal in soil bacterial communities, then deviations in species abundance relative to original states might be adjusted in time, facilitating stable coexistence of species at local scales [41, 42]. Thus, given a constant species pool (gamma diversity), stabilizing processes can reduce community divergence (beta diversity) by promoting stable coexistence of species in local communities (alpha diversity). In contrast, deviations in species abundance would persist if stabilizing processes were lacking or even be amplified by destabilizing processes (also known as positive frequency dependence). As a result, the coexistence of species in local communities would be disrupted; moreover, community divergence might be promoted by the reduced local diversity.

Here, we address the possibility that local biotic interactions drive the divergence of soil bacterial communities closed to microbial migration, with an abundance-manipulation experiment in Inner Mongolia’s typical steppe. The abundances of six resident bacterial species were individually increased in closed communities. By monitoring the change in community composition during the subsequent year of recovery in situ, we addressed whether a single episode of change in species abundance could drive community divergence via species interactions and whether the divergence could be maintained in the absence of migration.

Materials and methods

Study site and six distinct species manipulated in soil bacterial communities

The study site (N 43°32′30.5″, E 116°33′19.2″, H 1194 m) was located in a semiarid grassland of Xilin Gol, Inner Mongolia (Northern China), far away from anthropogenic disturbance, where the average annual air temperature was 2.0 °C and the average annual precipitation was 350 mm, falling primarily during the wet season from June to September. Chestnut soil was the major soil type, and Stipa grandis was the dominant species of vegetation. Soil samples from the top 15 cm of the four corners and the center in a 1 m × 1 m square were collected on 20 June 2012. The top 5 cm of soil was removed, and the remaining soil was then sieved with 2-mm filters and homogenized.

To obtain culturable bacterial colonies, a soil suspension was prepared by grinding and vortexing 5 g of homogenized soil sample in 15 mL of sterile water. The 10−3-diluted suspension was plated on beef-protein agar plates (3 g L−1 beef extract, 10 g L−1 peptone, 5 g L−1 NaCl and 20 g L−1 agar) and then incubated at 20 °C for 6 days. In total, 355 bacterial colonies were randomly chosen, inoculated into 24-well microplates containing 1.5 mL of fresh beef-protein medium, and then cultured at 20 °C in an orbital shaker at 200 rpm for 4 days. Subsequently, a 0.5-mL culture of each colony was stored at −80 °C, and bacterial DNA was extracted from the remaining culture with the Easy DNA Kit (Invitrogen, Paisley, UK). Each bacterial colony was identified by sequencing the 27–534 bp portion of the 16S rRNA gene via Sanger sequencing (Supplementary Methods S1). Six species varying in relative abundance were randomly chosen to serve as the manipulated species in the subsequent field experiment.

Establishing in situ microcosms with increased abundance of manipulated species

For the establishment of soil microcosms, we collected another 5 kg of soil on 9 July 2012 at the same site and depth as those used for the soil sample collection described above. After being sieved to <2 mm and homogenized, 15 g of soil was accurately dispensed in each 50-mL centrifuge tube as a soil microcosm (‘tube’ microcosm). At the same time, the six chosen species were grown at 20 °C in an orbital shaker at 200 rpm. After 4 days of incubation, 1 mL of each bacterial culture was centrifuged at 4000 rpm for 5 min, and the supernatant was discarded. The cell pellet was washed three times with 1 mL of sterile demineralized water by pipetting up and down, centrifuging at 4000 rpm for another 5 min, and discarding the supernatant again. Finally, the bacterial cells were resuspended in sterile demineralized water and starved for 4 h at 20 °C.

Approximately 109 starved bacterial cells of each chosen species in 1 mL of sterile demineralized water were added into each prepared soil microcosm and mixed well. There were five replicate soil microcosms for each of the six chosen species, in addition to ten replicates for the control group, in which only 1 mL of sterile demineralized water was added into the soil microcosms, yielding 40 soil microcosms in total. To stabilize the physical conditions of the soil after artificial disturbance, the soil microcosms were incubated statically at 20 °C for 10 days and aerated every day. Subsequently, 10 g of mixed soil from each ‘tube’ microcosm was transferred into a new ‘cage’ microcosm, which was a tube with a 0.22-μm-pore-size filter protected by a plastic grid at the bottom (Grace Catalysts Technologies, Columbia, Maryland, USA; Fig. S1). The ‘cage’ microcosms allowed molecules (water, nutrients, etc.) to pass through but blocked the migration of bacteria and any other larger organisms. The remaining soil from each ‘tube’ microcosm was stored at −80 °C, as a sample of ‘microcosms before recovery’.

The ‘cage’ microcosms were buried at the same site and depth as the soil had been sampled and covered with 5 cm of soil on 23 July 2012. One year later, the soil microcosms were retrieved and stored at −80 °C, as ‘microcosms after recovery’.

16S rRNA gene sequencing of bacterial communities in soil microcosms

16S rRNA gene sequencing was employed to detect bacterial community composition changes in soil microcosms during a whole year of in situ recovery after species abundance increased. DNA was extracted from each soil microcosm at the beginning and end of in situ recovery with the FastDNATM SPIN Kit for Soil (MP Biomedicals, Santa Ana, California, USA), yielding 80 DNA samples from soil communities in total. The 27–534 bp portion of the 16S rRNA genes was sequenced by 454 pyrosequencing (Supplementary Methods S2). 454 pyrosequencing data are archived at NCBI SRA under BioProject PRJNA491623. The raw read data were analyzed via mothur (v. 1.36.1) [44], following the mothur 454 standard operating procedure unless otherwise noted (Supplementary Methods S3). Operational taxonomic units (OTUs) were grouped at 97% similarity, defined as ‘species’ in this study.

Statistical analysis

To examine the change in relative abundance of the six chosen species during in situ recovery, we first matched the sequences of the six chosen species with 454 pyrosequencing data at 100% similarity to identify their relative abundance (the number of sequences of manipulated species divided by the total number of sequences in each community). Subsequently, the relative abundance of the six chosen species was analyzed separately using analysis of variance (ANOVA), with relative abundance in the corresponding treatment and control group as a response variable and treatment (two levels: abundance increased and control) and time as explanatory variables.

To eliminate the influence of artificially increased species abundance on the dissimilarities between communities, we removed the exact sequence variants of six manipulated species from the communities and rarefied sequences to 12497 per sample in subsequent analyses (Supplementary Methods S3). The phylogenetic distance between species was obtained by constructing a phylogenetic tree in QIIME (v. 1.2) [45] using FastTree (v. 2.1.3) [46]. Weighted UniFrac, which integrates species relative abundance and phylogenetic distance, was used as an index to quantify pairwise community dissimilarities. To examine whether communities responded differently to distinct species abundance increases, principal coordinate analysis (PCoA) in a two-dimensional space based on a dissimilarity matrix was performed to display the divergence of communities. We tested for species-specific effects of an increase in abundance on community composition across time using permutational multivariate analysis of variance (PERMANOVA) with the ‘adonis’ function in the ‘vegan’ package in R. In this model, treatment and time were treated as explanatory variables. To examine the extent to which the communities converged through time or diverged among treatments, Mann–Whitney tests were performed to compare pairwise dissimilarities of communities between time points or treatments.

Alpha diversity was measured as the number of observed OTUs and the Shannon index (\(- \mathop {\sum }\nolimits_{i = 1}^{S_{{\mathrm{obs}}}} \frac{{n_{i}}}{N}\ln \frac{{n_{i}}}{N}\), where Sobs: the number of observed OTUs, N: the total number of sequences in the community, and ni: the number of sequences of the ith OTU). Two-sample t-tests were used to compare alpha diversity between the control and each treatment.

To illustrate how manipulated species affected individual OTUs and shaped community composition, the differences of relative abundance of 188 common OTUs between treatment and control were tested by two-sample t-tests. The extent of each OTU’s relative abundance change was calculated as \(\frac{{R_{{\mathrm{iT}}} - {R}_{{\mathrm{iC}}}}}{{R_{{\mathrm{iC}}}}}\), where RiT: the average of ith OTU’s relative abundance in treatment, RiC: the average of ith OTU’s relative abundance in control. Both the statistical significance and extent of relative abundance change were visualized on the phylogenetic tree with the online tool iTOL [47]. To measure the overall phylogenetic shift in community composition, we calculated the weighted phylogenetic distance between manipulated species and the remaining species in the same communities (WPDj = \(\mathop {\sum }\nolimits_{i = 1}^{S_{{\mathrm{obs}}}} \frac{{n_{i}}}{N}{\mathrm{PD}}_{ij}\), where PDij: the phylogenetic distance between the ith OTU in the community and the manipulated species j). A larger WPDj in the corresponding treatment than in the control group indicates that the manipulated species caused a decrease in the relative abundance of its closely related species, and a smaller WPDj indicates that the manipulated species caused a decrease in the relative abundance of its distantly related species. The comparison of WPDj between the corresponding treatment and control group was conducted by a two-sample t-test. All statistical analyses were performed with the ‘vegan’ and ‘phyloseq’ packages in R [48], unless otherwise noted.

Results

Profile of the six manipulated bacterial species

The six chosen bacterial species were classified into six distinct genera: Chryseobacterium (species A), Acinetobacter (species B), Arthrobacter (species C), Microbacterium (species D), Pseudomonas (species E), and Bacillus (species F). According to NCBI database, their closest relatives are Chryseobacterium elymi, Acinetobacter calcoaceticus, Arthrobacter humicola, Microbacterium foliorum, Pseudomonas silesiensis, and Bacillus proteolyticus, respectively (Table S1). The GenBank accession numbers for 16S rRNA gene nucleotide sequences of these six manipulated species are MK063785, MK063786, MK063787, MK063788, MK063789, and MK063790. Five of the six bacterial species (species A, B, C, E and F) were detected in the soil bacterial community via 454 pyrosequencing (Fig. 1). The range of these five manipulated species’ relative abundance was from 0.088 to 10.48% in treatment microcosms after 10 days of the inoculations (Fig. 1), and the estimated survival rate ranged from 0.16 to 24.3% (Supplementary Methods S4 and S5, Table S2). Species D was undetected in all microcosm communities, possibly due to its extremely low relative abundance and survival rate in soil [49].

Fig. 1
figure 1

The relative abundance (mean ± se) of manipulated bacterial species in corresponding treatments (open circle) and the control (solid circle) before and after recovery (a–f). The relative abundance was calculated as the number of manipulated species’ sequence reads divided by the total number of sequences in each community. Greek letters denote significant differences in the change in relative abundance (ANOVA, p < 0.05) between the treatment and control group. Note that species D was not detected by 454 pyrosequencing in the soil bacterial communities

Full size image

Among the species detected, four species (species A, B, C and E) showed significant increases in the relative abundance in treatment communities compared with the control communities both before and after in situ recovery (Fig. 1, Table S3; ANOVA, treatment; FA1,26 = 114.83, FB1,26 = 25.64, FC1,26 =  1562.99, FE1,26 = 62.12; all p < 0.001). However, the relative abundance of species F in the treatment microcosms was not significantly different from that in the control (Fig. 1).

The species-specific divergence of soil bacterial communities

The treatment of species abundance changes was a primary factor driving the divergence of soil bacterial communities, which explained 36.2% of the variation in community composition (Fig. 2, Table S4; PERMANOVA; treatment, R2 = 0.362, F6,66 = 10.638, p < 0.001, permutations = 9999). This pattern was described as ‘species-specific divergence’, because each of the manually increased resident bacterial species exerted specific effects on community composition (Tables S4 and S5), and the divergence pattern was further strengthened after a year of recovery (Fig. 2, Table S4, S5, and S6; PERMANOVA; treatment × time, R2 = 0.091, F6,66 = 2.687, p < 0.001, permutations = 9999). In addition to divergence among treatments, a single episode of species abundance change also gave rise to greater pairwise dissimilarities within each treatment than in the control group (Fig. 3).

Fig. 2
figure 2

Metric multidimensional scaling ordination of soil bacterial communities before and after recovery by principal coordinate analysis (PCoA). Each point represents an experimental soil bacterial community. Distances among communities represent weighted UniFrac dissimilarities. The communities are colored based on the treatments, and the shapes of symbols distinguish the communities at different times (before or after recovery)

Full size image
Fig. 3
figure 3

Pairwise dissimilarities (weighted UniFrac index) between communities. The lower and upper ends of the boxes represent 25% and 75% of the range, respectively; lines in the boxes indicate medians; triangles indicate maximums and minimums; squares indicate mean values; and whiskers represent ±1.5 × the interquartile range (IQR, defined as the upper quartile minus the lower quartile). The dissimilarities are separated by time (black box: before recovery, red box: after recovery) and grouped as follows: between communities within specific treatments (A, B, C, D, E, F and control, where N represents the control) and among all communities in our experiment (All). Comparisons between two time points and between each treatment and the control group were conducted by Mann–Whitney tests, and significance is shown in the table below the figure (***p < 0.001, **p < 0.01, and *p < 0.05)

Full size image

There was also a trend of community convergence over time, shown as a significant decline in pairwise dissimilarities both between all communities and between replicates in four out of six treatment groups after a year of recovery (Fig. 3). However, the convergent trend within groups B, D, E, and F was much stronger than the overall trend (Fig. 3), which contributed to a greater divergence between treatments after a year of recovery (Table S4, S5, and S6).

The relative-abundance-change patterns of manipulated species in communities during a year of in situ recovery

The manually increased species showed diverse relative-abundance-change patterns during the year of recovery (Fig. 1). The relative abundance of species A, B, and C significantly declined throughout the year (Fig. 1, Table S3; ANOVA, time; FA1,26 = 56.82, pA < 0.001; FB1,26 = 11.44, pB = 0.002; FC1,26 = 18.03, pC < 0.001), and the trend of change in relative abundance in the treatments was significantly steeper than that in the control (Fig. 1, Table S3; ANOVA, time × treatment; FA1,26 = 113.54, FB1,26 = 22.78, FC1,26 = 22.53; all p < 0.001). The relative abundance of species E significantly increased throughout the year (Fig. 1, Table S3; ANOVA, time; FE1,26 = 5.05, pE = 0.033), and the trend of change in relative abundance in the treatment E was also significantly steeper than that in the control (Fig. 1, Table S3; ANOVA, time × treatment; FE1,26 = 9.77, pE = 0.004). Note that the relative abundance of species F in the corresponding treatment was not significantly different from that in the control before or after in situ recovery (Fig. 1 and Table S3).

The patterns of local community compositional change underlying divergence

At the OTU level, increase in the abundance of single species affected a wide range of OTUs spreading across the entire tree of common OTUs (Fig. 4 and Figs. S2–S6). However, the OTUs whose relative abundances were significantly different between the treatment and control microcosms (p < 0.000264, Bonferroni-corrected p value threshold) were few, ranged from 0 to 15 species (Fig. 4 and Figs. S2–S6). Considering that abundance change of rare species was very vulnerable to stochasticity, only common species were included in this analysis. The common species were chosen with the criterion that their relative abundance was >0.1%. Pairwise compositional dissimilarities (weighted UniFrac) calculated based on full dataset and the dataset with common species were highly correlated (Supplementary Methods S6, Fig. S7; Mantel test; spearman ρ = 0.953, p = 0.001). And the dataset consisting of 188 common species can adequately capture the overall compositional variation among bacterial communities (Supplementary Methods S6, Fig. S8, and Table S7).

Fig. 4
figure 4

Comparisons of the relative abundances of 194 bacterial OTUs between treatment A and control, shown with the phylogenetic relationship among those bacterial OTUs (188 common species and the six manipulated species). a, b show the observations before and after the one-year in situ recovery, respectively. Phyla are indicated by branch colors and ring 1, and the manipulated species A is indicated by the branch with red tip. Ring 2 indicates each common OTU’s average relative abundance in control group. The red strips in ring 3 indicate significance of two-sample t-tests of relative abundance between control and treatment group, 0.000264 is the Bonferroni-corrected p value threshold for multiple comparisons. Ring 4 indicates the extent of each common OTU’s relative abundance change, calculated as \(\frac{{R_{{\mathrm{iT}}} - {R}_{{\mathrm{iC}}}}}{{R_{{\mathrm{iC}}}}}\), where RiT: the average of ith OTU’s relative abundance in treatment A, RiC: the average of ith OTU’s relative abundance in control. Phylogeny and relative abundance change in other treatments are provided in Figs. S2–S6

Full size image

At the community level, increases of abundance of a single species resulted in a significant decrease in alpha diversity (Fig. 5). However, underlying the consistent pattern of alpha diversity decrease were different shifts in community composition among treatments. The  WPD index was employed to summarize the overall trends in abundance change of species. An increase in the abundance of species A, B, E, and F led to a higher WPD in the corresponding treatment than in the control (Fig. 6), indicating that manually increasing these species led to a decrease in the relative abundance of their closely related species. In contrast, an increase in the abundance of species C resulted in a lower WPD compared with the control, indicating a decrease in the relative abundance of its distantly related species in the corresponding treatment (Fig. 6).

Fig. 5
figure 5

Alpha diversity (mean ± se), number of observed OTUs (top panels) and Shannon index (bottom panels) in each treatment and the control before (left panels) and after recovery (right panels). The alpha diversities in the treatments were all significantly lower than those in the control (two-sample t-test, all p < 0.05)

Full size image
Fig. 6
figure 6

Weighted phylogenetic distances (WPD, mean ± se) between manipulated species and the remaining species within the same communities in each corresponding treatment and control (a–f). A value larger than that in the control suggests that the manipulated species caused a decrease in the relative abundance of closely related species and vice versa. The comparisons between treatments and the control were conducted by two-sample t tests (*p < 0.05)

Full size image

Discussion

As a supplement to major studies focusing on the roles of environmental heterogeneity [13, 25] and dispersal [26, 50] in shaping soil bacterial communities, our study found that local biotic interactions allowed community divergence presumably via interspecific competition. Specifically, the increase in abundance of one species could result in a decrease in the relative abundance of its closely related species, thus driving a species-specific compositional shift in the soil bacterial community. In addition, local processes did not universally lead to effective population regulations that force species relative abundance to revert to the original level. As a consequence, despite an overall trend of community convergence, there was still substantial variation between communities in different treatments after 1 year of in situ recovery.

We found that a single episode of abundance-increasing manipulation induced immediate disturbance in community structure (Fig. 2), accompanied by a significant decrease in alpha diversity (Fig. 5). Competition exerted by the manipulated species was the most likely reason for this decrease in alpha diversity, which further resulted in substantial variation in community composition between treatments. Although the manipulated species led to the loss of ~100–300 rare OTUs (Fig. 5), the compositional changes of communities in treatments were primarily rooted in the relative abundance changes in the common OTUs (Supplementary Methods S6, Figs. S7 and S8). Unexpectedly, only a few number (between 0 and 15 out of 188) of common species’ relative abundance in treatments were significantly affected by the inoculation of manipulated species compared with that in control (Fig. 4, Figs. S2–S6). Therefore the observed community divergence was resulted from collective subtle relative abundance changes of many community members. These results implied that influence of the manipulated species was not constrained within a few specific OTUs, but covered extensive range of OTUs. Such a pattern was consistent with a suggestion that most bacterial species share, and compete for, a number of common resources [51]. Another possible explanation is that the inoculated species affected community composition through complex biotic interaction network, as suggested by previous experimental evidence that antagonistic interactions among bacteria exhibited hierarchical structure [52]. Unlike the other manipulated species, Microbacterium sp. was undetected in all microcosms. Nevertheless, the inoculation of Microbacterium sp. (treatment D) still significantly affected community composition (Tables S4–S6). Such ‘unsuccessful but effective invasion’ events were also reported by an earlier experiment where Escherichia coli failed to invade the resident soil bacterial communities but caused shifts in the composition of resident communities presumably due to the transient exploitive competition [53].

The influence of species abundance manipulation on resident species could be described as ‘small changes in many OTUs’. The WPD index was calculated to describe the overall phylogenetic shift of bacterial communities in our experiment, we found that four of the six manipulated species, namely, Chryseobacterium sp., Acinetobacter sp., Pseudomonas sp., and Bacillus sp., caused a significant decrease in the relative abundance of their closely related species compared with the control (Fig. 6). We may cautiously suggest that interspecific competition was the primary local process driving the species-specific divergence in the soil bacterial community. This is consistent with previous suggestions that competition dominated bacterial interactions [54], and that closely related species were more likely to have similar ecological characteristics [36,37,38] thereby resulting in more intense competition for resources and interference competition [39, 40]. One notable exception is that an increase in the abundance of Arthrobacter sp. led to a decrease in the relative abundance of its distantly related species compared with the control (Fig. 6). One possible explanation is that there were synergistic interactions between Arthrobacter sp. and its closely related species, such as sharing public goods [55], which indirectly influenced community composition via competition exerted by species closely related to Arthrobacter sp.

While abundance increases of particular species triggered species-specific community divergence, there was a trend of community convergence both among all communities and within treatments over time (Fig. 3). This trend ruled out the possibility that stochastic processes predominately controlled community assembly, as stochastic processes would only result in further divergence over time [10]. Nevertheless, this convergence did not suggest that there was a single stable state for soil bacterial communities in our experiment. Indeed, the effects of species abundance change on community divergence were stronger after 1 year of recovery, and the convergent trend within treatments was much stronger than that among treatments (Tables S4–S6). This result implied that there might be alternative stable states, or more likely alternative transient states in soil bacterial communities [56]. The occurrence of alternative community states across homogenous habitat patches reduces environment-community matching, and may in turn exaggerate the unexplained variation in community composition, which could be falsely attributed to stochastic processes.

The lack of effective population regulations forcing species relative abundance to revert to the original level in our experiment suggested that the community divergence would persist. Stabilizing processes (negative frequency dependence) are thought to be the foundation of stable coexistence at local scales [41, 42]. The relative abundance of Chryseobacterium sp., Acinetobacter sp. and Arthrobacter sp. significantly decreased in the treatment during recovery compared with the control, which is consistent with negative frequency dependence; however, their relative abundance was still significantly higher than that in control. Meanwhile, the performance of Pseudomonas sp. clearly opposed the prediction of negative frequency dependence, which even became more abundant after the 1 year in situ recovery. Pseudomonads have the features of biofilm formation, toxin secretion, and siderophore production [34, 57, 58], allowing them to gain more fitness when their abundance is greater. This result is consistent with a scenario of positive frequency dependence (or positive feedback), which is a destabilizing process that can promote ASS [18, 24, 59,60,61]. Taken together, our results suggest that stabilizing processes are far from universal in soil bacterial community. Combined with the possibility of alternative community states in soil bacterial communities we noted above, it is reasonable to suggest that the community divergence observed in our experiment was not transitory.

The roles of biotic interactions in driving communities to diverge have been extensively studied within the theoretical framework of priority effects, in which early immigrants affect the fitness of later immigrants; variation in immigration order in the early stage of community assembly may result in different community compositions [62,63,64]. Most previous explorations of priority effects were restricted to the initial stage of community assembly [62,63,64], but the underlying ecological processes may extend to established communities. Our study provided solid evidence that established soil bacterial communities might diverge in response to changes in the abundance of a resident species due to interspecific interactions. Therefore, the complex biotic interactions between bacteria can be a continuous force shaping soil bacterial community structure.

The manipulated species’ survival cells per gram of soil after 10 days of preincubation were estimated to be 105–107 (Supplementary Methods S4 and S5, Table S2), except Microbacterium sp. undetected. Compared with the previous invasion experiment of soil with E. coli [65,66,67], the survival rates of the inoculated species were relatively higher in our experiment, especially for Chryseobacterium sp. A possible explanation was that species inoculated in our experiment were all resident species and preadapted to the local environmental conditions, unlike the alien invasion scenario in previous studies. The abrupt abundance increase scenario might be rare in natural soil bacterial community, but not impossible. For example, a field study found a maize-induced boom of Pseudomonas spp. in the rhizosphere microbiome, such that the relative abundance of three Pseudomonas OTUs increased from ~2.7 to 44.5% during 1 week [68]. Although the abundance manipulation treatment in this study did not ideally mimic natural events, our results suggest influences of soil bacterial abundance changes due to dispersal, environmental change, and/or disturbance can transmit through the biotic interaction network, and produce profound and lasting community-wide impacts. Note that mass immigration, as in our manipulated microcosms, inevitably results in substrate provision (due to decomposition of the introduced organisms), of which the contribution to community composition changes was not investigated here but deserves attention in future studies.

Conclusions

Our study demonstrates that local biotic interactions, especially interspecific competitions, might give rise to divergence in soil bacterial communities, and that such divergence will probably persist due to a lack of effective population regulations in soil bacterial communities. These results extend knowledge on the processes that can result in variation in community composition, highlight the importance of local biotic interactions, and might help explain the substantial variation in composition of the soil bacterial community at small scales, which was previously thought to be the consequence of stochastic processes or dispersal limitation. Further studies including the synthesis of environmental filtering, local biotic interactions and dispersal will shed new light on the mechanisms shaping soil bacterial communities.