Article pubs.acs.org/est
Computational Evidence for the Detoxifying Mechanism of Epsilon Class Glutathione Transferase Toward the Insecticide DDT Yanwei Li, Xiangli Shi, Qingzhu Zhang,* Jingtian Hu, Jianmin Chen, and Wenxing Wang Environment Research Institute, Shandong University, Jinan 250100, P. R. China S Supporting Information *
ABSTRACT: A combined quantum mechanics/molecular mechanics (QM/MM) computation of the detoxifying mechanism of an epsilon class glutathione transferases (GSTs) toward organochlorine insecticide DDT, 1,1,1-trichloro-2,2-bis(p-chlorophenyl)ethane, has been carried out. The exponential average barrier of the proton transfer mechanism is 15.2 kcal/mol, which is 27.6 kcal/mol lower than that of the GS-DDT conjugant mechanism. It suggests that the detoxifying reaction proceeds via a proton transfer mechanism where GSH acts as a cofactor rather than a conjugate. The study reveals that the protein environment has a strong effect on the reaction barrier. The experimentally proposed residues Arg112, Glu116 and Phe120 were found to have a strong influence on the detoxifying reaction. The influence of residues Pro13, Cys15, His53, Ile55, Glu67, Ser68, Phe115, and Leu119 was detected as well. It is worth noticing that Ile55 facilitates the detoxifying reaction most. On the basis of the structure of DDT, structure 2, (BrC6H4)2CHCCl3, is the best candidate among all the tested structures in resisting the detoxification of enzyme agGSTe2.
1. INTRODUCTION Glutathione transferases (GSTs) are a superfamily of multifunctional enzymes present ubiquitously in aerobic organisms, plants, and animals.1−13 Their primary function is to detoxify endogenous and xenobiotic electrophiles by catalyzing the conjugation of these substances to glutathione (GSH), the tripeptide γ-Glu-Cys-Gly.11,14−16 The resultant products are more water-soluble and excretable than the non-GSH conjugated substrates.17 In addition, GSTs are involved in a wide range of biological processes. They play an important role in protecting the cells from the harmful effects of oxidative stress and have been implicated in various biosynthetic pathways.18−23 GSTs in insects are of particular interest because of their potential role in insecticide metabolism producing resistance.24−27 Organochlorine insecticide DDT, 1,1,1-trichloro-2,2-bis(p-chlorophenyl)ethane, plays a prominent role in the malaria (the most severe insect transmitted disease with at least 300 million cases per year and 1.4 million deaths) eradication by controlling the population of malaria vector mosquitoes since 1950s.28−30 However, resistance to DDT has developed in some mosquito species (e.g., Anopheles Gambiae), raising the threat of malaria to humans.31,32 A mechanistic understanding of the detoxifying of GSTs toward DDT is critically warranted to understand the impact of DDT resistance and to develop more effective novel insecticides.33 At least six classes of insect GSTs have been identified: delta, epsilon, omega, sigma, theta, and zeta, being possibly the existence of novel GST classes.34,35 The delta and epsilon classes are the major GSTs in insects. The GSTs from the © 2014 American Chemical Society
malaria vector Anopheles Gambiae have been extensively studied.30−32,36−43 Twenty eight GSTs from Anopheles Gambiae have been identified, in which 12 belong to the deta class and 8 to the epsilon class.38 Although the DDT-detoxifying ability of the deta class GSTs was proved, they are not considered important in conferring DDT resistance based on the observations that none of them are overexpressed in DDT resistant strains.39,40 However, five of the eight epsilon members have elevated expression level. Among them, agGSTe2 is overexpressed about 8-fold and represents up to 92% DDT metabolism.41 The cytosolic GSTs are homo or heterodimeric proteins, that is, they are formed by two subunits or polypeptide chains of approximately 25 kDa in size each. Each subunit has a G-site where glutathione (GSH) substrate binds and an H-site pocket for electrophilic substrates. When GSH binds to the G-site, it is activated to a negatively charged thiolate anion (GS−) under the help of a water molecule.11,44−49 GS− is stabilized by a hydrogen-bond-network which distributes the negative charge to the surrounding networked atoms.43 The nucleophilic thiolate anion is capable of attacking the electrophilic center of the lipophilic compounds. Generally, the GSTs catalyze the detoxification reactions of halogenated hydrocarbons via proton transfer and/or glutathione conjugation.21,50,51 The proton Received: Revised: Accepted: Published: 5008
November 25, 2013 March 2, 2014 March 28, 2014 March 28, 2014 dx.doi.org/10.1021/es405230j | Environ. Sci. Technol. 2014, 48, 5008−5016
Environmental Science & Technology
Article
Scheme 1. Reaction Pathways of the Proton Transfer Mechanism and the GS-DDT Conjugation Mechanism of agGSTe2 with DDTa
a
The boundary between the QM and MM regions is indicated by wavy lines.
procedure,53 and the missing hydrogen atoms were complemented through the HBUILD facility in the CHARMM package. 54−56 DDT was built and docked with the agGSTe2−GS− and agGST1−6−GS− complexes through a grid-based receptor-flexible docking module (CDOCKER) by using Material Studio 4.4 and Discovery Studio 2.1 (Accelrys Software Inc.), and the best docking pose was identified by the lowest binding energy.57,58 Then, the agGSTe2−GS−−DDT and agGST1−6−GS−−DDT ternary complexes were placed in a water sphere (TIP3P model59) with a diameter of 65 Å, which guarantees that the whole protein was completely solvated. Water molecules overlapping within 2.5 Å of the ternary complex were deleted. The whole system was neutralized with seven sodium ions. After that, the system was heated from absolute zero to 298.15 K in 50 ps, and a trajectory of 500 ps was computed to reach an equilibration state (1 fs/step). A 9 ns stochastic boundary molecular dynamics (SBMD) simulation with canonical ensemble (NVT, 298.15 K) was performed to mimic the aqueous environment.60 During the SBMD simulations, the coordinates of GS− were restrained to keep its positions consistent with the crystal structure. A stochastic buffer region was used and a deformable potential was imposed at the edge of the water sphere to mimic the effect of the water solvent. Leap-frog algorithm and Langevin dynamics implemented in the CHARMM package were applied during the simulation, and the friction coefficients of Langevin thermostat for protein and water molecules were set to 200.0 and 62.0 (corresponding unit is 1/ps). The snapshots of the system were saved every picosecond resulting in 9000 conformations. For the system of agGSTe2-GS−-DDT, seven snapshots were taken every 0.5 ns from 6 to 9 ns of the molecular dynamics run. They are labeled as SN-6.0, SN-6.5, SN-7.0, SN-7.5, SN-8.0, SN-8.5, and SN-9.0. For the system of agGST1−6−GS−− DDT, three snapshots were taken from the MD run after 8.0, 8.5, and 9.0 ns (labeled as SS-8.0, SS-8.5, and SS-9.0). These structures serve as starting configurations for the following geometry optimizations and transition-state search. 2.2. QM/MM Calculations. The QM/MM calculations were performed with the aid of ChemShell61 integrating
transfer is considered to be one of the possible DDT resistance mechanism in mosquitoes based on the direct detection of DDE by mixing agGSTe2 and DDT in the presence of GSH, coupled with the fact that the GS-DDT conjugant has never been identified.51 However, there is a possibility that the GSDDT conjugant is unstable and it may quickly decompose to produce DDE. Thus, more evidence is required to clarify the DDT detoxifying mechanism. The crystal structures of agGSTe2 and agGSTe2−GSH binary complex were determined with a resolution up to 1.4 Å.43 However, attempts to cocrystallize the substrate DDT with agGSTe2 as well as the metabolism product DDE with agGSTe2−GSH complex were not successful.43 The comparison with the less active delta GST, agGST1−6, from the same mosquito vector indicates that in the structure of agGSTe2− GSH complex, residues Arg112, Glu116 and Phe120 is closer to the G-site resulting in a more efficient GS−-stabilizing hydrogenbond-network and higher DDT-binding affinity.43,52 However, the influence of residues Arg112, Glu116, and Phe120 as well as other key residues surrounding the active site on the DDT detoxifying reaction is still unknown. In the present work, the detoxifying mechanism of agGSTe2 toward DDT has been carried out with the aid of a combined quantum mechanics/ molecular mechanics (QM/MM) computations. QM/MM computations of the enzyme-catalyzed reaction can provide the atomistic details of the enzyme mechanism and is therefore becoming an increasingly important tool to complement experimental enzyme chemistry. In order to compare with the detoxifying ability of agGSTe2, the detoxifying mechanism of agGST1−6 toward DDT was also investigated in this study.
2. MATERIALS AND METHODS 2.1. System Setup and Molecular Dynamics. The initial models for the present simulation were built on the basis of Xray crystal structure of agGSTe2-GSH complex (PDB code: 2IMI) and agGST1−6−GSH complex (PDB code: 1PN9) obtained from the Protein Data Bank (www.rcsb.org). The protonation states of the ionizable residues were determined on the basis of the pKa values obtained through the PROPKA 5009
dx.doi.org/10.1021/es405230j | Environ. Sci. Technol. 2014, 48, 5008−5016
Environmental Science & Technology
Article
Figure 1. a, root-mean-square deviation for the system of agGSTe2-GS−-DDT during 9,000 ps MD simulation. b, B-factors for the system of agGSTe2-GS−−DDT (hydrogen atoms are not included). c, radial distribution of water oxygen atoms relative to sulfur atom in the systems of agGSTe2-GS− and agGSTe2-GS−-DDT. d, distance variations of S−Cβ, S−Hα, S−Hβ, S−Hδ, and S−Oβ during 9,000 ps MD simulation.
force field66 was chosen for the MM-region. Covalent bonds between the QM-region and MM-region were truncated and modified by adding hydrogen link atoms to the QM side. A
Turbomole62 and DL-POLY programs.63 Hybrid delocalized internal coordinate (HDLC)64 was adopted. The QM-region was treated by the DFT65 method, while the CHARMM22 5010
dx.doi.org/10.1021/es405230j | Environ. Sci. Technol. 2014, 48, 5008−5016
Environmental Science & Technology
Article
Table 1. Energy Barriers and Selected QM/MM Bond Distances in the Reactant (R), Transition State (TS), and Product (P) Involved in the Proton Transfer Mechanism of agGSTe2 in Seven Pathways distances (Å) Cα−Hα SN-6.0 SN-6.5 SN-7.0 SN-7.5 SN-8.0 SN-8.5 SN-9.0
Cα−Cβ
S−Hα
Cβ−Clα
barriers (kcal/mol)
R
TS
P
R
TS
P
R
TS
P
R
TS
P
28.2 39.9 20.5 32.4 21.3 14.1 21.1
1.09 1.09 1.09 1.09 1.09 1.10 1.09
1.76 1.67 1.67 1.83 1.68 1.74 1.75
4.08 3.73 4.21 4.02 4.27 3.35 3.15
4.23 4.25 4.07 4.06 4.28 3.92 4.14
1.48 1.55 1.42 1.50 1.53 1.54 1.49
1.34 1.35 1.35 1.35 1.35 1.34 1.34
1.56 1.56 1.55 1.55 1.55 1.55 1.55
1.42 1.44 1.53 1.43 1.46 1.44 1.43
1.34 1.34 1.34 1.35 1.35 1.34 1.34
1.79 1.80 1.80 1.80 1.80 1.80 1.80
2.17 2.09 2.21 1.96 1.91 2.04 2.10
2.90 2.88 2.99 2.59 3.52 2.73 2.92
charge shift model67 was used to avoid over polarization of the QM density of the QM-region. The QM-region contains cofactor GS−, DDT, a water molecule and part of Ser12. The truncation of Ser12 was done by cutting a C−C bond. This resulted in 71 QM atoms including a hydrogen link atom. The QM-region was optimized by the B3LYP/6-31G(d,p) method68,69 with a charge of −3 and a spin multiplicity 1. The potential energy profile from the reactant to the product was scanned to determine the transition state structure. The corresponding structure of the highest energy point along the reaction path was extracted and further optimized by combining partitioned rational function optimizer (P-RFO) algorithm70 and the low-memory Broyden-Fletcher-GoldfarbShanno (L-BFGS) algorithm.71 Harmonic vibrational frequency calculations were performed to validate its transition state character. A larger basis set of cc-pVTZ, which is comparable to basis set 6-311++G(2df,2pd), was employed in single point energy calculations. The energies include ZPE corrections with unscaled frequencies obtained at the B3LYP/6-31G(d,p) level. The natural bond orbital (NBO) analysis was performed at the B3LYP/6-311++G(d,p)//CHARMM22 level to investigate the key atom charge variations during the reaction process.
9,000 ps MD simulation. The corresponding root-mean-square deviation (RMSD) and B-factor analysis were provided in Figure 1a and 1b. The radial distribution function (RDF) analysis for the systems of agGSTe2−GS− and agGSTe2-GS−− DDT is given in Figure 1c. The water coordinate numbers for the systems of agGSTe2−GS− and agGSTe2−GS−-DDT are about 1.6 and 5.4, which indicates that about four water molecules are excluded from the active site after DDT binding. Five distance variations (S−Cβ, S−Hα, S−Hβ, S−Hδ, and S− Oβ) along the 9,000 ps trajectory were depicted in Figure 1d. The average S−Oβ distance (2.9 Å) and its small variation during the MD simulation indicate that a water molecule is involved in stabilizing the sulfur anion. Further analysis of S− Hβ or S−Hδ (∼2 Å) suggests that the hydrogen bond is formed through either Hβ or Hδ. S−Hα and S−Cβ distances are ∼4.0 Å and ∼5.5 Å, which indicates that ionized GS− can attack either DDT α-hydrogen or β-carbon, resulting in two possible reaction mechanisms, the proton transfer mechanism and the GS-DDT conjugation mechanism, as shown in Scheme 1. For the proton transfer mechanism, the DDT α-hydrogen transfers from DDT to GS−. A substantial energy barrier spread from 14.1 to 39.9 kcal/mol is found among different snapshots (Table 1). This suggests that the protein environment has a strong effect on the reaction barrier. To analyze the computed results, average energy barriers of the seven snapshots were calculated by exponential average:72,73
3. RESULTS AND DISCUSSION For convenience of description, several atoms in the QMregion are numbered, as presented in Scheme 1. The first step of this work is to identify the calculation method which not only is able to produce accurate theoretical results, but also is computationally feasible and economical for currently available hardware and software. Due to the absence of the X-ray crystal structure of the agGSTe2-GS−−DDT ternary complex, it is difficult to make a direct comparison between the calculated results and the experimental data. Thus, we took two snapshots from the MD simulation as starting configurations of the agGSTe2−GS− binary complex and optimized them at the B3LYP/6-31G(d,p)//CHARMM22 level. The calculated results agree well with the available experimental values. For example, the spatial distances between S and Oα (Scheme 1) are 3.24 and 3.26 Å, in agreement with the experimental value of 3.48 Å. The S−C bonds in GS− are 1.82 and 1.82 Å, which match perfectly with the X-ray value of 1.82 Å. A larger basis set 6-311++G(d,p) was also used to evaluate the accuracy of the B3LYP/6-31G(d,p) level, and the structural parameter differences are within 2% (Supporting Information (SI) Table S1). So, it might be inferred that the B3LYP/6-31G(d,p) level is appropriate for the QM-region geometrical optimization. 3.1. Reaction Mechanism. The agGSTe2-GS−−DDT ternary complex was extracted per picosecond during the
{
ΔEea = −RT ln
n ⎛ −ΔEi ⎞⎫ 1 ⎟⎬ ∑ exp⎜⎝ n i=1 RT ⎠⎭ ⎪
⎪
Where, ΔEea is the average barrier, R is gas constant, n is the number of snapshots, ΔEi is the energy barrier of path i, and T is the temperature. The calculated exponential average barrier of the proton transfer mechanism is 15.2 kcal/mol. The feasibility of an alternative GS-DDT conjugation mechanism was also investigated. The structural details are provided in Figure S1. In this mechanism, GS− attacks the βcarbon atom of DDT and results in a GS-DDT conjugant. The seven energy barriers calculated at the B3LYP/cc-pVTZ// CHARMM22 level vary from 41.6 to 67.9 kcal/mol with an exponential average barrier of 42.8 kcal/mol, which is 27.6 kcal/mol higher than that of the proton transfer mechanism (Figure S2). This shows that the proton transfer mechanism is more energetically feasible than the GS-DDT conjugation mechanism. This is supported indirectly by a fact that the experimentally determined specific activity of agGSTe2 toward DDT corresponds to an energy barrier ∼17 kcal/mol,41,42,74,75 which is in good agreement with the exponential average barrier of the proton transfer mechanism (15.2 kcal/mol) rather than 5011
dx.doi.org/10.1021/es405230j | Environ. Sci. Technol. 2014, 48, 5008−5016
Environmental Science & Technology
Article
that of the conjugation mechanism. As a result, only the proton transfer mechanism will be minutely investigated in the following study. 3.2. Catalytic Itinerary and Structural Details. For the proton transfer mechanism, lengths of bonds Cα-Hα, S−Hα, CαCβ and Cβ-Clα in the reactant, transition state, and product obtained at the B3LYP/6-31G(d,p)//CHARMM22 level were provided in Table 1. Since a majority of the catalytic reactions proceed through the lowest energy barrier of 14.1 kcal/mol, the following investigations will mainly focus on the pathway SN8.5. Figure 2 shows the active site structures of the reactant, transition state and product for the pathway SN-8.5. In the reactant, negatively charged sulfur is stabilized by Ser12 and a water molecule (S−Hγ, 2.04 Å; S−Hβ, 2.02 Å). Besides Ser12, other key residues which are spatially close to either DDT or GS− were also depicted in SI Figure S3 and Figure S4. The effect of these residues on the reaction mechanism will be discussed below. For the transition state structure, its character is confirmed by the vibrational mode and the corresponding imaginary frequency (1024i cm−1). In the product, the length of double bond CαCβ (1.34 Å) and the angles Clβ-Cβ-Clγ (114.8°), Cα-Cβ-Clβ (120.1°), and Cα-Cβ-Clγ (123.6°) suggest the formation of DDE. It is worth mentioning that the chlorphenyl groups and the CαCβ bond are not in a same plane. The distance of Cβ-Clα (2.73 Å) suggests the formation of a chloride anion. Figure 3 shows several selected bond length and atom charge variations along the pathway SN-8.5. The changes of the bond lengths of S−Hα, Cα-Hα, Cα-Cβ, and CβClα from the reactant to the product through the transition state indicate that the hydrogen transfer, the formation of the CαCβ double bond and the dechlorination process occur simultaneously. The anion character of Clα in the product was further confirmed by its negative charge (Clα, −0.87). 3.3. Hydrogen-Bond-Network and Individual Residue Influence. Experimental study about the hydrogen-bondnetwork is only focused on the interactions between GS− and agGSTe2 in the reactant structure.41 However, the hydrogenbond-network in the transition state should also be investigated because of its significant influence on the reaction barrier. Here variations of six hydrogen bonds (≤2.5 Å) between GS− and enzyme agGSTe2 along the reaction pathway of SN-8.5 are investigated (Figure 4). Variations of these hydrogen bonds fall into two groups: the hydrogen bond in the reactant is stronger than in the transition state, and the hydrogen bond in the reactant is weaker than in the transition state. For instance, the hydrogen bond variations along the reaction pathway SN-8.5 show that GS− is more stabilized by Ile55 in the transition state than in the reactant, whereas it is less stabilized by Arg112 in the transition state than in the reactant. The result is in excellent agreement with our following electrostatic influence study where it is shown that Ile55 facilitate the reaction by decreasing the reaction barrier of about 2.3 kcal/mol (ΔEi‑0 = 2.3 kcal/mol), whereas Arg112 suppresses the catalytic reaction by increasing the reaction barrier of about 1.9 kcal/mol (ΔEi‑0 = −1.9 kcal/mol). The electrostatic influence of 19 key individual residues on the energy barrier of the pathway SN-8.5 was estimated. The electrostatic influence of Ser12 is not studied because it is included in the QM-region. The activation energy difference caused by amino acid i can be described as
Figure 2. Structures of the reactant (R), transition state (TS), and product (P) involved in the pathway SN-8.5 of the proton transfer mechanism for agGSTe2 with DDT. QM atoms including a link hydrogen atom were shown through ball and stick representation. The unit for bond distances and imaginary frequency are in Å and cm−1.
Where ΔEi‑0 is the changes of the barrier, ΔEi is the energy barrier with charges on residue i set to 0, and ΔE0 is the original
ΔEi − 0 = ΔEi − ΔE 0 5012
dx.doi.org/10.1021/es405230j | Environ. Sci. Technol. 2014, 48, 5008−5016
Environmental Science & Technology
Article
energy barrier. In other words, the ith residue lowers the energy barrier and facilities the catalytic reaction. The ΔEi‑0 values of 19 residues were schematically represented in Figure 5. Besides
Figure 5. ΔEi‑0 values of individual residues toward the proton transfer mechanism of agGSTe2 with DDT.
electrostatic influence of residues Ile55 and Arg112 discussed above, we have also shown that residues Cys15, His53, Glu67, and Glu116 suppress the detoxifying reaction (ΔEi‑0 < −1 kcal/mol), whereas residues Pro13, Ser68, Phe115, and Phe120 facilitate the detoxifying reaction (ΔEi‑0 > 1 kcal/mol). There are nine other residues that are found to have a weaker contribution (−1 kcal/mol
