Climate warming enhances microbial network complexity and stability

Unravelling the relationships between network complexity and stability under changing climate is a challenging topic in theoretical ecology that remains understudied in the field of microbial ecology. Here, we examined the effects of long-term experimental warming on the complexity and stability of molecular ecological networks in grassland soil microbial communities. Warming significantly increased network complexity, including network size, connectivity, connectance, average clustering coefficient, relative modularity and number of keystone species, as compared with the ambient control. Molecular ecological networks under warming became significantly more robust, with network stability strongly correlated with network complexity, supporting the central ecological belief that complexity begets stability. Furthermore, warming significantly strengthened the relationships of network structure to community functional potentials and key ecosystem functioning. These results indicate that preserving microbial ‘interactions’ is critical for ecosystem management and for projecting ecological consequences of future climate warming.

In ecosystems, different species are interconnected via exchanges of materials, energy and information to form complex interactions¹, including predation, competition, mutualism, commensalism, parasitism, neutralism and amensalism², which are also affected by a diverse array of abiotic and biotic factors. Such complicated ecological relationships can be represented as networks, that is, the ecological networks with species as nodes and their relationships as links³, which is fundamental for characterizing species interactions and dynamics of ecosystems⁴. One fundamental yet hotly debated question is whether and how the complexity of the ecological network affects ecosystem stability^1,5,6,7. In the past, numerous insights have been revealed into the structure, function and dynamics of ecological systems on food web-based antagonistic networks and pollination-based mutualistic networks, but the results on the relationships between network complexity and stability remain controversial^6,8,9,10. Although it is well documented that changes in network structure affect ecosystem functioning and stability^11,12, little is known about whether and how the networks in ecology, particularly microbial ecology, will change under future climate change scenarios^8,13.

Global climate change is one of the most profound anthropogenic disturbances to our planet, which arises from increased greenhouse gases because of fossil fuel combustion and land use change¹⁴. Global surface temperature has risen 0.78 °C since industrialization and is predicted to increase an additional 1.5–4.5 °C by the end of the twenty-first century¹⁴. Climate warming can alter soil microbial community diversity, structure and activities¹⁵, but it remains uncertain whether and how it impacts network complexity and its relationships to stability in microbial communities.

Because temperature affects all levels of the biological hierarchy¹⁶, it should have profound influences on network complexity and stability. On one hand, as predicted by the metabolic theory of ecology (MTE)¹⁶, rising temperature would stimulate various biotic interactions (predation, parasitism, competition and symbiosis) owing to more active individual metabolic processes and faster growth at higher temperature. As a result, we expect that experimental warming would increase the complexity of species associations, resulting in more complex ecological networks. On the other hand, rising temperature and associated environmental changes (decreased moisture) would act as strong filtering factors against existing microbial species. Therefore, contradictory to the expectations of MTE, the microbial networks could diminish or even collapse, as recently witnessed in a marine food web¹⁰.

To understand whether and how climate warming affects the complexity and stability of ecological networks in soil microbial communities, we examined temporal dynamics of soil microbial communities in a long-term in situ warming experiment, which was carried out in a tallgrass prairie ecosystem of the US Great Plains in Central Oklahoma (34̊° 59ʹ N, 97̊° 31ʹ W)¹⁷. Several previous studies from this experiment revealed that warming significantly increased the rain-use efficiency of net primary productivity and belowground net primary productivity, mainly because proportionally more plant biomass was allocated to roots in response to warming-induced dry conditions¹⁸. Warming also shifted microbial community structure by acting as a deterministic filtering factor, resulting in divergent succession with reduced stochasticity^17,19. However, it is not clear whether warming will affect networks in microbial communities. Because the association networks in microbial ecology are typically reconstructed on the basis of molecular markers, here, we refer to them as molecular ecological networks (MENs)²⁰. Our main questions are the following: (1) whether and how experimental warming affects the complexity of soil MENs over time; (2) whether and how warming alters the relationships between the complexity and stability of soil MENs; (3) whether and how warming-induced changes of soil MENs shape ecosystem functioning.

To address these questions, the composition and structure of soil microbial (bacterial and archaeal) communities under warming and ambient temperature (control) in six consecutive years from 2009 to 2014 were determined with 16S ribosomal RNA gene (V3–V4 region). We then constructed 11 time-series MENs (Fig. 1a) on the basis of Pearson correlations of log-transformed operational taxonomic unit (OTU) abundances, followed by a random matrix theory- (RMT-) based approach²⁰, which is capable of providing a mathematically defined reference point as a threshold for network construction (see Methods for details). The empirical MENs were significantly different from random MENs (Supplementary Table 1). Similar to typical molecular biology and technology networks²¹, all empirical MENs exhibited scale-free features (Supplementary Table 1 and Supplementary Text A), with the distribution of network links independent of network size. In addition, they exhibited small-world characteristics with short geodesic distances (the average shortest path between two nodes) of 3.3–8.3 (Supplementary Table 1), which could allow the effects of a perturbation to distribute rapidly through the entire network, rendering the whole system efficient.

a, Visualization of constructed MENs in six years (Y) from 2009 (Y0) to 2014 (Y5). Large modules with ≥5 nodes are shown in different colours, and smaller modules are shown in grey. Details of network topological attributes are listed in Supplementary Table 1. b–h, Temporal changes of network topology, including n (b), L (c), average K (d), Con (e), average CC (f), RM (g) and the number of keystone nodes (h). In each panel, filled red symbols represent networks under warming, and open blue symbols represent networks under non-warming control condition. Slopes (b) and adjusted r² and P values from linear regressions are shown. *0.01 < P ≤ 0.05; **0.001 < P ≤ 0.01; ***P ≤ 0.001.

Source data

Full size image

Since the preceding networks were reconstructed using a correlation-based approach with taxa occurrence and abundance data, caution is needed in interpreting the mechanisms underlying these networks. Theoretically, species co-occurrence patterns could be driven mainly by three ecological processes: biotic interaction, environmental filtering and dispersal limitation^22,23. Multiple analyses, including multiple regression on distance matrices (MRM), Mantel test, variation portioning analysis (VPA), the link test for environmental filtering or dispersal limitation (Supplementary Fig. 1) and Goderma’s method²⁴, all indicated that these MENs were less likely due to environmental filtering or dispersal limitation, and hence, biotic interactions could be a major driver of the reconstructed networks (Extended Data Fig. 1, Supplementary Figs. 1 and 2, and Supplementary Text B), though roles of unmeasured environmental variables could not be assessed. Due to unmeasured environmental variables and the nature of the correlation-based network approach, the occurrence-estimated links in MENs should at best be treated as putative interactions²⁵, and thus, we use the term ‘biotic interactions’ throughout for simplicity.

