Molecular ecological network analysis of the response of soil microbial communities to depth gradients in farmland soils

Abstract

Soil microorganisms are considered to be important indicators of soil fertility and soil quality. Most previous studies have focused solely on surface soil, but there were numerous active cells in deeper soil layers. However, studies regarding microbial communities in deeper soil layers were not comprehensive and sufficient. In this study, phylogenetic molecular ecological networks (pMENs) based on the 16S rRNA Miseq sequencing technique were applied to study the response of soil microbial communities to depth gradients and the changes of key genera along 3 meter depth gradients (0–0.2 m, 0.2–0.4 m 0.4–0.6 m, 0.6–0.8 m, 0.8–1.0 m, 1.0–1.3 m, 1.3–1.6 m, 1.6–2.0 m, 2.0–2.5 m, and 2.5–3.0 m). The results showed that the modularity of microbial communities was consistently high in all soil layers and each layer was similar, which indicated that microbial communities were more resistant to depth changes. The pMENs further demonstrated that microbial community interactions were stable as the depth increased and they cooperated well to adapt to changes in different soil gradients. This was evidenced by similar positive links, average degree, and average clustering coefficient. In addition, key genera were obtained by analyzing module hubs in the pMENs. There may be at least one dominant genus in each layer that adapted to and resisted changes in the soil environment. It seems microbial communities demonstrate a stable and strong adaptability to depth gradients in farmland soils.

INTRODUCTION

Most natural and artificial habitats of terrestrial ecosystems are covered by soils (Bai, Wang, Deng, and He ()). Since precipitation, illumination, and temperature within an area are similar, high spatial heterogeneity and physiochemical properties are responsible for the diverse range of microorganisms found in soil ecosystems (Schimel & Schaeffer, ). Soil microorganisms are considered to be important indicators of soil fertility and quality (Paula, Rodrigues, Zhou, & Wu, ). Recent studies based on culture-independent techniques have shown that large-scale soil microbial diversity and community composition were predominantly driven by soil pH (Griffiths, Thomson, James, & Bell, ) and other soil properties, such as organic matter and salinity (Campbell & Kirchman, ; Tripathi, Kim, Tateno, & Kim, ; Zheng, Xue, Li, & Deng, ). It was known that biogeochemical processes were mainly driven by microorganisms throughout the soil profile; however, previous studies on the structure and diversity of soil microbial communities had solely focused on the top 15 cm of the soil column (Eilers, Debenport, Anderson, & Fierer, ). Numerous studies have suggested that microorganisms living in deeper soil layers were generally considered to be unimportant due to their low biomass density and low activity levels (Blume, Bischoff, Reichert, & Moorman, ; Fierer, Schimel, & Holden, ; Hartmann, Lee, Hallam, & Mohn, ). However, Fierer et al. (Fierer, Schimel, et al., ) noted that not only was a large number of microorganisms present in the underground soil but also that the potential activity of these underground microorganisms might be higher than that of surface microorganisms. Although microbial biomass generally exponentially decreased with depth (Blume et al., ; Fierer, Schimel, et al., ; Hartmann et al., ), in deeper soil layers, there were numerous active cells (Buss, Bruns, Schultz, & Moore, ; Richter & Markewitz, ), which may have a greater impact on the soil formation process than their surface counterparts because of their proximity to parent material (Buss et al., ). Fierer, Allen, Schimel, and Holden () also demonstrated that deep microbial communities differed in composition from those in the surface layer. On a depth-weighted basis, deep microbial communities have also been found to be diverse and abundant (Li, Yan, Tang, & Jia, ; Will, Thürmer, Wollherr, & Nacke, ). Previous studies have reported that in forest soil, soil depth affected the diversity and composition of the archaea community (Thoms, Gattinger, & Jacob, ; Too, Keller, Sickel, & Lee, ). Specific bacterial communities in deep wetland soils have been detected, and the interaction between these communities was more noticeable than that of the surface layer communities (Steinmuller, Dittmer, White, & Chambers, ). However, how a microbial community adapts to soil depth gradients through community-level adjustment of compositions and interactions has not yet been completely elucidated.

