A quantitative framework reveals ecological drivers of grassland microbial community assembly in response to warming

1Microorganisms are the most diverse groups of life presently known, inhabiting almost every imaginable environment on the Earth, and they typically form complex communities whose structure, functions, interactions, and dynamics are critical to ecosystem functioning and service. However, characterizing such complex communities, quantifying the accompanying ecological processes, and dissecting the mechanisms controlling biodiversity and community composition are extremely challenging. With the rapid development of high-throughput metagenomic technologies, large experimental data on community structure can be rapidly obtained. However, analyzing such massive data to address fundamental ecological questions, such as community assembly mechanisms is challenging due to various issues associated with detection specificity, sensitivity, quantification, reproducibility, and taxonomic resolution.

2Understanding community assembly rules is a longstanding issue of ecologists,. Niche-based theory asserts that deterministic processes, including environmental filtering (e.g., pH, temperature, moisture, and salinity) and various biological interactions (e.g., competition, facilitation, mutualisms, and predation), largely control the patterns of species composition, abundance, and distributions,. By contrast, neutral theory assumes that all species are ecologically equivalent, and species dynamics are largely controlled by stochastic processes of birth/death, speciation/extinction, and immigration,. After intensive debates in 2000s–, it is generally accepted that both deterministic and stochastic processes operate simultaneously in the assembly of local communities–, and the key question becomes how to define their relative importance in controlling community diversity, distribution, and succession,–.

3To unify niche and neutral perspectives on governing community structure, Vellend, proposed a conceptual framework that community diversity and dynamics are controlled by four high-level general ecological processes: selection, dispersal, speciation or diversification, and ecological drift,. Hereafter, we use the term ‘ecological processes’ particularly to represent these community assembly processes. Although the framework has recently received a great attention in microbial ecology,–, translating this conceptual framework into a quantitative operational model is even more challenging–,. Due to the lack of quantitative approaches, most analyses with respect to the relative importance of the four processes across different types of natural communities are qualitative and subjective, and replete with great uncertainty. As an exploratory effort, a null modeling-based operational approach was developed to obtain quantitative information on community assembly processes from the statistical perspective,, which is abbreviated as QPEN (Quantifying assembly Processes based on Entire-community Null model analysis) hereafter. QPEN uses phylogenetic metrics to infer selection since phylogenetic distance could reflect niche difference (so-called phylogenetic signal) within some threshold,.

4This statistical approach represents a significant advance in microbial ecology that enables microbial ecologists to obtain quantitative information on community assembly processes. It has provided valuable insights into the importance of various ecological processes in microbial ecology,,– and plant ecology. However, a major limitation is that various ecological processes are estimated based on the pairwise turnovers of the whole communities,. This may not be appropriate because it is well known that the actions of various ecological processes (e.g. natural selection) are typically on the finer biological organization levels, such as genotypes and populations rather than whole communities,,,,. Within a single microbial community, certain populations are under strong selection, whereas others could be under strong drift. This type of difference cannot be discerned using whole community level metrics. Also, various groups of organisms differ greatly in their responses to environmental changes. Similarly, the dispersal ability, diversification rates, and susceptibility to drift are substantially different among various microbial groups. Thus, it would be meaningful to consider selection and other ecological processes at the level of individual taxa/lineages rather than the entire community,. To this end, we developed a general framework to quantitatively infer Community Assembly Mechanisms by Phylogenetic-bin-based null model analysis, abbreviated as iCAMP, based on the turnovers of individual bins across communities (samples). We apply this approach to investigate whether and how experimental warming affects various ecological processes in the assembly of grassland soil microbial communities. Our results indicate that iCAMP provides a robust, reliable tool for quantifying the relative importance of ecological processes in controlling microbial community diversity and succession.

5To quantify various ecological processes, the observed taxa are first divided into different groups (‘bins’) based on their phylogenetic relationships (Fig. 1, Supplementary Fig. 1a). Then, the process governing each bin is identified based on null model analysis of the phylogenetic diversity using beta Net Relatedness Index (βNRI), and taxonomic β-diversities using modified Raup–Crick metric (RC) (Fig. 1, Supplementary Fig. 1b). For each bin, the fraction of pairwise comparisons with βNRI < −1.96 is considered as the percentages of homogeneous selection, whereas those with βNRI > +1.96 as the percentages of heterogeneous selection based on the threshold applied previously,. Next, taxonomic diversity metric RC is used to partition the remaining pairwise comparisons with |βNRI| ≤ 1.96. The fraction of pairwise comparisons with RC < −0.95 is treated as the percentages of homogenizing dispersal, while those with RC > + 0.95 as dispersal limitation,. The remains with |βNRI| ≤ 1.96 and |RC| ≤ 0.95 represent the percentages of drift, diversification, weak selection, and/or weak dispersal, hereafter, simply designated as ‘drift’ for convenience. The above analysis is repeated for every bin. Subsequently, the fractions of individual processes across all bins are weighted by the relative abundance of each bin, and summarized to estimate the relative importance of individual processes at the whole community level (Fig. 1, Supplementary Fig. 1e). Furthermore, various statistical analyses, e.g. Mantel test, multiple regression on distance matrix (MRM), variation partitioning, etc., are used to reveal the linkages of individual processes to different environmental factors for obtaining detailed insights into community assembly mechanisms (Fig. 1). Besides βNRI and RC, iCAMP also incorporated direct test based on null model distribution, which provided highly similar results as βNRI and RC in this study (Supplementary Fig. 17), and should be a preferred choice when the null model simulated values do not follow normal distribution (see Supplementary Note 1 for details).

6Due to lack of a gold-standard experiment to establish the true community assembly processes, using simulated communities is the predominant approach for assessing performances of various computational methods, as shown in many other computational studies,. Thus, we used a simulation model to generate communities with pre-defined relative importance of each ecological process (Supplementary Fig. 2). In this simulation model, four plots (LA, LB, HA, HB) in two islands (A, B) under two types of environments (L, H) are considered, with six local communities sampled from each plot (Supplementary Fig. 2a). Each local community consists of different types of species controlled by drift, selection, or dispersal (Supplementary Fig. 2c–e).