The microbial MENs under warming and control followed different successional trajectories over time on the basis of 22 network topological parameters (Supplementary Fig. 3). Although the numbers of OTUs used for network construction were on average 12% fewer under warming than under control, the resulting networks without isolated nodes were on average 14% larger in size (total nodes) under warming than under control (Supplementary Table 1), suggesting that these microbial taxa could associate more closely with each other, and hence, neutral process might be less prevalent under warming¹⁷. In addition, the composition of the networked communities (assemblages of microbial taxa detected in the networks) was significantly different between warming and control as revealed by detrended correspondence analysis (Extended Data Fig. 1a) and three non-parametric dissimilarity analyses (multiple response permutation procedure (MRPP), analysis of similarity (ANOSIM) and permutational multivariate analysis of variance (Adonis), Table 1). Although soil temperature, moisture and plant biomass played statistically significant roles in controlling the structures of the networked communities as revealed by canonical correspondence analysis (CCA, Extended Data Fig. 1b), all measured soil and plant variables were able to explain only minor portions (<13%) of the networked community variations (Extended Data Fig. 1c). In addition, using the recently developed method for alleviating compositional data bias²⁶, 58% of the phyla (11 out of 19) and 46% of the classes (22 out of 48) showed absolute values of differentials of larger than 0.5 in their abundances between warming and control (Supplementary Fig. 4 and Supplementary Table 2), with significant increase of some fast-growing bacteria (Actinobacteria) but decrease of some slow-growing bacteria (Negativicutes in Firmicutes) under warming (Supplementary Tables 3 and 4). Similar trends of warming’s impacts on different phyla and classes were observed with the Mann–Whitney U test (Supplementary Table 2). Finally, there were 43% more positive and 18% more negative associations under warming than under control (Supplementary Table 1). If these associations do in fact indicate biological interactions, positive associations (positive correlations between nodes) could represent mostly cooperative behaviours such as cross-feeding, syntrophic interactions, mutualistic interactions and commensalism as well as shared environmental requirements and common dispersal barriers. By contrast, negative associations (negative correlations between nodes) could reflect mainly competition for limiting resources as well as distinctive environmental niches and spatial isolation^4,27, although the exact specific mechanisms underlying these associations are unknown with correlation-based network analyses.

To determine how warming affected microbial network complexity, changes of various network topological parameters were regressed against time under warming and control. Network size (total number of nodes, n) strongly increased over time under warming (r² = 0.82, P = 0.009), so did network connectivity (total number of links, L; r² = 0.90, P = 0.003), average connectivity (average links per node, average K, or 2 L n^–1; r² = 0.92, P = 0.002), average strength (average weighted connectivity, average S, r² = 0.92, P = 0.002), average clustering coefficient (the extent to which nodes are clustered, average CC, r² = 0.94, P = 0.001), connectance (the proportion of realized links in all possible ones, Con, r² = 0.82, P = 0.008), number of nodes in the largest connected component (r² = 0.85, P = 0.005) and largest module (r² = 0.96, P < 0.001, Fig. 1b–f and Supplementary Fig. 5a,b). By contrast, no significant trend was observed under control across time. In addition, we quantified the degree of community complexity using a recently developed metric, named cohesion²⁸. Positive cohesion, derived from positive pairwise correlations, could reflect the degree of cooperative behaviours in a sample, whereas negative cohesion could indicate the magnitude of competitive behaviours among community members (see Methods for details). Positive cohesion increased significantly over time under warming (P = 0.013) but not under control (Extended Data Fig. 2). Strong correlations were observed between positive cohesion and various network complexity indices under warming (r > 0.83, P ≤ 0.05) rather than control (Extended Data Fig. 2c), suggesting that warming could enhance taxon connectedness and hence possible biotic interactions²⁸. These results collectively indicated that the MENs were substantially different between warming and control and that warming significantly increased network complexity over time.

Changes in network structure could further affect network organizational principles, such as modularity and nestedness. The MENs were highly modular (Supplementary Table 1) but not nested (Supplementary Text A). A total of 100 large modules (modules with ≥5 nodes) accounted for 56–72% of the nodes in the MENs under control, while 98 large modules accounted for 45–86% of the nodes in the MENs under warming (Supplementary Table 1). If two modules in different networks had significantly large proportions of shared nodes, they were considered as preserved module pairs (Supplementary Text C). There were only 15 module pairs (18% of all the pairs tested) preserved between warming and control networks of the same year, suggesting that warming drastically changed module identity (Extended Data Fig. 3, Supplementary Table 5 and Supplementary Text C). To visualize the changes in the higher-order organization of the MENs, module eigengene analysis was also performed²⁹. A module eigengene is a synthetic relative abundance profile, collectively representing the relative abundance changes of all individual OTUs within the module, and the analysis based on module eigengene can provide information on the relationships among different modules (Supplementary Text C). The eigengene analysis showed that the relationships among different modules were substantially different between warming and control (Supplementary Fig. 6 and Supplementary Text C), indicating that the higher-order network organization was not preserved under warming. In addition, the MENs had 58% higher relative modularity under warming than under control, although all of them were significantly higher than those from their corresponding random MENs (Supplementary Table 1). The relative modularity significantly increased over time under both warming (r² = 0.84; P = 0.007) and control (r² = 0.60, P= 0.043, Fig. 1g), with a larger slope under warming (slope = 0.065) than under control (slope = 0.021). Interestingly, the relative modularity (RM) under warming increased significantly with various network metrics such as network size, connectivity, average connectivity, average clustering coefficient and connectance, but not under control (Fig. 2). Altogether, those results indicated that the MENs between warming and control were poorly conserved at the higher-level organization, and that warming significantly enhanced relative modularity.

