The sampling area of this study extends about 1,092 km from approximately 30°N to 40°N and 109°E to 122°E across the North China Plain, with an average annual temperature varying between 8 and 15°C, and the average annual precipitation varies between 500 and 1,000 mm. This is the most important grain producing region of China, mostly planted with wheat and maize in a rotation, and supported by irrigation, fertilization and high yielding crop varieties (Chen et al., 2004). The soils from the sampling sites were classified as Ochric Aquic Cambosols (within the Chinese soil taxonomy) in this area (Zhu et al., 2005). A total of 90 samples were collected from 24 sites in North China Plain during the winter season (the 20th–30th of November 2014) (Supplementary Figure 1A), when winter wheat was at the tillering stage and fertilizer was not yet applied to the field, which effectively avoided the strong effect of fertilizer on microbial community. In order to ensure the integrity and uniformity of soil pH gradient, we sampled a varying number of plots (each plot is 100 m × 100 m) for each site (within 100 km) (Supplementary Figure 1B and Supplementary Table 1), and collected 12 cores within each plot at a depth of 0–15 cm and mixed as one sample. All the samples were stored in ice boxes and brought back to the laboratory within several hours to minimize temperature changes on the way before they could be stored in the freezer at −20°C. Soil samples were sieved using a 2 mm mesh to remove roots and stones, homogenized in the laboratory.

For each sample, a total of 25 edaphic factors were tested, while only those with Spearman correlation coefficient less than 0.6 were retained for subsequent analysis (Supplementary Figure 2). Finally, environmental factors include soil pH, SM (soil moisture), OC (organic carbon), DOC (dissolved organic carbon), DON (dissolved organic nitrogen), NH₄⁺ (ammonium), NO₃^– (nitrate), TP (total phosphorous), TK (total potassium), AP (available phosphorous), AK (available potassium), EC (electrical conductivity), K (potassium), Cr (chromium), Mn (manganese), Cu (copper), Zn (zinc), Cd (cadmium), Pb (plumbum), and As (arsenic). Specific measurement methods for each factor have been presented in detail in the earlier published article (Shi et al., 2018). The regional map of North China Plain was obtained from the Resource and Environment Science and Data Center¹.

DNA was extracted from 0.5 g soil using a MoBio Power Soil DNA extraction kit (MoBio Laboratories, Carlsbad, CA, United States) by following the manufacturer’s instructions, purified with an Ultra Clean 15 DNA purification kit (MO BIO), and stored at −20°C. Bacterial community analysis was carried out using a 16s rRNA genes primer pairs 515F (5′-GTGCCAGCMGC CGCGGTAA-3′)/907R (5′-CCGTCAATTCCTTTGAGT TT-3′) for the V4 hypervariable regions (Biddle et al., 2008), the fungal ITS2 region was amplified by primer sets ITS3 (5′-GCATCGATGAAGAACGCAGC-3′)/ITS4 (5′-TCCTCCGCTTATTGATATGC-3′) (Gade et al., 2013). DNA concentration was measured on a Nano-Drop ND-1000 spectrophotometer (Thermo Scientific, United States). After sequencing, sequences were analyzed using the QIIME pipeline² (Caporaso et al., 2010). The low-quality sequences that had a quality score < 20, contained ambiguous nucleotides, or did not match the primer and barcode, were removed. Operational taxonomic units (OTUs) were generated at 97% similarity cutoff using the UCLUST method in QIIME (Edgar, 2010). Using Greengenes database³ to annotate taxonomic information for each bacterial sequence, UNITE database (Kõljalg et al., 2005) to identify fungal taxonomy.

Meta-community co-occurrence network was constructed using the “WGCNA” package (Langfelder and Horvath, 2012). All core bacterial and fungal taxa in the top 20% relative abundance and presented in more than 60% of all soil samples were included in the network analysis, while those taxa having no robust correlation relationships (Spearman’s correlation coefficients of > 0.6 and false-discovery-rate-corrected P values of < 0.01) with other taxa were lost during the generation of the networks. Sub-networks for each soil sample were also generated by subgraph function in the “igraph” package in R (Csardi and Nepusz, 2006). These network images were visualized with Gephi⁴ (Bastian et al., 2009). To describe the topology structure of the network, we calculated a series of topological features (Supplementary Table 2). Node-level topological features describe the ecological location information of each node in the network, nodes with higher node-level topological feature values occupy the core position in the network, whereas lower value nodes are in the peripheral position (Ma et al., 2016). Network-level topological features with a high value (such as number of edges, average degree, clustering coefficient and density) indicate a more complex and connected network, whereas those with lower values (such as average path length and modularity) suggest closer connections and more concentrated within the network (Barberán et al., 2012; Ma et al., 2016).

