Statistical analyses
Two-ways analyses of variance (ANOVAs) were conducted to figure out the main and interactive effects of afforestation and season on ecosystem functioning. To test whether significant relationships and covariations of ecosystem functioning with microbial diversity existed, we used the Pearson correlation test and scatter plots to examine and describe significant BEF relationships. To test whether the sensitivity of BEF relationship varied with different land use types, a standardized major axis regression was performed to examine the slope with the “smatr” package. Variance analyses of abundance of microbes and microbial functioning were performed to assess the effects of land use types and season dynamics. Principal coordinates analyses (PCoA) were used to evaluate the previous treatments effect on the Euclidean distance of soil microbial communities and Bray-Curtis distances of functioning genes with the “vegan” package. Given these distance matrices, we computed Mantel correlation between microbial communities (keystone species, non-keystone species and total communities) with environmental data with “vegan” package, and visualized these relationships with the “ggcor” package.
The co-occurrence patterns of microbial communities were conducted based on the correlation across different land use types. Spearman’s correlation between any two OTUs was estimated, only the robust correlation with absolute value of correlation coefficients > 0.6 and FDR-correctedP -value < 0.01 were used to form networks in “igraph” R package (Benjamini et al. 2006), and visualized network images with Gephi (http://gephi.github.io/). Moreover, 1, 000 Erdös–Rényi random networks with the same numbers of nodes and edges as the corresponding real networks were generated with “erdos.renyi.game” function in “igraph” package, and the mean value of each network metric was compared to real networks with Z-test (Hu et al. 2017). In each network each node represents one OTU, and same color nodes represent the same modules. The network and sub-network properties were calculated in the “igraph” package. The more detail topological features were estimated by a set of metrics and present in Table S1. The algorithm of fast greedy modularity optimization was applied to isolate modules. Keystone species were defined based on their roles in the network structure with the index (among and within module connectivity; Pi and Zi score). The node with either a high value was defined as keystone species (Zi >2.5 or Pi > 0.62) (Deng et al. 2012).
Structural equation modeling (SEM) was further used to evaluate the direct and indirect effect of environmental factor and microbial and keystone species richness on ecosystem functioning. Environmental and microbial communities were performed with principal component analyses. The first two components were chosen to represent the environments and community composition to eliminate potential collinearity. Models were conducted via the R package “lavaan”. All statistical were performed in R and visualization with “ggplot2” and “ggpubr” packages (http://www.r-project.org/) unless otherwise indicated.