A total of two soil samples were collected from the Langfang research base, Hebei Province (LF), and the Jiansanjiang reclamation area, Heilongjiang Province (JSJ). Based on soil particle diameter, soils from LF and JSJ were classified as silty loam soil and silty clay soil, respectively. In the LF soil, the content of organic matter, available P, available K, and pH were 25.8 g kg⁻¹, 51.7 mg kg⁻¹, 289 mg kg⁻¹, and 7.24 mg kg⁻¹, respectively. In JSJ soil, the content of organic matter, available P, available K, and pH were 18.0 g kg⁻¹, 74.9 mg kg⁻¹, 289.8 mg kg⁻¹, and 7.07 mg kg⁻¹, respectively. A 2-mm mesh was used to sieve the soils, and the soils were preincubated for 2 weeks (Trabue et al., 2006). Trifluralin purity was 98%, and it was dissolved in acetone (analytical grade; Beijing Chemical Company). In total, three concentrations of trifluralin in soils were used. The active ingredients of trifluralin in 1 kg of dry soil were 0.84, 8.4, and 84 mg, and they were corresponding to 1 (L), 10 (M), and 100 (H) times of the recommended application rate, respectively. The M level represents excessive use of pesticides in the field, and the H level represents extremely contaminated soil by pesticides (e.g., soil near a pesticide factory). The procedure of pesticide exposure is as follows: 50 g of soil and 100 μL pesticide solution were added to the dark brown bottles, and then thoroughly mixed for 15 min; after that, 200 g of soil was transferred to each bottle and thoroughly mixed for 15 min. A control treatment was also needed and this consisted of treatment with a solution lacking trifluralin. Each treatment was repeated three times. The soil concentration was 1.5 g cm⁻³ (GB/T31270.1-2014, 2014). Soil moister was kept at 50% using deionized water. The weighing method was used to determine the loss of water every 2 days, and the soil moisture was kept constant according to the loss of weight. The laboratory experiment was carried out for 3 months in an artificial climate box at 25°C. Sampling times were 7, 15, 30, 60, and 90 days after experimental establishment. The samples were stored at −80°C until analysis.

A PowerSoil Isolation Kit (Mo Bio Laboratories, Carlsbad, CA, USA) was used to extract soil microbial DNA. An ND-1000 spectrophotometer (NanoDrop Technologies) was used to analyze microbial DNA quality. The forward primer ITS3_KYO2 (5′-GATGAAGAACGYAGYRAA-3′) and reverse primer ITS4 (5′-TCCTCCGCTTATTGATATGC-3′) were used to amplify internally transcribed spacer (ITS; Tian et al., 2017). Microbial DNA was amplified using a PCR amplifier, and each amplification system contained 1.5 μL of each 10 μM primer, 100–300 ng DNA template, 5 μL of 2 mM dNTPs, 1 μL 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. The temperature change steps were as follows: 94°C for 2 min, followed by 98°C for 10 s (for 35 cycles), 62°C for 30 s, and 68°C for 30 s, and the final extension temperature was 68°C for 10 min. Negative control with DNA solution was also settled. A PCR Purification Kit (Qiagen, Hilden, Germany) was used to purify PCR products after 1.5% agarose gel electrophoresis. An Illumina platform (Santiago, CA, United States) was used to sequence the purified PCR products using a 2 × 250 bp kit. USEARCH was used to process amplicon sequencing data (Edgar, 2010, 2013). The following rules were used to filter raw reads: (1) adaptors were cut, (2) reads which included more than 10% of unknown nucleotides were removed, and (3) reads that included < 80% of bases with quality (Q-value) > 20 were removed. Tags were assembled with clean reads according to more than 10 bp overlaps and <2% mismatch between paired-end reads. Clean data were clustered into operational taxonomic units (OTUs) with a similarity of 97%.

All co-occurrence networks were established on the basis of Pearson correlations of OTU abundances and performed on the Cytoscape platform using the CoNet plugin (Faust and Raes, 2016). Pearson's correlation was used to analyze the association of pairwise fungal OTUs with an absolute value of correlation coefficient (r) higher than 0.7. The topological indices were calculated using Gephi software include total nodes, total edges, average degrees, clustering coefficient, network density, and path length. Nodes in the network represent the OTUs in the network. The degrees of each node represents the connections of a node to others, and the average degree represents the main value of all degrees. Modularity based on the connections of nodes represents the level of a network divided into different modules. A network diagram was established on Gephi software.

The node's topological role was evaluated by among-module connectivity (Pi) and within-module connectivity (Zi; Guimerà and Nunes Amaral, 2005). The network nodes were classified as module hubs with Zi ≥ 2.5 and Pi < 0.62, connectors with Zi < 2.5 and Pi ≥ 0.62, and network hubs with Zi ≥ 2.5 and Pi ≥ 0.62 (Olesen et al., 2007; Chen et al., 2019). These three categories of nodes are referred to as keystone nodes (Banerjee et al., 2019; Röttjers and Faust, 2019). Upset plots were used to visualize shared keystone nodes between control and trifluralin treatment (Lex et al., 2014).

