Relationship between the binding free energy and PCBs’ migration, persistence, toxicity and bioaccumulation using a combination of the molecular docking method and 3D-QSAR

Zhao et al. Chemistry Central Journal Relationship between the binding free energy and PCBs' migration, persistence, toxicity and bioaccumulation using a combination of the molecular docking method and 3D-QSAR Xiao‑Hui Zhao 0 1 Xiao‑Lei Wang 0 1 Yu Li 0 1 0 College of Environmental Science and Engineering, North China Electric Power University , No. 2, Beinong Road, Beijing 102206 , China 1 The Moe Key Laboratory of Resources and Evironmental Systems Optimization, North China Electric Power University , Beijing 102206 , China The molecular docking method was used to calculate the binding free energies between biphenyl dioxygenase and 209 polychlorinated biphenyl (PCB) congeners. The relationships between the calculated binding free energies and migration (octanol-air partition coefficients, KOA), persistence (half‑ life, t1/2), toxicity (half maximal inhibitory concentration, IC50), and bioaccumulation (bioconcentration factor, BCF) values for the PCBs were used to gain insight into the degradation of PCBs in the presence of biphenyl dioxygenase. The relationships between the calculated binding free energies and the molecular weights, KOA, BCF, and t1/2 values for the PCBs were statistically significant (P < 0.01), whereas the relationship between the calculated binding free energies and the IC50 for the PCBs was not statistically significant (P > 0.05). The electrostatic field, derived from three‑ dimensional quantitative structure-activity relationship studies, was a primary factor governing the binding free energy, which agreed with literature findings for KOA, t1/2, and BCF. Comparative molecular field analysis and comparative molecular similarity indices analysis contour maps showed that the binding free energies, KOA, t1/2, and BCF values for the PCBs decreased simultaneously when substituents with electropositive groups at the 3‑ position or electronegative groups at the 3′‑ position were introduced. This indicated the binding free energy was correlated with the persistent organic pollutant characteristics of PCBs. Furthermore, low binding free energies improved the degradation of the PCBs and simultaneously decreased the KOA, t1/2, and BCF values, thereby reducing the persistent organic pollutant characteristics of PCBs in the environment. These results are expected to be beneficial in providing a theoretical foundation for further elucidation of the degradation and molecular modification of PCBs. Polychlorinated biphenyl; Molecular docking; Biphenyl dioxygenase; Pearson correlation; Three‑ dimensional quantitative structure-activity relationship Introduction Polychlorinated biphenyls (PCBs) are a typical persistent organic pollutant with an aromatic biphenyl skeleton containing one to ten chlorine atoms that can theoretically yield 209 different congeners [ 1 ]. PCBS are chemically and thermally stable and have accumulated in soil, sediments, and the atmosphere where they can harm human health and the environment because of their toxicity. From when they were first produced in the 1930s to when they were banned in the 1990s, about 1.3 million tons of PCBs were produced, and tens of thousands of tons are known to have been released into the environment causing widespread pollution [ 2 ]. PCBs are of great concern because of their persistence, high lipophilicity, and adverse effects on organisms. Consequently, numerous studies have investigated the levels of PCBs in various environmental samples, such as sediments [ 3, 4 ], indoor air [5], residential carpet dust [ 6 ], chicken egg yolks [ 7 ], fish and meat [ 8 ], and human breast milk [ 9, 10 ], blood [11], and adipose tissue [ 12 ]. Because PCBs tend to accumulate in the food chain, they present a high risk to human and animal health [ 13 ], and complete degradation of PCBs in the environment is essential to avoid this. Biodegradation of PCBs by microbes has been studied extensively. This process mainly occurs via oxidation by enzymatic catalysis. Since 1973, many aerobic bacteria capable of degrading PCB congeners have been isolated, including Burkholderia LB400 [ 14 ], Rhodococcus sp. strain RHA1 [ 15 ], Alcaligenes eutrophus H850 [ 16 ], and Enterobacter sp. LY402 [ 17 ]. Biphenyl dioxygenase (BphA) plays a major role in the degradation of PCBs [ 18 ]. BphA uses O2 and electrons to catalyze the dihydroxylation of an aromatic ring as the initial PCB degradation step and determine the substrate specificity for the overall degradation pathway. Cao et al. [ 19 ] analyzed the PCB degradation abilities of BphA1 from Enterobacter sp. LY402 experimentally and with molecular simulation. They found that the binding free energies of the PCBs were well matched with the degradation rate constants (k) for PCBs with different numbers of chlorine substituents. In other words, the binding free energies of the PCBs decrease as k increases. Wu et  al. [ 20 ] compared the binding free energies of pollutants and receptors, and found that a lower binding free energy indicated higher affinity. Therefore, the binding free energy can be used to evaluate the degradation abilities and analyze the biodegradation of PCBs. Liu et al. [ 21 ] studied the relationship between the molecular characteristics and degradation rates of substrates degraded by Enterobacter sp. LY402. They determined that the dipole moment of a substrate was positively correlated with its degradation rate, and the stretch-bend energy was negatively correlated with the degradation rates. However, the effectiveness of microbial degradation of highly chlorinated PCBs and specific substituted PCBs is limited [ 22 ]. To gain insight into the molecular basis of degradation, the key enzymes involved in the PCB biodegradation process, specifically BphA, have been studied intensively [ 23–25 ]. Three-dimensional quantitative structure–activity relationship (3D-QSAR) studies can be used to evaluate the structure–activity relationships of a set of molecules using comparative molecular field analysis (CoMFA) and comparative molecular similarity indices analysis (CoMSIA). These two methods can avoid the insufficiency of the traditional two-dimensional model that is used to characterize the relationships between properties and structures, and have clear physical significance and provide rich information about the molecular field energy [ 26 ]. The objective of this study was to reveal the relationships between the binding free energies of BphA with PCBs and the PCBs’ migration, persistence, toxicity, and bioaccumulation. First, molecular docking was used to explore the interactions between BphA and 209 PCB congeners. Next, Pearson correlation analysis was performed to study the relationships between the binding free energies and the molecular weights, migration abilities (octanol–air partition coefficients, KOA), bioaccumulation (bioconcentration factors, BCF), environmental persistence (half-life, t1/2), and toxicities (half maximal inhibitory concentration, IC50) of the PCBs. Finally, QSAR CoMFA and CoMSIA were used to investigate the relationships between the binding free energies and the structures of the PCBs. These results will provide a theoretical foundation for further elucidation of the degradation and molecular modification of PCBs. Materials and methods Protein structures Crystal structures of nine types of BphA (PDB IDs: 1ULJ [ 27 ], 1WQL [ 28 ], 2YFJ [ 29 ], 2YFL [ 30 ], 2GBX [ 31 ], 2XSH [ 32 ], 2E4P [ 33 ], 3GZX [ 34 ], and 3GZY [ 34 ]) used in Surflex-Dock were obtained from the Protein Data Bank (http://www.rcsb.org/pdb/). The catalytic activities these BphA were different. The active sites were determined by the amino acid residues in the region where catalytic activity occurred. As different types of BphA may contain multiple active sites, all active sites were docked into the 209 PCBs to ensure the molecular docking was accurate. Molecular docking The molecular docking was performed using the Sybyl-x 2.0 molecular modeling package (Tripos Inc., St. Louis, MO) running on a Windows 7 32-bit workstation. The Surflex-Dock method in the Sybyl package was used to carry out molecular docking simulations to dock the ligands into a receptor’s ligand binding site and represent the interaction strength. 209 PCBs were minimized under the Tripos force field with MMFF94 [ 35 ] atomic partial charges by the Powell method, with a maximum iteration of 10,000 to reach a convergence gradient value of 0.001 kcal mol−1 Å. Before docking, the natural ligand and structural water molecules were removed from the crystal structure. Using the Biopolymer module, polar hydrogen atoms were added into the standard geometry. Kollman all-atom charges were assigned to protein atoms [ 36 ]. In this process, ligands were automatically docked into the binding site of the protein using a ProtoMol-based approach [ 37 ] with a patented search engine and an empirical scoring function. We applied the automatic docking method. Based on the software protocol of Surflex-Dock, ProtoMol was produced. “Proto Bloat” and “Proto Thresh” are two important factors that greatly affect the size and extent of the ProtoMol. “Proto Bloat” determines how far ProtoMol extends outside of the concavity, while “Proto Thresh” represents how far ProtoMol extends into the concavity of the target site [ 38 ]. In this process, for all ProtoMols produced, “Proto Thresh” was set to 0.5 and “Proto Bloat” was set to 1. All other parameters were at the default settings. Then, a binding pocket was generated. Successful molecular docking is indicated by a root mean square deviation of less than 2 Å [ 39 ]. According to the total docking scores, the first twenty conformations of every ligand were saved. The best conformations of the ligands were analyzed for their binding interactions [ 40 ]. In the molecular docking process, the receptor protein was considered to be rigid, and the ligand compounds were regarded as being flexible [ 36 ]. Surflex-Dock’s scoring function, which contains hydrophobic, polar, repulsive, entropic, and solvation terms, was trained to estimate the dissociation constant (Kd) expressed as − log(Kd) unit [ 41 ]. The docking results were ranked in a molecular spreadsheet. The conformer with the highest score was taken as the docking result [ 20 ]. For a better comparison between the binding free energies of the PCBs and BphA, the binding scores were converted to binding free energies (kcal/mol). The binding free energy was calculated as follows: where RT = 0.59 kcal/ mol [ 38 ]: Free energy of binding = RT × ln10−pKd (1) RT = 0.59 kcal/mol. 3D‑QSAR analysis The calculated binding free energies for 209 PCBs and BphA were used for modeling. CoMFA and CoMSIA methods were applied to build 3D-QSAR models using the structural parameters as independent variables and the binding free energy as the dependent variable to obtain the relationship between the binding free energy and molecular structure. In construction of the 3D-QSAR model, molecular alignment is an important step [ 42 ]. In this approach, decachlorobiphenyl, which had the largest binding free energy, was used as the template to align the remaining PCBs. The 209 PCBs were well aligned. CoMFA and CoMSIA analysis The CoMFA model includes steric and electrostatic descriptor fields. A 3D cubic lattice with a grid spacing of 2 Å and a sp3 carbon probe atom was created to calculate the CoMFA descriptor fields. The distance dependent dielectric constant was set at 1.0. The default energy cutoff value was set at 30 kcal mol−1. The CoMSIA model was based on the CoMFA model, which was composed of hydrophobic, hydrogen bond donor, and hydrogen bond acceptor fields. A sp3 carbon carrying a charge of +  1.0 was used as the probe atom to generate the CoMSIA descriptor fields. The default value of 0.3 was used as the attenuation factor α [ 43 ]. CoMSIA uses a Gaussian-type distance dependence to measure the relative attenuation of the field position of each atom in the lattice, and leads to much smoother sampling of the fields around the molecules than CoMFA. Compared with the CoMFA model, CoMSIA can avoid inherent defects, but the results are not necessarily good [ 44 ]. Hence, both methods (CoMFA and CoMSIA) can be used to verify and complement each other to obtain a reliable prediction model. Partial least squares regression Partial least squares regression analysis was applied to establish linear correlations between the binding free energies and the CoMFA and CoMSIA descriptors. Cross-validation with the leave-one-out method was carried out to yield the square of the cross-validation coefficient (q2) and the optimum number of components (N). Then, non-cross-validation analysis was performed with N and a column filter of 2.0 kcal mol−1 to produce regression models for CoMFA and CoMSIA and accelerate the analysis and decrease the noise. The 3D colored contour maps represent the relationship between the binding free energies and each molecular field. Pearson correlation analysis Pearson correlation analysis was performed using SPSS software to study the relationship between the binding free energies and the molecular weights, KOA, BCF, t1/2, and IC50 values of the PCBs. The KOA [ 45 ], BCF [ 46 ], t1/2 [ 47 ], and IC50 [ 48 ] values were obtained directly from the literature. The binding free energies were significantly correlated to the molecular weight, KOA, BCF, t1/2, and IC50 values when the P value was less than 0.01. The binding free energy was significantly correlated to the molecular weight, KOA, BCF, t1/2, and IC50 when P was more than 0.01 but less than 0.05. The binding free energy was not significantly correlated to the molecular weight, KOA, BCF, t1/2, or IC50 when P was more than 0.05. Results and discussion Comparative analysis of affinities between PCBs and BphA The root mean square deviation range was 0.0154– 0.9753, with all values less than 2  Å, indicating that Surflex-Dock was reliable in causing reappearance of the binding pattern of the ligand, and the parameter setting for molecular docking was appropriate. The docking results are shown in Fig.  1. The binding free energy was calculated using Eq. (1) and was used to evaluate the affinity between the ligand and receptor. Low binding free energies indicate high affinity or catalytic activity [ 20 ]. The degradation rate constants of the PCBs decreased as binding free energy increased (i.e., the binding affinity between the ligand and receptor decreased) [ 19 ]. As different types of BphA may contain multiple active sites, we selected a group from every BphA that had the highest comprehensive scoring function (Fig. 2). The binding free energies of the PCBs were well matched with the degradation rate constants (k) for the different numbers of chlorine substituents. A lower binding free energy represents stronger binding affinity or catalytic activity [ 20 ]. Both 1ULJ and 1WQL could degrade the isomers from chlorobiphenyl to tetrachlorobiphenyl (Fig.  2). The PCB degradation rates of 1ULJ were higher than those of 1WQL for all isomers except tetrachlorobiphenyl, which had a higher degradation rate with 1WQL than 1ULJ. Both 3GZX and 2YFL could degrade the isomers from chlorobiphenyl to hexachlorobiphenyl. The PCB degradation rate of 3GZX was higher than that of 2YFL. 2GBX, 2XSH, and 2YFJ degraded the isomers from chlorobiphenyl to pentachlorobiphenyl with PCB degradation rates decreasing in the following order: 2GBX > 2YFJ > 2XSH. For hexachlorobiphenyl and heptachlorobiphenyl, the degradation rates with 2XSH were higher than those with 2GBX and 2YFJ. 2E4P was able to degrade the isomers from chlorobiphenyl to octachlorobiphenyl, and 3GZY degraded the isomers from chlorobiphenyl to nonachlorobiphenyl. Although all of the types of BphA could degrade the less-chlorinated PCBs (one to four Cl atoms), the degradation rates were different. The PCB degradation rates of the different types of BphA increased in the following order: 2E4P < 3GZY < 1WQL < 1ULJ < 2XSH < 2YFL < 2 YFJ  <  3GZX  <  2GBX. 3GZY was able to degrade highly chlorinated PCBs (five to nine Cl atoms), but the degradation rates were not high. For pentachlorobiphenyl and hexachlorobiphenyl, the degradation rates of seven types of BphA were in the following order: 2GBX > 2XSH > 2Y FJ > 2YFL > 3GZX > 3GZY > 2E4P. For heptachlorobiphenyl, the degradation rates of five types of BphA increased in the following order: 2YFJ < 2GBX < 2XSH < 2E4P < 3 GZY. The degradation rate for octachlorobiphenyl was higher with 3GZY than with 2E4P. Our results correspond with those from previous studies that found that low binding free energies were an indicator of high affinity [ 20 ]. Chen et  al. [ 45 ] studied the application of PCB industrial products in practical production. Their results showed that isomers from trichlorobiphenyl to heptachlorobiphenyl accounted for the majority of the commercial mixtures of PCBs, and chlorobiphenyl and dichlorobiphenyl also contributed a substantial proportion [ 45 ]. In the environment, the isomers from trichlorobiphenyl to heptachlorobiphenyl are the primary PCBs that undergo degradation. In this study, 2GBX, 2XSH, 2YFJ, 2E4P, and 3GZY had the ability to degrade the isomers from trichlorodiphenyl to heptachlorobiphenyl, and the PCB degradation rates were in the following order: 2GBX  >  2YFJ  >  2XSH  >  3GZY  >  2E4P. For single PCBs, 2GBX had the highest degradation rates for the isomers from chlorobiphenyl to pentachlorobiphenyl, 2XSH had the best ability to degrade hexachlorobiphenyl, and 3GZY had the highest degradation rates for isomers from heptachlorobiphenyl to nonachlorobiphenyl. In the microbial degradation system for PCBs, the biodegradation ratios of the PCBs were maximized with addition of 2GBX, 2XSH, and 3GZY. Correlation analysis of the relationships between the binding free energies and the persistent organic pollutant characteristics of the PCBs In this study, Pearson correlation analysis was used to study the relationships between the binding free energies and the molecular weight, KOA, BCF, t1/2, and IC50 values for the PCBs (Table 1). The Pearson correlation analysis revealed that the binding free energies were significantly correlated to the molecular weight, KOA, BCF, and t1/2 values (P = 0 < 0.01), but not significantly correlated to the IC50 values (P < 0.05). The binding free energies between the PCBs and BphA gradually increased as the number of Cl atoms increased (Fig. 2), that is, the binding free energy gradually increased with increasing molecular weight. This indicted a positive correlation between the binding free energy and molecular weight. Large KOA values for the PCBs imply they have stronger migration abilities [ 45 ]. Strong migration abilities may suggest that PCBs are not easily degraded in the environment because the atmospheric conditions are not suitable for biodegradation. Because the binding free energy was negatively correlated with the degradation ability, low degradation ability implied the binding free energy was strong. Therefore, the binding free energy was positively correlated with the KOA values of the PCBs. Large BCF values implied that the PCBs strongly bioaccumulated in the body [ 44 ]. With gradual accumulation in the food chain, the levels of PCBs in the body will gradually increase, and this will be accompanied by a decrease in the amount of PCBs in the environment. This is not beneficial for PCB degradation in the environment by microorganisms. Hence, PCBs are more difficult to degrade following an increase in the binding free energy, which shows that the binding free energy is positively correlated with the BCF. Large t1/2 values for the PCBs imply that they are persistent in “r” represents “Pearson correlation coefficient, the bigger the r, the stronger the correlation”; “sig.” means “significant, if sig < 0.05, indicating the significant correlation”; “N” represents “sample size”;“**” means “the coefficient statistically significant was at P = 0.01” (n = 209) N, optimum number of components; r2, correlation coefficient; q2, cross-validated val; SEE, standard error of the estimate; F, Fischer’s test value; S, steric; E, electrostatic; H, hydrophobic; D, hydrogen bond donor; A, hydrogen bond acceptor the environment [ 45 ]. Hence, PCBs are more difficult to degrade following an increase in the binding free energy, which agrees with earlier studies. CoMFA and CoMSIA contour analysis The q2 from the CoMFA models had a range of 0.585– 0.908 (>  0.5), the r2 range was 0.717–0.948 (>  0.6) [ 49 ], the SEE range was 0.357–1.380, and the F range was 32.668–269.003 (Table  2). Therefore, the CoMFA model was robust. For the CoMSIA models, the q2 range was 0.521–0.897 (> 0.5), the r2 range was 0.630–0.930 (> 0.6) [ 49 ], the SEE range was 0.394–1.750, and the F range was 20.393–197.637. Therefore, the CoMFA model was robust. The contribution rates of the electrostatic fields to the CoMFA and CoMSIA models were in the range of 0.609–0.848, which showed that the electrostatic field was important in each of the two models. Chen et  al. [ 45 ] found that the electrostatic field played an important role in determining the logKOA values of PCBs, Liu et  al. [ 44 ] determined that the electrostatic descriptor was a primary factor governing the logBCF, and Xu et al. [ 45 ] found that electrostatic interaction was crucial in determining the t1/2 values. These results are consistent with the positive correlations between the binding free energies and the KOA, BCF, and t1/2 values for the PCBs, which can be indirectly explained by the correlation between the binding free energies and the KOA, BCF, and t1/2 values for the PCBs from the molecular field perspective. The binding free energies between the nine types of BphA and 209 PCBs are shown in Table S1 (shown in the Additional file 1 available to the readers). The most interesting feature of CoMFA and CoMSIA models is the visualization of the results as 3D contour plots. We can observe the contour of the 3D space around the molecule, in which changes in the physicochemical properties are predicted to increase or decrease the effectiveness of the contour [ 50 ]. The results for the CoMFA and CoMSIA models are illustrated in field contribution graphs using the standard deviation times coefficient field. The contour maps of the two models (Figs. 3, 4) show the electrostatic fields with decachlorobiphenyl as a template. These contour maps showed 80 and 20% contributions for favorable and unfavorable regions, respectively. The electrostatic field is represented by blue and red contours, with blue regions indicating the electropositive groups close to these regions are favorable to the binding free energy, and the red regions indicate the electronegative groups near these regions can enhance the binding free energy. By contrast, the binding free energy decreased when electropositive groups were introduced in the red region or electronegative groups in the blue region. Although the electrostatic field was identified as a major factor affecting the binding free energies between BphA and the PCBs, there were some differences in the electrostatic contour maps (Figs. 3, 4, Table 3). The binding free energies between the PCBs and the nine types of BphA increased simultaneously with substituents possessing electropositive groups at the 3′-position of the B ring or electronegative groups at the 2- and 3-positions of the A ring. In the electrostatic fields of 1WQL, 2YFJ, 2YFL, 2XSH, and 2E4P, electropositive groups were desirable at the 3′- and 5-positions to increase the binding free energies. In the electrostatic fields of 2GBX and 3GZX, the electropositive groups close to the regions of the 2- and 3′-positions increased the binding free energies between the PCBs and BphA. Introducing electropositive groups at the 2-, 3′-, and 5-positions increased the binding free energies between the PCBs and 1ULJ. The binding free energies between the PCBs and 3GZY increased with substituents possessing electropositive groups at the 3′-position of the B ring. In the electrostatic fields of 1ULJ and 3GZX, introducing electronegative groups to the 2-, 3-, and 6-positions of the A ring increased the binding free energies. The binding free energies between PCB-24 and PCB-5 could be deciphered using this contour. In the electrostatic fields of 2YFJ and 2GBX, electronegative groups near the 2-, 3-, 3′-, 4-, 6-, and 6′-positions were favorable for enhancing the binding free energies, as shown by the higher binding free energy for PCB-109 than for PCB-55 or PCB-59. The binding free energies between the PCBs and 1WQL were enhanced with electronegative groups at the 2-, 3-, and 5-positions of the A ring. For example, the binding free energy of PCB-23 was higher than those of PCB-6 and PCB-9. Electronegative groups were desirable at the 2-, 3-, and 4-positions of the A ring and increased the binding free energies between the PCBs and 3GZY, as shown by the higher binding free energy for PCB-21 than for PCB-7 or PCB-12. Introducing electronegative groups to the 2- and 3-positions of the A ring increased the binding free energies between the PCBs and 2XSH. The same applied for the binding free energies between PCB-5 and PCB-2. The binding free energies between the PCBs and 2E4P increased on introduction of substituents with electronegative groups at the 2-, 3-, 6-, and 6′-positions. An example of this was the high binding free energy of PCB-5 compared with those of PCB-1 and PCB-2. Electronegative groups close to the regions of the 2-, 4-, 5-, 6-, and 6′-positions increased the binding free energies between the PCBs and 2YFJ, as shown by the higher binding free energy of PCB-30 than for PCB-7 or PCB-10. The logKOA increased on introduction of substituents with electropositive groups at the 3′- and 6-positions or electronegative groups at the 2-, 3-, 3′- and 5-positions [ 45 ]. The logKOA and the binding free energies between the PCBs and the nine types of BphA increased simultaneously on introduction of substituents with electropositive groups at the 3′-position of the B ring or electronegative groups at the 2-position of the A ring. In the electrostatic fields of 1ULJ, 2XSH, 2E4P, 3GZX, 3GZY, 1WQL, 2YFL, 2GBX, and KOA, introducing electronegative groups at the 3-position of the A ring increased the binding free energies and the KOA simultaneously. In the electrostatic fields of 2YFJ, 1WQL, and KOA, electronegative groups near the 5-position were favorable for enhancing the binding free energies and the KOA simultaneously. Therefore, using molecular modeling, we proved that the binding free energies were positively correlated with the migration abilities of the PCBs. The BCF values increased on introduction of substituents with electropositive groups at the 2-, 3′-, 5-, and 6-positions or electronegative groups at the 3-, 4-, and 5-positions of the A ring [ 44 ]. The BCF values of the PCBs and the binding free energies between the PCBs and BphA increased simultaneously on introduction of substituents with electropositive groups at the 3′-position of the B ring or electronegative groups at the 3-position of the A ring. In the electrostatic fields of 1ULJ, 3GZX, 3GZY, 2GBX, and BCF, electropositive groups were desirable at the 3′-positions of the B ring for increases in the binding free energies and the BCF values. In the electrostatic fields of 1WQL, 2YFJ, 2YFL, 2XSH, 2E4P and BCF, introducing electropositive groups to the 3′- and 5-positions increased the binding free energies and BCF values simultaneously. In the electrostatic fields of 1ULJ, 3GZX, 2XSH, 2E4P, and BCF, electronegative groups close to the regions of the 3-position of the A ring increased the binding free energies and the BCF values. In the electrostatic fields of 2YFJ, 2GBX, 3GZY, and BCF, electronegative groups near the 3- and 4-positions of the A ring were favorable for enhancing the binding free energies and the BCF values simultaneously. The introduction of electronegative groups to the 3- and 5-positions of the A ring increased the BCF values and the binding free energies between PCBs and 1WQL. Electronegative groups were desirable at the 4- and 5-positions of the A ring to increase the BCF values of the PCBs and the binding free energies between the PCBs and 2YFJ. Therefore, the binding free energies were positively correlated with the bioaccumulation of PCBs, which was proven by molecular modeling. It was determined that electronegative groups close to the regions of the 3-, 3′-, 5-, 6-, and 6′-positions increased the t1/2 values of the PCBs [ 45 ]. The t1/2 values of the PCBs and the binding free energies between the PCBs and the nine types of BphA increased simultaneously on introduction of substituents with electronegative groups at the 3-position of the A ring. In the electrostatic fields of 2YFL, 2GBX, and t1/2, electronegative groups were desirable at 3-, 3′-, 6-, 6′-positions to increase the binding free energies and t1/2 values. Electronegative groups close to the regions of the 3- and 5-positions of the A ring increased the t1/2 values and binding free energies between the PCBs and 1WQL. Electronegative groups near the 3-, 5-, 6-, and 6′-positions were favorable for simultaneously enhancing the t1/2 values and the binding free energies between PCBs and 2YFJ. Introducing electronegative groups at the 3- and 6′-positions increased the t1/2 values and the binding free energies between the PCBs and 2E4P. Electronegative groups were desirable at the 3-, 6-, and 6′-positions to increase the t1/2 values of the PCBs and the binding free energies between the PCBs and 3GZX. Molecular modeling showed that the binding free energies were positively correlated with the persistence of the PCBs. Conclusions The main research conclusions are as follows: 1. Among the nine types of BphA, 2GBX gives the highest PCB degradation rate. 2. The binding free energy is highly significantly correlated with the molecular weight, KOA, BCF, and t1/2 values, but is not significantly correlated with the IC50 values. 3. The electrostatic descriptors play a more significant role than steric descriptors and hydrophobic descriptors. This is consistent with literature findings of the electrostatic field as a major factor governing the KOA, t1/2, and BCF values. Therefore, that the influence of molecular structure on the binding free energies, KOA, t1/2, and BCF values is consistent. 4. In an electrostatic field, introducing electropositive groups in the 3-position or electronegative groups in the 3′-position of a PCB significantly reduces the binding free energies, molecular weights, KOA, BCF, and t1/2 values. Therefore, the binding free energies are correlated with the molecular weights, KOA, BCF, and t1/2 values of the PCBs because of structural features of the molecules. Additional file Additional file 1: Table S1. The binding free energies between nine types of BphA and 209 PCBs. Authors’ contributions All authors read and approved the final manuscript. Acknowledgements Not applicable. Competing interests Not applicable. Availability of data and materials Not applicable. Ethics approval and consent to participate Not applicable. Funding Fundamental Research Funds for the Central Universities in 2013 (JB2013146) and the Key Projects in the National Science & Technology Pillar Program in the Eleventh 5‑ Year Plan Period (2008BAC43B01). Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in pub‑ lished maps and institutional affiliations. 1. Lallas PL ( 2001 ) The stockholm convention on persistent organic pollutants . Am J Int Law 95 : 692 - 708 2. Birgül A , KurtKarakus PB , Alegria H , Gungormus E , Celik H , Cicek T , Güven EC ( 2017 ) Polyurethane foam (PUF) disk passive samplers derived polychlorinated biphenyls (PCBs) concentrations in the ambient air of Bursa‑ Turkey: spatial and temporal variations and health risk assessment . Chemosphere 168 : 1345 - 1355 3. Barakat AO , Mostafa A , Wade TL , Sweet ST , EI Sayed NB ( 2013 ) Distribution and ecological risk of organochlorine pesticides and polychlorinated biphenyls in sediments from the Mediterranean coastal environment of Egypt . Chemosphere 93 : 545 - 554 4. Saba T , Su S ( 2013 ) Tracking polychlorinated biphenyls (PCBs) congener patterns in Newark Bay surface sediment using principal component analysis (PCA) and positive matrix factorization (PMF) . J Hazard Mater 260 : 634 - 643 5. Frederiksen M , Meyer HW , Ebbehøj NE , Gunnarsen L ( 2012 ) Polychlorinated biphenyls (PCBs) in indoor air originating from sealants in contaminated and uncontaminated apartments within the same housing estate . Chemosphere 89 : 473 - 479 6. DellaValle CT , Wheeler DC , Deziel NC , De Roos AJ , Cerhan JR , Cozen W , Severson RK , Flory AR , Locke SJ , Colt JS , Hartge P , Ward MH ( 2013 ) Environmental determinants of polychlorinated biphenyl concentrations in residential carpet dust . Environ Sci Technol 47 : 10405 - 10414 7. Rawn DFK , Sadler AR , Quade SC , Sun WF , Kosarac I , Hayward S , Ryan JJ ( 2012 ) The impact of production type and region on polychlorinated biphenyl (PCB), polychlorinated dibenzo‑p ‑ dioxin and dibenzofuran (PCDD/F) concentrations in Canadian chicken egg yolks . Chemosphere 89 : 929 - 935 8. Su GY , Liu XH , Gao ZS , Xian QM , Feng JF , Zhang XW , Giesy JP , Wei S , Liu HL , Yu HX ( 2012 ) Dietary intake of polybrominated diphenyl ethers (PBDEs) and polychlorinated biphenyls (PCBs) from fish and meat by residents of Nanjing, China . Environ Int 42 : 138 - 143 9. Hassine SB , Ameur WB , Gandoura N , Driss MR ( 2012 ) Determination of chlorinated pesticides, polychlorinated biphenyls, and polybrominated diphenyl ethers in human milk from Bizerte (Tunisia ) in 2010. Chemosphere 89 : 369 - 377 10. Shen HT , Ding GQ , Wu YN , Pan GS , Zhou XP , Han JL , Li JG , Wen S ( 2012 ) Polychlorinated dibenzo‑p ‑ dioxins/furans (PCDD/Fs), polychlorinated biphenyls (PCBs), and polybrominated diphenyl ethers (PBDEs) in breast milk from Zhejiang, China . Environ Int 42 : 84 - 90 11. Jotaki T , Fukata H , Mori C ( 2011 ) Confirmation of polychlorinated biphenyl (PCB) distribution in the blood and verification of simple quantitative method for PCBs based on specific congeners . Chemosphere 82 : 107 - 113 12. Arrebola JP , Fernandez MF , Porta M , Rosell J , Ossa RMDL , Olea N , MartinOlmedo P ( 2010 ) Multivariate models to predict human adipose tissue PCB concentrations in Southern Spain . Environ Int 36 : 705 - 713 13. Brazova T , Hanzelova V , Miklisova D ( 2012 ) Bioaccumulation of six PCB indicator congeners in a heavily polluted water reservoir in Eastern Slovakia: tissue specific distribution in fish and their parasites . Parasitol Res 111 : 779 - 786 14. Erickson BD , Mondello FJ ( 1992 ) Nucleotide sequencing and transcriptional mapping of the genes encoding biphenyl dioxygenase, a multicomponent polychlorinated‑biphenyl‑ degrading enzyme in Pseudomonas strain LB400 . J Bacteriol 174 : 2903 - 2912 15. Kitagawa W , Miyauchi K , Masai E , Fukuda M ( 2001 ) Cloning and characterization of benzoate catabolic genes in the gram‑positive polychlorinated biphenyl degrader Rhodococcus sp . strain RHA1. J Bacteriol 183 : 6598 - 6606 16. Bedard DL , Haberl ML , May RJ , Brennan MJ ( 1987 ) Evidence for novel mechanisms of polychlorinated biphenyl metabolism in Alcaligenes eutrophus H850 . Appl Environ Microbiol 53 : 1103 - 1112 17. Jia LY , Jia LY , Zheng AP , Xu L , Huang XD , Zhang Q , Yang FL ( 2008 ) Isolation and characterization of comprehensive polychlorinated biphenyl degrading bacterium, Enterobacter sp . LY402. J Microbiol Biotechnol 18 : 952 - 957 18. Furukawa K , Suenaga H , Goto M ( 2004 ) Biphenyl dioxygenases: functional versatilities and directed evolution . J Bacteriol 186 : 5189 - 5196 19. Cao YM , Xu L , Jia LY ( 2011 ) Analysis of PCBs degradation abilities of biphenyl dioxygenase derived from Enterobacter sp . LY402 by molecular simulation. New Biotechnol 29 : 90 - 98 20. Wu B , Zhang Y , Kong J , Zhang XX , Cheng SP ( 2009 ) In silico predication of nuclear hormone receptors for organic pollutants by homology modeling and molecular docking . Toxicol Lett 191 : 69 - 73 21. Liu ZQ ( 2012 ) Expression of biphenyl dioxygense and the binding properties with substrates . Dalian university of Technology, Da Lian 22. Pieper DH ( 2005 ) Aerobic degradation of polychlorinated biphenyls . Appl Microbiol Biotechnol 67 : 170 - 191 23. Kumamaru T , Suenaga H , Mitsuo M , Furukawa K ( 1998 ) Enhanced degradation of polychlorinated biphenyls by directed evolution of biphenyl dioxygenase . Nat Biotechnol 16 : 663 - 666 24. Bruhlmann F , Chen W ( 1999 ) Tuning biphenyl dioxygenase for extended substrate specificity . Biotechnol Bioeng 63 : 544 - 551 25. Suenaga H , Goto M , Furukawa K ( 2006 ) Active‑site engineering of biphenyl dioxygenase: effect of substituted amino acids on substrate specificity and regiospecificity . Appl Microbiol Biotechnol 71 : 168 - 176 26. Gu WW , Chen Y , Zhang L , Li Y ( 2016 ) Prediction of octanol‑ water partition coefficient for polychlorinated naphthalenes through three‑ dimensional QSAR models . Hum Ecol Risk Assess 1-16 27. Furusawa Y , Nagarajan V , Tanokura M , Masai E , Fukuda F , Senda T ( 2004 ) Crystal structure of the terminal oxygenase component of biphenyl dioxygenase derived from Rhodococcus sp. strain RHA1 . J Mol Biol 342 : 1041 - 1052 28. Dong XS , Fushinobu S , Fukuda E , Terada T , Nakamura S , Shimizu K , Nojiri H , Omori T , Shoun H , Wakagi T ( 2005 ) Crystal structure of the terminal oxygenase component of cumene dioxygenase from Pseudomonas fluorescens IP01 . J Bacteriol 187 : 2483 - 2490 29. Mohammadi M , Viger JF , Kumar P , Barriault D , Bolin JT , Sylvestre M ( 2011 ) Retuning rieske‑type oxygenases to expand substrate range . J Biol Chem 286 : 27612 - 27621 30. Kumar P , Mohammadi M , Dhindwal S , Pham TTM , Bolin JT , Sylvestre M ( 2012 ) Structural insights into the metabolism of 2‑ chlorodibenzofuran by an evolved biphenyl dioxygenase . Biochem Bioph Res Co 421 : 757 - 762 31. Ferraro DJ , Brown EN , Yu CL , Parales RE , Gibson DT , Ramaswamy S ( 2007 ) Structural investigations of the ferredoxin and terminal oxygenase components of the biphenyl 2,3‑ dioxygenase from Sphingobium yanoikuyae B1 . BMC Struct Biol 7 : 1 - 14 32. Kumar P , Mohammadi M , Viger JF , Barriault D , Leticia GG , Lindsay DE , Jeffrey TB , Michel S ( 2011 ) Structural insight into the expanded PCB‑ degrading abilities of a biphenyl dioxygenase obtained by directed evolution . J Mol Biol 405 : 531 - 547 33. Senda M , Kishigami S , Kimura S , Fukuda M , Ishida T , Senda T ( 2007 ) Molecular mechanism of the redox‑ dependent interaction between NADH‑ dependent ferredoxin reductase and rieske‑type [2Fe ‑2S] ferre ‑ doxin . J Mol Biol 373 : 382 - 400 34. Christopher LC , Nathalie YRA , Pravindra K , Mathew NC , Sangita CS , Justin BP , Lindsay DE , Jeffrey TB ( 2013 ) Structural characterization of pandoraeapnomenusa B‑356 biphenyl dioxygenase reveals features of potent polychlorinated biphenyl‑ degrading enzymes . PLoS ONE . https://doi. org/10.1371/journal.pone.0052550 35. Halgren TA ( 1996 ) Merck molecular force field.1. Basis, form, scope, parameterization, and performance of MMFF94 . J Comput Chem 17 : 490 - 519 36. Li XL , Ye L , Wang XX , Wang XZ , Liu HL , Zhu YL , Yu HX ( 2012 ) Combined 3D‑ QSAR, molecular docking and molecular dynamics study on thyroidhormone activity of hydroxylated polybrominated diphenyl ethers to thyroidreceptors β . Toxicol Appl Pharm 265 : 300 - 307 37. Zhang J , Shan Y , Pan X , Wang C , Xu W , He L ( 2011 ) Molecular docking, 3D‑ QSAR studies, and in silico ADME prediction of p‑aminosalicylic acid derivatives as neuraminidase inhibitors . Chem Biol Drug Des 78 : 709 - 717 38. Holt PA , Chaires JB , Trent JO ( 2008 ) Molecular docking of intercalators and groove‑binders to nucleic acids using Autodock and Surflex . J Chem Inf Model 48 : 1602 - 1615 39. Ewing TJA , Makino S , Skillman AG et al ( 2001 ) DOCK 4.0: research strategies for automated molecular docking of flexible molecule database . J Comput Aided Mol Des 15 : 411 - 428 40. Li WL , Si HZ , Li Y , Ge CZ , Song FC , Ma XT , Duan YB , Zhai HL ( 2016 ) 3D‑ QSAR and molecular docking studies on designing inhibitors of thehepatitis C virus NS5B polymerase . J Mol Struct 1117 : 227 - 239 41. Jain AN ( 2007 ) Surflex‑Dock 2.1: robust performance from ligand ener ‑ getic modeling, ring flexibility, and knowledge‑based search . J Comput Aided Mol Des 21 : 281 - 306 42. Ai Y , Wang ST , Tang C , Sun PH , Song FJ ( 2011 ) 3D‑ QSAR and docking studies on pyridopyrazinones as BRAF inhibitors . Med Chem Res 20 : 1298 - 1317 43. Yang W , Shen S , Mu L , Yu H ( 2011 ) Structure-activity relationship study on the binding of PBDEs with thyroxine transport proteins . Environ Toxicol Chem 30 : 2431 - 2439 44. Wang XL , Gu WW , Guo EM , Cui CY , Li Y ( 2017 ) Assessment of long‑range transport potential of polychlorinated Naphthalenes based on threedimensional QSAR models . Environ Sci Pollut Res 24 : 14802 - 14818 45. Chen Y , Cai XY , Jiang L , Li Y ( 2016 ) Prediction of octanol‑air partition coefficients for polychlorinatedbiphenyls (PCBs) using 3D‑ QSAR models . Ecotoxcol Environ Safe 124 : 202 - 212 46. Liu H , Liu HL , Sun P , Wang ZY ( 2014 ) QSAR studies of bioconcentration factors of polychlorinated biphenyls (PCBs) using DFT. PCS and CoMFA . Chemosphere 114 : 101 - 105 47. Xu Z , Chen Y , Qiu YL , Gu WW , Li Y ( 2016 ) Prediction of stability for polychlorinated biphenyls (PCBs) in transformer insulation oil through 3DQSAR pharmacophore model assisted with full factor experimental design . Chem Res Chin Univ 32 : 348 - 356 48. Li XL , Ye L , Wang XX , Shi W , Liu HL , Qian XP , Zhu YL , Yu HX ( 2013 ) In silico investigations of anti‑androgen activity of polychlorinated biphenyls . Chemosphere 92 : 795 - 802 49. Krishna KM , Inturi B , Gurubasavaraj VP , Madhusudan NP , Vijaykumar GS ( 2014 ) Design, synthesis and 3D‑ QSAR studies of new diphenylamine containing 1,2,4‑triazoles as potential antitubercular agents . Eur J Med Chem 84 : 516 - 529 50. Mao YT , Li Y , Hao M , Zhang SW , Ai CZ ( 2012 ) Docking, molecular dynamics and quantitative structure‑activity relationship studies for HEPTs and DABOs as HIV‑1 reverse transcriptase inhibitors . J Mol Model 18 : 2185 - 2198