7Three scenarios were simulated with three different levels of phylogenetic signal in the regional species pool: low (Blomberg’s K = 0.15), medium (K = 0.9), and high (K = 5.5). Each scenario has 15 simulated situations, where the expected (‘true’) relative importance of each process (selection, dispersal, drift) was set from 0% to 100% (with an interval of 25%, Supplementary Fig. 2b and Supplementary Table 1). Then, the relative importance of different ecological processes was assessed by iCAMP. The performance was evaluated with six quantitative and qualitative indexes: quantitative accuracy and precision, and qualitative accuracy, precision, sensitivity, and specificity (detailed in the Methods section, Eqs. (16)–(21)). The quantitative performance indexes are based on the difference between the expected and estimated relative importance of each process, while the qualitative performance indexes are calculated from the true or false identification of dominant process. According to the performances with simulated communities, we optimized the binning algorithms (Supplementary Figs. 3 and 4), metrics (Supplementary Fig. 5), and null model algorithms (Supplementary Fig. 6), and explored the impact of randomization times and resampling taxa on the performance (Supplementary Figs. 7 and 8; detailed in Supplementary Note 1).

8After appropriate parameter settings were determined, iCAMP and several previously reported approaches were compared for their performances with the simulated communities (Fig. 2, Supplementary Figs. 9–12). First, the ecological stochasticity was quantified by five approaches, including abundance-weighted neutral taxa percentage (NP) based on neutral-theory model, normalized stochasticity ratios (NST) based on taxonomic (tNST) or phylogenetic metrics (pNST), and the relative importance of stochastic processes (homogenizing dispersal, dispersal limitation, and drift) based on QPEN or iCAMP. Under high-phylogenetic-signal and medium-phylogenetic-signal scenarios, iCAMP consistently showed the highest quantitative accuracy (0.978–0.997) and precision (0.903–0.930), while pNST exhibited similar accuracy (0.924–0.954) but lower precision (0.658–0.722, p < 0.001, Fig. 2a–c). Under low-phylogenetic-signal scenario, iCAMP continued to show the highest precision (0.807) and the second-high accuracy (0.770), while pNST showed similar precision (0.723) and the highest accuracy (0.947). By contrast, tNST, NP, and QPEN showed lower precision (<0.57, down to −0.75, p < 0.0001) than iCAMP in simulated scenarios (Fig. 2a–c).

9Only QPEN and iCAMP can quantify relative importance of different ecological processes, so we compared their performances. On average, iCAMP had higher accuracy (0.93–0.99 against 0.81–0.97, 9.9% higher), precision (0.82–0.94 against 0.33–0.52, 120.2% higher), sensitivity (0.83–0.94 against 0.54–0.58, 61.1% higher), and specificity (0.96–0.98 against 0.87–0.88, 10.6% higher) than QPEN (Fig. 2i, Supplementary Fig. 9k, l). The setting of phylogenetic signal also had significant impacts on iCAMP performance. When the phylogenetic signal increased from low/medium to high (Fig. 2i, Supplementary Fig. 9k, l), the accuracy and specificity of iCAMP remained high (>0.92) without significant changes (p > 0.20), but the precision and sensitivity of iCAMP increased from 0.80–0.82 to 0.90–0.94. By contrast, the overall performance of QPEN was improved by higher phylogenetic signal. In addition, iCAMP demonstrated good robustness to the uncertainty in bin determination (Supplementary Figs. 10 and 11, see Supplementary Note 2).

10The performance varied considerably among different ecological processes (Fig. 2d–h, Supplementary Fig. 12). In the simulated communities under medium and high phylogenetic signals, all performance indices were higher than 0.78 for iCAMP (Supplementary Fig. 12), indicating considerable improvement from QPEN, particularly in estimating homogeneous and heterogeneous selections. However, with low phylogenetic signal, iCAMP had low sensitivity (down to 0.17) for homogeneous selection (Supplementary Fig. 12), albeit still higher (p < 0.05) than QPEN (sensitivity < 0.04). These results confirmed that low phylogenetic signal of niche preference can limit the capability of phylogenetic metrics to infer selection, which can be partly but not completely overcome by iCAMP. Nevertheless, the quantitative performance of iCAMP remained relatively high for all processes under all scenarios, with quantitative accuracy and precision 0.71–1.00 (averagely 129% higher than QPEN), indicating that iCAMP can substantially improve the quantitative estimation of community assembly processes.

11To determine the effectiveness of iCAMP in real-world studies, iCAMP was applied to an empirical data of soil bacterial communities in a grassland under experimental warming, with focus on within-treatment spatial turnovers. Based on iCAMP analysis, homogeneous selection and drift were more important than other processes in bacterial community assembly, with average relative importance of 37.0–38.5% and 58.3–59.9% (Fig. 3a, b), respectively. Warming significantly altered the relative importance of different processes (p < 0.01, permutational ANOVA). Since other processes had quite low estimated relative importance (<3.4%), we primarily focused on the effects of warming on homogeneous selection and drift in subsequent analyses. Overall, warming decreased the relative importance of drift and increased homogeneous selection. Significant year-to-year variations were observed (Fig. 3c, d). In the first year, the communities under warming showed significantly higher ratio of drift (Cohen’s d = 2.9, p = 0.001), but lower ratio of homogeneous selection (Cohen’s d = −2.7, p < 0.001) than those under control, suggesting the bacterial community assembly was even more stochastic under warming than control in the beginning. In the second year, the difference between warming and control became insignificant. In the third to fifth years, the communities under warming had significantly higher ratio of homogeneous selection (Cohen’s d = 0.6–1.7) and lower ratio of drift (Cohen’s d = −0.8 to −1.3), suggesting that the selection pressure imposed by warming on the soil bacteria gradually increased with time.