Many network analysis approaches have been developed, such as the equation-based network (Gardner, Di, Lorenz, & Collins, ; Yeung, Tegner, & Collins, ), the related/co-expression network (Horvath, Zhang, Carlson, & Lu, ; Oldham, Horvath, & Geschwind, ), and Bayesian network (Gerstung, Baudis, Moch, & Beerenwinkel, ). One frequently studied and widely used method is the phylogenetic molecular ecological networks (pMENs) proposed by Deng, Jiang, Yang, and He (), who were the first to apply metagenomic techniques (such as sequencing and microarrays) (Deng, Zhang, Qin, & Tu, ; Faust & Raes, ). In the entire network, various parameters, such as topology, node modules, topological roles, and network composition, better reflect the relationship between microbial communities and their associated niche functions (Hahn, Konwar, Louca, & Hanson, ; Zhou, Deng, Luo, & He, , ). Thus, in the present study, the pMENs were applied to investigate the responses of soil microbial communities to soil depth gradients and the underlying microbial interactions, including the following: (a) phylogenetic diversity and structure of microbial community shift; and (b) changes of key genera in microbial network interactions along soil depth gradients.

MATERIAL AND METHODS Study sites and soil sample collection

The study sites were located in Tianjin, a municipality close to Bohai Bay, which is influenced by a continental monsoon climate. Both hydrology and climate played the significant roles on the physical and chemical properties of soil. In November 2013, soil samples were collected from eight farmland sites in five regions of Tianjin, China (S1: 39°36′54.49″N, 116°58′04.20″E; S2: 39°32′09.31″N, 116°59′32.31″E; S3: 39°40′37.50″N, 117°21′03.07″E; S4: 39°40′32.66″N, 117°20′39.53″E; S5: 40°04′14.87″N, 117°20′01.90″E; S6: 40°05′25.91″N, 117°38′10.21″E; S7: 38°43′09.78″N, 117°26′44.00″E; S8:38°43′06.33″N, 117°26′20.41″E; Figure ). Soil samples were continuously collected using a 5.5-cm-diameter hollow-stem hand auger in a vertical profile of 10 layers (from the ground surface to a depth of 3 m), including layer A: 0–0.2 m, layer B: 0.2–0.4 m, layer C: 0.4–0.6 m, layer D: 0.6–0.8 m, layer E: 0.8–1.0 m, layer F: 1.0–1.3 m, layer G: 1.3–1.6 m, layer H: 1.6–2.0 m, layer I: 2.0–2.5 m, and layer J: 2.5–3.0 m (the soil profile was over 3 m, as the shallow groundwater table is close to this depth). The 80 collected samples were stored in polyethylene bags. The physicochemical parameters of soil water extracts were determined by adding 200 ml Milli-Q water to 100 g soil sample, followed by shaking (45 min at room temperature), centrifuging (2,057 g), and filtering. The pH and salinity were determined by a portable analyzer (Orion Star A329, Thermo, USA). The cation concentrations of K⁺, Na⁺, Ca²⁺, and Mg²⁺ were determined by an inductively coupled plasma atomic emission spectrometer (Optima 8300, PE, USA), and anion concentrations of ${\text{SO}}_4^{2-}$ and Cl⁻ were determined by ion chromatography (ICS-2100, Dionex, USA). ${\text{NO}}_3^{-}$, ${\text{NO}}_2^{-}$, and ${\text{NH}}_4^{+}$ were determined by a continuous flow analyzer (Auto Analyzer 3, Seal, Germany). Part of the soil samples were freeze-dried and stored at −20°C for genomic DNA extraction.

Study area sampling point diagram

DNA extraction, Illumina MiSeq sequencing and analysis