To test the significance and importance of the environmental variables for network-level topological features, first, multiple regression model (MRM) on Euclidean distance matrices with the R “ecodist” package (Martiny et al., 2011) were used to assess the relative contribution of each non-collinear edaphic variable in shaping the change pattern of overall network-level topological features. Then, kriging interpolation maps were performed in ArcGIS 10.4⁵ to estimate the geographic patterns of the co-occurrence network-level topological features. Pearson correlation coefficients and P values of the predicted values of topological attributes and observed values at the point calculated by cor.test function in “stats” package in R (Field et al., 2012) were used as the result of cross validation (“CV”), and a simple linear regression between pH and the network level topological characteristic was shown in each map. Additionally, to test the potential roles of soil pH in bacterial and fungal community assemblages, we clustered the core taxa into two ecological preference groups (High pH and Low pH) based on the significant Spearman correlations (P < 0.05) between OTUs abundance and environmental factors (Chen et al., 2021). Random forest analysis was performed to identify the major environmental factors contributing to the relative abundance of dominant taxa in high-pH and low-pH cluster for soil microbiota. The analysis was performed using the randomForest function in the “randomForest” package in R (Liaw and Wiener, 2002).

In order to further illustrate the effect of pH on the network structure, we identified the ratio of link to node as network complexity and compared the network stability of high and low pH clusters. A regression analysis between network complexity and pH based on Pearson correlation was further conducted. Robustness, as a measure of network stability, was calculated by removing nodes in the static network to estimate how easily robustness degraded. Edge information of the high and low clusters pH were obtained and then nodes were randomly and repetitively removed, and natural connectivity of the nodes was used to assess the network robustness (Peng and Wu, 2016). Moreover, to quantify habitat specialization for high and low pH clusters communities, Levins’ niche breadth (B) index was calculated using the formula:

where Bi represents the habitat niche breadth of species i; N represents the total number of samples; Pij is the proportion of species i in sample j (Levins, 1968; Pandit et al., 2009; Wu et al., 2018). The habitat niche breadth was calculated using the “niche.width” function in “spaa” package of R. Boxplot were used to illustrate the habitat niche breadth of high and low pH clusters, with t-test were used to reveal the significant difference between them.

The within-module connectivity (Zi) and among-module connectivity (Pi) of each node are calculated based on Markov clustering algorithm with the “rJava” package in R to reflect the topological role of each node, the Zi reflects how close a node is connected to other nodes within its own module, and Pi describes how close a node contacts with different modules. The topological roles of different nodes can be categorized into four types: peripherals (Zi ≤ 2.5, Pi ≤ 0.62), connectors (Zi ≤ 2.5, Pi > 0.62), module hubs (Zi > 2.5, Pi ≤ 0.62) and network hubs (Zi > 2.5, Pi > 0.62) (Deng et al., 2012; Cong et al., 2015). Generally, connectors, module hubs and network hubs are considered as putative keystone taxa of ecological network (Olesen et al., 2006). These relatively rare but highly connected taxa play an important role in improving soil functional potential in agricultural ecosystems (Shi et al., 2020). Spearman correlation and multiple regression model were performed to explore the determinants of the identified keystone taxa.

Across 90 soil samples collected from the North China Plain (Supplementary Figure 1), we obtained a total of 1,800,450 and 996,840 high-quality bacterial and fungal sequences, respectively, which were clustered into 65,761 bacterial operational taxonomic units (OTUs) and 4,033 fungal OTUs based on 97% sequence similarity. In this study, we selected OTUs that accounted for the top 20% in terms of relative abundance and occurring in more than 60% of all samples to construct co-occurrence networks (Jiao et al., 2019). In total, we identified 1,129 bacterial and 191 fungal dominant taxa (Supplementary Table 3), which occupied a small part of the individual total taxa (bacteria: 1.7%; fungi: 4.7%), but accounted for the majority of the total number of sequences, respectively (bacteria: 70.9%; fungi: 86.3%), which can maximize the representation of the whole community without data redundancy. The majority of bacterial sequences belonged to the phyla Proteobacteria (33.8%), Actinobacteria (21.6%), Acidobacteria (16.2%), Chloroflexi (7.2%), Gemmatimonadetes (6.3%) and Planctomycetes (3.9%). And the fungal sequences were mainly from the phyla Ascomycota (91.1%) and Basidiomycota (5.2%).

