There were two soils from the Jiansanjiang reclamation area (JSJ) and the Langfang research base of the Chinese Academy of Agricultural Sciences (LF). According to soil particle diameter, the soil form JSJ was identified as silty clay, and the soil form LF was identified as silty loam. The silty loam had 18 g organic matter kg⁻¹, 74.9 mg available P kg⁻¹, 289.8 mg available K kg⁻¹, and a pH of 7.07; the silty clay had 25.8 g organic matter kg⁻¹, 51.7 mg available P kg⁻¹, 289.8 mg available K kg⁻¹, and a pH of 7.24. The soils were sieved with 2-mm mesh and preincubated for 2 weeks (Trabue et al., 2006). The concentration transfer of this study was based on soil depth of 10 cm with a bulk density of 1.5 g cm⁻³ (GB/T31270.1-2014, 2014). The purity of clomazone is 98.4% and purchased from Beijing Qinchengyixin Technology Development Co., Ltd. (Beijing, China). Three clomazone treatments were prepared in brown bottles: 0.8 mg kg⁻¹ (active ingredients per soil dry weight; L), 8 mg kg⁻¹ (a.i./dw; M), and 80 mg kg⁻¹ (a.i./dw; H). The L level represents recommended application rate in the field; the M level represents excessive use of clomazone in the field, and the H level represents extremely polluted soil (e.g., soil near a pesticide factory). In addition, a control treatment was also needed. These treatments were all prepared in triplicate. Soil moisture was adjusted daily by deionized water to 50% of the maximum water-holding capacity. The soil samples were kept for 90 days in an artificial climate box, and the temperature was maintained at 25°C. The samples were taken on days 7, 15, 30, 60, and 90 and kept in a refrigerator at -80°C.

Soil microbial DNA was extracted using a PowerSoil Isolation Kit (Mo Bio Laboratories, Carlsbad, CA, USA), according to the instructions, and the DNA quality was evaluated using an ND-1000 spectrophotometer (NanoDrop Technologies). 16S rRNA gene was amplified by the primer sets of 341F (5′-CCTAYGGGRBGCASCAG-3′) and 806R (5′-GGACTACNNGGGTATCTAAT-3′) (Yu et al., 2005). Microbial DNA was amplified in 50 μl reactions per sample, and each PCR solution contained 100–300 ng of DNA template, 1.5 μl of each 10 μM primer, 5 μl of 2 mM dNTPs, 1 μl of KOD-Plus-Neo enzyme (Toyobo, Shanghai, China), 5 μl of 10 × PCR Buffer for KOD-Plus-Neo, 3 μl of 25 mM MgSO₄, and water to 50 μl. Reaction procedures were as follows: an initial step was 94°C and kept for 2 min, followed by 35 cycles of 98°C for 10 s, 62°C for 30 s, and 68°C for 30 s, and the final extension temperature was 68°C for 10 min. Negative-control reactions were also needed. PCR products were analyzed by 1.5% agarose gel electrophoresis and purified with a PCR Purification Kit from QIAGEN (Hilden, Germany). Purified PCR products were sequenced using Illumina equipment (Santiago, CA, USA). Amplicon sequencing data were processed through the USEARCH pipeline (Edgar, 2010, 2013), and clean data were clustered into operational taxonomic units (OTUs) with 97% similarity.

All networks were established on the basis of Pearson's correlations and performed on Cytoscape (Faust and Raes, 2016). The correlation coefficient was set as 0.9. The network topological indices were also based on Cytoscape. For each network, the topological roles of each node were classified and calculated by nodes' within-module connectivity (Zi) and among-module connectivity (Pi) (Guimerà and Nunes Amaral, 2005). An adopted criterion in previous studies (Olesen et al., 2007; Zhou et al., 2011; Shi et al., 2016; Yuan et al., 2021) was used in this study to identify module hubs (Zi ≥ 2.5, Pi < 0.62), connectors (Zi < 2.5, Pi ≥ 0.62), and network hubs (Zi ≥ 2.5, Pi ≥ 0.62). Module hubs referred to nodes that were highly connected to other members in a module, connectors referred to nodes that linked different modules, and network hubs were nodes that were both a module hub and a connector. These nodes were referred to as keystone nodes. Other nodes were classified as peripherals (Banerjee et al., 2019; Röttjers and Faust, 2019).