12QPEN was also applied to quantify the ecological processes. The results from QPEN indicated that homogeneous selection predominated (>73%) bacterial assembly with higher relative importance under warming (83.3%) than control (73.3%, Supplementary Fig. 13a, b) although not significant (p = 0.174). QPEN suggested 0.0% of drift under warming, 0.0% of heterogeneous selection, and homogenizing dispersal across all years, and 100% of homogeneous selection in some years (Supplementary Fig. 13c, d). This appears not reasonable, considering important roles of stochastic processes have been widely reported across various ecosystems,,,,, including desert and acidic soils.

13Based on the principle of the null models employed by iCAMP and QPEN, the fractions of dispersal limitation, homogenizing dispersal, and drift are largely considered stochastic. Thus the sum of their estimated relative importance can be used to estimate stochasticity of community assembly. Based on iCAMP results, the relative importance of stochastic processes was 62.6% under control and 61.3% under warming (Supplementary Fig. 14). By contrast, QPEN estimated the relative importance of stochastic processes was 26.7% under control and 16.7% under warming, which were much lower than those estimated by other approaches (Supplementary Fig. 14). For instance, variation partitioning analysis (VPA) revealed that substantial portions of the community variations (68.4%) could not be explained by all measured environmental variables. The tNST and pNST were on average 48.8% under warming and 52.3% under control, and NP ranged from 74% to 79% in different years for both warming and control (Supplementary Fig. 14). It appears that VPA, NST, and NP showed more consistent results with iCAMP than QPEN.

14All approaches did not reveal significant (p > 0.10) differences of the 5-year mean stochasticity between warming and control, except tNST with medium effect size (p < 0.05, Supplementary Fig. 14). But in the third to fifth year (Fig. 3e), both tNST and iCAMP revealed that warming had significant (p < 0.05) decrease in stochasticity, and there was slight decrease in stochasticity with pNST and NP under warming (small effect size) though it was insignificant (p > 0.10). On the contrary, QPEN showed a slight but insignificant increase in stochasticity. Collectively, consistent with our previous analysis, various approaches supported that stochastic processes could play more important roles in grassland soil bacterial assembly and that warming decreased the stochasticity after 3 years.

15In contrast to QPEN and other approaches, iCAMP can provide information on the relative importance of different ecological processes in individual lineages (bins). For this purpose, the observed 18,123 OTUs were divided into 658 phylogenetic bins, each of which was then analyzed separately as outlined in Fig. 1. Our results revealed that homogeneous selection dominated 59 bins (9% of bin numbers and 33% of relative abundance, Fig. 4a). Two of the major bins were Bacillales (Bin 1, 26.7% in total abundance of bins controlled by homogeneous selection) in Firmicutes and Spartobacteria (Bin 2, 18.8%) in Verrucomicrobia (Fig. 4b, Supplementary Fig. 15a). By contrast, drift dominated 598 bins (91% of bin numbers and 67% of relative abundance, Fig. 4a), which mainly belonged to Class Alphaproteobacteria (22.2% in total abundance of drift-controlled bins) and Phylum Actinobacteria (23.5%, Supplementary Fig. 15a).

16To understand how different lineages respond to warming, we further determined the bacterial groups contributing to the warming-induced changes of homogeneous selection and drift in the third to fifth years (Fig. 4c, d). Our results revealed that Firmicutes contributed 58.2% of the warming-induced increases in homogeneous selection (Supplementary Fig. 15b). The most abundant Firmicutes bin (Bin 1, Bacillales, average 74.8% in Firmicutes) was governed by homogeneous selection (Supplementary Fig. 15c). After Year 3, which had severe drought, Firmicutes were significantly more abundant under warming than control (Supplementary Fig. 15c). 1 by contrast, the decrease of drift under warming was due to similar negative responses of many bins in five phyla (Proteobacteria, Verrucomicrobia, Bacteroidetes, Planctomycetes, and Acidobacteria, Fig. 4c, Supplementary Fig. 15b). For instance, Bin 4 of Rhizobiales in Alphaproteobacteria had lower relative abundance and reduced relative importance of drift under warming, especially in later 3 years (Supplementary Fig. 15d). From Year 3 to Year 5, 27 bins showed significantly higher relative abundances under warming than control (Wilcoxon p < 0.05, Cohen’s d > 0.5). However, 81% of them (22 bins) were mainly governed by drift, while two of them switched from governed by drift in controls to homogeneous selection under warming. These results demonstrated complex assembly mechanisms of different bins in response to warming, and possible misinterpretation could occur without bin-level mechanistic information.

17The influences of environmental factors on various ecological processes were further determined with Mantel test and MRM. Since the relative importance of other processes (heterogeneous selection, homogenizing dispersal, and dispersal limitation) was small, we primarily focus on homogeneous selection and drift. Overall, these two processes were significantly linked to the environmental factors related to water, temperature, and plants (Fig. 5, Supplementary Fig. 16, Supplementary Tables 2 and 3). Such linkages were greatly changed by warming (Fig. 5, Supplementary Fig. 16).