The constructed co-occurrence network consisted of 1,320 nodes (OTUs), including 1,129 bacterial OTUs and 191 fungal OTUs. A total of 33,032 edges (associations between OTUs) were inferred for the consensus microbial network (Figure 1A). Supplementary Figure 3 showed that the curves of network connectivity distribution fitted well with the power-law model, which was indicative of scale-free networks and significantly different from random networks. The associations were mainly observed among phylum Proteobacteria, Actinobacteria, Acidobacteria and Chloroflexi (Supplementary Figure 4 and Supplementary Table 4). The presence of positive edges (65.3%) is much greater than negative edges (34.7%) in the network. Bacterial nodes had the most internal connections, fewer edges between bacteria and fungi, and the fungal nodes had the fewest internal connections (Figure 1A). Moreover, values for the node-level topological features, including degree, eigenvector and closeness centrality were significantly higher (P < 0.001) in bacterial OTUs than in fungal OTUs (Figure 1C), which suggested that bacterial taxa played a more critical role in the network than fungal nodes.

Additionally, to explore the ecological preferences of soil microorganisms in the North China Plain, we divided the dominant bacterial and fungal taxa into two ecological clusters: high and low pH based on the Spearman’s correlations (P < 0.05) with soil pH, which had different habitat preferences (Supplementary Table 5). Co-occurrence network captured 730 nodes belonging to the high pH cluster and 416 nodes for low pH cluster, these taxa within the cluster that shared the same habitat and environmental preferences tended to co-occur (Figure 1B). Each of the ecological clusters consisted of multiple orders of soil microbiota. iii1–15, Actinomycetales and Rhizobiales tended to aggregate in soils with high pH, Gaiellales, Solibacterales and Burkholderiales prederred low soil pH, while Xanthomonadales and Rhodospirillales were relatively abundant in both high-pH and low-pH soils (Supplementary Table 5). Random forest analysis revealed that pH was the dominant environmental factor affecting the relative abundance of high and low pH clusters for bacterial and fungal communities (Supplementary Figure 5). Regarding the node-level topological features, the values of all parameters (including degree and betweenness, eigenvector and closeness centrality) in the low pH cluster were significantly higher than those in the high pH cluster (Figure 1D). Similar trends were seen in subnetworks constructed for bacterial community (Supplementary Figure 6A). These results showed that taxa of the low pH cluster, in comparison to those of the high pH cluster, were more often occupied the central positions in the network. Putative keystone species composition information below also confirms this conclusion (Supplementary Table 6). However, there was no significant difference in the node-level topological characteristics between high and low pH clusters in the fungal subnetwork (Supplementary Figure 6B). Moreover, when comparing the subnetwork stability of the high and low pH clusters, the natural connectivity at low pH was greater than that at high pH after removing the same proportion of nodes, which indicates the subnetwork of low pH cluster was more stable (Figure 2A).

To examine the relative contribution of multiple environmental factors to the co-occurrence patterns on the North China Plain, we singled out 20 relatively independent (Spearman’s correlations: ρ² < 0.6) environmental factors based on correlation (Supplementary Figure 2), and calculated a series of network-level topology characteristics for each soil sample. Based on correlation and hierarchical clustering analysis, the network-level topological features could be divided into two clusters (Figure 3A). The first cluster included number of edges, clustering coefficient, density, average degree, and degree centralization. The second cluster contained betweenness centralization, number of nodes, modularity, average path length, and eigenvector centralization. The first cluster of network-level topological features was positively correlated with soil pH, TP, EC, Cd, and As, and negatively correlated with NH₄⁺, AP, and Pb, whereas the second clusters were on the contrary. Multiple regression on distance matrices (MRM; Figure 3B) showed that the contribution of pH (R² = 34.2%, P < 0.001) to the network-level topological features overwhelmed than of other environmental factors. And the variance interpretation of each environment variable was shown in Table 1. The variation of each network-level topology characteristics could be well explained by the environmental variables. However, the explanatory quantities were not the same among the parameters, soil pH contributed the most toward explaining the variation in the topological features.