a–e, Network topological properties include n (a), L (b), average K (c), Con (d) and average CC (e). In each panel, filled red circles with the red solid line (least square best fit) represent networks under warming, and open blue circles with the blue dashed line represent networks under non-warming control condition. Pearson correlation coefficients (r) are shown for warming and control networks in corresponding colours. **0.001 < P ≤ 0.01; ***P ≤ 0.001.

Source data

Full size image

The altered network complexity could lead to changes in the role of individual members within the network. On the basis of their within-module connectivity (Z_i) and among-module connectivity (P_i)²⁰, a total of 128 module hubs (nodes highly connected to other members in a module), 39 connectors (nodes linking different modules) and 2 network hubs (nodes being both a module hub and a connector) were detected across all the MENs (Extended Data Fig. 4, Supplementary Table 6 and Supplementary Text D), all of which could be regarded as keystone nodes (Methods) that play key roles in shaping network structure³⁰. Many (62 of 169) of the keystone nodes-affiliated taxonomic groups appeared to be important in carbon compound degradation, nitrification, denitrification and phosphorus utilization, as well as adaptation to elevated temperatures (Supplementary Table 6). Importantly, there were 27% more keystone nodes under warming than under control, and the number of keystone nodes increased significantly over time under warming (r² = 0.94, P = 0.001, Fig. 1h) but not under control. In addition, phylogenetic analysis showed that only 9.8% of all the keystone nodes were shared between warming and control (Extended Data Fig. 4g and Supplementary Text D). Therefore, keystone nodes were also not preserved between warming and control, and warming altered network structure at the keystone node level.

All the MENs detected are scale free, modular and small world, which may have profound implications on microbially mediated ecosystem functions and stability²¹ (see Supplementary Text E for detailed possible reasons). Therefore, we simulated species extinction and calculated robustness (the resistance to node loss)⁹ of the MENs. On the basis of either random species loss or targeted removal of module hubs, the MENs had significantly higher (P < 0.001) robustness under warming than under control after 2011 (Fig. 3a,b). The network vulnerability (the maximum decrease in network efficiency when a single node is deleted from the network)³⁰ (see Methods for details) was also on average lower under warming (0.10 ± 0.07) than under control (0.18 ± 0.13) and decreased marginally significantly over time under warming (r² = 0.62, P = 0.060) but not under control (Fig. 3c), showing higher robustness in warming networks.

a, Robustness measured as the proportion of taxa remained with 50% of the taxa randomly removed from each of the empirical MENs. b, Robustness measured as the proportion of taxa remained with five module hubs removed from each of the empirical MENs. In a and b, each error bar corresponds to the standard deviation of 100 repetitions of the simulation. Significant comparisons (two-sided t-test) between warming (W) and control (Ctl) are indicated by ***P < 0.001. c, Network vulnerability measured by maximum node vulnerability in each network. The adjusted r² and P values from linear regressions are shown. d, Compositional stability of the networked community over time from every two adjacent years. The slopes (b from Y = a + bX), adjusted r² and P values are shown. e, Pearson correlation of compositional stability and node persistence. The line indicates least square best fit under warming (red) or control (blue) condition. In d and e, each error bars represent the standard deviation in 24 plots. f, Pearson correlations between network complexity and stability indices under warming (boxed in red) or control (boxed in blue). Significant (P ≤ 0.05) correlations are shown here, with orange for positive correlations and green for negative correlations. Numbers inside cells are corresponding correlation coefficients. Correlations with P > 0.05 are marked in grey.

Source data

Full size image

Multiple stability indices calculated from empirical data also supported that warming enhanced MEN stability. First, multi-order compositional stability (the temporal invariability of community composition)³¹ for microbial communities in the MENs increased significantly over time (r² = 0.86, P = 0.015, Fig. 3d and Supplementary Fig. 7a–e) under warming, indicating that microbial community composition became more stable under warming than under control. The constancy of nodes and links (the inverse of their temporal variations) and node persistence (the proportion of nodes persisting over time) of the MENs were higher under warming, revealing smaller temporal changes of the MENs under warming than under control (P ≤ 0.001, Extended Data Fig. 5 and Supplementary Fig. 7f–j). In addition, node persistence strongly correlated with the compositional stability under warming (r = 0.95, P = 0.014) rather than under control (Fig. 3e). These results suggested that the MENs appeared more stable under warming than under control.

More important, significant correlations were observed between various stability measurements and network complexity. Network robustness, compositional stability or node persistence of the MENs was positively correlated with various indices of network complexity under warming, such as node and link numbers, average connectivity, average clustering coefficient, geodesic distance, harmonic geodesic distance, relative modularity and connectance (Pearson r = 0.83–0.99, P ≤ 0.036, Fig. 3f), but with very few under control. Consistently, network vulnerability significantly decreased with connectance and the number of keystone nodes under warming, but not under control (Fig. 3f). Similarly, positive cohesion was significantly correlated with robustness and persistence under warming, while negative cohesion was negatively correlated with persistence (Extended Data Fig. 2c). No significant correlations between cohesion and network stability were observed under control (Extended Data Fig. 2c). Collectively, the preceding results indicated that warming enhanced the stability of the soil MENs, probably due to warming-enhanced network complexity associated with connectivity and relative modularity.

An intriguing question is whether the warming-enhanced network complexity and stability affect microbial community functional structure and associated ecosystem processes. We thus address it with three analyses. First, the networked communities exhibited generally stronger correlations with various variables of ecosystem functioning, particularly GPP, ecosystem respiration (ER) and net ecosystem exchange (NEE), under warming than under control (Fig. 4a and Supplementary Fig. 8a). Second, microbial community functional structure was assessed by using a comprehensive functional gene microarray, GeoChip 5. Interestingly, there were 57% more significant correlation pairs between network indices and carbon (C) degradation genes under warming than under control (Fig. 4b and Supplementary Fig. 8b), suggesting that the changes in network structure could be more tightly associated with the observed stimulated C degradation functional potentials and heterotrophic respiration under warming. CCA showed that GeoChip-based functional traits were tightly linked with various ecosystem functional processes (Supplementary Table 7). In addition, the potential metabolic functions of individual nodes were predicted on the basis of their 16S rRNA gene sequences using PICRUSt2³², and warming significantly shifted the community structures of all PICRUSt2-derived individual metabolic pathways harboured in the networked communities, including carbohydrate metabolism, energy metabolism, and glycan biosynthesis and metabolism (Supplementary Table 8). Altogether, these results indicate that the warming-induced network complexity and stability might have strong effects on microbial community functional structure and hence ecosystem functional processes.