18Homogeneous selection under control showed the strongest correlations with the C-4 biomass difference and plant richness (Mantel R2 > 0.33, p < 0.01, Fig. 5a, Supplementary Table 2), and had slightly lower correlations with drought, precipitation, and moisture (R2 > 0.3, p < 0.01). By contrast, homogeneous selection had the highest correlations with drought and precipitation (R2 = 0.52–0.57, p < 0.1) under warming, followed by C-4, total plant biomass, and soil temperature (R2 > 0.32, p < 0.1) (Fig. 5a). Interestingly, the total plant biomass under warming and the difference of C-4 biomass under control had strong correlations with homogeneous selection when the effects of drought or any other factors were controlled (partial Mantel, Supplementary Table 3). Considering potential significant correlations among multiple factors, MRM was further used to determine the contributions of different environmental factors to homogeneous selection. Our result showed that the MRM models were able to explain a large portion of the plot-wise variations of homogeneous selection under warming (R2 = 0.94, p < 0.001, Fig. 5b) and control (R2 = 0.86, p < 0.001, Fig. 5c). The most important variables explaining homogeneous selection were soil temperature in sampling month under warming (partial regression coefficient b = 0.92, p < 0.001, Fig. 5b), and the between-plot difference of annual soil temperature under controls (b = −1.2, p < 0.001, Fig. 5c). Meanwhile, some water-related and plant-related variables also showed significant, although smaller, effects (|b| > 0.16, p < 0.1) on homogeneous selection, including total plant biomass, drought, and moisture under warming (Fig. 5b), and plant richness, C-4 biomass, precipitation, and moisture in controls (Fig. 5c). These results indicated soil temperature, plant productivity, and drought were more effective in explaining the plot-wise variation of homogeneous selection under warming.

19Drift had similar relationships to various environmental factors as homogeneous selection, but with opposite direction, based on Mantel (Supplementary Fig. 16a, Supplementary Table 2) and partial Mantel (Supplementary Table 3) tests. Under control condition, MRM models showed that soil total carbon and nitrogen, various plant biomass differences, and precipitation, had relatively strong effects on drift (|b| > 0.36, p < 0.06, Supplementary Fig. 16c). Under warming, the three most important variables shaping drift included total plant biomass, soil temperature, and drought (|b| = 0.34–0.66, p < 0.005, Supplementary Fig. 16b). Collectively, these results suggested that the environmental factors shaping homogeneous selection and drift were more similar under warming than control.

20Disentangling ecological drivers controlling community assembly is crucial but difficult in ecology, especially in microbial ecology — mainly due to the huge diversity of microorganisms and the difficulties in their detection and quantification,,,,. Although metagenomics and associated technologies have revolutionized microbial ecology research, a great challenge is how to use such massive data to address compelling ecological questions such as community assembly mechanisms. Thus, in this study, we developed a framework, iCAMP, to quantify the relative importance of different ecological processes underlying community diversity and dynamics based on individual phylogenetic groups (bins) rather than the entire community. Various analyses demonstrated that iCAMP improved performance substantially with higher precision, sensitivity, specificity, and accuracy compared with previous approaches. Also, the results from iCAMP indicated important roles of stochastic processes in shaping the grassland soil microbial community with an average of ~60% stochasticity, which are consistent with those from various previous studies,–. The developed framework would provide an effective and robust tool to quantify community assembly processes in microbial ecology towards mechanistic understanding of community diversity and succession. Given that our framework and simulation are general to high diversity communities rather than specific to microbes, it should also be applicable to plant and animal communities with high diversity.

21To quantify selection, phylogeny-based approaches–,,,– require that the phylogenetic distances among taxa reflect their niche difference, i.e. there is phylogenetic signal or niche conservatism,. Although phylogenetic niche conservatism of microbial traits was reported,, the signals were mostly at medium or low levels (i.e. close to or lower than Brownian Motion expectation). Fortunately, significant phylogenetic signals were frequently found within a short phylogenetic distance,,, which was employed by recent microbial studies,,,, particularly in QPEN,. However, QPEN did not always perform well in inferring selection (e.g. in our simulated communities), possibly because it does not distinguish the differential influences of selection on distinct phylogenetic groups. By contrast, iCAMP partitions selection based on individual phylogenetic groups with essential phylogenetic signal embedded, and hence it can greatly improve quantitative performance (accuracy and precision > 0.7) with robustness to low phylogenetic signal across the tree and the complex assembly within individual bins. However, if competition predominates, closely related organisms may lead to strong competitive exclusion, which could disrupt phylogenetic conservatism, and thus significantly decrease the performance of iCAMP as showed in our simulation (Supplementary Fig. 11a–f). But given that the significant phylogenetic signal within short phylogenetic distances is widely observed,,,, the disruption may not happen often under natural settings as theoretically predicted.

22Unraveling the drivers controlling the responses of ecological communities to climate change is a critical topic in ecology and global change biology. Several previous studies demonstrated that climate warming have significant impact on microbial diversity, structure, functional gene composition, and activities,, but the underlying community assembly processes were rarely examined. In our grassland site, experimental warming increased the soil temperature by ~3 °C (ref. ), thus it may gradually impose selective pressure as a deterministic force to decrease stochasticity as evident by our previous studies,. Here, iCAMP further revealed that warming gradually enhanced homogeneous selection and decreased drift in bacterial community assembly, and that the warming-induced selection was enhanced by drought and lower plant biomass, besides the soil temperature. These results are consistent with several previous studies showing that the effects of warming on microbial communities intertwines with precipitation/drought– and plant variables. In addition, our results showed that the warming-enhanced homogeneous selection was mainly attributed to the positive responses of a group of Bacillales in Firmicutes, which are ubiquitous Gram-positive bacteria in nature. Bacillales are endospore-forming and can remain in the dormant state for years. These traits could offer their competitive advantages under selective pressure from increased temperature and drying.

23Although iCAMP has better performance over traditional approaches and provided valuable insights into the ecological processes governing the responses of grassland soil microbial communities to climate warming, there are still some limitations. For instance, one of the fundamental processes, diversification, is important to govern community assembly,,, but it is not explicitly accounted for in iCAMP. Diversification is still embedded with drift, weak selection, and/or weak dispersal in the ‘drift’ part of this framework. Given the lower performance of iCAMP when competition predominates selection, environmental filtering and biotic interactions need to be differentiated from each other in selection. In addition, although built on within-bin beta diversity, iCAMP should generally be able to capture important cross-bin selection (Supplementary Note 3). However, iCAMP might underestimate selection when cross-bin selection does not lead to detectable within-bin difference. Thus, further developments are needed by incorporating null model of evolution to infer the relative importance of diversification, and by integrating functional traits (genes) and network approaches with iCAMP to disentangle biotic interactions from abiotic filtering and capture special cross-bin selection.