Furthermore, in order to better visualize the relationship between pH and network-level topological characteristics, we performed Kriging interpolation on pH and each network-level topological feature to obtain the spatial distribution maps (Figure 4). The predicted spatial patterns showed that the network-level topological features in the first cluster had higher values in the low pH regions than high pH regions. In contrast, the features in the second cluster were higher in high-pH soil than in low-pH soil. These observations were confirmed by multiple regression analysis. Significant and negative linear regressions were found between soil pH and the network-level topological features in the first cluster. However, the mostly features in the second cluster strongly increased with increasing pH, while the number of nodes peaked at neutral pH (Figure 4). Altogether, these results indicated that the ecological network of dominant bacterial and fungal taxa was more complex in low pH regions than that of high pH regions. In order to avoid the effect of microbial abundance to the ecological network when illustrate the change of the network pattern along pH gradient, the ratio of links to nodes were calculated as the network complexity. The result of Pearson correlation between network complexity and pH reconfirmed that the network of dominant bacterial and fungal taxa was significantly affected by pH, with an R² of 0.7057 (P < 0.001) and was more complex in low pH than that of high pH (Figure 2B). For the sub-network generated for bacteria, the same network-level topological characteristics and pH variation patterns were observed, but not for fungal subnetworks (Supplementary Figures 7, 8).

In addition, to determine how the number of connections between subnetworks changed, we applied regression analysis and found there were significantly negative relationships between soil pH and the number of positive/negative edges within the bacterial community. Moreover, the number of positive edges between the two kingdoms decreased with increasing pH, whereas the number of negative edges peaked at neutral pH. For the fungal community, there was no significant correlation between the number of positive and negative edges and pH (Supplementary Figure 9). The links between the number of positive and negative associations within and between communities and each environmental factor were detected by Mantel tests (Figure 5). Number of connections in the bacterial network were significantly correlated with pH, TP, AP, EC, and Cd. The number of positive edges between kingdoms was significantly affected by pH and Cd, and AK was correlated to the number of negative edges between kingdoms. However, no significant relationships among fungal community connections number and soil factors were observed.

Keystone nodes (network hubs, module hubs, and connectors) were identified by analyzing the topological roles that each node played in the network of dominant bacterial and fungal taxa (Figure 6A). In this study, a total of 2 network hubs, 26 module hubs, and 17 connectors were detected in the microbial network. Both network hubs belonged to the phylum Proteobacteria. The detected module hubs were composed of six taxa within the Proteobacteria, five taxa within Ascomycota, two taxa within Firmicutes, and one taxon within each of the phyla Acidobacteria, Chloroflexi, Gemmatimonadetes and Zygomycota. The connectors consisted of eight taxa within the Actinobacteria, six taxa within Proteobacteria, three taxa within Acidobacteria, two taxa within Planctomycetes, Gemmatimonadetes, Chloroflexi, one taxon within Verrucomicrobia, Ascomycota, and Bacteroidates (Supplementary Table 6). Additionally, significant and negative linear correlations were found between pH and the relative abundance of keystone taxa (Supplementary Figure 10), suggesting that increasing soil pH reduced the number of keystone taxa significantly correlated with network complexity.

To explore the linkages between environmental parameters and the relative abundance of keystone taxa, we then correlated the keystone taxa to soil properties and identified the major drivers for each keystone taxon (Figure 6B). Spearman correlation analysis revealed that there were distinct environmental preferences between keystone taxa. Compared with network hubs and module hubs, connectors had less correlation with environmental variables, suggesting that they were more resistant to environmental changes. Moreover, the results of multiple regression model showed that pH was the most important predictor for the relative abundance of keystone taxa, TP, SM, OC, Cu, DOC, and NO₃^– also played important roles.

