Dibutyl phthalate alters the metabolic pathways of microbes in black soils

Dibutyl phthalate (DBP) is well known as a high-priority pollutant. This study explored the impacts of DBP on the metabolic pathways of microbes in black soils in the short term (20 days). The results showed that the microbial communities were changed in black soils with DBP. In nitrogen cycling, the abundances of the genes were elevated by DBP. DBP contamination facilitated 3′-phosphoadenosine-5′-phosphosulfate (PAPS) formation, and the gene flux of sulfate metabolism was increased. The total abundances of ABC transporters and the gene abundances of the monosaccharide-transporting ATPases MalK and MsmK were increased by DBP. The total abundance of two-component system (TCS) genes and the gene abundances of malate dehydrogenase, histidine kinase and citryl-CoA lyase were increased after DBP contamination. The total abundance of phosphotransferase system (PTS) genes and the gene abundances of phosphotransferase, Crr and BglF were raised by DBP. The increased gene abundances of ABC transporters, TCS and PTS could be the reasons for the acceleration of nitrogen, carbon and sulfate metabolism. The degrading-genes of DBP were increased markedly in soil exposed to DBP. In summary, DBP contamination altered the microbial community and enhanced the gene abundances of the carbon, nitrogen and sulfur metabolism in black soils in the short term.

Microbial community structure. The taxonomic profiles of the metagenome analysis, which contained domain, kingdom, phylum, class, order, family, genus and species, were constructed. With increasing DBP concentration, the percentages of Eukaryota were decreased from 2% to 0.9%, and the percentage of total bacteria was increased from 81% to 92%. We focused on the genus level to explain the changes in the microbial community. Based on the metagenome analysis, the results were shown in Fig. 2. The relative populations of the genera Arthrobacter and Nocardioides were increased along with the increase of DBP concentration. The relative abundances of Arthrobacter in CK, DBP1, DBP2, DBP3 and DBP4 were 1.67%, 6.82%, 10.56%, 22.82% and 39.86%, respectively. The relative abundance of Nocardioides in CK, DBP1, DBP2, DBP3 and DBP4 were 1.75%, 1.91%, 2.05%, 3.07% and 6.09%, respectively. On the other hand, the relative levels of some bacteria were decreased by DBP in all samples and showed a negative correlation with DBP concentration; these bacteria included Nitrososphaera, Lactococcus, Streptomyces, Frankia, Solirubrobacter, Conexibacter, Mycobacterium, Rubrobacter, Tetrasphaera, Pseudomonas, Intrasporangium, Actinoplanes, Actinomadura, Enterococcus, Amycolatopsis, Rhodococcus, Pseudonocardia, Acidothermus, Nocardiopsis, and others.  Annotation of metabolic pathways by KEGG analysis. Metabolic pathways, including nitrogen cycling, carbon metabolism (glycolysis, TCA cycle and pentose phosphate), sulfur metabolism and signal regulatory pathways, were identified in DBP-contaminated black soils based on the KEGG database. The original data are shown in Supplementary Tables 1 and 2. As shown in Fig. 3, DBP contamination changed the flux of nitrogen cycling. The total abundance of nitrogen cycling genes rose significantly in response to the increase of DBP concentration (F = 25.71, P < 0.01) ( Fig. 3A and Supplementary Table 7). Nitrite reductase (NirBD) (F = 234.76, P < 0.01), nitrate reductase (NarGHIJ, NasAB and NxrAB) (F = 136.06, P < 0.01), ferredoxin-nitrite reductase (NirA) (F = 188.95, P < 0.01) and nitric oxide reductase (NorBC) genes (F = 43.89, P < 0.01) are the key enzyme genes in nitrogen cycling, and their copy numbers were significantly increased after DBP contamination (Fig. 3B,C and Supplementary Table 7). The copy numbers of related enzyme genes in nitrogen cycling, such as the formamidase (F = 32.49, P < 0.01), nitrilase (F = 7.64, P < 0.05) and nitronate monooxygenase (F = 177.61, P < 0.01) genes, were increased significantly by DBP contamination (Fig. 3C and Supplementary Table 7).
The total abundance of sulfur metabolism genes was higher in the contaminated samples than in CK (F = 21.69, P < 0.01) ( Fig. 7A and Supplementary Table 11). The gene abundance of sulfate adenylyl-transferase (F = 22.91, P < 0.01), which included Sat, CysND and CysC, was promoted to various extents with increasing DBP concentration (Fig. 7B). The gene abundances of sulfite reductase (F = 18.78, P < 0.01) was significantly enlarged in DBP1 and DBP3, but it was decreased in DBP2 and DBP4 ( Fig. 7B and Supplementary Table 11). The gene Nonmetric multidimensional scaling (NMDS) analysis. To explain the overall differences among the CK and DBP contamination samples, NMDS analysis was performed at the levels of genes, taxonomy and metabolic pathways. The results indicated that DBP contamination altered the distribution of genes, taxons and metabolic genes (Fig. 9). Furthermore, the difference became larger with increasing DBP concentrations.
DBP degradation genes. The main DBP degradation genes are pcmA, pehA, phtAb, phtB and phtC. The copy numbers of pcmA, pehA, phtAb, phtB and phtC were significantly increased under DBP contamination (Fig. 10). When the concentration of DBP increased, the copy numbers increased as well, and the range of increase varied. The lowest degree of change was that of pcmA, which increased from 3.06 × 10 4 to 1.05 × 10 5 . The greatest change was that of phtB, which increased from 1.55 × 10 3 to 4.39 × 10 5 . Additionally, the copy numbers of pehA and phtAb genes were elevated along with increasing DBP concentration, increasing from 3.29 × 10 3 to 2.42 × 10 5 and from 3.13 × 10 4 to 1.62 × 10 5 , respectively. After DBP contamination, the copy number of phtC was higher compared with CK, increasing from 94.49 to 5.35 × 10 4 . The results were consistent with the change of DBP concentration in black soil (Fig. 1).