24Conceptually, selection under homogeneous abiotic and biotic conditions in space and time is referred to as constant selection or homogeneous selection, by which low phylogenetic compositional variations or turnovers are expected. By contrast, selection under heterogeneous conditions leads to high phylogenetic compositional variations, which is referred to as variable selection, or heterogeneous selection. Similarly, dispersal is also divided into two categories, — homogenizing dispersal and dispersal limitation. The former refers to the situation that high dispersal rate can homogenize communities and hence lead to little taxonomic compositional variations, whereas the later signifies the circumstance that low dispersal rates could increase community taxonomic variations. When neither selection nor dispersal is dominated, community assembly is governed by drift, diversification, weak selection and/or weak dispersal, which is referred to be ‘undominated’ or simply designated as ‘drift’.

25To quantify these processes, iCAMP includes three major steps (Fig. 1). The first step is phylogenetic binning (Supplementary Figs. 1a and 3). Three binning algorithms were compared. One is based on the distance to abundant taxa (Supplementary Fig. 3a). The most abundant (i.e. the highest mean relative abundance in the regional pool) taxon is designated as the centroid taxon of the first bin. All taxa with distances to the centroid taxon less than the phylogenetic signal threshold, ds, are assigned to this bin. The next bin is generated from the rest taxa in the same way. Consequently, a series of bins are generated with strict radiuses less than ds, so-called strict bins. However, some strict bins may have too few taxa to provide enough statistical power for further analysis. Thus, each small bin is merged into its nearest-neighbor bin until all bins reach the minimal size requirement, nmin. The second algorithm is based on pairwise distances (Supplementary Fig. 3b). The first bin consists the most abundant taxon, and all other taxa among which all pairwise distances are lower than ds. The second bin includes the next most abundant taxon among the remaining taxa. This procedure continues until all taxa are classified into different bins. To ensure each bin have enough size (≥nmin), a small bin less than nmin is merged into the nearest neighbor until all bins reach the minimal requirement nmin. The third algorithm is based on phylogenetic tree (Supplementary Fig. 3c). The phylogenetic tree is truncated at a certain phylogenetic distance (as short as necessary) to the root, by which all the rest connections between tips (taxa) are lower than the threshold ds. The taxa derived from the same ancestor after the truncating point are grouped to the same strict bin. Then, each small bin is merged into the bin with the nearest relatives. This procedure is repeated until all merged bins have enough taxa (≥nmin). Although not used in this study, another option is also provided in our tool to omit small bins when they are negligible. However, all binning algorithms require a reliable phylogenetic tree, which might be difficult to construct for highly divergent marker genes such as ITS. In this case, certain special phylogenetic tree construction approaches (e.g. hybrid-gene or constrained phylogenetic tree construction) should be considered.

26The objective of phylogenetic binning is to obtain adequate within-bin phylogenetic signal. To evaluate phylogenetic signal within each bin, the correlation between the pairwise phylogenetic distances and niche preference differences were analyzed by Mantel tests, where niche preference means the niche leading to optimum fitness (or relative fitness reflected by relative abundance) of a taxon. The bins with Pearson correlation coefficient R > 0.1 and p < 0.05 (one tail) are considered as bins with significant phylogenetic signal. In simulated communities, the niche preference difference between two taxa is treated as the key trait value difference. For empirical data, a practical index, i.e., niche value, is estimated as the relative-abundance-weighted mean of an environmental factor for each taxon. For instance, if OTU1 has relative abundances of 10%, 20%, and 10% in three samples under the temperature of 10, 20, and 30 °C, respectively, the temperature niche value of OTU1 is (10 × 10% + 20 × 20% + 3 × 10%)/(10% + 20% + 10%) = 20 °C. Then, the difference of niche values between taxa reflects niche difference, which are used for phylogenetic signal estimation. An optimized nmin should give the highest number of bins with significant phylogenetic signal and relatively high average correlation coefficient (average R) within bins. Accordingly, the optimized nmin is identified as 24 for simulated datasets and 12 for the empirical data. The index ‘niche value’ works under the assumption that relative abundances can represent the fitness of a taxon and the available environmental factors can measure the niche profile. Otherwise, an alternative way should be considered, i.e. choose the nmin value that makes the estimated relative importance of stochastic processes close to the stochasticity assessed by referable approaches (e.g. pNST).

27The second step is the null model analysis within each bin shown in Supplementary Fig. 1b. Accordingly, an operator is defined to count whether a bin is governed by a process. The operator can be calculated from βNRI and RC as Eqs. (1)–(10). Another significance testing index directly based on null model distribution is also provided in Supplementary Note 1.where ‘’ is operator for heterogeneous selection, to count whether the turnover of the kth phylogenetic bin (Bin k) between community u and v governed by heterogeneous selection.,,, and are analogous operators for homogeneous selection, dispersal limitation, homogenizing dispersal, or ‘drift’, respectively. ‘’ is bNRI of Bin k between community u and v. ‘’ is beta mean pairwise distance of Bin k between communities u and v, and ‘’ is the βMPD of the null communities randomized according to a null model. ‘’ is standard deviation. ‘’ is modified RC. ‘Nr’ is total randomization times, usually 1000 times. ‘’ is an operator to calculate RC value. ‘’ is Bray–Curtis dissimilarity index. ‘’ is Bray–Curtis dissimilarity of Bin k between null communities u and v of the nrth time randomization according to a null model. ‘’ and ‘’ represent relative abundance of taxon i in community u or taxon j in community v, respectively. ‘’ represents taxa number in Bin k. ‘’ and ‘’ are abundance of taxon i in communities u and v, respectively. For microbial data from sequencing, it is usually difficult to get accurate estimation of absolute abundances of taxa in a community, thus relative abundances can be used to calculate Bray–Curtis index as a common practice.