Network stability can be used to evaluate ecological system stability to disturbance (Thebault and Fontaine, 2010). Generally, it is usually evaluated by network robustness and vulnerability (Wu et al., 2021; Yuan et al., 2021). Robustness is the remaining proportion of nodes in the network after removing some nodes (Montesinos-Navarro et al., 2017). In this study, 0.05% of nodes in the network have been removed to simulate random species removal each time. Vulnerability is also an index used to evaluate network stability based on node removal (Yuan et al., 2021).

Network dissimilarity was used to evaluate the dissimilarity of two fungal networks in this study (Poisot et al., 2012; Mo et al., 2021). It was evaluated by shared nodes and edges of two different networks, and the shared nodes and edges were also used to evaluate dissimilarity (Poisot et al., 2012; Mo et al., 2021).

FUNGuild is a database of fungal functions, and it clusters almost all published studies on fungal functions (Nguyen et al., 2016). Based on fungal amplicon sequencing data and taxonomy, FUNGuild can be used to predict fungal functions. There are three categories of trophic modes of fungi, namely, pathotroph, saprotroph, and symbiotroph. Furthermore, these three categories can be divided for better evaluation of fungal functions. Saprotroph was divided into dung saprotroph, leaf saprotroph, plant saprotroph, soil saprotroph, and wood saprotroph; pathotroph was divided into animal pathogen, plant pathogen, fungal parasite, lichen parasite, bryophyte parasite, and endophyte; and symbiotroph was divided into ectomycorrhizal, ericoid mycorrhizal, and endophyte. Fungal network OTUs abundance was used to analyze the correlation of fungal community structure with fungal functions based on the mantal test (Duan et al., 2020).

A total of eight fungal networks were established for each treatment based on Pearson's correlation coefficients of fungal OTUs (Figure 1, Table 1) presents each network's topological indexes. In the networks, the nodes were assigned to four fungal phyla in LF soil and three fungal phyla in JSJ soil. Among these, the phyla Ascomycota and Basidiomycota were most abundant in both soils. Compared with the control, the total nodes, total links, and average degree were all increased by trifluralin. In LF soils, the total nodes, total links, and average degrees were increased by 11–45, 252–392, and 1.151–1.468 in trifluralin treatments, respectively; in JSJ soil, the total nodes, total links, and average degrees were increased by 6–29, 134–234, and 0.169–1.227 in trifluralin treatments, respectively. The average path length was decreased by 0.558, 0.304, and 0.70 in L, M, and H treatments in LF soils and by 0.562, 0.373, and 0.415 in L, M, and H treatments in JSJ soils.

Based on nodes' Zi and Pi, there were 93, 111, 86, and 96 keystone nodes in control, L, M, and H treatments in LF soil, respectively; there were 134, 118, 105, and 117 keystone nodes in control, L, M, and H treatments in JSJ soil, respectively. In LF soil, the shared keystone nodes were 35, 23, and 30 for the comparison of control-L, control-M, and control-H, respectively; in JSJ soil, they were 45, 37, and 36 for the comparison of control-L, control-M, and control-H, respectively (Figure 2).

Network dissimilarity is an effective tool to evaluate network similarity. According to Table 2, for the comparison of control-L, control-M, and control-H, the shared nodes were 231, 227, and 219 in LF soil and 285, 272, and 279 in JSJ soil, separately. The shared nodes accounted for a significant part of the total nodes in each network, but the shared links were significantly low for each comparison. Specifically, for the comparison of control-L, control-M, and control-H, shared links were 20, 16, and 27 in LF soil and 19, 26, and 23 in JSJ soil. According to the shared nodes and links, the network composition was significantly influenced by trifluralin. The dissimilarity was 0.97–0.99 between control and trifluralin treatments in both soils.

On the basis of random species loss, the network robustness was increased by trifluralin in the two soils (Figure 3). In LF soil, it was increased by 0.009, 0.003, and 0.002 in L, M, and H, respectively. In JSJ soil, it was increased by 0.003, 0.006, and 0.007 in L, M, and H, separately. For vulnerability, it was decreased by trifluralin in both soils (Figure 3). In LF soil, it was decreased by 0.00032, 0.00022, and 0.00014 in L, M, and H, separately. In JSJ soil, it was decreased by 0.0001, 0.00015, and 0.00019 in L, M, and H, separately.

There were 11 guilds identified in this study, and they were dung saprotroph, endophyte, lichen parasite, fungal parasite, plant pathogen, animal pathogen, plant saprotroph, soil saprotroph, wood saprotroph, ectomycorrhizal, and mycorrhizal. The correlations of the fungal network community with functions are shown in Figure 4. In LF soil, the fungal community was significantly correlated with dung saprotroph in the control treatment; in H treatments, the fungal community was significantly correlated with dung saprotroph and wood saprotroph. In JSJ soil, the fungal community was significantly correlated with dung saprotroph in L treatment; in M treatments, the fungal community was significantly correlated with endophyte, lichen parasite, plant pathogen, animal pathogen, soil saprotroph, and wood saprotroph; in H treatment, the fungal community was significantly correlated with dung saprotroph and endophyte.