Modeling Based Structural Insights into Biodegradation of the Herbicide Diuron by Laccase-1 from Ceriporiopsis subvermispora

The herbicide diuron (3-(3,4-dichlorophenyl)-1,1-dimethylurea) is used in many agricultural crops and non-crop areas worldwide, leading to the pollution of the aquatic environment by soil leaching. White rot fungi and its lignin modifying enzymes, peroxidases and laccases, are responsible for its degradation. Therefore, it is of interest to explore the potential use of Ceriporiopsis subvermispora laccase (CersuLac1) in the biotransformation of this herbicide by using its enzyme laccase. However, the structure of laccase from Ceriporiopsis subvermispora is still unknown. Hence, a model of laccase was constructed using homology modeling. The model was further used to dock p-methylbenzoate in the presence of four copper ions to analyze molecular basis of its binding and interaction. The ligand-protein interaction is stereo-chemically favorable in nature. The presence of the single protonated Lys457 was necessary for catalysis, being coordinated by a cupper ion. The best pose of diuron on CersuLac1 has a theoretical Ki of 2.91 mM. This is comparable to the KM values for laccases from other organisms with similar compounds. Thus, we document the insights for the potential use of laccase from Ceriporiopsis subvermispora in the biotransfrormation of diuron.


Background:
The presence of pesticides in the environment is a matter of particular concern for the conservation of ecosystems and for human health. The natural process of transformation of such substances in the environment, as well as their elimination is called bioremediation. Therefore, the understanding of the biochemical reactions involved in their metabolism is the basis for identifying the time of persistence of such compounds in nature. Diuron (IUPAC name: 3-(3,4-dichlorophenyl)-1,1dimethylurea; CAS number: 330-54-1) is a systemic herbicide, largely used in agriculture, belonging to the phenylamide family and the subclass of phenylurea. This substituted urea herbicide inhibits photosynthesis by preventing oxygen production [1] and blocks the electron transfer at the level of photosystem II of photosynthetic microorganisms and plants [2].
Diuron is classified as moderately toxic to the environment and contaminant of groundwater [2]. The time for its degradation in the soil depends on biological and nonbiological factors that may interact with it. In the absence of biological processes, the degradation rates vary between 55 to 70% in a period of three weeks. With the participation of microorganisms, this rate increases between 80 to 85%. Under anaerobic conditions, the rate reaches 85% after 18 days [3]. Diuron hydrolyses very slowly in neutral solution at 25 o C. The abiotic degradation in water solution is an irreversible reaction giving 3,4-DCA (3,4-dichloroaniline) as the only product. 3,4-DCA is considerably more toxic than diuron itself and has a higher water solubility, so it can leach out from treated agricultural land. It cans eventually condensate into chloroazobenzene, which is highly toxic [4].
Biodegradation is a major way of diuron decomposition in the environment. Many studies were conducted and they showed that degradation happens under both aerobic and anaerobic conditions and it is promoted by both Gram positive and Gram negative bacteria as well as by some fungi [2]. The capability of ligninolytic or white-rot fungi to degrade diuron has also been considered. White rot fungi (WRF) produce several extracellular lignin modifying enzymes that include lignin peroxidase, manganese peroxidase and laccase [5]. The degradation of diuron by WRF is partially attributed to the action of these enzymes [6].
Ceriporiopsis subvermispora is a well-known WRF used in the bio-pulping industry [7], but its role in biodegradation and potential application in bioremediation remains untapped. Therefore, this work aimed to investigate the potential application of the fungus C. subvermispora in degrading the compound diuron by means of molecular docking simulations, using the modeled structure of enzyme laccase-1 produced by this fungus as a target. In the field of molecular modeling, docking is a prediction method to find the preferred orientation of one molecule to a target when bound to each other to form a stable complex. The knowledge of the preferred orientation may in turn be used to predict the strength of association, or binding affinity, between two molecules by means of scoring functions [8]. These studies are important for characterizing the structure of this enzyme, which has not yet been solved by experimental methods of X-ray diffraction, as well as for a better understanding of the reaction mechanism, which is still not fully understood.