28The null model algorithm for phylogenetic metrics is ‘taxa shuffle’,, which randomizes the taxa across the tips of the phylogenetic tree, and thus it randomizes the phylogenetic relationship among taxa. The null model algorithm for taxonomic metric is the one constraining occurrence frequency of each taxon proportional to observed and richness in each sample fixed to observed,. The null model algorithm results heavily depend on the selection of the regional pool, within which randomization is implemented. Thus, the algorithms randomizing taxa within each bin and across all bins were compared in iCAMP analysis for the simulated communities. No matter whether the randomization is within or across bins, the beta diversity metrics are calculated in each bin as defined in Eqs. (6)–(10).

29Null model analysis is most computational resource — and time-consuming, which largely depends on the times of randomization and taxa number. But decreasing randomization times or taxa number can reduce reproducibility of the null model analysis. Considering that most reported null model analyses used 1000-time randomization, iCAMP were performed for simulated data with randomization times ranging from 25 to 5000 and repeated 12 times with each number of randomization times. The results from 60,000-time randomization served as a standard for evaluation. In addition, three methods for reducing taxa number were tested. The method ‘rarefaction’ means to randomly draw the same number of individuals (sequences) from each sample and reduce the taxa number. The method ‘average abundance trimming’ ranks all taxa from abundant to rare according to their average relative abundances across all samples and only keeps the taxa before a certain rank. The method ‘cumulative abundance trimming’ ranks taxa in each sample from abundant to rare, then only keeps the abundant taxa in each sample so that every sample has the same cumulative abundance. The iCAMP results from the three methods were compared to that from the original simulated communities.

30The third step of iCAMP is to integrate the results of different bins to assess the relative importance of each process (Supplementary Fig. 1c–f). Defining neutrality at individual level has been proved a key to successfully develop the unified neutral theory. Therefore, the relative importance of a process can be quantitatively measured as abundance-weighted percentage for each bin (Eq. (11)) or the entire communities (Eqs. (12) and (13)). Qualitatively, for each pairwise comparison between communities (samples), the process with higher relative importance than other processes is regarded as the dominant process.where ‘’ is relative importance of the th ecological process in governing the turnovers of Bin k among a group of communities (e.g. samples within a treatment, a region, etc.; Supplementary Fig. 1d) or between a pair of groups (e.g. between treatment and control, which can be enabled by set ‘between.group’ as TRUE for functions ‘icamp.bins’ and ‘icamp.boot’ in iCAMP package). ‘’ is relative importance of the th ecological process in governing the turnover between communities u and v (Supplementary Fig. 1c). ‘’ is relative importance of the th ecological process in governing the turnovers among a group of communities (Supplementary Fig. 1c) or between a pair of groups. Thus, can be,,,, or for heterogeneous selection, homogeneous selection, dispersal limitation, homogenizing dispersal, or ‘drift’, respectively. ‘’ and ‘’ are total relative abundance of Bin k in community u and community v, respectively. ‘’ is operator counting whether the kth bin is governed by the th ecological process, including,,,, and (Eqs. (1)–(5)). ‘’ is total number of bins. ‘’ is number of pairwise comparisons in a group of communities (e.g. within a treatment) or between a pair of groups (e.g. between treatments). ‘’ is average relative abundance of Bin k in the group of communities.

31As shown in Eq. (13), the relative importance of each process is the sum of the terms, by which we can define the contribution of different bins to (Eqs. (14) and (15)).where ‘’ is Bin contribution to Process, measuring the contribution of Bin k to the relative importance of th ecological process in the assembly of a group of communities (Supplementary Fig. 1e). ‘’ is Bin Relative contribution to Process, measuring the relative contribution of Bin k to the th ecological process (Supplementary Fig. 1f).

32In the simulation model (Supplementary Fig. 2), all samples are from the same region sharing the same metacommunity (the regional species pool) with 20 million individuals. The relative abundances of species in metacommunity are simulated using metacommunity zero-sum multinomial distribution model (mZSM) derived from Hubbell’s Unified Neutral Theory Model, using R package ‘sads’ (version 0.4.2) with J = 2 × 107 and θ = 5000. The whole region has two separated islands of A and B (Supplementary Fig. 2a). For species controlled by dispersal, migration is unlimited within each island but nearly impossible between islands. Each island has two plots: plot LA and HA at island A, and plot LB and HB at island B. The two plots at the same island are under distinct environments. The environment variable is as low as 0.05 in the north plots at each island (LA and LB), but as high as 0.95 in the south plots (HA and HB), which is a critical setting for species under niche selection. At each plot, six local communities are simulated and sampled as biological replicates. Each local community contains 20,000 individuals of 100 species.

33A phylogenetic tree was retrieved from a previous publication, which simulated evolution from a single ancestor to the equilibrium between speciation and extinction and generated a tree with 1140 species. A trait defining the optimal environment of each species (Ei) evolves along the phylogenetic tree with a certain phylogenetic signal. We simulated three pools of species as three scenarios to explore the performance of iCAMP under distinct levels of phylogenetic signals. (i) The low-phylogenetic-signal pool was generated using Stegen’s evolution model. The Blomberg’s K value is as low as 0.15, close to the mean K value of 91 continuous prokaryotic traits. The phylogenetic signal is low if counting the phylogenetic distance across the whole tree. However, the trait still shows significant phylogenetic signal within a short phylogenetic distance, in accordance with general observations in microbial communities in various environments,. (ii) The medium-phylogenetic-signal pool was generated by simulating the trait according to Brownian motion model, using the function ‘fastBM’ in R package ‘phytools’ (version 0.6–99) with an ancestral state of 0.5, an instantaneous variance of Brownian process of 0.25, and the boundary from 0 to 1. The final K value is 0.9, close to the mean phylogenetic signal level of 899 prokaryotic binary traits. (iii) The high-phylogenetic-signal pool was simulated according to Blomberg’s ACDC model with a g value of 2000. The final K value is as high as 5.5, close to the highest phylogenetic signal of prokaryotic traits to date.

