Network construction
To understand whether and how the co-occurrence networks including fungal and bacterial community varies across the Tibetan Plateau in relation to environmental and plant richness gradients, two different kinds of networks were constructed basing on the Spearman correlation matrix by “WGCNA” R package (Langfelder & Horvath 2012), viz., molecular ecological network (MEN) including fungi-only and bacteria-only network, and plant–microbiota interkingdom ecological networks (IDEN) including plant-fungi network and plant-bacteria networks.
Given the generally observed relationships between broad habitat zones and microbiota and floristic composition, we also divided the 60 samples according to the three main types of vegetation in our samples (desert steppe, alpine meadow, alpine steppe), and constructed plant-fungi and plant-bacteria networks for each vegetation type.
To avoid the bias of the correlation matrix causing by rare taxa, only OTUs with average relative abundances > 0.01% of each subgroup were retained. Since the number of plant species amongst these 60 sites is low, we kept all the plant species for the plant-microbiota network construction. The Spearman correlations between OTUs and plants were filtered by the thresholds r > 0.6 and false discovery rate adjusted p < 0.05(Huang et al. 2019). The OTUs and plants presented at each site were retained and generated subnetworks for each soil sample from the combined interkingdom ecological networks by the “igraph” R package. Only the correlations between plants and fungi (or bacteria) in each site were kept by the “startswith” function in python, and were chosen as the adjacent matrix of the bipartite graph. The obtained adjacent matrix associated with the bipartite graph consisted of 1 or 0, showing the presence/absence of corresponding plant–microbiota associations (Fenget al. 2019). The plant–fungi and plant-bacteria network architecture of each group was visualized based on the “ForceAtlas2” layout algorithm (Jacomy et al. 2014) using the program Gephi (Bastian et al. 2009). We then examined the number of edges, plant and microbial species richness in the observed IDEN in 60 sites. The observed IDEN topological features (Table S1) was evaluated at both network and group (plants or microbiota) levels using “bipartite” v.2.08 package of R v.3.1.1 (Dormann et al. 2009). Note that low Nestedness values indicate nestedness, while high Nestedness values (0 means cold, i.e. high nestedness, 100 means hot, i.e. chaos) indicate antinestedness. In a nested network, specialists (that is, species with narrow partner ranges) interact with subsets of the partners of generalists (that is, species with broad partner ranges) (Toju et al. 2015). To further determine the compartmentalization of the observed IDEN, modularity was calculated by module detection algorithm for example simulated annealing (Guimeraet al. 2005), and high value indicates modular structure. Modularity is a measure of the extent to which the network is structured as cohesive subgroups of nodes (modules), in which the density of interactions is higher within subgroups than among subgroups (Olesen et al. 2007).
We then conducted the microbial intrakingdom ecological networks analysis using the same thresholds for OTU and correlations mentioned above. To caculated the network-level topological features (Table S2), 60 subnetworks were generated by retaining the OTUs and associated edges for each site using the “subgraph” function in “igraph” R package. Network-level topological features with a high value (such as edge density, degree centralization and betweenness centralization) indicate closer connections within the network, whereas those with lower values (such as average path length and modularity) suggest a more aggregated network (Barberán et al. 2012; Ma et al. 2016). We then calculated the absolute value of negative/positive cohesion to explore the stability of microbial networks along gradient(Yuan et al.2021).