Increased OCP levels decreased soil bacterial diversity and enriched functional gene abundance.
Soil samples were divided into three groups based on the OCP levels (i.e., clean group with OCP below the limit of detection, low pollution group with OCP at 14.53-79.69 mg kg-1 and high pollution group with OCP at 193.58-339.40mg kg-1) (Figs. S1-S2). As expected, different soil groups exhibited significant variations in bacterial compositions and diversity (Figs. 1A-B, Figs. S3-S8) [36]. Specifically, with increasing OCP levels, the Simpson index of iDNA decreased from 0.67 for clean soil to 0.56 for the highly polluted soil, and the Simpson index of eDNA decreased from 0.79 for clean soil to 0.72 for the high pollution soil. The Shannon index exhibited the similar trends as the Simpson index, indicating that species diversity of e-iDNA declines significantly with increasing OCP levels (p<0.01) (Fig. S3). Notably, the extent of diversity reduction was greater in the iDNA (Decreasing slope of bacterial species abundance>5%) relative to the eDNA (Decreasing slope of bacterial species abundance<5%), suggesting the bacterial species information in iDNA was more sensitive to OCP-induced environmental stress(Fig. 1A, Figs. S4-S6) [39-40].
Functional genes associated with microbial metabolism and environmental resilience were analyzed to gain insight into microbial response to increased OCP-induced stress(Figs. 1C-D, Figs. S9-S14) [36, 39]. Based on annotation using the KEGG database, the e-iDNA functional genes were classified to five groups (i.e., G1: OCPs degradation, G2: C/N/P/S metabolism, G3: Fe transport, G4: DNA repairand recombination, and G5: biofilm formation)(Figs. S9-S11) [2]. As the OCP levels increased, the intracellular genes types related to the five functional groups increased in e-iDNA annotations (Fig. 2, Figs. S9-S11). G1genes increased from 161 types (clean soil) to 174 types (high-OCP soil), G2 genes increased from 362 types (clean soil) to 402 types (high-OCP soil), G4 in iDNA increased from 21 types in clean group to 29 types in high pollution groups, while G3 genes (20~21 types) and G5 genes (92~95 types) showed a slight increase as OCP levels increased (Fig. 2, Fig. S11). In particular, the proportion of the relative abundance of G1, G2, G3, G4 and G5 functional genes increased significantly in iDNA annotations (p<0.05)(Figs. S9-S11). As OCP levels increased, G1~ G5 functional genes in high pollution group showed an increase by 5.09%~16.6% compared with clean groups (Increasing slope of functional genes abundance>5%). Consistently, as the OCP levels increased, the extracellular annotated genes related to the five functional groups showed a similar increasing trend with intracellular annotations. Specifically, in eDNA annotations, G1 genes increased from 80types (clean soil) to 90 types (high-OCP soil). Moreover, G2 genes increased from 171 types (clean soil) to 204 types (high-OCP soil), while G3 (11~12 types), G4 (7~11 types) and G5 (53~59 types) functional genes showed a slight increase as OCP levels increased (Fig. 2, Fig. S11). Correspondingly, as OCP levels increased, the proportion of the relative abundance of G1 genes in low and high pollution soil groups increased by 4.93% and 5.47% compared with the clean group, respectively. Functional genes related to G2~G5 in high pollution group increased by 0.01%~4.19% compared with clean groups (Increasing slope of functional genes abundance<5%)(Fig. S9-S11), suggesting the functional genes information in iDNA was more sensitive to OCP-induced environmental stress.
Increased OCP levels enhancedintracellular and extracellular gene associations.
With increasing OCP levels, the overlapping genes types in e-iDNA annotations increased from 41.70% in clean soil to 45.59% in highOCP pollution soils (p<0.05, and Increasing slope of functional genes abundance>5%). Among these overlapping genes, the proportion of the relative abundance of OCPs degradation genes increased from 3.02% to 4.32% (p<0.05, Increasing slope of functional genes abundance>5%).The proportion of the relative abundance of the remaining functional genes groups increased from 3.28% to 4.85% (p>0.05, Increasing slope of functional genes abundance>5%).In eDNA, the proportion of the relative abundance of these overlapping genes exhibited no significant change (p>0.05, Increasing slope of functional genes abundance<5%)(Fig. 1C, Figs. S12-S14).It indicated that the overallintracellular and extracellular gene information transfer was maybe gradually enhanced with the increasing OCP-induced environmental stress (Fig. 1C).At the meantime, these overlapping genes in e-iDNA annotations were divided into four types according to their relative abundance variation to explain the gene information transfer in e-iDNA (Type Ⅰ, Type Ⅱ, Type Ⅲ and Type Ⅳ, Method and Fig. 1D).These four types of genes accounted for 19.07-19.12% (eDNA) and 17.00-17.25% (iDNA) of all annotated genes, respectively. As increasing OCP levels, the proportion orders of the four types of intracellular genes exhibited a significant variation (e.g., Type Ⅲ>Type Ⅱ>Type Ⅳ>Type Ⅰ (clean group), Type Ⅱ>Type Ⅲ>Type Ⅰ>Type Ⅳ (low-OCP pollution groups), and Type Ⅱ>Type Ⅰ>Type Ⅲ>Type Ⅳ (high-OCP pollution groups)(Fig. 1D, Figs. S12-S15). The proportion orders of the extracellular genes remained stable as OCP levels increased (e.g.,Type Ⅱ>Type Ⅳ>Type Ⅰ>Type Ⅲ)(Fig. 1D, Figs. S12-S15), indicating that the intracellular overlapping genes were more sensitive to OCP pollution than those in extracellular.As expected, the proportion of Type Ⅱ genes was the highest in e-iDNA among the three soil groups, suggesting that intracellular and extracellular genes exhibitedthe opposite variations with increasing OCPs pollution. Moreover,as OCP levels rise, the proportion orders of the G1 functional genes in four types remained the same in e-iDNA (Type Ⅱ>Type Ⅰ>Type Ⅲ>Type Ⅳ), except for the clean groups. Similarly, the proportion orders of the G2, G3, G4 and G5 genes remained constant as follows, Type Ⅱ>Type Ⅰ>Type Ⅲ>Type Ⅳ (Fig. 2, Fig. S14). It further demonstrated that the relative abundance of G1~G5 functional genes exhibited an opposite variation in both intracellular and extracellular fractions, proving the strongbuffer capacityof eDNA fraction and information transfer within e-iDNA annotations.
Furthermore, co-occurrence networks were constructed to illustrate gene interactions and evaluate the potential variation ranges of microorganism structures and functions among the OCP-induced soil groups (Fig. 2, Fig. S12) [40-41].In the co-occurrence networks formed by iDNA annotated genes, the relative abundance of G1~G3 functional genes increased significantly with rising OCP-levels (p<0.01). While in eDNA fractions, the proportions of the relative abundance of G2 gene types in the networks increased significantly as OCP concentrations increased (p<0.01, Figs. 2A-D). Within the co-occurrence network of clean soil group, the interaction types of the copresence and mutual exclusion accounted for 52.40% and 47.60% (eDNA) and 39.60% and 60.40% (iDNA), respectively. In the other twosoil groups, there was 100% copresence interaction both in extracellular and intracellular gene networks, demonstrating that the OCP-induced environmental stress could significantly decrease mutual exclusion correlations among genes in co-occurrence networks (Fig. 2 and Figs. S12-S15).
OCP-induced stress enriched keystone genes associated with metabolism and gene transfer.
According to the connection degree of individual node within the co-occurrence networks, keystone genes (KGs) were calculated and detected as module hubs(Figs. S16-S17) [2-3]. In total, 37, 78 and 29 types of keystone genes were obtained in the three soil groups, respectively (e.g.,KGs-Clean, KGs-Low and KGs-High, Fig. 3, Fig. 4, Figs. S16-S19).Interestingly, there were no overlapping KGs among different soil groups, indicating their types variations with increasing OCP-levels(Fig. 3A, Fig. S16). Intheclean soil group, the top five KGs were corresponding to “ABC transporters (K02030)”, “Metabolic pathways (K01209)”, “Oxidoreductases (K11209)”, “curved DNA-binding (K00516)” and “Sulfur metabolism (K04091)”. In the low-OCP pollution groups, genes associated with “Transferases”, “DNA repair and recombination”, “Metabolic pathways”, “Secretion system” were identified as KGs (e.g.,K12132, K01971, K05349, K06131, K02652). Meanwhile, gene aggregate (e.g.,K01996, K03654, K03296,K00342, K00185) corresponding to “ABC transporters”, “DNA repair and recombination”, “Oxidative phosphorylation”, and “Sulfur metabolism” were observed as the top five KGs in high-OCP pollution soils. As OCP-induced stress increased, the orders of KGs associated with DNA repair, ABC transporters, and Metabolic pathways were advanced, proving that the microbial metabolic and gene transfer potentials were promoted.
Furthermore, functional KGs (G1~G5) were analysed to gain their metabolic and gene transfer variations with increasing OCP-exposing environments (Fig. 3A, Fig. 4 & Fig. S19). Respectively,the relative abundance ratios of functional KGs increased by 2.14 times in intracellular fractions (p<0.05), but decreased by 3.46%in extracellular annotations (Fig. S19) as OCP-levels increased. Correspondingly, in the low- and high-OCP pollution groups, the total relative abundance ratios of functional KGs related to G1, G2 and G3 types in iDNA annotations were higher than those in clean soil group and increased significantly with increasing OCP pollution (Fig. 4A, p<0.01). While in eDNAfractions, the total proportions of G1~ G3 functional KGs were higher in clean and low-OCP pollution soil groups than those in high-OCP groups and decreased significantly with increasing OCP levels (Figs. 4A-4B& Fig. S19, p<0.01). It indicated that the functional potentials associated with OCPs degradation and metabolic pathways were promoted in intracellular microbiome and restrained in extracellular fractions with increasing OCP levels [36, 39].
Keystone genes had niche specificity devoting for microbiome interactions in OCP-induced soils.
Despite of their lowtotal contents in the soil environment (ranging from 0.90% to 5.94%), the functions and stability of these KGs were conductive to exploring the structure maintenance of ecosystems (Fig. S18).Among these KGs, the four specific gene types (e.g., Type Ⅰ~Type Ⅳ) in KGs ranged from 9.59% to 14.57%.Specifically,Type Ⅱ (iDNA) and Type Ⅳ (eDNA)genes were observed with highest proportionsin the clean soil group(Figs. S20). In low-OCP pollution soils, the highest proportions of KGs were evaluated as Type Ⅱ (iDNA) and Type Ⅳ (eDNA). While in high-OCP pollution soils, Type Ⅰand Type ⅢKGswere the most abundantgene types in intracellular and extracellular annotations, respectively. As we expected, KGs also exhibited an opposite variation in both intracellular and extracellular fractions as the total gene annotations (Fig. 2, Fig. S15).In addition, the total relative abundance of KGs obtained from the three soil groups (KGs-Clean, KGs-Low and KGs-High) ranged from 0.90-0.99%, 5.32-5.94%, 1.94-2.05% in iDNA fractions and 3.25-3.31%, 4.40-4.43%, 2.59-2.61% in eDNA fractions, showing no significant difference with OCP-levels increasing (ANOVA analysis, p>0.05).
Subsequently, the relative abundance ratios of KGs-Clean, KGs-Low and KGs-High in the co-occurrence networks were analysed to evaluate their interactions frequency (Figs. 3D-3E). In intracellular annotations, the variation in the relative abundances of KGs-Clean, KGs-Low and KGs-High ranged from 0.37-0.99%, 1.48-5.73% and 0.80-1.83%, respectively. Similarly, in eDNA fractions, the proportions of the relative abundance of KGs-Clean, KGs-Low and KGs-High ranged from 1.02-3.25%, 2.15-3.94% and 1.57-2.50%, showing a significant variation among different soil groups of co-occurrence networks (p<0.01). Meanwhile, both in e-iDNA fractions, KGs-Clean, KGs-Low and KGs-High only accounted for the highest relative abundance ratios in clean, low- and high-OCP pollution soil groups, respectively and decreased dramatically in the other two soil types, supporting that KGs had niche specificity in microbiome interaction networks in OCP-induced soils(Figs. 3D-3E) [15, 42-43].
eDNA harboured keystone genes to facilitate intercellular exchanges with and after OCP stress
Simplified cultivation medium was applied to verify the niche specificity of KGs and intercellular gene exchanges with and after OCP-induced stress (e.g.,p-nitrochlorobenzene (p-NCB): 0,10,100 mg L-1, Figs. 5A-5B, Figs.S21-S23). Plasmid pUC19 containing KGs fragments was expressed by E. coli DH5α, while the types of the KGs included KGs-Clean (e.g., livH,parB, ssuD) (Fig. 5A-1, Fig. 5B-1), KGs-Low (e.g., gst, accC, fadA) (Fig. 5A-2, Fig. 5B-2) and KGs-High (e.g., yfcG, nuoD, coxA) (Fig. 5A-3, Fig. 5B-3,Fig. S22).As p-NCB-exposing increased,the KGs-Clean, KGs-Low, KGs-High reached the highest gene copy numbers in 0, 10 and 100 mg L-1concentrations, respectively, around4.04×109-1.13×1012 copies μL-1.In the other two concentrations, the absolute abundance of these genes declined significantly, which was consistent with the results in Fig. 3D (Fig. 3). Except for the genes parB and coxA, the copy numbers of the KGs-Clean, KGs-Low and KGs-High in eDNA fractions were the lowest at 0 mg L-1p-NCB-levels compared with the other two concentrations, ranging from 0-7.24×1011 copies μL-1. As p-NCB levels increased, their absolute abundance increased significantly and reached around9.52×107-6.50×1012 copies μL-1 (Fig. 5), further verifying that the KGs had niche specificity to combat with p-NCB-induced stress and the genes exchanges between intracellular and extracellular fractions were enhanced. Copy numbers of the non-KGs (includingprkC, fusA, livK, uvrA, asnB, yfbK, dnaE, ftsH and topA) in iDNA declined significantly and increased dramatically in eDNA fractions as p-NCB-levels increased. At the meantime, the copy numbers of the non-KGs and KGs differed by 1.11~2.50 orders of magnitude(Fig. S23).
As p-NCB-exposing eliminated (e.g., 100→10mg L-1:p-NCB-levels declined from 100 mg L-1 to 10 mg L-1; 10→0mg L-1:p-NCB-levels declined from 10mg L-1 to 0 mg L-1;), the absolute abundance of the KGs-Clean, KGs-Low, KGs-High ranged from 4.28×109-8.42×1011 copies μL-1, which fell between the best and poorest growth states compared with the p-NCB-increasing systems, except for the genesparB and yfcG. While in eDNA fractions, the gene copy numbers of the KGs fell the range of 5.15×105-7.55×1011 copies μL-1, which decreased by 1.23~2.37 orders of magnitude compared with the 100 mg L-1systems and increased by 1.00~5.41 orders of magnitude compared with the 0 mg L-1systems, respectively (Fig. 5B& Fig. S22). Correspondingly, the copy numbers of these KGs both in intra-and extracellular fractions of the p-NCB-eliminating systemswere not comparable to the optimal growth states of the p-NCB-increasing systems. Additionally, the variations of KGs-Clean(e.g., livH,ssuD)and KGs-Low(gst,accC,fadA) both in e-iDNA fractions ofp-NCB-eliminating systems were significantly higher than those in eDNA fractions of p-NCB-increasing systems (p<0.01,Fig. 5B& Fig. S22), indicating that as p-NCB-levels declined, the increasing variations of KGs in iDNA and decreasing variations in eDNA were promoted and might offset the environmental stressof the microbiome. As expected, variations of the non-KGs in iDNA fractions showed the similar promotion in the p-NCB-eliminating medium, butdiffered by 0~1.44 orders of magnitude from the KGs (Fig. S23).
Intercellular exchanges of keystone genes were observed in soil microbial adaptation after OCPs stress
According to the gene exchanges of the KGs in the simplified cultivation, a 100-day field observationstudywas carried out in ten different OCP-polluted sites, in China, to determine the recovery potential of KGs with OCPs-induced environmental stress (Figs. 5C-5E, Method and Figs. S24-S26) [14, 18, 44]. The OCP levels were in the range of below the limit of detection and as high as 6656.04mg kg-1, meanwhile the site in Jiangsu (Nanjing) was taken as the background reference. During the 100-day observation, the gene copy numbers of the KGs among different sites showed a similar increase, which were of the same orders of magnitudes at the same observation time, except for gene livHin Jiangsu (Suzhou). It demonstrated that the total contents of KGs remained stability, regardless of the soilphysicochemical properties and the OCP-levels, which was consistent with the results in Figs. 3B-3C.Simultaneously, the copy numbers of the KGs ranged from 0-7.74×105 copies μL-1and 0-2.14×106 copies μL-1, while the non-KGs ranged from 0-1.77×107 copies μL-1 and 1.14×102-5.92×108 copies μL-1ine-andi-DNA fractions, respectively, indicating that the absolute abundance of KGs was not comparative with that of non-KGs (except livH) (Figs. 5D-5E & Figs. S24-26).
In addition, during the observation period, aside from the gene livH,the copy numbers of the KGs showed 0~1.50 orders of magnitudeat day 20 and day 100 in Jiangsu (Nanjing) (Fig. 5D). While inthe other sites (e.g., Zhejiang (Hangzhou), Jiangxi (Nanchang), Shanxi (Yuncheng) and Heilongjiang (Suihua)), copy numbers of the KGs at day 20 and day 100 differed by 0~7.4 (iDNA) and 0~4.2 (eDNA) orders of magnitude(Figs. 5D-5E & Figs. S24-26), reflecting that the increasing slopes of the KGs in OCP-induced sites were significantly higher than those in the background reference. Moreover, the gene variations were more sensitive in iDNA fractions than those in eDNA fractions, which were consistent with the metagenomic results in Fig. 1. As we expected, the increasing slopes of the non-KGs showed a high similarity with the KGs, but differed by 1.51~3.74 orders of magnitude from the KGs (Figs. S24-26).
Furthermore, the absolute abundance of KGs differed by 1.12~3.60 orders of magnitude in intra- and extracellular DNA fractions, except for Jiangsu (Nanjing) (Fig. 5D). As the soil properties and OCP-levels varied, the difference magnitudes of the gene copy numbers of KGs between iDNA and eDNA fractions increased, implying that the intracellular gene variations of KGs were enhanced (Fig. 5E& Figs. S24-26). Notably, some KGs (e.g., ssuD in Shapingba, livH in Nanchang, coxAand gst in Beichen,livH and yfcGin Yuncheng)showed opposite variations in iDNA and eDNA (Figs. 5D-5E & Figs. S24-26), which further indicated that the gene correlations and exchanges between e-iDNA were might promoted in these observation sites.