DNA of the microorganisms in 80 soil samples was extracted by biotechnological methods using the Ezup Genomic DNA Extraction Kit (Sangon Biotech, China, Cat# SK8264). After extraction, a NanoDrop Spectrophotometer was used to detect the concentration and mass of DNA. The DNA extraction was diluted to 10 ng/μl.

The operational taxonomic unit (OTU) data used to construct the pMENs were generated by the 16S rRNA gene. Miseq sequencing of the V4 hypervariable region of the 16S rRNA gene was performed using the universal primers 515F (5'-GTGCCAGCMGCCGCGGTAA-3') and 806R (5'-GGACTACHVGGGTWTCTAAT-3'): Amplification was then measured using a MiSeq sequencer (Caporaso, Lauber, & Walters, ; Lin, De, Li, & Li, ; Zheng et al., ). Ten nanograms of soil genomic DNA, 0.4 mM of each deoxynucleotide triphosphate, 1.0 μM of each primer, 0.5 U of Ex Taq (TaKaRa, Dalian), 1× PCR buffer, and 1.5 mM MgCl₂ were completely mixed in a 25 μl mixture for performing PCR (Zheng et al., ). The PCR amplification was conducted as described by Li, Rui, Mao, and Yannarel (). Under the reading conditions with the unique identification bar code of the sample, sequencing was performed on the Illumina Miseq platform of the Environmental Genome Platform of the Chengdu Institute of Biology. Eighty samples were sequenced using the Reagent Kit v2 2 × 250 bp. The results were processed by QIIME pipeline version 1.7.0 (http://qiime.org/). The original data were first screened, and an average base quality score of 20 was the lowest; when the low-quality scoring sequence was removed, the sequence length was at least 300 bp and without an ambiguous base “N.” These sequences were used for downstream analysis (Zheng et al., ). Next, the chimeric sequence was removed using the UCHIME algorithm (Edgar, Haas, Clemente, & Quince, ). The processed sequence was clustered by the complete-linkage clustering method of the Ribosomal Database project pyro pipeline. A 97% nucleotide sequence similarity cutoff was used to classify OTUs. The RDP classifier was used to assign taxonomy (Wang, Garrity, Tiedje, & Cole, ). Re-sampling was performed to a total of 2,590 reads of the same sequence depth. Chao1 estimator of richness and Shannon index were calculated using the RDP classifier. The original sequence data were deposited at the European Nucleotide Archive by accession PRJEB21751 (http://www.ebi.ac.uk/ena/data/view/PRJEB21751) (Zheng et al., ).

Statistical analysis

Microbial community diversity was compared using IBM SPSS 20 with a univariate ANOVA based on Chao1 richness and Shannon diversity indices. Correlation analysis between microbial communities and environmental variables was performed by the Mantel test, conducted by PCORD 5.0. In addition, based on the study of Fukuyama (Fukuyama, Mcmurdie, Dethlefsen, Relman, & Holmes, ), double principal coordinate analysis (DPCoA) was performed on the relative abundance data of the OTUs.

Establishment and analysis of pMENs

Construction of pMENs is based on the random matrix theory (RMT). The network was built using the online pipeline provided by the Institute of Environmental Genomics at the University of Oklahoma (http://ieg2.ou.edu/MENA) (Deng et al., ; Zhou, Deng, Luo, & He, , ). The establishment of pMENs was divided into nine key steps, namely collection of metagenomic sequences, data standardization, Pearson's correlation estimation, adjacency matrix determination by an RMT-based approach, network characterization, module detection, eigengene network analysis, network comparison, and correlation based on the RMT method (Zhou, Deng, Luo, & He, ). High-throughput metagenomic data were collected and transformed to construct Pearson's correlation matrix (Zhou, Deng, Luo, & He, ), which was further converted into a similarity matrix. The method of RMT is to automatically select the optimal threshold of the network according to the data to build an optimal network. The automatically determined threshold was used to derive the adjacency matrix from the similarity matrix. Using the adjacency matrix to encode the strength of the connection between each pair of nodes, module analysis and network characterization were performed to generate different network topology attributes. The procedure was detailed in a previous MENA study (Deng et al., ). Different topological roles were described by the following two parameters: the within-module connectivity (Zi) and among-module connectivity (Pi). The Zi reflected the degree to which nodes were connected to other nodes within their own module, while the Pi reflected the degree to which nodes were connected to different modules. Since we were primarily interested in the impact of soil depth on network interactions, pMENs were constructed based on the relative abundance of samples obtained from each layer of the 8 sampling sites.

RESULTS

A total of 834,485 high-quality 16S rRNA gene sequences were obtained for 80 samples. They were resampled to 2,590 sequences per sample, which were clustered into 7,403 OTUs. The rarefaction curves showed that the number of samples was reasonable (Figure A1). Both Shannon and Chao1 indices showed no significant differences between the depth layers (Table A1). To understand the effects of soil depth on the linkages between microorganisms, OTU data obtained from MiSeq sequencing was used to construct the microbial pMENs based on 10 layers of the 8 sampling sites (Table ). One node represented an OTU, with each link representing the correlation between two connected nodes. The results presented in Table showed that microbial connectivity had similar thresholds (0.930–0.960). Modularity acted as an indicator of system resistance (Carpenter, Arrow, Barrett, & Biggs, ). In our study, the modularity values were similar (0.771–0.810) and higher than 0.4, which indicated the modular structure of the network (Newman, ) and the strong resistance to environmental changes (Bai et al., ). The numbers of module hubs varied from 1 to 5, which indicated that the degree of soil generalization differed from that in the 10 layers, with smaller values representing lower degrees of generalization (Wang, Zhang, Zheng, & Deng, ). The average clustering coefficient (avgCC) described the relationship between a particular node and its adjacent node (Shen, Huang, Zeng, & Yu, ). The avgCC values varied over a similar range from 0.240 to 0.365 in the 10 layers, showing connections similar to that of the neighbors (Deng et al., ). The avgKK values varied from 3.973 to 5.559, which indicated that each layer of the network had different degrees of complexity. The R-squared values were all greater than 0.700, implying that the RMT-based ecological network should be scale-free (Zhou et al., ). Scale-free implied that most nodes in the network had limited connections with other nodes, while a few nodes had many connections (Sun, Wang, Lin, & Zhou, ). As shown in Table , the percentages of positive links (70%–86%) were much higher than those of the negative links (14%–30%), which indicated that the relationship between microbial communities was based more on mutual cooperation than competition (Newman, ). Alternatively, positive links have also been found to be related to mutualism among genes during coevolution processes (Faust & Raes, ).

Topological properties of the phylogenetic molecular ecological networks of microbial communities in ten layers (A to J)

Depth	No. of original OTUs	Similarity threshold	Modularity	Module	Numbers of module hubs	Average clustering coefficient (avgCC)	Average degree (avgKK)	Total nodes	Total links	R square of power-low	Percentage of positive links	Percentage of negative links
A	7,188	0.93	0.796	39	3	0.240	3.973	375	745	0.857	70	30
B	7,188	0.93	0.771	40	4	0.244	4.105	362	743	0.860	75	25
C	7,188	0.95	0.781	38	3	0.279	5.088	351	893	0.791	77	23
D	7,188	0.93	0.810	33	2	0.328	5.559	354	984	0.841	86	14
E	7,188	0.95	0.790	26	4	0.343	5.190	294	763	0.856	79	21
F	7,188	0.94	0.782	38	4	0.302	4.361	366	798	0.808	76	24
G	7,188	0.94	0.795	35	5	0.267	4.133	331	684	0.902	79	21
H	7,188	0.96	0.796	26	1	0.347	4.070	257	523	0.822	82	18
I	7,188	0.94	0.789	29	5	0.365	4.476	319	714	0.856	78	22
J	7,188	0.93	0.778	31	1	0.341	5.288	351	928	0.795	82	18

Abbreviation: OTU, operational taxonomic unit.

Results of the Pi and Zi analysis, as depicted in the Z-P diagram (Figure ), revealed that different nodes played different roles in the pMENs. From an ecological point of view, each module in the pMENs indicated one niche (Wang et al., ). Peripherals (Zi < 2.5 and Pi < 0.625) referred to specialists, while module hubs (Zi > 2.5 and Pi < 0.625) and connectors (Zi < 2.5 and Pi > 0.625) showed more generalists, and network hubs (Zi > 2.5 and Pi > 0.625) embodied supergeneralists (Deng et al., ; Lewinsohn, Prado, Jordano, & Bascompte, ). No network hubs were found in the Z-P diagram in all ten layers. The vast majority of OTUs in each layer occurred in peripherals, which represented most of their links inside their own modules (Deng et al., ). The module hubs and connectors, which occurred in each layer of the community, as shown in Figure , indicated that the microbial communities in all ten layers were stable. According to the module hubs in the Z-P diagram, corresponding OTUs can be used to determine the key genus of each layer (see Table ). It can be observed that the key genus in each layer was basically different and will be further discussed in the discussion section.

The network of the pMENs of ten layers. The Z-P plot showing the topological roles of each OTU based on Zi (within-module connectivity) and Pi (among-module connectivity) of microbial communities in ten soil layers. According to values of Zi (2.5) and Pi (0.625), the roles of nodes were classified into four categories

Operational taxonomic units (OTUs) and corresponding keystone genus in ten layers (A to J)

Depth	Numbers of module hubs	OTU	Phylum	Class	Genera
A	3	OTU832	Gemmatimonadetes	Gemmatimonadetes	Gemmatimonas
		OTU2510	Acidobacteria	Acidobacteria_Gp25	Gp25
		OTU4954	Bacteroidetes	Bacteroidetes_incertae_sedis	Ohtaekwangia
B	4	OTU71	Spirochetes	Spirochaetia	Treponema
		OTU498	Thaumarchaeota	–	Nitrosopumilus
		OTU774	Proteobacteria	Deltaproteobacteria	–
		OTU1532	Actinobacteria	Actinobacteria	–
C	3	OTU650	Actinobacteria	Actinobacteria	–
		OTU1466	Acidobacteria	Acidobacteria	Gp4
		OTU1457	Proteobacteria	Proteobacteria	Cystobacter
D	2	OTU723	Actinobacteria	Actinobacteria	–
		OTU967	Chloroplast	Chloroplast	–
E	4	OTU631	Proteobacteria	Gammaproteobacteria	–
		OTU1510	Proteobacteria	Betaproteobacteria	–
		OTU1853	Verrucomicrobia	Spartobacteria	Spartobacteria_genera_incertae_sedis
		OTU4767	Chloroflexi	–	–
F	4	OTU396	Spirochetes	Spirochaetia	Treponema
		OTU1362	Acidobacteria	Acidobacteria_Gp7	Gp7
		OTU1366	Actinobacteria	Actinobacteria	Gaiella
		OTU2089	Actinobacteria	Actinobacteria	Gaiella
G	5	OTU316	Proteobacteria	Gammaproteobacteria	Steroidobacter
		OTU1664	Actinobacteria	Actinobacteria
		OTU2980	Gemmatimonadetes	Gemmatimonadetes	Gemmatimonas
		OTU1853	Verrucomicrobia	Spartobacteria	Spartobacteria_genera_incertae_sedis
		OTU2885	Firmicutes	Clostridia	Ruminococcus2
H	1	OTU679	–	–	–
I	5	OTU15	–	–	–
		OTU156	Bacteroidetes	Flavobacteriia	Croceibacter
		OTU570	Proteobacteria	Betaproteobacteria	–
		OTU408	Proteobacteria	Betaproteobacteria	–
		OTU2365	Proteobacteria	Alphaproteobacteria	Sphingomonas
J	1	OTU2930	Acidobacteria	Acidobacteria_Gp18	Gp18

The Mantel test was applied to determine the relationship between microbial communities and environmental factors (K⁺, Ca²⁺, Na⁺, Mg²⁺, ${\text{NH}}_4^{+}$, ${\text{NO}}_2^{-}$, ${\text{NO}}_3^{-}$, Cl⁻, ${\text{SO}}_4^{2-}$, total organic carbon (TOC), total nitrogen (TN), pH and salinity). As shown in Table , most of the environmental factors did not have significant effects on microbial communities in each layer, excluding Mg²⁺ in layer A, Ca²⁺ in layer C, ${\text{NO}}_2^{-}$ in layer F, TOC and TN in layer B, and pH in layers E and I. It was obvious that the microbial communities in all ten layers were adapted to depth gradients and demonstrated stable resistance to environmental changes in oil profiles (Table ). In addition, DPCoA was performed on the microbial communities in each layer (Figure ). The DPCoA results revealed an independent relationship among the ten layers of microbial communities.

Correlations between environmental variables and the microbial community composition by the Mantel test

Environmental variable	A	B	C	D	E	F	G	H	I	J
K+	−0.2526	0.2784	0.0069	0.2629	0.0323	0.1102	0.2133	−0.1142	−0.2386	−0.1356
Ca2+	0.1232	0.1828	−0.3250 *	−0.0835	0.4148	0.3644	0.3378	0.0620	0.2396	0.2697
Na+	0.1031	0.0392	0.1531	−0.0123	0.1002	0.1490	0.2646	0.6197	0.0087	0.0022
Mg2+	0.4541 **	−0.0086	0.0099	−0.0211	0.0160	0.1873	0.0447	−0.0481	0.0403	0.0567
${\text{NH}}_4^{+}$	−0.1450	−0.0292	−0.1802	−0.2909	0.0031	−0.2291	−0.1221	–	−0.2189	0.0495
${\text{NO}}_2^{-}$	−0.0023	−0.2641	−0.2089	−0.1498	−0.3799	0.0259 **	−0.0240	0.0773	0.3007	0.3038
${\text{NO}}_3^{-}$	−0.0503	0.0296	–	0.0295	−0.0436	0.0573	0.4301	0.0077	0.3232	0.1084
Cl−	−0.0296	−0.0215	0.0336	0.0990	0.1781	0.0962	0.2342	0.4263	−0.0937	0.1597
${\text{SO}}_4^{2-}$	−0.2695	−0.1983	0.0338	0.0487	−0.1010	−0.0365	0.0196	0.0316	−0.0016	0.1273
TOC	0.2977	−0.3063 **	−0.1422	0.2033	−0.2898	−0.2273	−0.3077	−0.2347	−0.0410	0.2256
TN	0.2379	−0.3035 *	0.0121	0.2963	−0.2296	−0.2503	–	0.0370	0.1700	0.0865
pH	0.2063	0.3800	−0.1289	−0.0567	0.5459 *	0.4617	0.0454	0.1930	0.4508 *	−0.0622
salinity	−0.3322	−0.0446	0.3782	0.1613	0.1321	0.1311	0.3876	0.1579	−0.0108	0.1887

Bold: P-value of the correlation between physical and chemical factors and microorganisms in each layer.

^*p < .05.

^**p < .01.

Double principle coordinate analysis of microbial communities based on the OTU relative abundances

DISCUSSION

According to the results of the pMEN analysis, the modularity values of each layer in our study were relatively high and close, reflecting habitat heterogeneity, different selection mechanisms, and phylogenetic clustering of closely related species (Lewinsohn et al., ), which might lead to nonrandom interaction patterns and ecological network complexity (Olesen, Bascompte, Dupont, & Jordano, ). This can be explained as follows: (1) modules with closely linked species may be the key units of coevolution, in which reciprocal selection leads to trait convergence of unrelated species (Olesen et al., ); (2) species converge on a combination of related traits shaped by similar interaction patterns (Olesen et al., ); and (3) such a process may result in an interactive heterogeneous network, with taxonomically or functionally related taxa packed into different modules (Deng et al., ). In this study, modularity values were high, which indicated that the ecological network was complex, and the microbial community can fully adapt to the environment. Wang et al. () suggested that the lack of module hubs would lead to fragmentation of the module, and essentially there was no cascading effect on other modules. Furthermore, the modular structure exists to dampen the rapid spread of disturbance in the community (Wang et al., ). Module hubs were detected in each layer in our study, indicating cascade effects between modules and strong resistance of the microbial community to the environment. In addition, positive links can be explained by the mutualism among genes in coevolution (Zhang, Zhao, & Dai, ), or to have similar niches (Zheng et al., ) or cross-feeding (Sun et al., ), while negative relationships are due to competition (Faust & Raes, ). Strong positive links were detected in all networks in this study, which implied that microbial communities were potentially more inclined to function together to adapt to depth-changing environments. In contrast, some microbes that did not possess the ability to compete with other microbes were filtered out (Pointing, Chan, Lacap, & Lau, ). Although modularity changes in our study were not as apparent as the changes reported in the study conducted by Bai et al. () (Figure ), most molecular ecological network parameters varied over a relatively narrow range, which showed that the microbial communities in our study were more stable and adaptable to environmental changes. This was related to many reasons, such as microbial diversity, soil parent material, availability of resources along soil depth gradients, and environmental conditions, which required further research in the future.

The phylogenetic molecular ecological network parameters change along the depth gradient. (a) Results of Bai's study (A: 0 – 0.05 m, B: 0.05 – 0.2 m, C: 0.2 – 0.4 m, D: 0.4 – 0.6 m), (b) Results of this study (A: 0 – 0.2 m, B: 0.2 – 0.4 m, C: 0.4 – 0.6 m, D: 0.6 – 0.8 m, E: 0.8 – 1.0 m, F: 1.0 – 1.3 m, G: 1.3 – 1.6 m, H: 1.6 – 2.0 m, I: 2.0 – 2.5 m, J: 2.5 – 3.0 m)

Since module hubs represented key species in the network (Bai et al., ), we obtained key genera that corresponded to key nodes at different soil layers (Table ). It was evident that multiple module hubs occurred in most layers excluding those of D, H and J with only 1-2 module hubs. A total of 16 module hubs were found from layers A to E, which was consistent with other studies that showed microbial communities were more active within one meter (Blume et al., ; Eilers et al., ; Fierer, Schimel, et al., ; Hartmann et al., ). However, a total number of 16 module hubs were obtained from layers F to J indicating potential activity of subsurface microorganisms. Approximately, 56% of the key genera belonged to the phyla Acidobacteria, Actinobacteria and Proteobacteria. Acidobacteria appeared in layers A, C, F, and G with the dominant species, being Gp25, Gp4, Gp7, and Gp18, respectively. These microbiota are particularly sensitive to pH variations (Pointing et al., ). Although Acidobacteria was discovered from layers A to J, it was not the invariably key phylum in module hubs. The key genus in each layer was potentially related to microorganism diversity, their related network, and environmental conditions, resulting in different module hubs in each layer. The classes Gammaproteobacteria, Deltaproteobacteria, Alphaproteobacteria and Betaproteobacteria belong to the phylum Proteobacteria. Gammaproteobacteria was found in layers E and G and was the aerobic nitrogen-fixing bacteria involved in soil nitrogen cycling (Tsoy, Ravcheev, & Cuklina, ). Deltaproteobacteria, which was also present in layer B, was considered a major choline user. These results provided new insights into our understanding of the use of choline by soil microorganisms (Eleanor, Jason, & Helen, ). Moreover, choline may be involved in the carbon mineralization process (Mussmann, Ishii, & Rabus, ). Betaproteobacteria, which preferred relatively oxygen-rich conditions (Lüdemann, Arth, & Liesack, ), appeared in both layers E and I. Similarly, Sphingomonas was likely to occur under aerobic conditions (White, Sutton, & Ringelberg, ), but was found in layer I and requires further study. The dominant genera found in layers B, C, D, F, and G were members of the phylum Actinobacteria. Most Actinobacteria appeared in the upper soil, which was related to aerobic degree (Eilers et al., ), and its presence may accelerate the decay of animal and plant remains in soil (Stackebrandt, Rainey, & Ward-Rainey, ). Firmicutes and Chloroflexi appeared in the layers G and E, which can be explained by their adaptation to low-nutrient environments (Holanda & Hedrich, ). Nitrosopumilus was found in layer B and was involved in ammonia oxidation in the nitrogen cycle of the ecosystem (Nakagawa & Stahl, ). The phylum Verrucomicrobia appeared in the layers E and G, which was consistent with previous studies reported that they were prone to occur in relatively anoxic environments (Eilers et al., ). Therefore, there should be at least one dominant genus in each layer that must adapt to and resist changes in the soil environment. Factors such as carbon source, microbial decomposition, and species interactions also affect the distribution of microorganisms in the soil depth gradients as well (Marilley & Aragno, ). Notably, module hubs were extremely high at 2 m soil depth. Only a few studies exist on deep microbial communities; thus, more attention should be given to the study of subsurface microbial communities in future.

We performed the Mantel test on key genera obtained by module hubs and environmental factors (K⁺, Ca²⁺, Na⁺, Mg²⁺, ${\text{NH}}_4^{+}$, ${\text{NO}}_2^{-}$, ${\text{NO}}_3^{-}$, Cl⁻, ${\text{SO}}_4^{2-}$, TOC, TN, pH, and salinity) along depth gradients (Table A2). Most environmental factors did not have significant effects on key genera obtained by module hubs along depth gradients, excluding K⁺ and Mg²⁺ in layer A, K⁺ in layer C, ${\text{NO}}_3^{-}$ in layers G and H, Na⁺, Cl⁻, ${\text{SO}}_4^{2-}$, salinity in layer I, and K⁺, Na⁺, and ${\text{SO}}_4^{2-}$ in layer J. It was evident that the key genera obtained by module hubs in all ten layers had no significant difference in pH, TOC, and TN. Soil layers dominant with Acidobacteria were mostly associated with potassium. Ding, Jiang, and Ma () noted that inorganic fertilizers may lower the pH of soil and lead to its acidification. Members of the phylum Proteobacteria had a variety of metabolic types (Song, Liu, & Liang, ). This may affect the sensitivity of Proteobacteria to environmental factors. Bacteria within the phylum Bacteroidetes can be distributed throughout the ecological niches (Garrity & Holt, ), some of which were likely to occur in deep soil layers. During the formation and development of soil, each layer has its own physical and chemical properties due to the migration of energy and matter (Agnelli, Ascher, & Corti, ). In general, depth-related differences in the physicochemical and structural characteristics of soil profiles can include many microenvironmentally complex microbial populations that can evolve (Ranjard & Richaume, ). In contrast, by aiming to determine the correlation between each physical and chemical factor index in each layer for whole microorganisms in our study, it can be concluded that the microorganisms were generally stable. Therefore, some mechanisms affected the adaptation of microbial communities to the environment. Our findings revealed that the adaptability of microbial communities to depth gradients may result in an approximate ten-layer network analysis. Key genera play an important role in maintaining community stability. Notably, salinity has been shown to be a key driving factor for microbial communities in the same research region, as described by Zheng et al. (), but in our study, soil depth seems to also create important effects in the microbial communities.

CONCLUSIONS

In this study, based on the ecological network, we found that microbial communities along the depth gradient had a strong overall resistance. The change in the environmental conditions with soil depth represents an ecological filter; however, the microbial communities fully adapted to the depth gradient, which indicated that there are certain mechanisms that affect the adaptation of microbial communities to the environment.

ACKNOWLEDGMENTS

We gratefully acknowledge Linzhen Guo and Mengfan Yang for soil sampling and Key Laboratory of Environmental and Applied Microbiology, Chinese Academy of Sciences for 16S rRNA Miseq sequencing analysis. This work was financially supported by Ministry of Science and Technology of the People's Republic of China (2018YFD0800400), the National Natural Science Foundation of China (41973017), the Natural Science Foundation of Tianjin (19JCZDJC40700), and the Innovation Team Training Plan of the Tianjin Education Committee (TD12-5037).

CONFLICT OF INTEREST

None declared.

AUTHOR CONTRIBUTIONS

Hang Yu: Conceptualization; Formal analysis; Writing-original draft; Writing-review & editing. Dongmei Xue: Conceptualization; Funding acquisition; Writing-review & editing. Yidong Wang: Funding acquisition; Writing-review & editing. Wei Zheng: Conceptualization. Guilong Zhang: Funding acquisition. Zhongliang Wang: Funding acquisition.

ETHICS STATEMENT

None required.

DATA AVAILABILITY STATEMENT

The datasets used and/or analyzed during the current study are available from the corresponding author upon reasonable request. The original sequence data are deposited at the European Nucleotide Archive: https://www.ebi.ac.uk/ena/data/view/PRJEB21751.

Appendix

Fig. A1. Rarefaction curves in 10 layers

A1 The univariate ANOVA based on Shannon and Chao1 indices

Depth	Shannon	Chao1
A	7.98 ± 0.26	1015.48 ± 96.92
B	7.67 ± 0.78	986.46 ± 199.05
C	7.83 ± 0.79	1034.79 ± 91.86
D	7.53 ± 0.83	1026.91 ± 109.15
E	7.63 ± 0.79	1025.48 ± 171.69
F	7.60 ± 0.60	1023.90 ± 108.78
G	7.52 ± 0.89	984.15 ± 123.27
H	7.82 ± 0.93	1040.34 ± 181.42
I	7.58 ± 0.93	1035.21 ± 171.58
J	7.81 ± 0.39	1013.16 ± 78.05

The values in the table were the mean ± standard deviation (SD) of eight replicates. No significant difference in depth gradient (p > .05).

A2 Correlations between environmental variables and key genus obtained by the module hubs along depth gradient by the Mantel test

Environmental Variable	A OTU832 OTU2510 OTU4954	B OTU71 OTU498 OTU774 OTU1532	C OTU650 OTU1466 OTU1457	D OTU723 OTU967	E OTU631 OTU1510 OTU1853 OTU4767	F OTU396 OTU1362 OTU1366 OTU2089	G OTU316 OTU1664 OTU2980 OTU1853 OTU2885	H OTU679	I OTU15 OTU156 OTU570 OTU408 OTU2365	J OTU2930
K+	0.4474 *	−0.1732	−0.2505 *	−0.2347	0.0621	−0.0162	0.1845	−0.1926	−0.2084	−0.2024 *
Ca2+	−0.1757	−0.1326	−0.2785	0.1966	0.3680	0.0851	0.4014	−0.0175	0.0853	0.1758
Na+	0.0859	0.1264	0.3348	−0.1997	−0.0563	0.1755	−0.0877	−0.0631	−0.5248 *	−0.1960 *
Mg2+	−0.3225 *	−0.1878	−0.0016	0.2329	0.1755	0.1049	−0.0454	0.0425	0.2017	−0.2367
${\text{NH}}_4^{+}$	−0.1156	−0.0810	−0.0673	−0.2864	−0.2065	0.0586	0.2035	−	0.0185	0.2037
${\text{NO}}_2^{-}$	−0.1066	0.2282	−0.0259	0.0508	0.0905	0.1427	−	0.0182	−0.0544	−0.1956
${\text{NO}}_3^{-}$	−0.0288	−0.1500	−	0.3382	0.1822	−0.0630	0.3516 *	0.3547 *	−0.0001	−0.1995
Cl−	0.0520	0.0577	0.3111	−0.1237	0.1017	0.0973	−0.1182	−0.1553	−0.5618 *	−0.1395
${\text{SO}}_4^{2-}$	−0.0939	0.1519	0.1277	0.1251	0.2083	−0.0863	−0.0776	−0.0644	−0.3232 *	−0.2408 *
TOC	0.0596	0.2642	0.0393	0.1997	0.1807	−0.0481	−0.4587	−0.2181	−0.0410	−0.1434
TN	0.07764	0.2453	0.1979	0.1101	0.0962	0.0132	−	−0.149500	−0.0903	−0.0952
pH	0.1327	−0.0847	−0.0094	−0.2744	0.1294	0.0678	−0.1846	−0.0379	−0.0903	−0.0141
Salinity	−0.0306	−0.0493	0.1228	−0.0614	0.2043	−0.0623	0.0147	0.1850	−0.5251 *	−0.1903

^*p < .05

^**p < .01.