a, Correlations of the networked community structures (Bray–Curtis distance) with time, soil and plant variables and ecosystem C fluxes. Edge width corresponds to the Mantel’s r value, and the edge colour denotes the statistical significance. Pairwise correlations of these variables are shown with a colour gradient denoting Pearson’s correlation coefficients. Soil variables include soil nitrate (NO₃^–), ammonium (NH₄⁺), total nitrogen (TN), total organic C (TOC), pH, moisutre and temperature; plant variables include C₃ and C₄ aboveground biomass, total biomass and plant richness; ecosystem C fluxes include ecosystem respiration (ER; total respiration from below- and aboveground), gross primary productivity (GPP; photosynthetic rate), net ecosystem exchange (NEE; the differece between GPP and ER), autotrophic respiration (R_a; those from root), heterotrophic respiration (R_h; mostly from soil microbes) and total soil respiration (R_t; the sum of R_a and R_h). b, Pearson correlations of MENs’ topological properties and relative abundances of C degradation genes. The annotations of the C degradation genes are listed in Supplementary Table 9. For significant (P ≤ 0.05) correlations, positive correlation coefficients are highlighted in green, while negative correlation coefficients are highlighted in orange. Among the detected C degradation genes, only those with significant correlations to MEN properties are shown here.

Source data

Full size image

In summary, by examining the temporal network dynamics of grassland soil microbial communities in response to warming, this study provides several important insights into the roles of experimental warming in mediating the relationships between microbial network complexity and stability. First, as predicted by the MTE¹⁶, experimental warming significantly accelerated the dynamic changes of network complexity over time. This is probably because warming acts as a deterministic filtering factor¹⁷ to select for certain fast-growing bacteria, but against some slow-growing bacteria (Supplementary Tables 2 and 3), which could further trigger a series of dynamic changes among different microorganisms. Consequently, new network structures that could quickly adapt to the increased temperature and other associated changes were established. By contrast, neutral taxa that are less responsive to warming could prevail in the microbial communities under control^17,19, resulting in less complex network structures. Second, although climate warming is known to trigger complex interactive effects on ecosystem structure and functioning^10,17,19, it is not clear whether and how it affects network stability. This is new documentation that climate warming enhances network stability over time. In addition, our study provides explicit evidence that network stability increases with network complexity, particularly relative modularity, which is consistent with that observed in macroecology¹. These results also lend support to MacArthur’s argument that the complexity of ecosystems begets their stability³³, but contradict the theoretical analysis showing that higher complexity destabilizes ecological systems³⁴. Finally, several previous studies documented that long-term warming significantly alters microbial biodiversity, and thus various ecosystem processes such as soil carbon fluxes and nutrient cycling^14,35,36,37, but it is not clear whether warming-induced network changes will affect ecosystem functional processes. Our study further suggested that the warming-enhanced network complexity and stability appeared to be more tightly associated with microbial community functional structure and ecosystem functional processes.

Our findings have several important implications for projecting ecological consequences of future climate warming and ecosystem management. First, similar to microbial biodiversity, which is dependent on both space and time, the network features of MENs are also temporally dynamic, particularly under warming. Thus, future studies on microbial networks need to integrate the temporal dynamics along with their spatial turnovers³⁸. Further, since the networked communities have strong linkages with ecosystem functioning, preserving network structure is important for future ecosystem conservation³⁹. In addition, on one hand, as previously demonstrated, due to accelerated dynamic responses, biodiversity could change more quickly under future climate change scenarios¹⁹, which could result in potential biodiversity loss. Consequently, its linked ecosystem functions and services may become more vulnerable^17,19. On the other hand, warming significantly stimulates the dynamic responses of network complexity, particularly relative modularity, which leads to higher community stability¹², and hence, the linked ecosystem functions could be less vulnerable in a warmer world. Therefore, preserving microbial ‘interactions’ could be important to mitigate the detrimental effects of warming-induced biodiversity loss on ecosystem functions⁴⁰ if the warming-enhanced network complexity and stability observed in this study are applicable to other ecosystems.

Field experiment description

The field experiment was conducted at the Kessler Atmospheric and Ecological Field Station of the US Great Plain in McClain County, Oklahoma, USA (34° 59ʹ N, 97̊° 31ʹ W). As described previously¹⁸, the average air temperature was 16.3 °C, and the average annual precipitation was 914 mm, according to the Oklahoma Climatological Survey, from 1948 to 1999. Categorized by the type of photosynthesis and plant functional group, the experimental site was dominated by C₃ forbs (Solanum carolinense, Ambrosia trifida and Euphorbia dentate), C₃ grasses (Bromus spp) and C₄ grasses (Tridens flavus and Sorghum halapense). The soil type was Port–Pulaski–Keokuk complex, with a neutral pH, a high available water-holding capacity (37%) and a deep (~70 cm), moderately penetrable root zone⁴¹. The soil texture is loam with 51% sand, 35% silt and 13% clay. The experiment used a blocked split-plot design, in which warming (continuous heating at a target of +3 °C above ambient) and precipitation alteration (targets of −50% and +100%) were primary factors nested with annual removal of aboveground biomass in peak growth season as the secondary factor. Altogether, there were 24 plots under warming and 24 plots under ambient temperature as control. The experimental site was initiated in the fall of 2009.

Field measurements, soil geochemistry and microbial community characterization

Soil temperature, moisture (volumetric water content), total C, total nitrogen (N), nitrate (NO₃^–) and ammonia (NH₄⁺), plant biomass and richness (separated into C₃ and C₄ species), and ecosystem C flux were measured and analysed as described previously^17,19. In 2009, 24 pre-warmed surface (0–15 cm) soil samples were collected using a soil core (2.5 cm in diameter, 15 cm deep). In 2010–2014, a total of 240 samples were collected, one sample per plot per year. Although microscale soils are highly heterogeneous^42,43,44, centimetre-scale sampling is demonstrated to be appropriate to study ecosystem-level responses^45,46. Soil microbial DNA was extracted from 1 g of well-mixed soil for each sample by grinding, freeze-thawing and sodium dodecyl sulfate (SDS)-based cell lysis⁴⁷, and purified with the PowerMax Soil DNA Isolation Kit (MO BIO). The DNA quality was assessed on the basis of 260/230 nm and 260/280 nm absorbance ratios using a NanoDrop ND-1000 Spectrophotometer (NanoDrop Technologies) to ensure that the 260/230 ratios were larger than 1.7 for all samples, and 260/280 nm were larger than 1.8. Although this extraction method is more labour intensive and time consuming, it is capable of recovering high-molecular-weight DNA with high yield and high quality from diverse representative soil samples^47,48.