Discussion
Using metagenomic analysis, this study found that the microbial community structure was changed in black soil under DBP contamination in the short term. The abundances of some genera, such as Arthrobacter and Nocardioides, were significantly enriched in response to DBP treatment. Previous studies showed that Arthrobacter and Nocardioides could degrade organic contamination 19,20 . This phenomenon could be further demonstrated by comparing the copy numbers of DBP degradation genes (Fig. 10). However, the growth for some genera that are indispensable for nutrient cycling 21 33 and Tetrasphaera 34 are beneficial to soil health and plant growth. Amycolatopsis and Intrasporangium were also decreased even though they could remediate heavy metal pollution 35,36 . The changes in the microbial community indicated that DBP contamination disturbed the soil ecosystem in the short term.
Nitrogen cycling, carbon metabolism and sulfur metabolism are the basis of nutrient cycling in soil ecosystems [37][38][39] . The turnover of the nitrogen cycle, the glycolysis pathway, the TCA cycle, the pentose phosphate pathway and sulfur metabolism were accelerated significantly when soil microorganisms were exposed to DBP contamination. In nitrogen cycling, the genes NarGHIJ, NasAB, NxrAB, NirBD, and NorBC are the key genes to the dissimilatory nitrate reduction, assimilatory nitrate reduction, denitrification and nitrification, and their abundances were elevated due to DBP contamination. Meanwhile, the gene abundances of NorBC and NarGHIJ were increased. Because NorBC and NarGHIJ catalyze the generation of nitrous oxide (N 2 O) and nitrite 40 , nitrite accumulation and N 2 O emission could be increased. In this study, it was found that the turnover of the carbon metabolic pathways was accelerated after DBP treatment such that more carbons in the substrate were converted to CO 2 , which meant that more carbon could be consumed in the soils. It was also found that the gene abundance of sulfate adenylyltransferase was increased in the soil microorganisms. The phenomenon implied that DBP contamination facilitated 3′-phosphoadenosine-5′-phosphosulfate (PAPS) formation. This result could indicate that the flux of sulfate metabolism was increased and that more sulfate turned toward the direction of oxidative metabolism after DBP contamination 39,41,42 . For all of these results, the metabolic processes of soil microorganisms were accelerated including nitrogen cycling, carbon metabolism and sulfur metabolism, which could affect the nutrient transformation and the quality of black soils in the short term.
ABC transporters are found in all domains of organisms and transfer a remarkably wide range of substrates into and out of living cells 43,44 . ABC transporters are useful for transporting substrates including sugars, amino acids, peptides, ions and other molecules as well as for absorbing nutrients into cells 45,46 . Increases in the total abundance of ABC transporters and the relative gene abundances of the monosaccharide-transporting ATPases MalK and MsmK could be the reason for the accelerated metabolism of nitrogen, carbon and sulfate, which may accelerate the depletion of nutrients in black soils. In prokaryotes, TCS is a sophisticated signal transduction strategy and senses any alteration in the environment 47 . It regulates downstream gene expression, thus orchestrating several cellular responses, and participates in biochemical processes, acting as an energy carrier 48 . The total abundance of TCS genes and the abundances of the malate dehydrogenase, histidine kinase and citryl-CoA lyase genes were increased after DBP contamination, suggesting that it could lead to accelerated metabolism in black soils, including nitrogen, carbon and sulfate metabolism. PTS is the key signal transduction pathway involved in the regulation of central carbon and nitrogen metabolism in bacteria 49,50 . The total abundance of PTS genes and the gene abundances of phosphotransferase, Crr and BglF were raised by DBP. These manifestations indicated that