34For each scenario, we simulated 15 situations with different levels of expected relative importance of various processes (Supplementary Fig. 2b). The situations can be classified into two types. In the first type, all species under each situation are governed by the same kind of processes, i.e. pure selection, or dispersal, or drift. In each of the other situations, species in the regional pool are assigned to different types controlled by various processes. Once a species is assigned to be controlled by selection or dispersal rather than drift, its nearest relatives within ds will also be assigned to the same type of processes considering the phylogenetic signal of traits. Species controlled by each type of processes are simulated as below. (i) To simulate strong selection due to abiotic filtering without stochasticity, the relative abundance of each species is determined by the difference between the environment variable and their trait values (optimal environment), following a Gaussian function (Eq. (16), Supplementary Fig. 2d).where ‘’ is abundance of species i in local community j. ‘’ is the value of the key environmental variable in local community j, which is 0.05 in the north plots (LA and LB) and 0.95 in the south plots (HA and HB). ‘’ is the optimum environment of species i. ‘’ is the standard deviation, which is 0.015. Consequently, the turnovers of these species under the same environment (i.e. within north plots, or within south plots) are solely governed by homogeneous selection, and those between distinct environments (i.e. between north and south plots) are governed by heterogeneous selection.

35(i) To simulate competition without stochasticity, a geometric series model was modified to consider stronger competition between species with similar niche preference. Competitive species in a local community are ranked from the strongest competitor to the weakest with their relative abundances proportional to 0.5, 0.52, 0.53, …, 0.5h, …. The strongest competitor is randomly selected from species with the best fitness, i.e. from the top 10 species with the lowest |EVu–Ei|. Then, the next competitor is the one with the largest niche difference with prior competitor(s) in the rank, based on abundance-weighted Euclidean trait distance to previous competitor(s) (Eq. (17)). The total relative abundance of species controlled by competition is determined as the designated ratio of competition in selection multiplied by the designated ratio of selection in a simulated situation. The turnovers of these species are defined as governed by selection, without distinguishing between homogeneous and heterogenous selections.where ‘’ is the index to assess niche difference between species i and (h−1) prior competitors in the rank. The species with the highest ndhi will be the hth competitor in the rank, and assigned relative abundance proportional to 0.5h. ‘’ is the optimum environment of species i which is not included in the (h−1) prior competitors. ‘’ is the optimum environment of species j which is the jth prior competitor with relative abundance proportional to 0.5j.

36(ii) To simulate extreme dispersal without selection, we modified Sloan’s simulation model which was derived from Hubbell’s neutral theory model (Supplementary Fig. 2e). Each island has a unique species pool, simulated as a local community under the regional metacommunity following neutral theory model but with a relatively low dispersal rate (m1 = 0.01). However, the unique species pools of the two islands are constrained to have no overlapped species, regarding extreme dispersal limitation between the two islands. Then, the local communities in each island are simulated as governed by neutral dispersal from both the regional metacommunity with a low rate (m1 = 0.01) and the unique species pool of the island with a high rate (m2 = 0.99). It means 99% of dead individuals in a local community are replaced by species from the small island–unique species pool at each time step. Therefore, all the turnovers within an island are governed by homogenizing dispersal, and those between islands are controlled by dispersal limitation.

37(iii) Drift is simulated as neutral stochastic processes at a moderate dispersal rate rather than limited or strong dispersal. To simulate drift, all local communities are generated under neutral dispersal from the regional metacommunity with a medium rate (m1 = 0.5, Supplementary Fig. 2c). Since 50% of dead individuals are replaced by species randomly drawing from a relatively large regional pool, all the turnovers among local communities are neither affected by homogenizing dispersal nor under dispersal limitation.

38Under each situation, the dataset of the 24 local communities is simulated as a combination of species governed by different ecological processes, with ratios defined by the situation setting (Supplementary Table 1, Supplementary Fig. 2b). To simulate complex assembly of bins, the species pool is divided into bins with different bin size limitation (nmin = 3, 6, 12, 24, 48) and phylogenetic distance cutoff (ds = 0.1, 0.2, 0.4), and each bin is simulated as controlled by a certain process. Then, as iCAMP analysis still uses nmin = 24 and ds = 0.2, some estimated bins can have members governed by different processes in the same bin. For each turnover between a pair of local communities, the mean relative abundance of species governed by a process defines the expected relative importance of the process (Eq. (18)). The process with the highest relative importance is the expected dominant process of the turnover. Since dispersal and drift are simulated as pure stochastic processes, the expected stochasticity is defined as the sum of expected relative importance of homogenizing dispersal, dispersal limitation, and drift (Supplementary Table 1).where ‘’ is the expected relative importance of the th ecological process in community turnover between samples u and v. ‘’ is total relative abundance of Bin k in community u. ‘’ is total relative abundance of Bin k in community v. ‘’ is operator, equal to 1 if the turnover of the kth bin between communities u and v is governed by the th ecological process, and equal to 0 if not.

39We simulated three scenarios with different levels of phylogenetic signal, 15 situations per scenario with 1 dataset per situation, thus a total of 45 datasets. In each dataset, we applied both QPEN and iCAMP to estimate the relative importance of different processes (quantitative estimation) and the dominant process (qualitative estimation). QPEN cannot assess relative importance of processes for each turnover, but can estimate their relative importance as the percentage of turnovers governed by the process in all turnovers within a plot (e.g. plot HA) or between a pair of plots (e.g. plot HA vs. HB). Then, the ecological stochasticity of community assembly can be quantified as the relative importance of stochastic processes (i.e. homogenizing dispersal, dispersal limitation, and drift) based on QPEN and iCAMP, respectively. For comparison, the ecological stochasticity in each dataset is also estimated with NP, tNST, and pNST,.