The 16S rRNA gene was PCR-amplified with the primer set 515F (5’-GTGCCAGCMGCCGCGGTAA-3’) and 806R (5’-GGACTACHVGGGTWTCTAAT-3’), targeting the V4 hypervariable regions⁴⁹, using a two-step PCR protocol with phasing priming technique⁵⁰. Briefly, the DNA (10 ng) from each sample was first PCR-amplified for ten cycles in triplicate in 50 μl reaction with the primers without adaptors to avoid potential amplification bias, due to the long tail of adaptors and other added components. The PCR products were purified and resolved in 50 µl deionized water. Then, 15 µl of the PCR products from each sample was subjected to second-round amplification using the primers with all adaptors, barcodes and spacers in triplicate an additional 15 cycles to ensure the PCR amplification was not saturated so that semi-quantitative information could be achieved and amplification artifacts could be limited. Finally, all triplicate PCR-amplified products were combined, purified and quantified for subsequent sequencing. Various control experiments demonstrated the two-step PCR method with phasing primers in triplicate can help to reduce sequencing errors, minimize amplification bias and preserve the semi-quantitative nature of PCR amplification^50,51. Although amplicon sequencing has low reproducibility and poor quantification^52,53, amplification in triplicate in two steps with lower total cycle numbers (25–30 cycles) will help to improve sequence representation and quantification. This is critical for subsequent data analysis, data interpretation and biological inference⁵³.