The role of environmental pH shifts is a prominent topic in the study of interspecies interactions. Following on from the study by Shi et al. (2021) which explored high and low pH as discontinuous and discrete categories, we found that in fact the topological characteristics of the co-occurrence network influenced by soil pH vary along a continuous gradient, and that the influence of pH was far greater than any other measured environmental factors. This is broadly consistent with recent studies of soybean fields showing that soil pH plays a key role in rhizosphere microbial interactions (Zhang et al., 2018).

Soil pH shapes microbial metabolism in different ways. Firstly, as an integrating index to measure soil conditions (Lauber et al., 2009), soil pH is highly correlated with a wide range of biogeochemical conditions [e.g., NH₄⁺, TP, AP and many metal ions (Figure 4)], which significantly affect the environmental conditions related to microbial growth and survival (Pan et al., 2014; Ma et al., 2018). The pH also affects the activity of extracellular enzymes and the reactivity of natural organic matter (Paul et al., 2006). Secondly, soil pH directly imposes a physiological constraint on soil microorganisms (e.g., in affecting homeostasis of intracellular pH), and any slight change in pH value unit will have a significant impact on microbial growth and metabolic activity (Fernandez-Calvino and Baath, 2010). Additionally, pH may also affect the rate of energy expenditure in microbial respiration and hence microbial community structures by regulating the thermodynamics and kinetics of redox reactions (Jin and Kirk, 2018). Therefore, it is reasonable that one would find a close connection between pH and microbial interactions at the community level.

In general, taxa of microorganisms (especially prokaryotes) show preferences for ranges of soil pH, with the greatest species diversity and abundance in neutral soils, “acidity specialists” exist in acidic soils as well as “alkalinity specialists” in alkaline soils (Barberán et al., 2012; Jones and Bennett, 2017; Vieira et al., 2020). Among the low-pH cluster of microbial taxa that we detected, the abundance of Xanthomonadales has been found to be significantly negatively correlated with pH (Zhang et al., 2020), while Nocardioidaceae and Sordariales from the high-pH cluster, favors high-pH environments (Huang et al., 2017). In addition, we found a unique bacterial class Gammaproteobacteria, which proved to be widespread in wheat fields with low pH and high NO₃^–-N concentrations (Hamamoto et al., 2018). The main families within this class were Chromatiaceae, Sinobacteraceae, and Xanthomonadaceae, all of which have been reported as only being found in agricultural soils (Kuramae et al., 2012), which could be used as indicator groups in farmland ecosystems.

Co-occurrence networks represent strong correlations between OTUs, indicating potential species-species interaction in communities (Ma et al., 2018). In the present study, we found that bacteria and fungi were more closely related within their respective community, suggesting a wide range of interactions between species such as the exchange of metabolites (Woyke et al., 2006) and cooperation in biofilm construction (Rodríguez-Martínez and Pascual, 2006). We also observed some correlations between microbes in the two kingdoms. It is no surprise that several bacterial taxa, such as Burkholderia, Sphingomonas and Pseudomonas, have been found to interact with fungi in soil (Gamalero et al., 2008; Warmink and van Elsas, 2008). In general, the products (such as water-soluble sugars and phenolic compounds) released by fungi during the degradation of refractory organic matter (lignin and cellulose) in soil can be utilized by bacteria (Boer et al., 2005). Although fungi can produce a portion of carbon, plant root exudates (i.e., amino acids, sugars and organic acids) are considered to be the main carbon source for bacteria (Philippot et al., 2013). Therefore, the nutrient dependence of bacteria on fungi in farmland ecosystems may be low, resulting in a low proportion of bacteria-fungus edges in the co-occurrence network. In addition, most of the network links were between bacteria: often linking Proteobacteria, Actinobacteria and Acidobacteria with other bacterial phyla. Amongst the fungi, there were many positive associations among Ascomycota, which is consistent with previous results, that is, phylogenetically related taxa are also ecologically related (Barberán et al., 2012). This suggests that these ascomycetes may either have synergies or share similar ecological niches in the farmland environments.

