General statistics of Pseudomonas putida YC-AE1 whole genome
The whole genome of Pseudomonas putida YC-AE1 consists of a circular chromosome and three plasmids with a total of 6,992,587 bp and GC content of 60.62%. The size of chromosome is 6,053,381 bp with GC content of 61.29% while, the 3 plasmids size was 504,084 bp, 388,915 bp and 46,207 bp with GC content of 56.28%, 56.09% and 57.73%, respectively. We identified Eighty-two copies of tRNA and twenty-two copies of rRNA in the bacterial genome.
RNA-seq analysis
The total clean reads of control samples (CK) were 16.1, 14.89, and 16.07 Mb obtained from CK1, CK2, and CK3, respectively, with an average of 15.68 Mb and 90% clean reads. On the other hand, the total clean reads of treated samples (T) were 15.90, 15.82, and 14.74 Mb for T1, T2, and T3, respectively, with an average of 15.48 Mb 88.92% clean reads ratio. The Q20 and Q30 of all samples were more than 98 and 95%, respectively, indicating the high quality of clean reads [41], (Table S2). The clean reads were mapped onto the genomic sequences of Pseudomonas putida YC-AE1, with high mapping rates, an average of 96% for control and 95.75% for treated groups. The new transcript prediction analysis yielded 1405 new transcripts, with 107 and 1298 for coding and non-coding transcript, respectively. The new transcript prediction analysis yielded 1405 new transcripts, most of them representing the non-coding region by 1298 reads, in which 710 for intergenic region and 588 for antisense to mRNA, coding region was represented by 107, calculated as 63 for intergenic region and 44 for antisense to mRNA (Table S3). Gene expression level was determined based on the expected number of FPKM attributed to gene numbers. As shown in Fig. S1, in FPKM < = 100, the average of gene numbers expressed in control and treated were 1523 and 1662, respectively, while, in FPKM > = 100, the expressed gene numbers were 1884 and 1635 on average for control and treated samples, respectively. The largest gene numbers were expressed in FPKM between 10–100, approximately 2558 and 2666 on average for control and treated samples, respectively.
To evaluate the degree of similarity between samples, the correlation and variance were calculated. In this study, a principal component analysis (PCA) and a heat map were performed to visualize this correlation. As depicted in the heat map and PCA (Fig. S2 and S3), the correlation between control and treated samples was high, which indicates the correctness of the sampling.
Identification and analysis of DEGs
The expression level of 6165 genes identified from RNA-seq between control and treated group were quantified by FPKM and folding change |Log2FC| with DEGseq2 into three groups. There was a total of 725 genes with an expression level at |Log2FC|≥ 1 and Padj ≤ 0.05 were identified as significantly up-regulated genes [41], while 504 genes at expression level |Log2FC|≤ -1 and Padj ≤ 0.05 were identified as significantly down-regulated, while, 4936 genes were regarded as non-significantly regulated at expression level abs|Log2FC|< 1 and Padj > 0.05. It is plausible that, BPA degradation related genes are up-regulated, while the down-regulated genes may reflex the response of Pseudomonas putida YC-AE1 to toxic effects upon exposure to BPA. [42], reported a total of 839 genes as DEGs in response to biodegradation of chlorimuron-ethyl by Klebsiella jilinsis 2N3 with 386 up-regulated and 453 down-regulated.
A volcano plot illustrated in Fig. 1a, which can quickly screen out genes with significant statistical differences, was used to visualize these data. The horizontal axis represents the value of the multiple of the difference after log2 conversion, and the vertical axis represents the significance value after -log 10 conversion. A higher position means that the gene has a smaller p-value in the differential pair. According to the color of the legend, it is divided into significantly up-regulated DEGs, significantly down-regulated DEGs, and non-regulated DEGs. A heat map was created to visualize the clusters of significantly DEGs between the control and treated groups, as shown in Fig. 1b. The heat map was adopted relying on FPKM of RNA-seq data and represented the genes expression levels. DEGs were classified in the bar according to the log10 (FPKM + 1) value.
Functional enrichment analysis of DEGs in Pseudomonas putida YC-AE1
Gene Ontology is an international classification system to unify and describe the gene and gene product attributes in organisms [43]. GO describes the molecular function, cellular component, and biological processes. All DEGs were categorized into the three terms in which biological process with 916 (30.1%), cellular component 112 (36.5%), and molecular function 1014 (33.3%) (Fig. 2a). DEGs in response to BPA degradation by YC-AE1 were associated mainly with the metabolic and cellular processes represented by 298 and 284 genes, respectively. In the context of molecular function, catalytic activity and binding were mapped with the highest number of DEGs, 437 and 357 genes, respectively. Altogether, these data indicated that the metabolic functions of Pseudomonas putida YC-AE1 were affected upon exposure to BPA, and the DEGs of these GO terms may be linked to the metabolizing process of BPA by YC-AE1 strain.
The pathway-based analysis is crucial to elucidate the biological functions of genes. To further analyze the potential role of DEGs, we performed a KEGG pathway analysis to identify the metabolic pathways enriched in the BPA degradation process [45]. Results showed that the 135 different metabolic pathways between CK and T were involved in five categories, cellular process (CP), environmental information processing (EIP), genetic information processing (GIP), metabolism (M), and organismal systems (OS). A total of 929 DEGs were involved in the pathways with 669 genes (72%) in M category, followed by EIP with 139 genes (~ 15%) (Fig. 2b).
The numbers of regulated genes were illustrated in Fig. 3; the highest number of up and down-regulated genes belong to microbial metabolism in diverse environments pathway with 63 and 41 DEGs, respectively, followed by two-component system pathway and carbon metabolism pathways with 86 and 43 DEGs. Efflux pump related genes have been reported to play a crucial role in microbial environmental adaptability, for instance, SrpABC from Pseudomonas putida strain B6-2 (DSM 28,064) enhanced the degradation capacity toward polycyclic aromatic hydrocarbon (PAH) via extruding toxic intermediates during the degradation [46]. In contrast, EmhABC pump reduced the degradation capability of Pseudomonas fluorescens strain LP6 via excreting out its substrate [47]. TtgABC efflux machinery enhances phenol tolerance of Pseudomonas putida strain DOT-T1E because this pump contributes to toluene tolerance [48, 49]. Adaptation of bacterial cell to environmental pollutants and orchestrating their gene response accordingly largely reside on two-component signal systems (TCS) [50]. Pseudomonas putida ColR (TCS) deleted strain showed attenuation in phenol tolerance compared with wild type [51]. In our study, thirty efflux pump related genes were regulated; in which 12 genes were positively regulated and 18 genes were negatively regulated. In the same context, 193 genes related to two-component system have regulated, 86 up regulated and 107 down regulated. However their specific roles in PBA degradation still unknown and required further investigation.
According to the intermediate compounds and the proposed pathway of BPA degradation by YC-AE1 in our previous study [21], accompanied by data generated from transcriptome analysis, we hypothesized that 15 genes could be involved in BPA metabolism. The up-regulated genes in the degradation process are probably associated with BPA metabolism or the improvement of BPA tolerance in YC-AE1 strain. Overall, the mRNA expression levels of the investigated 15 genes were significantly higher than the control group. The information about the regulated genes, including name, ID, and EC number for encoding enzymes, are summarized in Table S4 in the Supplementary Material. The degradation pathways of BPA by YC-AE1 obtained from transcriptome analysis and pathway enrichment are illustrated in Fig. 4. As shown in Fig. 4, CYP450 monooxygenase (Gene ID: YCAE1GL005919) was proposed to initiate the first step in BPA degradation in which BPA is transformed to 1,2-bis(4-hydroxyphenyl)-2-propanol and/or 2,2-bis(4-hydroxyphenyl)-1-propanol, this finding was consistent with the previous report by Sasaki et al. [52]. CYP450 monooxygenases are a superfamily of ubiquitous heme- monooxygenases that participate in the microbial degradation of diverse materials such as oils, fuel additives, and chlorinated hydrocarbons [53]. Cheng et al. [54] revealed the potential role of CYP450 in the biodegradation of Chlorimuron-ethyl by Rhodococcus erythropolis D310-1. Moreover, the CYP450 mRNA level was significantly up-regulated in the aminobenzoate degradation pathway, inconsistent with the production of the degradation products [42].
Our results also showed that the dehydratase enzymes were significantly up-regulated. These enzymes are involved in the transformation of 1,2-bis(4-hydroxyphenyl)-2-propanol to 4,4-dihydroxy-alpha-methylstilbene which is further converted to p-hydroxyacetophenone and p-hydroxybenzaldehyde. To detect whether YC-AE1 can further metabolize p-HBAL and p-HAP, a TEM medium supplemented with 100 mg l−1 p-HBAL and p-HAP as a sole carbon source was inoculated with YC-AE1 and incubated at 30 oC for 24 h. YC-AE1 strain showed 100% degradation capacity of p-HBAL and p-HAP in 24 h (data not shown); this is evidence that YC-AE1 strain could mineralize p-HBAL and p-HAP. Our finding can interpret the high mineralization rate for BPA by YC-AE1 strain in our previous study [21]. Also, our RNA seq data showes significant upregulation to oxidoreductase related genes.
The activity of oxidoreductase act on CH-CH group and CH-OH group [55]. Twenty-one up-regulated genes were detected in our transcriptome with an expression level at |Log2FC|≥ 1 and Padj ≤ 0.05. Our results indicated that Pseudomonas putida YC-AE1 over-expressed oxidoreductase encoding genes to improve the biodegradation of BPA oxidation and to form hydroxylated benzene rings of BPA. Our data was compatible with data reported by wang et al. [55] in BPA biodegradation by the green alga Desmodesmus sp.WR1. García-Rodríguez et al. [56] reported that phenolic compounds could be oxidized by oxidoreductases like polyphenol oxidase and peroxidase.
qRT-PCR Verification
Ten genes involved in the BPA degradation pathway were selected for further verification using qRT-PCR analysis. Genes were abbreviated as follows: CYP450 monooxygenase (cyp), 4-hydroxybenzaldehyde dehydrogenase (hbd), 4-hydroxybenzoate 3-monooxygenase (hbm), γ-carboxymucono-lactone hydrolase (clh), 4-carboxymuconolactone decarboxylase (cmd), 3-oxoadipate enol-lactonase (oel), protocatechuate 3,4-dioxygenase (pcd), 3-oxoadipate CoA-transferase (oct), 4-hydroxyacetophenone hydrolase (hah) and hydroxyquinol 1,2-dioxygenase (hyd). Figure 5 showed the relative transcription level of selected genes in both RNA-seq data and qRT-PCR. The result was consistent with RNA-seq data and confirmed the upregulation of the BPA degradation genes in the YC-AE1 strain.
CYP
450
mediates BPA degradation in YC-AE1
To confirm the involvement of CYP450 in BPA degradation by YC-AE1, we used ABT, which is known as an inhibitor for CYP. The addition of 0.1, 0.4, 0.7, 1, and 2 mmol of ABT significantly decreased the efficiency of YC-AE1 to degrade BPA to 83, 76.7, 62.6, 46, and 2%, respectively, compared to 95% BPA degradation in control, Fig. 6a. These findings indicated that CYP is involved in BPA catabolism, consistent with the previous studies that reported the involvement of CYP in the hydrolysis of many xenobiotic compounds [57, 58]. Mtibaà et al. [57] Reported a slight decrease in BPA removal rate in the presence of ABT as a CYP450 inhibitor. Further, Wei et al. [58] reported a 55% decrease in triphenyl phosphate (TPHP) degradation by Bacillus brevis in the presence of 1 mmol of piperonyl butoxide (PB) as a CYP450 inhibitor. Sasaki et al. [59] reported that CYP encoded by bisd in Sphingomonas bisphenolicum AO1 is responsible for the transformation of BPA. This finding indicates the indispensable role of CYP450 as an essential initiator during BPA degradation in YC-AE1 strain. E. coli genome encodes ferredoxins and ferredoxin reductases, indicating that E. coli cells bearing bisdB might degrade BPA [25]. We cloned the bisdA and bisdB genes from YC-AE1 in E. coli BL21 (DE3). As shown in Fig. 6b, the recombinant E.coli carrying pET-32a-bisdB could degrade 12 mg l−1 of BPA after 24 h of incubation, while E.coli harboring pET-32a-bisdAB possess the ability to degrade about 66 mg l−1 after 24 h incubation. Ferredoxin of Pseudomonas putida YC-AE1 was implicated in BPA degradation by CYPbisd, and that can be seen obviously from the high difference in BPA degradation rate in E.coli harboring only bisdB, and that contain bisdAB. Together, these results suggested that the ferredoxins of E. coli could not meet the requirement well enough for a high degradation rate of BPA [20].
Knockout of bisdB in YC-AE1
The previously mentioned data in this study provided evidence that CYPbisd is involved in BPA degradation by YC-AE1; however, this is not reliable evidence that no other genes in YC-AE1 can perform the same function. Thus, it is necessary to test whether bisdB deleted YC-AE1 could metabolize BPA. Therefore, bisdB YC-AE1 knockout strain was constructed and tested for BPA degradation. Interestingly, YC-AE1ΔbisdB completely lost its ability to degrade BPA and failed to grow on BPA as a sole source of carbon and energy, Fig. 6c. This implied that bisdB gene is indispensable for initiating BPA catabolism in our strain. In the same context, Zhang et al. [42] reported the role of CYP encoded by Kj-CysJ in the biodegradation of Chlorimuron-ethyl by Klebsiella jilinsis 2N3, 2N3 deleted Kj-CysJ showed a lower degradation rate compared to wild type and was no longer able to utilize Chlorimuron-ethyl as carbon source. This finding, coupled with those of qRT-PCR and transcriptomic analysis, confirmed that CYP450 was a key enzyme in the biodegradation process of BPA in YC-AE1 strain.
Phylogenetic analysis of P450 encoded by YC-AE1 bisdB and related bacterial CYP450s are depicted in Fig. 7. P450bisdB was clustered with P450bisdB from Sphingobium sp. YC-JY1 showed high similarity (99.5%) and 98.3% with P450 from Sphingomonas bisphenolicum AO1. P450bisdB showed little similarity with other P450s. The entire sequence of amino acids for BisdA and BisdB of YC-AE1 strain were shown in the supplementary data.