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

Sep 18, 2020

### Procedure of iCAMP

Conceptually, selection under homogeneous abiotic and biotic conditions in space and time is referred to as constant selection16 or homogeneous selection20, 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 selection16,20 or heterogeneous selection3. Similarly, dispersal is also divided into two categories19,20 — 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’20 or simply designated as ‘drift’19.

To 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-gene52 or constrained phylogenetic tree construction53) should be considered.

The 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 taxon21. 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).

The 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.

$$W_{{mathrm{HeS}}uvk} = left{ {begin{array}{*{20}{l}} 1 & {beta {mathrm{NRI}}_{uvk} > 1.96} \ 0 & {{mathrm{else}}} end{array}}, right.$$

(1)

$$W_{{mathrm{HoS}}uvk} = left{ {begin{array}{*{20}{c}} 1 & {beta {mathrm{NRI}}_{uvk} < – 1.96} \ 0 & {{mathrm{else}}} end{array}}, right.$$

(2)

$$W_{{mathrm{DL}}uvk} = left{ {begin{array}{*{20}{c}} 1 & {left| {beta {mathrm{NRI}}_{uvk}} right| le 1.96,{mathrm{and}},{mathrm{RC}}_{uvk} > 0.95} \ 0 & {{mathrm{else}}} end{array}}, right.$$

(3)

$$W_{{mathrm{HD}}uvk} = left{ {begin{array}{*{20}{c}} 1 & {left| {beta {mathrm{NRI}}_{uvk}} right| le 1.96,{mathrm{and}},{mathrm{RC}}_{uvk} < – 0.95} \ 0 & {{mathrm{else}}} end{array}}, right.$$

(4)

$$W_{{mathrm{DR}}uvk} = left{ {begin{array}{*{20}{c}} 1 & {left| {beta {mathrm{NRI}}_{uvk}} right| le 1.96,{mathrm{and}},left| {{mathrm{RC}}_{uvk}} right| le 0.95} \ 0 & {{mathrm{else}}} end{array}}, right.$$

(5)

$$beta {mathrm{NRI}}_{uvk} = frac{{beta {mathrm{MPD}}_{uvk} – overline {beta {mathrm{MPD}}_{{mathrm{null}}uvk}} }}{{{mathrm{Sd}}(beta {mathrm{MPD}}_{{mathrm{null}}uvk})}},$$

(6)

$${mathrm{RC}}_{uvk} = 2frac{{mathop {sum }nolimits_{n_{mathrm{r}} = 1}^{N_{mathrm{r}}} delta _{uvk}^{(n_{mathrm{r}})}}}{{N_{mathrm{r}}}} – 1,$$

(7)

$$delta _{uvk}^{(n_{mathrm{r}})} = left{ {begin{array}{*{20}{c}} 1 & {{mathrm{BC}}_{{mathrm{null}}uvk}^{(n_{mathrm{r}})} < {mathrm{BC}}_{uvk}} \ {0.5} & {{mathrm{BC}}_{{mathrm{null}}uvk}^{(n_{mathrm{r}})} = {mathrm{BC}}_{uvk}} \ 0 & {{mathrm{BC}}_{{mathrm{null}}uvk}^{(n_{mathrm{r}})} > {mathrm{BC}}_{uvk}} end{array}}, right.$$

(8)

$$beta {mathrm{MPD}}_{uvk} = frac{{mathop {sum }nolimits_i^{S_k} mathop {sum }nolimits_j^{S_k} f_{iu}f_{jv}d_{ij}}}{{mathop {sum }nolimits_i^{S_k} mathop {sum }nolimits_j^{S_k} f_{iu}f_{jv}}},$$

(9)

$${mathrm{BC}}_{uvk} = frac{{mathop {sum }nolimits_i^{S_k} left| {x_{iu} – x_{iv}} right|}}{{mathop {sum }nolimits_i^{S_k} left( {x_{iu} + x_{iv}} right)}},$$

(10)

where ‘(W_{{mathrm{HeS}}uvk})’ 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. (W_{{mathrm{HoS}}uvk}), (W_{{mathrm{DL}}uvk}), (W_{{mathrm{HD}}uvk}), and (W_{{mathrm{DR}}uvk}) are analogous operators for homogeneous selection, dispersal limitation, homogenizing dispersal, or ‘drift’, respectively. ‘(beta {mathrm{NRI}}_{uvk})’ is bNRI of Bin k between community u and v. ‘(beta {mathrm{MPD}}_{uvk})’ is beta mean pairwise distance of Bin k between communities u and v, and ‘(beta {mathrm{MPD}}_{{mathrm{null}}uvk})’ is the βMPD of the null communities randomized according to a null model. ‘({mathrm{Sd}})’ is standard deviation. ‘({mathrm{RC}}_{uvk})’ is modified RC. ‘Nr’ is total randomization times, usually 1000 times. ‘(delta _{uvk}^{(n_{mathrm{r}})})’ is an operator to calculate RC value. ‘({mathrm{BC}}_{uvk})’ is Bray–Curtis dissimilarity index. ‘({mathrm{BC}}_{{mathrm{null}}uvk}^{(n_{mathrm{r}})})’ is Bray–Curtis dissimilarity of Bin k between null communities u and v of the nrth time randomization according to a null model. ‘(f_{iu})’ and ‘(f_{jv})’ represent relative abundance of taxon i in community u or taxon j in community v, respectively. ‘(S_k)’ represents taxa number in Bin k. ‘(x_{iu})’ and ‘(x_{iv})’ 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.

The null model algorithm for phylogenetic metrics is ‘taxa shuffle’21,37, 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 observed19,54. The null model algorithm results heavily depend on the selection of the regional pool, within which randomization is implemented54. 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).

Null 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.

The 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 theory6. 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.

$$P_{tau k} = frac{{mathop {sum }nolimits_{uv}^m frac{{f_{uk} + f_{vk}}}{2}W_{tau uvk}}}{{mathop {sum }nolimits_{uv}^m frac{{f_{uk} + f_{vk}}}{2}}},$$

(11)

$$P_{tau uv} = mathop {sum }limits_{k = 1}^K frac{{f_{uk} + f_{vk}}}{2}W_{tau uvk},$$

(12)

$$P_tau = frac{{mathop {sum }nolimits_{uv}^m P_{tau uv}}}{m} = mathop {sum }limits_{k = 1}^K f_kP_{tau k},$$

(13)

where ‘(P_{tau k})’ is relative importance of the (tau)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). ‘(P_{tau uv})’ is relative importance of the (tau)th ecological process in governing the turnover between communities u and v (Supplementary Fig. 1c). ‘(P_{tau})’ is relative importance of the (tau)th ecological process in governing the turnovers among a group of communities (Supplementary Fig. 1c) or between a pair of groups. Thus, (P_{tau}) can be (P_{{mathrm{HeS}}}), (P_{{mathrm{HoS}}}), (P_{{mathrm{DL}}}), (P_{{mathrm{HD}}}), or (P_{{mathrm{DR}}}) for heterogeneous selection, homogeneous selection, dispersal limitation, homogenizing dispersal, or ‘drift’, respectively. ‘(f_{uk})’ and ‘(f_{vk})’ are total relative abundance of Bin k in community u and community v, respectively. ‘(W_{tau uvk})’ is operator counting whether the kth bin is governed by the (tau)th ecological process, including (W_{{mathrm{HeS}}uvk}), (W_{{mathrm{HoS}}uvk}), (W_{{mathrm{DL}}uvk}), (W_{{mathrm{HD}}uvk}), and (W_{{mathrm{DR}}uvk}) (Eqs. (1)–(5)). ‘(K)’ is total number of bins. ‘(m)’ 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). ‘(f_k)’ is average relative abundance of Bin k in the group of communities.

As shown in Eq. (13), the relative importance of each process (P_{tau}) is the sum of the terms (f_{k} P_{tau k}), by which we can define the contribution of different bins to (P_{tau}) (Eqs. (14) and (15)).

$${mathrm{BP}}_{tau k} = f_kP_{tau k} = frac{{mathop {sum }nolimits_{uv}^m frac{{f_{uk} + f_{vk}}}{2}W_{tau uvk}}}{m},$$

(14)

$${mathrm{BRP}}_{tau k} = frac{{{mathrm{BP}}_{tau k}}}{{P_tau }} = frac{{mathop {sum }nolimits_{uv}^m frac{{f_{uk} + f_{vk}}}{2}W_{tau uvk}}}{{mathop {sum }nolimits_{k = 1}^K mathop {sum }nolimits_{uv}^m frac{{f_{uk} + f_{vk}}}{2}W_{tau uvk}}},$$

(15)

where ‘({mathrm{BP}}_{tau k})’ is Bin contribution to Process, measuring the contribution of Bin k to the relative importance of (tau)th ecological process in the assembly of a group of communities (Supplementary Fig. 1e). ‘({mathrm{BRP}}_{tau k})’ is Bin Relative contribution to Process, measuring the relative contribution of Bin k to the (tau)th ecological process (Supplementary Fig. 1f).

### Simulation model

In 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 Model55, using R package ‘sads’ (version 0.4.2)56 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.

A phylogenetic tree was retrieved from a previous publication20, 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 model20. The Blomberg’s K value is as low as 0.15, close to the mean K value of 91 continuous prokaryotic traits42. 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 distance20, in accordance with general observations in microbial communities in various environments19,38. (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)57 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 traits42. (iii) The high-phylogenetic-signal pool was simulated according to Blomberg’s ACDC model58 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 date42.

For 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).

$$A_{ij} propto {mathrm{exp}}left[ { – frac{{({{mathrm{EV}}_j – E_i})^2}}{{2sigma _{mathrm{E}}^2}}} right],$$

(16)

where ‘(A_{ij})’ is abundance of species i in local community j. ‘({mathrm{EV}}_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). ‘(E_i)’ is the optimum environment of species i. ‘(sigma _{mathrm{E}})’ 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.

(i) To simulate competition without stochasticity, a geometric series model59 was modified to consider stronger competition between species with similar niche preference37. 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 |EVuEi|. 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 distance37 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.

$${mathrm{nd}}_{hi} = sqrt {mathop {sum }limits_{j = 1,,i > j}^{h – 1} 0.5^j({E_i – E_j})^2},$$

(17)

where ‘({mathrm{nd}}_{hi})’ 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. ‘(E_i)’ is the optimum environment of species i which is not included in the (h−1) prior competitors. ‘(E_j)’ is the optimum environment of species j which is the jth prior competitor with relative abundance proportional to 0.5j.

(ii) To simulate extreme dispersal without selection, we modified Sloan’s simulation model60 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.

(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.

Under 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).

$${mathrm{EP}}_{tau uv} = mathop {sum }limits_{i = 1}^K frac{{f_{uk} + f_{vk}}}{2}omega _{tau uvk},$$

(18)

where ‘({mathrm{EP}}_{tau uv})’ is the expected relative importance of the (tau)th ecological process in community turnover between samples u and v. ‘(f_{uk})’ is total relative abundance of Bin k in community u. ‘(f_{vk})’ is total relative abundance of Bin k in community v. ‘(omega _{tau uvk})’ is operator, equal to 1 if the turnover of the kth bin between communities u and v is governed by the (tau)th ecological process, and equal to 0 if not.

We 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 NP61, tNST33, and pNST33,34.

The performance of quantitative estimation is evaluated by accuracy (Eq. (19)) and precision coefficients (Eq. (20)) derived from concordance correlation coefficient (CCC)62. 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)).