In this study, we observed that the number of nodes in the co-occurrence network peaked at about pH 7, which is consistent with previous studies confirming the highest diversity of soil bacterial communities with near-neutral pHs (Lauber et al., 2009). However, microbial taxa were more closely correlated in the low pH area than in the high pH area, this may be due to the close contact between Acidobacteria and other microbial groups (Supplementary Table 4), which has been shown to increase toward lower pH (Jones et al., 2009). In addition, soil acidity supports basic soil properties and functions, such as the solubility of exchangeable ions and nutrients, or the use of internal and external acid gradients by bacteria to produce ATP and rotate flagella motors (Braus and Whitman, 2021). As far as is known, the large pH gradient between the inside and outside of the cell causes microbes living in an alkaline environment to greatly reduce the proton dynamics (PMF) used to produce ATP—much less so than for acidophiles (REF). This makes it energetically less favorable for alkaliphiles to generate energy through oxidative respiratory chains and oxidative phosphorylation (Padan et al., 2005), and results in a decrease in the frequency of microbial interaction in the community.

Understanding ecological niches is important for determining the mechanisms of community assembly (Jiao and Lu, 2020). Here, we further investigated the differences of environmental adaptive capabilities between different ecological clusters (high pH and low pH). We found that the niche breadth of low pH cluster was narrower than that of high pH cluster (Supplementary Figure 11), implying that interspecific interactions (such as competition, etc.) are stronger in low pH clusters, thus the co-occurrence network is more closely connected. Deterministic processes tend to have a stronger impact on the habitat specialists with a narrow niche breadth than on generalists with a wide niche breadth (Pandit et al., 2009). Therefore, it is perhaps unsurprising that the lower pH regions have more stable co-occurrence network structures. However, this observation is inconsistent with previous studies on the Tibetan Plateau, and we reason that this difference might be mainly since Chen et al. (2021) explored the associations of microbial communities in the pH range under the influence of different vegetation types (alpine steppe and alpine meadow) and did not study a unified microhabitat.

It is also interesting that in our NCP study, the pH-related trend is essentially dominated by bacteria. Fungi by themselves showed no significant trend, while bacteria alone—and fungi in combination with bacteria—showed the same clear trend. Generally, fungal community structure shows much weaker trends in relation to pH than bacteria, and fungi generally exhibit wider pH ranges for optimal growth. For instance, the fungal genus Mariannaea, has been found to grow between pH 5.5–8.5 in culture-based studies (Domsch et al., 1980), which may relate to the larger cells and much lower surface to volume ratio of fungi, which enables better cellular pH homeostasis (Beales, 2004).

Unlike the trend in network connectivity, the associations within communities were most stable at neutral pH (Supplementary Figure 9), where they were dominated by negative associations. This greater stability may reflect the greater available energy and resources for niche specialization by bacterial cells at around neutral pH, since high and low pH are each thought to act as a drain on resources due to the requirements for cellular homoeostasis, relative to neutral pH (Lauber et al., 2009). The high ratio of negative interactions may reflect narrow specialized niches and strong competitive exclusion at around neutral pH.

Through network analysis, keystone taxa can often be detected as network hubs, module hubs, and connectors (Olesen et al., 2007; Faust and Raes, 2012). Only a few hubs or connectors have a wide distribution across different plots, supporting the environmental dependency theory that keystone taxa do not dominate anywhere or at any time, but only play a key role in specific environments (Abdul Salam et al., 2013). We noted that the relative abundance of keystone taxa decreased as soil pH increased, consistent with the previous observation that the putative keystone taxa change with conditions (Lu et al., 2013). Although they were not abundant, these low abundance keystone taxa may have disproportionately large impacts on ecosystem functions and services (Lynch and Neufeld, 2015). Such a network hub, Oxalobacteraceae, has been previously found to be closely related to plant growth and nitrogen acquisition (Yu et al., 2021). Module hub Burkholderiales and connector Rhizobiles, are dominant members of the rhizosphere microbiome and are considered keystone taxa across different ecosystems (Banerjee et al., 2018). Another connector, Mortierella, has the ability to decompose complex organic substrates and with an important phosphate lysozyme that promotes soil activity and plant growth (Zhang et al., 2011). In addition, taxa within the phyla Gemmatimonadetes and Planctomycetes are often associated with the rhizosphere of plants and indirectly promote plant growth by participating in a variety of ecological processes, which have also been identified as kinless hubs in previous studies (Shi et al., 2020). The presence of these keystone taxa could thus be critical to maintaining soil health and crop productivity. Hence, manipulation of keystone taxa in microbial network structures, such as the addition or removal of functionally specific keystone taxa, may provide a promising approach for agricultural management to improve crop yields (Banerjee et al., 2018).