Conclusion
Based on the results of metagenome analysis and qPCR analysis in the short term (20 days), we propose that DBP contamination altered the structure of the microbial community; improved the gene abundances of nitrogen cycling, glycolysis, the TCA cycle, the pentose phosphate pathway and sulfate metabolism; and increased the abundance of DBP degradation genes. The increased abundances of signal regulatory pathways for ABC transporters, TCS and PTS could be the reason for the changes of gene abundances in carbon, nitrogen and sulfur metabolism. However, it remains to be seen whether structure and function are impacted over the long-term.

Materials and Methods
Sample collection and treatment. DBP (chromatographically pure grade) was purchased from the Information Center for Standard Samples (Beijing, China). Acetone (chromatographically pure grade) was acquired from Traditional Chinese Medicine (Beijing, China). DBP was dissolved in acetone at a ratio of 1 to 9 (m/m) and then stored in the dark at 4 °C.   Table 1. Soils were sieved by 0.3 mm mesh. After blending and sieving, the soils were divided into 15 samples of 1000 g each. Afterwards, the soil was moistened up to 30% and preincubated at 25 °C for 7 days. Varied amounts of DBP were added to the soil samples and mixed well, as follows: Control (CK), 0 mg DBP kg −1 soil; DBP sample 1 (DBP1), 5 mg DBP kg −1 soil; DBP sample 2 (DBP2), 10 mg DBP kg −1 soil; DBP sample 3 (DBP3), 20 mg DBP kg −1 soil; and DBP sample 4 (DBP4), 40 mg DBP kg −1 soil. The measured concentrations of DBP from the different treatments are shown in the Supplementary Figure 1. Each concentration stress was repeated 3 times. All the samples were exposed to the air for 3 h until the acetone was evaporated completely. The soil samples were placed in porcelain pots and incubated in the dark at 25 °C for 20 days. During the period of experiment, all samples were maintained at 30% moisture. DBP concentrations were measured in black soils incubated for 0 d and 20 days in the dark at 25 °C and 70% air humidity in constant temperature and humidity incubator.  Table 3). DNA concentration, purity, and integrity were detected by NanoDrop2000, TBS-380 and agarose gel electrophoresis, respectively. 1 μg DNA was used to build library. Briefly, the DNAs were sheared using an M220 Focused-ultrasonicator ™ (Covaris Inc., Woburn, MA, USA) and the parameters of instrument were set according to Supplementary Table 4. The PE library was composed of ~300 bp DNA fragments. The TruSeq ® DNA Sample Prep Kit was used to prepare DNA library based on the manufacturer's protocol (www. illumina.com). Primer hybridization sites were included in dual-index adapters and ligated to the blunt-end fragments. To improve quality of the template DNA, paired-end sequencing (2 × 151 bp) was conducted on an Illumina Genome Analyzer (Illumina Inc., USA) by Majorbio Bio-Pharm Technology Co., Ltd. (Shanghai, China). All the raw metagenomic datasets have been deposited into the NCBI Sequence Read Achieve (accession numbers: SRP094732, SRP096214).