In a network, the node composition in one network module differs from that in other modules. However, some modules may preserve with some shared nodes after environmental change (Deng et al., 2012). Fisher's exact test is used to evaluate an association of two categorical variables (Warner, 2013) and has been used to identify the preserved modules in previous studies (Horvath, 2011; Langfelder et al., 2011; Deng et al., 2012; Dong et al., 2021). This method has also been used to evaluate preserved modules between control and clomazone treatments in two soils. There were four categories for the nodes within the two networks during the evaluation of preserved modules. In the first case, members were included in two modules; in the second case, members were included in one module of the pair; in the third case, members were included in the other module of the pair; in the fourth case, members were not included in these two modules. To determine whether nodes in the two modules were independent or exclusive, the observed frequency of the four categories was placed into four cells of a contingency table for one-sided exact testing. Each p-value from the exact tests was adjusted through the Bonferroni procedure within each network. Through Fisher's exact test, two modules from different networks that both included a significant part of the same nodes were considered as preserved modules.

Network stability could evaluate ecological system stability to disturbance (Thébault and Fontaine, 2010), and it was usually evaluated by network robustness and vulnerability (Wu et al., 2021; Yuan et al., 2021). Robustness and vulnerability were used to evaluate network stability (Wu et al., 2021; Yuan et al., 2021). Robustness is defined as the remaining proportion of species after random removal in the network (Montesinos-Navarro et al., 2017). In this study, every 0.05% of nodes was randomly removed to simulate random species removal. Vulnerability is calculated as V = max [(E – Ei)/E], in which E is the global efficiency and Ei is the global efficiency after removing node i and its entire links (Deng et al., 2012). Global efficiency is calculated as E = ∑_j≠i[1/d(i,j)]/n(n – 1), in which d(i,j) is the number of edges in the shortest path of node i to j (Deng et al., 2012).

The metabolic function of each sample was predicted by tax4fun based on 16S rRNA gene data (Aßhauer et al., 2015). There was one cellular process, three genetic information processing, and 11 metabolism categories that have been used to analyze the correlation with the bacterial community using the Mantel test (Duan et al., 2020). The cellular process was cell growth and death (CGD); the genetic information processing categories were folding, sorting, and degradation (FSD), replication and repair (RR), and translation; and the metabolism categories were carbohydrate metabolism (CM), lipid metabolism (LM), amino acid metabolism (AAM), metabolism of cofactors and vitamins (MCV), xenobiotic biodegradation and metabolism (XBM), biosynthesis of other secondary metabolites (BOSM), energy metabolism (EM), metabolism of terpenoids and polyketides (MTP), metabolism of other amino acids (MOAA), nucleotide metabolism (NM), and glycan biosynthesis and metabolism (GBM). Some functions were fundamental for ecological balance, for example, XBM is important for chemical pollution cleaning (Thelusmond et al., 2019).

Analysis of similarity (ANOSIM) has been used to evaluate differences in bacterial community structure based on Bray–Curtis distance (Oksanen et al., 2012). Network dissimilarity is an effective tool to evaluate networks' differences, and it is based on network nodes and edges (Poisot et al., 2012; Mo et al., 2021). Shared nodes and edges of two networks are used to evaluate coexisting elements of different networks. Correlation coefficient “r” has been used to evaluate the correlation of functions to the bacterial community; Fisher's least significant difference test has been used to evaluate significant differences, and a 5% level was set (p < 0.05).

There were eight networks that have been established in Figure 1A. In all networks, the node comprised mostly of eight phyla: Acidobacteria, Actinobacteria, Bacteroidetes, Chloroflexi, Gemmatimonadetes, Planctomycetes, Proteobacteria, and Verrucomicrobia. In the JSJ soil, the node percentages in clomazone treatments were increased in Acidobacteria, Proteobacteria (except for H treatment), and Verrucomicrobia (except for L treatment) but decreased in Actinobacteria, Gemmatimonadetes, Bacteroidetes (except for L treatment), Chloroflexi (except for M treatment), and Planctomycetes (except for M treatment) (Table 1). In the LF soil, Acidobacteria (except for H treatment), Chloroflexi (except for L treatment), Planctomycetes, and Verrucomicrobia were increased, while others were decreased (except for Actinobacteria in H treatment and Proteobacteria in M treatment) (Table 1).