$${mathrm{qACC}} = frac{{2sigma _xsigma _y}}{{sigma _x^2 + sigma _y^2 + left( {mu _x – mu _y} right)^2}},$$

(19)

$${mathrm{qPRC}} = frac{{sigma _{yx}}}{{sigma _xsigma _y}},$$

(20)

where ‘qACC’ and ‘qPRC’ are quantitative accuracy and precision, respectively. ‘(sigma _{yx})’ 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. ‘(sigma _x^2)’ and ‘(sigma _y^2)’ are variance of x and y, respectively. ‘(mu _x)’ and ‘(mu _y)’ are mean of x and y, respectively.

$${mathrm{ACC}} = frac{{{mathrm{TP}} + {mathrm{TN}}}}{{{mathrm{TP}} + {mathrm{TN}} + {mathrm{FP}} + {mathrm{FN}}}},$$

(21)

$${mathrm{PRC}} = frac{{{mathrm{TP}}}}{{{mathrm{TP}} + {mathrm{FP}}}},$$

(22)

$${mathrm{SST}} = frac{{{mathrm{TP}}}}{{{mathrm{TP}} + {mathrm{FN}}}},$$

(23)

$${mathrm{SPC}} = frac{{{mathrm{TN}}}}{{{mathrm{TN}} + {mathrm{FP}}}}.$$

(24)

In the qualitative performance indexes, ‘({mathrm{ACC}})’ is accuracy; ‘({mathrm{PRC}})’ is precision measured as positive predictive value; ‘({mathrm{SST}})’ is sensitivity measured as true positive rate; ‘({mathrm{SPC}})’ is specificity measured as true negative rate. ‘({mathrm{TP}})’ 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. ‘({mathrm{TN}})’ 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. ‘({mathrm{FP}})’ 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. ‘({mathrm{FN}})’ 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.

For 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.

### Experimental data and analyses

We applied iCAMP to an empirical dataset from our previous study34, 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)34. 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 lysis63 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 MiSeq34, 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)34 built on the Galaxy platform (version 17.01)64 and OTUs were generated by UPARSE65 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)34. Plant biomass was measured with a modified pin-touch method and the plant richness was based on identification of all species in each plot34. The drought index is calculated as additive inverse of standardized precipitation–evapotranspiration index (SPEI) retrieved from SPEIbase66.

### Statistical analyses

The 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 d67. NST33, NP61 and QPEN19,20 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.

For correlation test between each process and various measured factors, we applied Mantel test, partial Mantel68, and MRM69 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)70. All significance tests are two-side unless specified.

### Reporting summary

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