Sequence analysis
The sequence from laccase-1 of Ceriporiopsis subvermispora (CersuLac1) was obtained from UniProt database (UNIPROT: O59944) and then subjected to InterProScan-5 and TargetP-1.1 servers in order to characterize the structural domains and predict their subcellular localization. After, the amino acid sequence on the functional protein domain was used to search for structural templates by BLASTp against Protein Data Bank. The presence of disulfide bonds was evaluated by Ensemblebased Connectivity Disulfide Bonding Pattern prediction server (EDBCP) and the prediction of secondary structure was performed by PsiPred server.

Model building and validation
The tertiary structure of CersuLac1 bonded to the pmethylbenzoate (ZINC: 330134) and four copper ions was obtained by program Modeller-9.11. First, the amino acid sequence of the target was aligned with the sequence of the template PDB: 2HRG to generate 1,000 models of the enzyme by using restrictions to make two disulfide bonds between residues 87 to 487 and 119 to 207. The top five models ranked by Modeller Dope Score were evaluated by their stereochemical quality through software Procheck. The best model was used to generate other 1,000 new models for corrections of residues positioned at unfavorable regions of Ramachandram plot, and the final model was chosen by applying again the evaluation method described above. For the docking studies, the final model was fully minimized with 20,000 steps of conjugate gradient through NAMD2-2.9 software. For this purpose, the protein was immersed in a periodic water box with limits at least 10 Å away from the outer surface of the protein, and also with sufficient sodium counter ions for the neutralization of the system charges. The force field Charmm-c35b2 was used to define the residues of the protein and the force field of the ligand was generated, in the same format, through Swissparam server.

Figure 1:
ClustalO sequence alignment of CersuLac1 with templates 2HRG from Trametes trogii and 3KW7 from Trametes sp. The first line on the alignment corresponds to the prediction of secondary structure by PsiPred server, having the alpha-helix in blue boxes and beta-strands in green boxes. Yellow boxes indicate the regions of CersuLac1 with no match in templates. Residues in red correspond to the glycosylation sites, and residues in blue correspond to copper-binding site as Uniprot notation for template 3KW7. Asterisks represent identical residues while points, and two points represent similar residues.

Molecular docking
To establish the optimal docking simulation protocol, the redocking of the ligand p-methylbenzoate on the proteinligand complex of the template 2HRG was performed, because it is assumed that the crystallized final structure of the complex represents the lowest energy conformation and therefore the most stable [9]. In this procedure, the ligand was taken from its original structure and redocked by means of the software AutoDock-4.2.3 in order to find the conditions that could reproduce the pose of the modeled/crystallographic ligand with a root mean square deviation (rmsd) lower than 2.0 Å, as well as 100% reproducibility. The docking protocol was established by keeping the single bonds of the ligands free rotatable, the enzyme treated as a rigid body, and applying the Lamarckian genetic algorithm for search, based at the already known active site location in the template structure. This protocol was then applied for the redocking of pmethylbenzoate into the minimized structure of CersuLac1 and in the docking of the compound diuron (ZINC57287) obtained from Zinc database in the template and minimized CersuLac1 structures.