The amplified products were then sequenced with the 2 × 250 bp kit on the Illumina MiSeq platform (Illumina). Amplicon sequencing data were processed through a pipeline (http://zhoulab5.rccc.ou.edu:8080/root) by the Institute for Environmental Genomics at the University of Oklahoma^17,19. In brief, after demultiplexing, primer trimming and quality score-based cleaning, forward and reverse sequences were combined and clustered into OTUs at 97% similarity. We used OTUs at 97% similarity instead of amplicon sequence variants because the latter generates more taxa, and the abundance matrix is sparser, which could result in more bias in estimating correlation for network construction. Then, singleton OTUs were removed, and sequences were resampled to the same sequence depth (30,000 sequences per sample) across all samples.

MEN construction

All MENs were constructed on the basis of Pearson correlations of log-transformed OTU abundances, followed by an RMT-based approach that determines the correlation cut-off threshold in an automatic fashion^20,54,55. Random matrix theory was initially proposed in the 1960s as a procedure to identify phase transitions associated with noise in physics and material science and was later adopted for studying the behaviours of many other complex systems, including gene co-expression network construction for predicting gene functions^54,55 and molecular ecological network construction^20,56,57. The core of RMT states that level (eigenvalue) fluctuations of real random matrices follow two different universal laws, depending on the correlation property of eigenvalues. The nearest-neighbour spacing distribution of eigenvalues (the distribution of the difference of two nearest-neighbour eigenvalues) follows Gaussian orthogonal ensemble (GOE) statistics if there is correlation between nearest-neighbour eigenvalues, while it follows Poisson statistics if there is no correlation⁵⁸. Deviations from GOE universal predictions can be used to distinguish system-specific, nonrandom properties of complex systems from random noise⁵⁹.

A RMT-based network framework^20,29,56 was developed to automatically identify molecular ecological networks in microbial communities by assuming that the two RMT predictions are applicable to ecological communities as complex systems. We then also hypothesize that there is a transition of nearest-neighbour spacing distribution of eigenvalues from GOE to Poisson distributions, which can be used as a reference point to disentangle random noise from system-specific, nonrandom properties embedded in high-throughput metagenomic data. In this study, since correlation matrix is used, this transition point (St) is then used as the correlation cut-off. This RMT-based approach avoids an arbitrary St determination commonly used in association-based network methods. Once the adjacency matrix is selected with a defined St threshold, an undirected network graph can be drawn.

The RMT-based network approach has four advantages^20,54,55,56. First, it has a sound theoretical foundation because it is based on the two universal laws of RMT⁵⁵. Second, its automatic threshold detection avoids arbitrary cut-off determination, which is a typical serious problem and a major source of uncertainties in association network construction. Third, it reliably removes noises from nonrandom, system-specific features. Last, it applies to diverse formats of ecological datasets (hybridizations, sequencing, geochemistry and physiology). Although no method is perfect, we believe that the RMT-based network approach is the most appropriate choice for this study due to its mathematically defined non-arbitrary correlation cut-off to minimize the uncertainty in network construction and comparison. This RMT-based network tool is named the Molecular Ecological Network Analysis Pipeline (MENAP) and is available at the Institute for Environmental Genomics, University of Oklahoma (http://ieg4.rccc.ou.edu/MENA/).

One MEN was constructed on the basis of the abundance profile of OTUs in the 24 samples collected in 2009 as the initial network. Then, an additional ten MENs were constructed, one for warming and one for control each year, using samples collected from 2010 to 2014. Each of the MENs was constructed independently on the basis of a unique set of biological samples, which captured average rather than transient effects of warming on the soil microbial communities. To ensure the reliability of correlation calculation, only OTUs present in 12 of the 24 samples were included for correlation calculation. The Pearson correlation coefficient was calculated for each pair of OTUs on the basis of the OTUs’ log-transformed relative abundances. The resulting correlation matrix was analysed with the RMT-based network approach to determine the correlation threshold for network construction. After St was determined for each correlation matrix, an adjacency matrix was generated, containing only the correlations whose absolute values of coefficient (correlation strengths) were larger or equal to St. Nodes in isolation after the St threshold (no correlation coefficient to other nodes ≥St) were removed from the network.

Due to large inherent variations of the sequence output among different samples, relative abundance is generally used for subsequent analysis. As a result, the data become compositional, which could result in spurious indirect correlations^26,60,61,62. To assess the potential influence of compositional data bias, we compared the results from the networks based on log-transformation and centred log-ratio transformation, which are expected to mitigate the bias induced by compositionality²⁵. There were very strong correlations of various topological properties between the networks based on these two transformed datasets (r = 0.96–0.99, P < 0.001, Supplementary Fig. 9), indicating that the networks are very similar. We also used SparCC^63,64, an algorithm that is compositionally robust in calculating correlations to construct networks. The topological properties of SparCC-based networks revealed a similar trend as the RMT-based MENs we presented in this study (Supplementary Fig. 10). All the preceding results suggested that the effects of the compositional bias on network structure of highly diverse microbial communities as included in this study could be negligible⁶⁵.

Network characterization

We tested the potential contribution of environmental filtering or dispersal limitation in shaping network topology (Supplementary Text B). To test whether soil and plant variables and spatial distance had impacts on the composition of the networked communities, MRM was carried out in R software⁶⁶ (version 3.5.3) with the MRM function in the package {ecodist} version 2.0.1⁶⁷. The Mantel tests between network topological parameters and plant and soil variables were performed in R with the mantel function in the package {vegan} version 1.4-2⁶⁸. To detect taxon–taxon–environment co-variation links⁶⁹, we developed a Python 3 script. To unveil links possibly caused by dispersal limitation, we developed an R script. We named this method the link test for environmental filtering or dispersal limitation (Supplementary Fig. 1), and the scripts and examples are publicly available on the GitHub repository⁷⁰. A tool recently reported by Goberna et al. was further used to quantify the relative contributions of community assembly processes (biotic interactions, environmental filtering and dispersal limitation) to the observed patterns of MENs²⁴. Although the preceding analyses could reveal the relative importance of biotic interactions, particularly for non-trophic biotic interactions (facilitation, competition), in shaping the MENs, we still do not have a way to prove that the links are truly due to biotic interactions. Interpretation of these results thus needs to consider this point, and they should be best used for relative comparison across different conditions or treatments^{29,52,71,72,73,74}. Thus, in this study, we are focusing mainly on the relative changes of the estimated network parameters between warming and control.

Various network topological indices were calculated in the MENAP interface to characterize the topological structure of the MENs⁵⁶: n, L, power-law fitting of node degrees, average K, average CC, average path distance (geodesic distance), geodesic efficiency, Con, modularity, harmonic geodesic distance, maximal degree, centralization of degree, maximal stress centrality, centralization of stress centrality, maximal betweenness, centralization of betweenness, maximal eigenvector centrality, centralization of eigenvector centrality, density, transitivity and efficiency. In addition, two topological indices for weighted networks, average strength and average weighted clustering coefficient, were calculated in R package {igraph}⁷⁵ to reflect potential effects of edge weight on the topology of the MENs. The correlation coefficients between pairs of nodes were used as the weights of edges. To test the significance of the constructed empirical MENs, 100 random networks were generated for each empirical network by randomly rewiring the links among the nodes while constraining n and L, following the Maslov–Sneppen procedure⁷⁶. With each randomization, the same suite of network topological properties was calculated. The means and standard deviations of these properties from the 100 randomizations were calculated and compared with those from the corresponding empirical MENs. Network randomization was performed in the MENAP.

Modularity measures the degree to which a network is compartmentalized into different modules, in which the nodes within a module are highly connected but have very few connections with nodes from other modules. Since the network size and connectivity varied substantially among the different MENs (Supplementary Table 1), RM, which measures how modular a network is as compared with the mean expected modularity¹², is more meaningful for comparing modular structure across different networks. In this study, RM was calculated as the ratio of the difference between the modularity of an empirical network and the mean of modularity from the random networks over the mean of modularity from the random networks¹²:

where M is the modularity of the empirical MEN, and \(\overline {M_r}\) is the mean of modularity from the random networks.

Nestedness describes the level to which the species interacting with specialists are subsets of those interacting with generalists. Nestedness ranges from zero to one⁷⁷. Zero means that the network is not nested (species interacting with specialists are not subsets of those interacting with generalists), and one means that the network is perfectly nested (species interacting with specialists are subsets of those interacting with generalists). As nestedness also depends on network size and connectivity, relative nestedness (RN) was calculated as

where N is the nestedness of the empirical MEN, and \(\overline {N_r}\) is the mean of nestedness values from the random networks.

For each node, its within-module connectivity (Z_i) and among-module connectivity (P_i)⁷⁸ were calculated and used for classification of its topological roles in the network. We adopted criteria used in previous studies^29,57,79 and identified module hubs (Z_i ≥ 2.5, P_i < 0.62), connectors (Z_i < 2.5, P_i ≥ 0.62) and network hubs (Z_i ≥ 2.5, P_i ≥ 0.62). Module hubs, connectors and network hubs are referred to as keystone nodes^80,81. All other nodes were categorized as peripherals. The calculations of Z_i and P_i were performed in the MENAP.

The relationships among different modules in a network represent the higher-order organization of the network. Based on the abundances of the nodes in these modules, a single-dimensional standardized relative abundance vector, termed ‘module eigengene’^56,82, was calculated for each module through single-value decomposition. Those module eigengenes were used to visualize the organization of modules in the MENs through a hierarchical clustering analysis. The eigengene analysis was performed in the MENAP.

To verify the findings from the network analysis, we also calculated a different metric named cohesion, which is an abundance-weighted, null model-corrected metric based on pairwise correlations across taxa²⁸:

where m is the total number of taxa in a community. Briefly, pairwise correlations are first obtained between all taxa on the basis of the data matrix of taxa relative abundance x samples, followed by subtracting the pairwise correlations based on multiple null model iterations to generate the null model-corrected correlations for individual taxa. For each sample, all positive and negative null model-corrected correlations are averaged, respectively, to obtain a connectedness matrix with average positive and negative correlations for different samples. Finally, positive and negative cohesions are calculated on the basis of the preceding formula for each sample. Thus, because positive and negative cohesions in a sample are the sum of abundance-weighted, null model-corrected positive and negative correlations, respectively, they could reflect the degree of cooperative behaviours or competitive interactions. It is argued that cohesion can be used as a measure of the strengths of biotic interactions if taxa are influenced to the same degree by environmental drivers, but differentially affected by species interactions²⁸. To be consistent with all other network analyses, we focused on the taxa present in the network. We used the author-recommended ‘taxa shuffle’ null model with provided R code to calculate both positive and negative cohesions for each microbial community over time. Then, the community’s cohesions were used for correlation analyses with various network complexity indices and stability measurements.

Network comparison

To evaluate the overall difference of MENs, 22 topological indices calculated for each empirical MEN were used as input for a principal component analysis, carried out in R with the function prcomp in the package {stats}. To understand how each network topological property changed over time, we fit a linear model between each individual property and time (in years). The linear models were constructed using R with the lm function in the package {stats}.

Since the node composition differed considerably among modules, Fisher’s exact test was performed to identify preserved module pairs in networks (1) under warming or control over time and (2) between warming and control in the same year⁵⁶. A total of 12,363 pairs of modules were tested. For each pair of modules to be tested, the pool of all the nodes in the corresponding two networks was divided into four categories: (1) nodes belonging to both modules, (2) nodes present only in one of the two modules, (3) nodes present only in the other one of the two modules and (4) nodes absent from both modules but present in the two networks. Frequency of observations in each of these four categories was placed in each of the four cells of a contingency table for a one-sided exact test, which determined whether members in the two modules were independent or exclusive. P values from the exact tests were adjusted through the Bonferroni procedure within each network. If these two modules consisted of a significant portion of the same nodes (OTUs) by Fisher’s exact test (adjusted P ≤ 0.05), they were considered as preserved modules from different networks. The exact tests were performed in R with the fisher.test function (alternative = ‘greater’) in the package {stats}. P-value adjustment was done with the p.adjust function in the package {stats}.

Stability analyses

To understand whether and how warming affects the stability of the constructed MENs, multiple indices were used to characterize the stability of the networks and their embedded members, including robustness, vulnerability, node and link constancies, node persistence and compositional stability. We assessed the stability of the networks and the microbial communities in the networks on the basis of both simulations of extinction and empirical observations.

Network stability on the basis of simulation

Robustness

Robustness of an MEN is defined as the proportion of the remaining species in this network after random or targeted node removal^9,11. For simulations of random species removal, a certain proportion of nodes was randomly removed. For simulations of targeted removal, certain numbers of module hubs were removed. To test the effects of species removal on the remaining species, we calculated the abundance-weighted mean interaction strength (wMIS) of node i as

where b_j is the relative abundance of species j and s_ij is the association strength between species i and j, which is measured by Pearson correlation coefficient. Thus, in this study, s_ij = s_ji. After removing the selected nodes from the MEN, if wMIS_i = 0 (all the links to species i have been removed) or wMIS_i < 0 (not enough mutualistic association between species i and other species for its survival), node i was considered extinct/isolated and thus removed from the network. This process continued until all species had positive wMISs. The proportion of the remaining nodes was reported as the network robustness. We measured the robustness when 50% of random nodes or five module hubs were removed.

Vulnerability

The vulnerability of each node measures the relative contribution of the node to the global efficiency. The vulnerability of a network is indicated by the maximal vulnerability of nodes in the network as

where E is the global efficiency and E_i is the global efficiency after removing node i and its entire links⁵⁶. The global efficiency of a graph was calculated as the average of the efficiencies over all pairs of nodes:

where d(i,j) is the number of edges in the shortest path between node i and j. Efficiency describes how fast information spread within the network. In ecological networks, it can provide information on how fast the consequence of biological/ecological events traverse to parts or the entire network.

Network stability on the basis of empirical data

Node constancy and overlap

Constancy measures temporal stability of species and is defined as μ/σ, where μ is the mean of abundance over time and σ is the standard deviation⁸³. Here, we calculated the constancy of node i as

where μ_i is the mean and σ_i is the standard deviation of abundances of species i in the networks from different time points. The abundance of species i at a certain time point was positive only if species i was in the MEN at that time point. Otherwise, the abundance of species i was considered zero for that time point. For σ_i = 0, the constancy was not finite and thus removed from subsequent analyses. Finally, the average of all the constancy values was reported as the constancy of network nodes under warming or control. The difference of link constancy between treatments was tested using a Mann–Whitney U test (R package {stats}).

Then we followed the method introduced by Hui et al.⁸⁴ and calculated the number of overlapping nodes among multiple networks. The time points included were referred to as ‘order’. Higher numbers of overlapping nodes among networks indicate slower turnover of species composition in the networks.

Link constancy

We adapted the concept of node constancy and calculated the link constancy in a similar way. We let l_ij+ = 1 if nodes i and j were positively linked in a network, l_ij– = 1 if nodes i and j were negatively linked in a network. Let l_ij+ = l_ij– = 0 if there was no link between i and j. The constancy of a link, for example, l_ij+ or l_ij- was calculated as

where \(\mu _{l_{ij + }}\) is the mean and \(\sigma _{l_{ij + }}\)is the standard deviation of l_ij+ in all the networks from different time points, and \(\mu _{l_{ij - }}\) and \(\sigma _{l_{ij - }}\)are the mean and standard deviation of l_ij– in all the networks from different time points. For \(\sigma _{l_{ij + }} = 0\), the constancy of l_ij+ was not finite and thus removed from subsequent analyses, and similarly for l_ij-. Finally, the average of all the constancy values was calculated as the constancy of network links under warming or control. The difference of link constancy between treatments was tested using a Mann–Whitney U test (R package {stats}).

Node persistence

The node persistence is defined as the proportion of coexisting species (over the total number of species) at an ecological regime⁶. We thus calculated the persistence of a node as the percentage of nodes present in the network in multiple consecutive years, measured by

where v is the number of samples taken from the same field plot at multiple consecutive time points. S is the total OTU number in the networks, and δ_i,k is a Dirac delta function with δ_i,k = 1 if the abundance of OTU k in sample i is larger than 0 and δ_i,k = 0 if OTU k is not present in sample i.

Compositional stability

The compositional stability evaluates the change in community structure over time³¹. We calculated the compositional stability for microbial community in the networks, using the sample × OTU matrix as

where v is the number of samples taken from the same field plot at multiple consecutive time points. S is the total OTU number in the networks. y_i,k is the abundance of OTU k in sample i. If community structure has no change, that is,\(y_{i,k} = \mathop {{\min }}\limits_i y_{i,k}\), this stability index is equal to 1; If community structure is totally different among time points, that is, \(\mathop {{\min }}\limits_i y_{i,k} = 0\), this stability index is 0. When v = 2, it is equal to the square root of the Bray–Curtis similarity between two adjacent time points. Different v (the number of time points) enables it to assess stability at different timescales, so it is called the multi-order stability index.

When v = 2,

where y_i,k and y_{i – 1,k} are abundances of OTU k in samples from Year i and Year i – 1, respectively.

Microbial functional potential assessment

We used two methods to assess functional potentials harboured in the networks. First, microbial community DNAs were hybridized with a comprehensive functional gene array named GeoChip 5.0 M^19,85. GeoChip 5.0M enables detection and quantification of the abundance of >365,000 functional genes from 1,447 gene families involved in biogeochemical cycling of C, N, sulfur, phosphorus and other functional categories. The hybridization signal intensities can provide quantitative measurements of the absolute abundance of the corresponding genes in the community. GeoChip has been widely used in various studies^15,20,52 to understand the functional potentials of complex microbial communities across different environments. Second, we used PICRUSt2³² software, which predicts metabolic pathways on the basis of taxonomy annotation from 16 S rRNA gene sequences, in its default settings to calculate the abundances of metabolic pathways. Since this information is predicted but not directly measured⁸⁶, we focused only on the overall change of the pathways between warming and control, but did not rely on detailed functions in these predicted genomes to interpret their specific influence on the ecosystem processes. To assess the warming effect on the functional potentials predicted by PICRUSt2, we used a general linear model with R package {stats} ‘Pathway abundance ~ Warming*Year + Block’ to assess warming and sampling year’s effect on the predicted pathway abundances.

Differences of the networked community composition under warming and control

To assess whether the networked communities are different under warming and control, we performed three non-parametric multivariate analyses of dissimilarity, including MRPP, ANOSIM and Adonis, on the basis of Bray–Curtis distance (R package {vegan}⁶⁸). Mantel tests were also performed between networked community structures and time, plant and soil variables, and ecosystem C fluxes in R with the mantel function in the package {vegan} version 1.4-2⁶⁸. We further used detrended correspondence analysis to visualize the differences of the networked communities across time and treatment conditions in a two-dimensional ordination space (R package {vegan}⁶⁸). To assess how the soil, plant variables and spatial distance might have shaped the structure of the networked communities, a CCA model was used, followed by a VPA to discern the contributions of these variables to the overall variations of the networked communities (R package {vegan}⁶⁸). In addition, we analysed the taxonomic composition of the networked communities under warming and control conditions at the phylum and class levels, used Mann–Whitney U tests to evaluate the abundance change of each taxa under warming (R package {stats}) and calculated the differentials of their relative abundances under warming compared with control using a method proposed to avoid compositional data bias with the proposed Songbird software²⁶.

Reporting Summary

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

16S rRNA gene sequences were deposited to the National Center for Biotechnology Information (NCBI) under the project accession number PRJNA331185. The OTU table and OTU representative sequences, soil physical and chemical attributes, and plant biomass and richness are downloadable online at http://www.ou.edu/ieg/publications/datasets. GeoChip signal intensity data can be accessed through the URL (https://www.ou.edu/ieg/publications/datasets). Source data are provided with this paper.

The R scripts and Python 3 scripts are publicly available on GitHub at https://github.com/Mengting-Maggie-Yuan/warming-network-complexity-stability with the identifier https://doi.org/10.5281/zenodo.4383469.

We thank numerous former and current members in the Institute for Environmental Genomics for their help in maintaining the long-term field experiment. This work was supported by the US Department of Energy, Office of Science, Genomic Science Program under award numbers DE-SC0004601 and DE-SC0010715, and the Office of the Vice President for Research at the University of Oklahoma. X.G. and X.Z. were generously supported by China Scholarship Council (CSC) to visit the University of Oklahoma. The statistical analyses performed by X.G. were also supported by the China Postdoctoral Science Foundation (2018M641327 and 2019T120101).

1. These authors contributed equally: Mengting Maggie Yuan, Xue Guo, Linwei Wu, Ya Zhang.

Affiliations

1. Department of Environmental Science, Policy, and Management, University of California, Berkeley, CA, USA

   Mengting Maggie Yuan

2. Institute for Environmental Genomics and Department of Microbiology and Plant Biology, University of Oklahoma, Norman, OK, USA

   Mengting Maggie Yuan, Xue Guo, Linwei Wu, Ya Zhang, Naijia Xiao, Daliang Ning, Zhou Shi, Xishu Zhou, Liyou Wu & Jizhong Zhou

3. State Key Joint Laboratory of Environment Simulation and Pollution Control, School of Environment, Tsinghua University, Beijing, China

   Xue Guo & Yunfeng Yang

4. Gladstone Institutes and Chan-Zuckerberg Biohub, San Francisco, CA, USA

   Zhou Shi

5. Key Laboratory of Biometallurgy of Ministry of Education, School of Minerals Processing and Bioengineering, Central South University, Changsha, China

   Xishu Zhou

6. Center for Microbial Ecology, Michigan State University, East Lansing, MI, USA

   James M. Tiedje

7. School of Civil Engineering and Environmental Sciences, University of Oklahoma, Norman, OK, USA

   Jizhong Zhou

8. Earth and Environmental Sciences, Lawrence Berkeley National Laboratory, Berkeley, CA, USA

   Jizhong Zhou

Contributions

All authors contributed intellectual input and assistance to this study. The original concepts were conceived by J.Z. and J.M.T. Field management was carried out by M.M.Y., X.G., Linwei W., Y.Z., Z.S., D.N. and Liyou W. Sampling collection, soil chemical and microbial characterization were carried out by M.M.Y., X.G. and X.Z. Data analyses were done by M.Y., X.G., Linwei W., Z.S. and N.X. with the assistance provided by D.N. and J.Z. All data analysis and integration were guided by J.Z. The manuscript was prepared by J.Z, M.M.Y., X.G., Linwei W. and Y.Z. with substantial input from J.M.T. and Y.Y. Considering their contributions in terms of site management, data collection, analyses and/or integration, M.Y., X.G., Linwei W. and Y.Z. were listed as co-first authors.

Corresponding author

Correspondence to Jizhong Zhou.

Competing interests

The authors declare no competing interests.

Peer review information Nature Climate Change thanks Johannes Bjork, Dongmei Xue and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Reprints and Permissions