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.