Discussion:
Analysis of the structure The analysis of the sequence of CersuLac1 suggested the presence of a signal peptide between residues 1 and 21, which signals to the secretory pathway, suggesting an extracellular enzymatic activity and showed 70 and 67% of identity with laccases from the fungi Trametes sp, (PDBid: 3KW7), and from Coriolopsis trogii, (PDBid: 2HRG) respectively. The template 3KW7 was not used in the modeling, only in the comparative analysis of results due its low resolution. The template 2HRG was chosen for the modeling because it has high sequence identity and high resolution (1.58 Å, Rvalue=0.175), and it was co-crystallized in the presence of the ligand p-methylbenzoate, a compound found in industrial wastes [10] with shares certain structural similarity to diuron (Tamimoto index: 0.11).
The modeling of a protein in the presence of a ligand has some advantages. It facilitates the location of the binding site, and the binder itself, gives the proper guidance in the selection of the best pose of an unknown ligand. There were three regions in the target sequence without match the templates (yellow boxes at Figure 1). Without this structural information, the modeling program tends to include them in the model as bulges. According to PsiPred server, these three unmatched regions in the templates would not be part of any element of secondary structure in CersuLac1. As a result, no restriction for modeling secondary structure was applied during the construction of the models. Nevertheless, these three unmatched regions in model sequence are placed in loops, away from the ligand binding site, without compromising the stability of the structure. However, it was observed that in the final model there were two regions with deletions in CersuLac1 (162 to 163) and (392 to 393), which should be in direct contact with the ligand p-methylbenzoate. This fact led to a significant change in the area of catalytic site that was reflected in a different affinity of the CersuLac1 by the ligands compared to the template. Based on the analysis of the Procheck software, the final CersuLac1 model presented better stereochemical quality than their templates, having all residues placed in the allowed regions of Ramachandram plot. The residue Asp205 from 2HRG template structure is located at the p-methylbenzoate binding site. According to the Procheck evaluation, this residue is located in a generously allowed region of the Ramachandram plot, which could compromise the redocking simulations of the p-methylbenzoate in this structure.
The final modeled structure of CersuLac1 bonded to the pmethylbenzoate and four copper ions is shown in Figure 2 and, its analysis by PDBsum server showed that its tertiary structure is not very different from the templates used in its construction. As a result, we can infer a high confidence level that the same CATH structural classification identified in the template must also be present in the CersuLac1 model. According to CATH server domain, the regions between residues 1 to 134, 135 to 310 and 311 to 498 would form three