The network's topological indices are shown in Table 2. Compared with the network of the JSJ control soil, network size (total nodes), the number of links, and the degree of clomazone treatments were decreased by 9–127, 648–829, and 2.249–2.429, respectively. While in the LF soil, network size (total nodes), the number of links, and the degree of clomazone treatments were decreased by 228–384, 667–753, and 0.703–1.371, respectively.

The influenced network could result in changing the roles of the networked members. Based on the criteria of node classification, keystone nodes were 473, 360, 466, and 392 for control, L, M, and H in the JSJ soil, respectively; network keystone nodes were 549, 360, 466, and 466 for control, L, M, and H in the LF soil, respectively (Figure 2). The shared keystone nodes were 173, 231, and 165 for the comparison of control and L, control and M, and control and H in the JSJ soil. The shared keystone nodes were 86, 119, and 119 for the comparison of control and L, control and M, and control and H in the LF soil (Figure 3).

ANOSIM is for the comparison of control and clomazone treatments. They were 0.617–0.879 in the JSJ soil and 0.985–0.997 in the LF soil (Table 3). The network dissimilarities between control and clomazone treatments were 0.954–0.983 in the JSJ soil and 0.979–0.990 in the LF soil (Table 3). These results revealed that the composition of each network was significantly impacted by clomazone in both soils. In addition, the shared nodes between control and clomazone treatments were used to evaluate the effects of clomazone in the JSJ soil and LF soil. Compared with the nodes in the JSJ soil control, the shared nodes were 196, 260, and 194 for L, M, and H treatments, respectively (Table 3). For the LF soil, there were 131, 151, and 171 shared nodes for L, M, and H treatments, compared with the nodes in the control treatments (Table 3). Compared with the links in the JSJ soil control, the shared links were 71, 100, and 36 for L, M, and H treatments, respectively (Table 3). For the LF soil, there were 17, 25, and 12 shared nodes for L, M, and H treatments, compared with the nodes in the control treatments (Table 3).

Based on random species loss, the robustness was decreased by 0.014–0.016 in clomazone treatments in the JSJ soil; for the LF soil, it was decreased by 0.01–0.013 in clomazone treatments (Figure 4). For vulnerability in the JSJ soil, it was increased by 0.00079–0.00147 in clomazone treatments; in the LF soil, it was decreased by 0.00023–0.00044 for clomazone treatments.

The influenced networks suggested that clomazone could alter network organization. The big modules (i.e., ≥5 nodes) were used to analyze preserved modules based on Fisher's exact test. In total, there were 25 preserved module pairs in the two soils (Table 4), and most of the preserved module individuals belong to the phyla of Proteobacteria, Acidobacteria, and Planctomycetes (Figures 1A, B). Specifically, in the JSJ soil, there were five module pairs accounted for 0.39% of total module pairs between CK and L, seven module pairs accounted for 0.55% of total module pairs between CK and M, and six module pairs accounted for 0.58% of total module pairs between CK and H (Table 4). In the LF soil, there were two module pairs accounted for 0.13% of total module pairs between CK and L, three module pairs accounted for 0.18% of total module pairs between CK and M, and two module pairs accounted for 0.09% of total module pairs between CK and H (Table 4).

An intriguing issue is whether the alterations in bacterial network composition as a result of clomazone treatment caused alterations in microbial community functions and their associated ecosystem processes. We used the Mantel test to address this issue, and the relationships are shown in Figure 5. In the JSJ soil, bacterial network community correlated with MOAA and MTP in the control treatment (r ≥ 0.4); in clomazone treatments, the correlations of MOAA (only in L treatment), BOSM, RR, MTP, and CM (only in H treatment) with the bacterial community are higher than 0.4. In the LF soil, bacterial network community correlated with BOSM, CM, AAM, and LM in the control treatment (r ≥ 0.4); in L treatment, it correlated with EM, RR, and MOAA; in M treatment, it correlated with EM; in H treatment, it correlated with translation, RR, MCV, MOAA, and FSD (r ≥ 0.4).