Sequence quality control, assembly and gene prediction.
To improve the quality and reliability of subsequent analysis, the raw data were processed using the following steps: (1) SeqPrep (https://github.com/jstjohn/SeqPrep) was used to remove the reads that contained N bases in the sequence and strip the 3′-adaptors and 5′-adaptors. (2) To retain high-quality pair-end reads and single-end reads, the reads of sequence length <50 bp and quality score <20 were excised with Sickle (https://github.com/najoshi/sickle). Clean reads were assembled with SOAPdenovo (http://soap.genomics.org.cn,Version 1.06) with a default k-mer length (39)(40)(41)(42)(43)(44)(45)(46)(47), and the de Bruijn graphs were used to build contigs. For the optimal assembly, scaffolds over 500 bp were divided into new contigs without gaps. The open reading frames (ORFs) of the contigs from the assembled results were predicted using MetaGene Annotator (http://metagene.cb.k.u-tokyo.ac.jp/). , and protein information resources from corresponding coding sequences (CDS) on RefSeqG and GenBank. The BLAST parameters were set as follows: E-value cutoff = 1*10 −5 , num_alignments = 250, and num_descriptions = 500. The taxonomy annotation was obtained based on the taxonomic information database corresponding to the NCBI NR library. The species abundance was calculated using the total gene abundance across all species, and the abundance profile was built based on taxonomic levels, which were the domain, kingdom, phylum, class, order, family, genus and species. The KEGG database (Kyoto Encyclopedia of Genes and Genomes, http://www.genome.jp/kegg/) is a huge library for systematic analysis of gene function, combining genome and function information. It provides sequence information on the genes and proteins identified in genome projects. Further, the database includes various pathways, such as metabolic, biosynthetic, membrane transport, signaling, cell cycle and disease pathways. In addition, it collects records of all kinds of enzymes and other molecules as well as enzymatic reactions and related information. The KEGG pathway annotation was done using KOBAS 2.0 (KEGG Orthology Based Annotation System, http://kobas.cbi.pku.edu.cn/home.do) according to the alignment results between the gene set sequence and the KEGG database using BLAST (BLAST Version 2.2.28+, http://blast.ncbi.nlm.nih.gov/Blast.cgi). The amassed read data for each sample are presented in Supplementary Table 5.
Nonmetric multidimensional scaling (NMDS). Nonmetric multidimensional scaling (NMDS) is a monotonic function that evaluates the similarities or differences of different samples based on distance. In this study, NMDS represented the overall difference between control and DBP (varied concentration) contamination samples in taxonomy, genes and KEGG pathways. NMDS was generated by R language release package (MASS package).
Real-time fluorescent quantitative PCR (qPCR). qPCR was performed to detect the genes associated with the degradation of DBP, including phthalate ester hydrolase (pehA), 3,4-dihydrophthalate dehydrogenase (phtB), phthalate dioxygenase small subunit (phtAb), protocatechuate 4,5-dioxygenase (pcmA) and 3,4-dihydrophthalate decarboxylase (phtC) 52 . The primer information for these genes is shown in Table 2. PCR was conducted in 20-μL reactions containing 10 μL of 2 × SYBR green qPCR master mix (premix of dNTPs, Taq DNA polymerase, PCR buffers and SYBR green), 0.4 μL of primer F and 0.4 μL of primer R of each gene-specific primer pair, 7.2 μL of ddH 2 O and 2 μL of template (DNA) based on the manufacturer's instructions (LightCycler480, Roche). DNA concentration was shown in Supplementary Table 6. The qPCR thermal cycling steps are shown in Table 3. The melting curves for the genes were obtained during the detection, representing the specificity of the amplification. A linear relationship was shown based on standard curves ranging from 10 3 to 10 10 copies. The amplification efficiencies ranged from 71.2% to 107.8%, and R 2 ranged from 0.9901 to 0.9984 (Supplementary Figure 2). Analyses of variance (ANOVA) and the least significant ranges test (Duncan's method) were performed with SPSS 22.0 to test the significance (p < 0.05) of the differences among the treatments.   Table 3. qPCR thermal cycling steps.