40The performance of quantitative estimation is evaluated by accuracy (Eq. (19)) and precision coefficients (Eq. (20)) derived from concordance correlation coefficient (CCC). The performance of qualitative estimation is assessed with respect to accuracy, precision, sensitivity, and specificity by counting the true and false positive/negative results (Eqs. (21)–(24)).where ‘qACC’ and ‘qPRC’ are quantitative accuracy and precision, respectively. ‘’ is covariance of x and y. In our study, x and y are the expected and estimated stochasticity or relative importance of a process, respectively. ‘’ and ‘’ are variance of x and y, respectively. ‘’ and ‘’ are mean of x and y, respectively.

41In the qualitative performance indexes, ‘’ is accuracy; ‘’ is precision measured as positive predictive value; ‘’ is sensitivity measured as true positive rate; ‘’ is specificity measured as true negative rate. ‘’ is true positive number. A true positive for a process means a turnover is correctly identified as dominated by this process. Overall true positive of a method is calculated as the sum of true positive numbers of all processes. ‘’ is true negative number. A true negative for a process means a turnover is correctly identified as not dominated by this process. Overall true negative is calculated as the sum of true negative numbers of all processes. ‘’ is false positive number. A false positive for a process means a turnover is incorrectly identified as dominated by this process. Overall false positive is calculated as the sum of false positive numbers of all processes. ‘’ is false negative number. A false negative for a process means a turnover is incorrectly identified as not dominated by this process. Overall false negative is calculated as the sum of false negative numbers of all processes.

42For example, a turnover is in fact dominated by drift. If the estimated dominating process is drift, this is a true positive for drift, and a true negative for other processes. If the estimated dominating process is dispersal limitation, this is a false positive for dispersal limitation and a false negative for drift, but a true negative for other processes.

43We applied iCAMP to an empirical dataset from our previous study, with sequencing data available in the NCBI Sequence Read Archive under project no. PRJNA331185. Briefly, the grassland site is located at the Kessler Atmospheric and Ecological Field Station (KAEFS) in the US Great Plains in McClain County, Oklahoma (34°59ʹN, 97°31ʹW). The field site experiment was established in July of 2009. Surface soil temperature in warming plots (2.5 m × 1.75 m each) is increased to 2–3 °C higher than the controls by utilizing infrared radiator (Kalglo Electronics, Bath, PA, USA). Surface (0–15 cm) soil samples were taken annually from four warming and four control plots. A total of 40 samples over 5 years after warming (2010–2014) were analyzed in this study. Soil DNA was extracted by from 1.5 g of soil by freeze-grinding and SDS-based lysis and purified with a MoBio PowerSoil DNA isolation kit (MoBio Laboratories). Then the V4 region of 16S rRNA gene was analyzed by amplicon sequencing on Illumina MiSeq, using the primers 515F (5ʹ-GTGCCAGCMGCCGCGGTAA-3ʹ) and 806R (5ʹ-GGACTACHVGGGTWTCTAAT-3ʹ). Sequencing results were analyzed with our pipeline (http://zhoulab5.rccc.ou.edu:8080) built on the Galaxy platform (version 17.01) and OTUs were generated by UPARSE at 97% identity. Soil properties were analyzed using a dry combustion C and N analyzer (LECO), a Lachat 8000 flow-injection analyzer (Lachat), pH meter, a portable time domain reflectometer (Soil Moisture Equipment Corp.), and constantan–copper thermocouples with CR10x data logger (Campbell Scientific). Plant biomass was measured with a modified pin-touch method and the plant richness was based on identification of all species in each plot. The drought index is calculated as additive inverse of standardized precipitation–evapotranspiration index (SPEI) retrieved from SPEIbase.

44The significance of difference for each evaluation index (e.g. qualitative accuracy, precision, sensitivity, etc.) between different methods was calculated by bootstrapping for 1000 times (one-side test). To assess the effects of warming on ecological processes, the standardized effect size (Cohen’s d) was calculated as the difference of means between warming and controls divided by the combined standard deviation, and the magnitude of effect is defined as large (|d| > 0.8), medium (0.5 < |d| ≤ 0.8), small (0.2 < |d| ≤ 0.5), and negligible (|d| ≤ 0.2) according to Cohen’s d. NST, NP and QPEN, were applied to the dataset. The significance of difference in stochasticity or relative importance of ecological processes between warming and control was calculated by permutational t test (1000 times). The empirical study only investigated within-treatment spatial turnovers at each time point.

45For correlation test between each process and various measured factors, we applied Mantel test, partial Mantel, and MRM with constrained permutation considering repeated measures design of the experiment. For Mantel test, both linear model and general linear model with a logit link function and a ‘quasibinomial’ distribution were tested, and the relative importance of each process and each factor were either log-transformed or not, to explore the best model. To log-transform a factor with zero or negative values, all its values were subtracted by the lowest value and the resulted zero values were replaced by 0.05 (i.e. −3.00 in natural log) of the minimum positive value before natural −log transformation. Partial Mantel was performed on factors with significant correlation with homogeneous selection or drift. For MRM, the factors were log-transformed and standardized to zero-mean and unit-variance, then the best model was forward selected based on Akaike information criterion (AIC). For each measurement (e.g. soil nitrate), both the variation (e.g. |Nitrateu–Nitratev|, where u and v represent samples) and the mean (e.g. [Nitrateu + Nitratev]/2) in each pair of samples were investigated for correlation with the relative importance of each process (e.g. PHoSuv). All statistical analyses were implemented by R (version 3.5.3). All significance tests are two-side unless specified.

46Further information on research design is available in the Nature Research Reporting Summary linked to this article.