Molecular docking
When a protein is crystallized in the presence of ligand, both the macromolecule and ligand assume a position that maximizes their intermolecular interactions. However, in a structure that was modeled in the presence of a ligand, this mutual interaction is not taken into account, since modeling programs are based on semi-empirical force fields to establish such interactions and parameters for small ligands are unavailable. As a result, the structure of the ligand is copied as it appears in the template structure. Energy minimization of the whole complex provides a way to strengthen the interactions between protein and the ligand in modeled complexes. Since accurate force fields for the ligands are parameterized, their use in minimization softwares can cause the complex interactions to be comparable to those found in a crystal structure.
Taking this fact into consideration, the docking protocol was validated by redocking the compound p-methylbenzoate in the structure of template 2HRG and in the minimized structures of CersuLac1. In both targets, the simulations were repeated ten times in order to verify the reproducibility of the data, thus avoiding false-positive results. The protocol applied uses the software AutoDock-4.2.3, Lamarckian GA, box size 15, 15, 35 and 0.375 grid spacing around the center of the ligand. The results obtained with redocking of p-methylbenzoate in 2HRG showed a good overlap of the ligand in the first three poses of lower energy (Figure 3A). All ten attempts provide an average ΔGbinding of -2.46 kcal.mol -1 (± 0.009), which corresponds to a theoretical Ki of 15.78 mM (± 0.218) and an rmsd of 1.57 Å (± 0.013). In the redocking of p-methylbenzoate in the CersuLac1, the first four poses of all ten attempts appeared to be well overlapped ( Figure 3B). The best pose had an average ΔGbinding of -2.08 kcal.mol -1 (± 0.013), which corresponds to a theoretical Ki of 30.0 mM (± 0.73) and an rmsd of 1.71 Å (± 0.034).
These results suggest that the 2HRG must have twice as much the affinity by the p-methylbenzoate as the CersuLac1. Anyway, both Ki values are relatively high, which indicates that this product should not be a good ligand for the laccase. This may justify why there were no references in the literature regarding experimental data for the reaction of laccase with pmethylbenzoate. Those rmsd values, although smaller than 2.0 Å, are relatively high, which would be expected in the docking of a small ligand within a wide cavity, providing a greater freedom of movement. This observation is compatible with the highest value of thermal vibration factor (Bfactor) of the binder in the crystallographic template 2HRG, which is 33.53 Å 2 , being relatively high compared to the mean from the whole structure, 19.70 Å 2 . Regarding the docking of the compound diuron in the 2HRG template structure, an average ΔGbinding of -4.53 kcal.mol -1 (± 0.027) corresponding to a theoretical Ki of 479.2 μM (± 21.11) was obtained. The docking of this compound in CersuLac1 structure provides an average ΔGbinding of -3.46 kcal.mol -1 (± 0.008), which corresponds to a theoretical Ki of 2.91 mM (± 0.043). The images of the best poses of diuron in the template and model structures are shown in Figure 4. The best pose in the docking of diuron on CersuLac1 is comparable to the position of the ligand p-methylbenzoate. This pose is also catalytic favorable for the oxidation of diuron [11], where His457 plays an important catalytic role, since the presence of copper ions causes the protonation of this residue to be moved from Nδ1 to the Nδ2 in order to make a hydrogen bond with the dimethylurea oxygen of diuron ( Figure 5). According to Claus [12], the various copper centers of laccase oxidize the substrate by transferring their electrons to molecular oxygen without the release of toxic peroxide intermediates. This is achieved through four monoeletronic oxidations of the substrate catalyzed by the type 1 copper (Cu1 which interacts with the residue His457 of modeled structure). The electrons would be transferred to a trinuclear copper cluster (Cu2, Cu3 and Cu4), where the reduction of molecular oxygen and the release of water occur.
The enzymatic oxidation of xenobiotic phenolic compounds by laccases produces free radicals which react among themselves to form dimers, oligomers or polymers covalently linked by type C-C, C-O and C-N bonds [13]. In other cases, the reaction may be accompanied by partial demethylation and/or dehalogenation. The ability of laccase to carry out such reactions is the basis for its potential application in bioremediation, because it allows these compounds to interact with the organic matrix of the soil [14]. It is well established that the aerobic biodegradation pathways for diuron proceed by successive demethylation steps to form DCPMU [1-(3,4dichlorophenyl)-3-methylurea], DCPU [1-(3,4-dichlorophenyl) urea], and finally, 3,4-DCA (3,4-dichloroaniline) [15]. We believe that the first two reactions are the most likely to be catalyzed by laccase. It is probable that 3,4-DCA is not produced by the laccase reactions, an interesting point if one takes into account that this compound is more toxic than diuron itself. It is important to note that diuron is a chlorophenol compound and laccases have a well-known dechlorination activity of chlorophenol compounds.  [20]. After careful analysis of these observations, we believe that the results of the simulations presented in this work suggest that the template structure 2HRG should have greater affinity to the compound diuron than the CersuLac1. However, the affinity of CersuLac1 for diuron determined in silico is also considerable, and suggests that the enzymatic extract of this fungus containing this enzyme could be a very promising alternative in the field of bioremediation.

Conclusion:
The structural model of the laccase-1 from Ceriporiopsis subvermispora (CersuLac1) solved in this study has a similar fold with lacases from Trametes trogii and Trametes sp (cath id: 2.60.40.420). This work gives evidences suggesting that Ceriporiopsis subvermispora, a fungus commonly used in bioremediation and degradation of lignocellulosic biomass may also be used for the biotransformation of the herbicide diuron, since there are no reports in literature regarding this purpose. Our results also show structural insights into the catalytic mechanism of fungal laccases, as the need of Histidin Nδ2 protonation is coordinated by a copper ion, in order to interact with ligand.