Seeking the Source of Catalytic Efficiency of Lindane Dehydrochlorinase, LinA.

Herein we present the results of an in-depth simulation study of LinA and its two variants. In our analysis, we combined the exploration of protein conformational dynamics with and without bound substrates (hexachlorocyclohexane (HCH) isomers) performed using molecular dynamics simulation followed by the extraction of the most frequently visited conformations and their characteristics with a detailed description of the interactions taking place in the active site between the respective HCH molecule and the first shell residues by using symmetry-adapted perturbation theory (SAPT) calculations. A detailed investigation of the conformational space of LinA substates has been accompanied by description of enzymatic catalytic steps carried out using a hybrid quantum mechanics/molecular mechanics (QM/MM) potential along with the computation of the potential of mean force (PMF) to estimate the free energy barriers for the studied transformations: dehydrochlorination of γ-, (-)-α-, and (+)-α-HCH by LinA-type I and -type II variants. The applied combination of computational techniques allowed us not only to characterize two LinA types but also to point to the most important differences between them and link their features to catalytic efficiency each of them possesses toward the respective ligand. More importantly it has been demonstrated that type I protein is more mobile, its active site has a larger volume, and the dehydrochlorination products are stabilized more strongly than in the case of type II enzyme, due to differences in the residues present in the active sites. Additionally, interaction energy calculations revealed very interesting patterns not predicted before but having the potential to be utilized in any attempts of improving LinA catalytic efficiency. On the basis of all these observations, LinA-type I protein seems to be more preorganized for the dehydrochlorination reaction it catalyzes than the type II variant.

Introduction

For many years now, quite a substantial amount of research has been devoted to understanding how enzymes work. Being essential for all living organisms and existing as ensembles of configurations, they are capable of creating life by performing chemistry with remarkable efficiency. In recent years experimental endeavors aiming at improving enzymes efficiency have been supported by modern computational tools. Although computational enzyme design alone is considered to be a relatively young field, it has growing possibilities coming from the development of the high-performance computing (HPC) technology as well as advances in computations acceleration that allow recognition of major factors contributing to biocatalysis.[1−4]

In the present work we focus our attention on lindane dehydrochlorinases, LinA, which although isolated from bacteria only in 1991[5] have been already demonstrated to be present in areas in which large amounts of various hexachlorocyclohexane isomers (HCH) had been deposited in the soil leading to severe contamination. HCH compounds have six chlorine substituents in the six-membered saturated carbon ring. The so-called technical HCH (t-HCH) is an industrially synthesized mixture of five stable isomers α, β, γ, δ, and ε-HCH. The composition of t-HCH is dominated by α-HCH in 60–70%, a chiral molecule found as a racemate, while γ-HCH (also known as lindane) is the only component of the mixture possessing the insecticidal properties. Altogether they constitute a noxious and potentially carcinogenic persistent pollutant.[6,7] Broad use of t-HCH as a pesticide in agriculture and forestry along with unregulated disposal after the industrial manufacturing of t-HCH has led to a global problem caused by HCH contamination.[8−10] Among HCH isomers, α-HCH has been assigned as group B2 possible human carcinogen[11] and β-HCH has been found to accumulate in human tissues.[12]

Due to their different physical and chemical properties governed by the orientation of chlorine substituents, the HCH isomers exhibit different susceptibility to various ways of degradation. The more chlorine atoms are at axial positions, the more accessible for biodegradation the compound is. The isomers with very low (one in δ-HCH) or zero (in β-HCH) axial substituents are more stable and persist in the environment for a longer time. α- and γ-isomers having four and three axial chlorine atoms, respectively, undergo readily microbial biodegradation. One of the already well-known and studied species responsible for the aerobic transformation of HCH isomers is Sphingobium japonicum UT26.[13]S. japonicum UT26 metabolizes γ-HCH to succinyl-CoA and acetyl-CoA by means of eight enzymes encoded by lin genes; first, sequential reactions catalyzed by LinA, LinB, and LinC lead to 2,5-dichlorohydroquinone, which is then further transformed by LinD, LinE, LinF, LinGH, and LinJ. In addition to the aforementioned enzymes, a putative ABC transporter system encoded by linKLMN was found to be an essential factor for γ-HCH utilization in the strain.[14] The first steps of dechlorination of α- and γ-HCH are metabolized by a dehydrochlorinase LinA via elimination mechanism (Scheme ).[15−19]

Scheme 1

Degradation Processes of γ/α-HCH Isomers Catalyzed by LinA

There are three well-characterized variants of LinA: LinA-type I, LinA-type II, and LinA1. LinA-type I has been isolated from the Sphingobium japonicum UT26 strain.[5] LinA-type II has been obtained from a soil metagenome containing HCH.[20] LinA1 has been found in other species: Sphingobium indicum B90A.[21] A number of previous studies have demonstrated the enantioselectivity of known LinA variants toward various HCH isomers.[12,22] LinA-type I has been shown to transform γ-HCH to a single γ-pentachlorocyclohexene enantiomer (γ-PCCH-2), and subsequently γ-PCCH-2 was shown to be metabolized faster than γ-PCCH-1.[12] The same LinA variant exhibits preference toward (−)-α-HCH enantiomer leading to the β-PCCH-1 formation.[22] In contrast, LinA-type II and LinA1 prefer the transformation of (+)-α-HCH.[23] Since the crystal structures of LinA-type I[23] and LinA-type II[24] variants have been resolved, it was possible to analyze these two enzymatic forms structurally and use their coordinates in docking studies.[22,23,25] It has been shown that type I variant differs from type II by 10 residues (K20C, A23G, F68Y, C71T, L96C, F113Y, D115N, R129L, A131G, and T133M). On the basis of the subsequent mutagenesis studies, it has been suggested that 4 out of these 10 residues vicinal to the active site of LinA might govern the enantioselectivity toward α-HCH enantiomer but do not necessarily play a similar role in the transformation of γ-HCH.[25]

In general, LinA variants are not very efficient enzymes. Their catalytic efficiency depends on the type of protein and the isomer to be metabolized; however, it does not exceed 104 M–1·s–1.[26] The only exception is the reaction of LinA-type I with (−)-α-HCH which has been very recently hypothesized to reach the diffusion-controlled limit.[27] Therefore, from this perspective, there is room and need for improvement, in particular, if any successful bioremediation strategy is of interest. So far, attempts in this direction were only based on the analysis of inter-residual distances within the active sites of the respective variants resulting from docking and mutagenesis studies. In 2001, Nagata et al.,[26] based on activities of mutant proteins, suggested that, apart from the catalytic dyad D25-H73, R129 is also directly involved in the catalysis. In the same article they proposed that the affinity of LinA for the substrate can be influenced by mutations of W42, L64, L96, and F113. Additionally, their analysis indicated that substrate binding and product release are not the rate-limiting steps.[28] Sharma et al. postulated that residues 110 and 111 influences γ/α isomer specificity of the enzyme, while sequence differences in positions 68, 71, 96, 133 affect enantioselectivity of LinA variants for the (±)-α-HCH.[29] Furthermore, it has been shown that a LinA-type I mutant that carries four changes, namely, K20Q, L96C, A131G, and T133M has reversed enantioselectivity compared to the wild type. However, the complete change of enantioselectivity has not been reached yet, so the molecular factors responsible for this behavior have not been fully described and understood.[25] All those results have contributed to disclosing the molecular mechanism of LinA enzymes and making a step forward toward improving their efficiency; however there are still many key aspects unaddressed. As long as the source of the differences between LinA variants remains elusive, any rational modification of the enzyme in order to improve it toward selected isomers will be prohibited.

In order to provide more insight at the molecular level into the recognition of the substrate and its subsequent interactions within the binding site of LinA, we have explored conformational dynamics of the protein with and without the respective ligands bound and computed and analyzed the reaction free energy profiles and the interaction energies between each isomer and the most important first shell active site residues.

Methods

Preparation of Models

Initial coordinates of LinA enzymes were taken from Protein Data Bank (PDB codes 3A76 and 3S5C for type I and type II, respectively). LinA is a homotrimeric protein, where monomers form cone-shaped α + β barrel folds, with identical active sites placed within each monomer.[23] As the crystal structures of LinA-type I and LinA-type II were not resolved with any ligand bound, docking of the substrates was conducted using Glide 5.0, as described previously,[18] for type I complexes and using PatchDock[30,31] for type II complexes. Top-ranked candidate structures of complexes according to a geometric shape complementarity score were analyzed, and the best candidate was used for further model preparation. The protonation states of titratable amino acids were determined by the PROPKA 3.1 program[32] at pH 7. The AM1-BCC charges[33,34] for γ- and α-HCH isomers were computed using the Antechamber module of the AMBER package.[35] LinA variants were described by the ff14SB[36] force field, whereas for the ligands GAFF[37] was used. Twenty-one sodium ions were added to neutralize each system. The ligand–enzyme complexes were subsequently solvated with a 10 Å radius layer of water molecules using the TIP3P water model.[38] The resulting systems consisted of approximately 51 000 atoms. Minimization was performed as previously described.[39] Then, the systems were heated from 100 to 300 K during 300 ps using the NVT ensemble with a time step of 1 fs and equilibrated using the NPT ensemble. Three replicas of classical MM MD simulations were performed for each complex. The production runs were performed using the NPT ensemble with a time step of 1 fs, and the length of the simulations was equal to 500 ns. Periodic boundary conditions with the particle mesh Ewald method[40] and 12 Å cutoff distance for nonbonding interactions were applied throughout the simulations.

Conformational Analysis

The trajectories from the production runs were analyzed using the principal component analysis (PCA) method as implemented in the PyEMMA package.[41] PCA is a technique for reducing the dimensionality of large data sets providing increased interpretability with minimum information loss. It has been quite successfully used in the analysis of trajectories resulting from MD simulations.[42−44] The analysis was performed with the coordinates of α-carbon atoms on every 50th frame of the original trajectories. Additionally, the root-mean-square fluctuation (RMSF) values and selected distances were calculated using the CPPTRAJ module[45] implemented in AMBER. Volume calculations were performed using the Computed Atlas of Surface Topography of proteins (CASTp) Web server with the radius probe of 1.4 Å.[46] Visualization of analyzed complexes was performed using VMD[47] and the Discovery Studio BIOVIA software.[48]

Symmetry-Adapted Perturbation Theory Calculations

In order to investigate interactions between the selected first shell residues and HCH, analysis using the symmetry-adapted perturbation theory (SAPT)[49] was conducted. The method allows computation of the noncovalent interaction between two molecules and provides a decomposition of the interaction energy into four physically motivated components: electrostatics, exchange–repulsion, induction, and dispersion. In order to reduce the influence of geometry differences between the analyzed models and to focus on the HCH spatial structure effect on its interaction with the enzyme, the analysis has been carried out for models composed of single structures obtained by short energy minimization of the complexes conducted using the AMBER package. We performed the interaction energy calculation using the high-order SAPT (SAPT2+) with correlation-consistent augmented double Dunning basis (aug-cc-pVDZ)[50] in the PSI4 software package.[51] Furthermore, in order to speed up SAPT calculations considerably (without significantly affecting overall accuracy), density fitting was used. Given the specificity of the method of complex preparation, only the total interaction energy values calculated using SAPT (ΔEint) were compared.

Quantum Mechanics/Molecular Mechanics Free Energy Simulations

A hybrid QM/MM scheme[52,53] with the DFTB3 Hamiltonian[54] to treat the QM part of each model as implemented in the AMBER package was used to obtain the minimum free energy paths (MFEPs) by means of the adaptive string method.[55] The QM region of each complex comprised the HCH ligand and the side chains of H73 and D25 residues, whereas the rest of the protein and the solvent were treated by the respective force fields as during the classical MD simulations. The hydrogen-link atoms were used as a boundary between the QM and the MM regions. MD simulations were performed with a Langevin thermostat at 300 K, velocity Verlet integrator,[56,57] and periodic boundary conditions with particle mesh Ewald[58,59] to treat long-range electrostatic interactions. Replica exchange was applied during both the string optimization and the umbrella sampling stages with exchange attempts performed every 100 fs. Three collective variables (CVs) were defined to monitor the elimination of the H/Cl pair from the HCH molecule by the imidazole ring: the Cx–Hx, Hx···Nε2, and Cy–Cly interatomic distances. Additionally, two variables were used to control the proton transfer between H73 and D25: Nδ1–Hδ1 and Hδ1···Oδ1 or Oδ2 depending on the protein variant and ligand preference. Fifty string nodes were used to follow the reaction progress to the products. These nodes are propagated according to the mean force and kept equidistant, converging to the MFEP. Then a path CV, denoted as s, is defined to measure the advance of the system along the MFEP from reactants to products and used as the reaction coordinate to trace the potential of mean force (PMF) associated with the reaction under analysis. At this stage simulations are performed until the statistical error of the free energy barriers drops below 1 kcal·mol–1. In order to keep the abstracted hydrogens on the path, we used a biasing potential perpendicular to the minimum free energy path Vb(z) and increased the mass of both protons to 2.0 amu in some of the studied complexes. Basic simulation parameters are given in Table S1.

Results and Discussion

Active Site Dynamics and Catalysis

We have initiated our study by exploring whether the LinA variants can adopt distinct conformations (open vs close) upon substrate (HCH) binding and whether this conformational change would have an impact on reactivity, as it has been postulated earlier.[23] If such a conformational change does exist, it should be related to the extensive movement of the C-terminal fragment from the neighboring chain as illustrated in Figure S1 and Figure . In order to test this hypothesis, we constructed enzyme–substrate complexes for each studied protein variant (LinA-type I and LinA-type II) and HCH isomer ((±)-α and γ). In this way we obtained six different complexes (holo conformations of LinA). Additionally, we also studied both protein types without any ligand bound (apo conformations of LinA). All those models were subjected to long-time classical MD simulations performed in three replicas of 500 ns each. In order to find out whether prevalent atomic motions were somehow connected to the motion of C-terminus in respective simulations, we resorted to one of the dimensionality reduction methods: principal component analysis (PCA). Two principal components were extracted (PC1 and PC2); in the case of type I they constitute almost 50% of overall variance whereas in type II 40% at maximum (Figure S2). First, we looked at the position and orientation of C-terminus as well as loop 7 (Figures S3–S6).

Figure 1

(A) LinA-HCH complex. HCH and catalytic dyad, D25 and H73, are represented using sticks-and-balls. C-terminus and loop 7 are highlighted in the darker shades of the color representing the chains they belong to. (B) Active site of LinA, where residues 64–66 have been omitted for clarity. (C) Visualization of the QM region of the complex, where atoms included in the QM part of the system are colored in red.

Our simulations showed large displacements of the C-terminal regions of the three chains that, in some cases, can get close to the active site of the neighboring chain. The analysis of our simulations showed some correlation between closing movement of C-terminus and the distance to loop 7 in the case of LinA-type II in such a way that if the C-terminus closes, these two regions get closer to each other (Figure , Figure S7). Importantly, changes in the position and orientation of the C-terminal region of the neighboring chain are not correlated with the presence of HCH, as indicated by the values of the distance between C-terminus and the catalytic H73, which ranges from 13.0 to 29.0 Å and from 13.3 to 23.5 Å for holo and apo conformations of LinA, respectively. Furthermore, the shortest distances were observed for enzyme–substrate complexes of α-HCH isomers with LinA-type I and for enzyme–substrate complexes with γ-HCH and for the apo conformation of LinA-type II. On the basis of these results, we can conclude that C-terminus and loop 7 dynamics are not substrate-related, since the movement of C-terminal region leading to close conformation occurs similarly in apo conformations of LinA as well as in its complexes with different HCH isomers. Furthermore, the range of this movement differs in individual replicas of the same complex, suggesting that conformational changes result from the intrinsic flexibility of the protein motifs rather than from any specific protein–ligand interactions (Figures S3–S6).

These results might indicate that the position of C-terminus is not an important factor when catalytic abilities of LinA are considered. In order to confirm that (open/close) conformation of the enzyme does not influence the reaction, we examined changes of the HCH positioning within the LinA active site (measured as a distance from the catalytic base, H73) in relation to changes of the C-terminus position. This distance, in particular between the proton to be abstracted in the HCH molecule and the Nε2 atom of H73, might be an indication of whether the HCH molecule is at the catalytically competent position, ready to be attacked by the base. Present analysis (Figure and Figures S8 and S9) indicates that changes in the HCH–H73 distance are unrelated to the C-terminus movements, which allows us to argue that the substrate binding and catalytic activity of LinA are not necessarily influenced by these conformational motions. This lack of correlation is especially evident in the case of LinA-type II in which movement of the substrate is minimal regardless of the wide fluctuations observed of the C-terminus position and HCH spatial structure differences. Although parallel changes of HCH–H73 and C-terminus-H73 distances can be observed for some LinA-type I complexes (i.e., γ-HCH replica 1 or (−)-α-HCH replica 1; Figure ), these trends do not occur in all three replicas of the same complex (Figures S8 and S9), which suggests that they are not substrate-specific, and therefore they cannot influence the enzyme enantioselectivity. Apart from that, we analyzed the position of the substrate in the case of a complex of single mer from LinA-type I with γ-HCH to further test whether movement of the C-terminus from the neighboring chain leading to close conformation affects substrate binding. The results obtained for the single monomer model (Figure ), where we did not observe any substantial movement of the substrate molecule within the binding pocket (e.g., like substrate release), confirm that the presence of the C-terminal region of the neighboring chain is not a prerequisite for catalysis. A similar conclusion can be reached comparing the similarity between the results obtained from theoretical analysis of the γ-HCH reaction mechanism in a single monomer model of LinA[60] and in a model that included the entire protein.[18]

Figure 2

Evolution of the HCH–H73 (red) and C-terminus-H73 (black) distances during the 500 ns MD simulations of one of the replicas of the studied complexes. The remaining replicas are shown in Figures S8 and S9.

Figure 3

Evolution of the HCH–H73 distance resulting from the 500 ns simulation of the γ-HCH complex with LinA-type I single monomer model (rep 1–3). The inset plot shows exactly the same results but using a smaller y-axis scale.

In the case of α-HCH isomers, although we could easily find examples of conformations in which the distance between the C-terminus and the catalytic base was getting shortened throughout the simulation (Figures S4 and S6), their presence did not affect the distance between H73 and the ligand molecule (Figure , Figures S8 and S9). The difference between the two protein variants in terms of changes in this key distance (HCH–H73) is also worth noting. In type I it tends to fluctuate quite a lot, whereas in type II much less pronounced changes were observed (Figure , Figures S8 and S9).

Modeling of the elimination reaction for α-HCH isomers in the two protein variants seems to additionally confirm these observations (Figure ). In all cases the same net mechanism comprising a concerted E2 pathway was found (Scheme ): the H1 proton is abstracted from the HCH molecules by H73, the chloride anion is eliminated, and the C1–C2 double bond is formed. The PMFs obtained for the reaction mechanisms (Figure ) agree with the experimental observations regarding the substrate preference of each type of enzyme: LinA-type I catalyzes more efficiently the reaction with (−)-α-HCH, while type II is more efficient against (+)-α-HCH.[61] The overall larger catalytic efficiency of type I is also reflected in our simulations, where smaller activation free energies are found for this enzyme. The free energies of activation obtained for all studied complexes are provided in Table S1. Decomposition of the total PMFs in the contributions of each CV to the path (Figures S10–S13, panels A) demonstrates very similar free energy changes along each of the CVs. In all cases the ones that contribute the most to the total free energy change are the carbon–chlorine bond cleavage and the H1–Nε2 bond formation. The difference is however spotted when the variation of the CVs is explored along the reaction coordinate (Figures S10–S13, panels B). In the case of the reaction of (−)-α-HCH with both LinA variants the Hδ1 proton is transferred from H73 to D25, whereas in the case of (+)-α-HCH it is not (Scheme ).

Figure 4

Free energy profile as a function of the path CV for the elimination of one H/Cl pair from (−)-α-HCH (black curves) and (+)-α-HCH (red curves) by LinA-type I (continuous lines) and LinA-type II (dashed lines).

Scheme 2

Postulated Reaction Mechanism for the α-HCH Isomers Dehydrochlorination by LinA Protein Variants

Atoms are numbered and named as in Figure .

Altogether the findings presented in this section do not support the open/close conformation hypothesis postulated previously based on the structural similarity to scytalone dehydratase.[23] This conclusion is also consistent with experimental results showing that mutation within C-terminus of the neighboring chain (i.e., F144L) does not decrease an activity level of the enzyme suggesting that contact of this region with the substrate does not influence the reaction.[28]

Plasticity of the Binding Site

Analysis of the HCH position within the binding pocket drew our attention to significant differences between the two types of LinA that had not been analyzed before; differences between the HCH movement patterns in complexes with LinA-type I and -type II suggest that the structure of the LinA-type II active site hampers this movement (the standard deviation of the HCH–H73 distance ranges from 1.55 to 3.76 Å and from 0.69 to 1.17 Å for LinA-type I and LinA-type II complexes, respectively). Since the smaller size and the plasticity of the type II binding pocket are the most plausible sources for the aforementioned differences, we compared the pocket volume based on the solvent-accessible surface of the two enzymes and measured the RMSF for their holo conformations. The results of these analyses are presented in Table and Figure .

Table 1

Binding Pocket Volume Values, V (in Å3), Calculated for Different Holo Conformations Using CASTp along with Their Standard Deviation, SD (in Å3)

protein: LinA-type I			protein: LinA-type II
ligand	V	SD	ligand	V	SD
γ-HCH	358	67	γ-HCH	214	20
	226			179
	314			179
(−)-α-HCH	305	49	(−)-α-HCH	148	14
	230			167
	213			176
(+)-α-HCH	235	51	(+)-α-HCH	155	18
	250			187
	329			155

Figure 5

Root-mean-square fluctuations (RMSFs) of all residues calculated for all studied holo conformations of LinA complexes with γ-, (−)-α-, and (+)-α-HCH (shown in left, middle, right plots, respectively). Light gray and red colored areas represent mers occupied by the ligand in type I and II variant, respectively.

As expected, LinA-type II has a noticeably smaller binding pocket than LinA-type I since volume values range from 147 to 213 and from 213 and 357 Å3, respectively. Additionally, type II is generally more rigid than type I as indicated by the smaller standard deviation (SD) of the volume measured for different replicas or different enzyme–substrate complexes (on average 17 vs 57 Å3). This finding is confirmed by inspection of the RMSF for the two enzyme types with the three substrates present in one of the active sites (Figure ). The RMSF presents high values for the C-terminal regions in each of the three units in all the combinations of substrates and enzymes, as expected from our previous discussion. Regarding the active site residues, type II presents lower values in all cases, regardless of the kind of substrate or even of the presence or not of substrate bound in the active site. These results are consistent with previously reported differences in the LinA variants thermostability, i.e., higher thermostability of LinA-type II due to larger intersubunit buried surface area and the presence of an additional salt bridge.[24]

The lack of distinction between RMSF values for the monomer of LinA-type II where HCH was bound and those that are empty indicates that the differences observed in the plasticity are not a result of substrate binding (Figure ). The ligand behavior in the active site was analyzed to obtain further information on the correlation between enzyme plasticity and ligand binding. In particular, conformational mobility of the ligand was analyzed as it is known to undergo ring flipping, which leads to the change of trans-1,2-diaxial moieties (H–C–C–Cl) within the cyclohexane ring that can be eliminated from the molecule. Therefore, this aspect is of special interest as it can directly affect catalysis by producing conformers that do not lead to products. This issue has been very recently explored using isotopic analysis by Schilling et al.[26] It has been shown that the C and H isotope effects associated with dehydrochlorination of γ-HCH are independent of the HCH conformational mobility.

However, our analysis shows that ring flipping occurs only solely in the type II variant (Figure S14), which is the more rigid variant. In order to analyze the consequences of this ring flipping motion for catalysis, we have analyzed the impact on the distances between the proton to be abstracted and the base and the angle formed between the proton donor, the proton, and the base (Figure S15). Interestingly, we found that in the case of LinA-type II complexes with (−)-α and γ-HCH, ring flipping leads to conformations that are less competent for catalysis (with larger distances or angles differing from linearity), while in the case of (+)-α-HCH ring flipping does not lead to substantial variations in these two parameters or the changes induced result in conformations more competent for catalysis.

This suggests that plasticity and size of the active site are not the key factors from the catalytic point of view and that steric restrictions resulting from the changes in the amino acid sequence can be compensated (to a certain extent) by HCH conformational changes. In other words, the protein environment is less/worse preorganized in type II for the reaction to occur, despite being more rigid and having smaller volume, and the ligand molecules have to undergo adaptation in order to find the catalytically competent geometry. As the environment does not make it any easier, this adaptation comes with a larger energy cost. It is also worth noting that in the reaction catalyzed by type II protein the products of the reaction are by far less stabilized than in the case of the reaction metabolized by type I, in particular during the reaction with the α-HCH isomers. In type I the chloride leaving group is stabilized by three residues: R129, K20, and W42 (Figure S16). In type II, though, K20 and R129 are replaced by Q and L, respectively, and chloride is eliminated in a different side of the active site where there are no residues capable of establishing similar stabilizing interactions. In type II there is only one Tyr residue nearby that might presumably play this stabilizing role (Figure S16). This evidently worse chloride anion stabilization is also reflected in the free energy of the reaction (Figures. and 6, Table S1).

Figure 6

Free energy profile as a function of the path CV for the elimination of H1/Cl2 pair from γ-HCH (black curves) and H5/Cl4 pair (red curves) by LinA-type II.

In order to understand the possible influence of substrate conformational changes on enzyme selectivity, we analyzed the preference of the carbon atom from which the proton is abstracted during the reaction for different enzyme–substrate complexes (Figures S17–S22). Our studies showed similar discrimination patterns for both LinA types: we either see a preference for one of the carbons (in the cases of γ-HCH and (+)-α-HCH) or no preference is observed ((−)-α-HCH). This discrimination among carbon atoms is most evident for γ-HCH and can be also illustrated by the free energy plots obtained for dehydrochlorination of γ-HCH by LinA-type II (Figure ), where the abstraction from position C5 is clearly favored versus position C1. The process is more exothermic and presents a smaller activation free energy when abstraction takes place from position C5.

Lack of larger distinctive differences between preferences patterns for LinA variants based on conformational studies suggests that the source of LinA variants selectivity should be sought in specific interactions between HCH and the active site amino acids rather than in conformational changes of the protein, especially since it has been shown that the two types differ in their ability to catalyze proton abstraction from different positions.[26]

Interactions of the HCH Ligands with the First Shell Residues

In order to analyze HCH interactions with selected first shell residues, the SAPT2+ method has been used. As mentioned before, minimizing the influence of geometry differences between the models on the results was crucial for this part of our study; therefore, the models of enzyme–substrate complexes used for these calculations were obtained using molecular mechanics. Consequently, the results cannot be compared quantitatively, but they can be used to select those residues that may play an important role in enzyme–substrate interactions since the only differences between complexes result from the HCH spatial structure. Only those interaction energies with absolute values of calculated interaction energy (|ΔEint|) greater than 1.0 were considered as significant due to the benchmarked method error.[50]

As expected, the analysis showed important interactions between the substrate and H73 regardless of the enzyme type and HCH isomer (Figure ). Interestingly, significant differences between type I and type II complexes were observed in terms of HCH isomers interactions with residues reported as relevant for enantioselectivity (K20Q, L96C, A131G, and T133M).[25] Our results indicate that interaction with residue 20 is among the strongest ones in every type I complex while much smaller in type II complexes. Similarly, L96 contributes to favorable enzyme–substrate interactions in all type I complexes, unlike C96 in LinA-type II that does not interact with γ- and (+)-α- isomers and causes unfavorable interaction with the (−)-α- isomer. On the other hand, (−)-α-HCH is the only isomer that interacts significantly with residue 131 regardless of the A/G modification from type I to type II. Residue 133 interacts weakly with the (+)-α- isomer in the type I complexes but significantly with the (+)-α- and γ- isomers when type II complexes are considered.

Figure 7

Interaction energies ΔEint (in kcal mol–1) for the selected first-shell residues in a complex with the HCH isomers. Gray background corresponds to |ΔEint| < 1.0 considered to be insignificant, and sequence differences between type I and type II are framed with a purple box.

Additionally, modifications at positions 68, 71, and 113 cause changes in the interactions, in contrast to those of amino acids 23, 115, and 129. The last two are especially interesting since structural changes at these positions cause charge alterations; therefore, a more significant role of these mutations could be anticipated. However, the results are also consistent with previously published data showing that mutation at position 115 does not influence activity of the enzyme.[28] Furthermore, it should be remembered that R129 plays a role in the stabilization of the leaving chloride ion but does not interact directly with the substrate at the reactants state.[60] On the other hand, K20 is part of the same positively charged region but mutation at this position drastically changes the interactions, indicating the important role of this residue for the specificity of the enzyme types. Influence of K20 on the enzyme–substrate interaction was additionally confirmed by ΔEint calculation conducted for the models of K20–HCH where the amino acid was unprotonated (Table ). Alteration of K20 charge not only weakens its interaction with the substrate but results in unfavorable interaction energies. Importantly, a less significant effect was observed in the case of (+)-α-HCH complex which points to differences in the interaction pattern with the protein between this substrate and the two others, more efficiently metabolized by LinA-type I.[24,25,27]

Table 2

Interaction Energies, ΔEint (in kcal mol–1), Calculated for Protonated and Unprotonated K20 Residue

substrate	protonated K20	unprotonated K20
γ-HCH	–7.0	3.5
(−)-α-HCH	–3.0	5.3
(+)-α-HCH	–6.8	–0.9

Interestingly, the results for I44 and I109 are unexpected since these amino acids seem to play an important role in enzyme–substrate interactions even though none of them was previously mentioned as significant from the catalytic point of view. As for the rest of amino acids, the two isoleucine residues were selected for this analysis due to their proximity to the substrate.[23] The results obtained here unexpectedly demonstrate that these two residues may establish important interactions with the substrate with significant differences among the different enzyme–substrate complexes. I44 stabilizes γ-HCH more efficiently in type II than in type I, while the opposite behavior is observed for (−)-α-HCH. I109 plays a stabilizing effect with the three ligands in type I that almost vanishes in type II. These results suggest that I44 and I109 could be taken into consideration as important objectives during studies directed to improve the catalytic efficiency of LinA.

Conclusions

In the present work we have explored possible sources of lindane dehydrochlorinase, LinA catalytic efficiency using 0.5 μs classical MD, QM(DFTB)/MM simulations along with the adaptive string method as well as extensive analysis of resulting trajectories and conformations of studied systems.

On the basis of the analysis of the conformational dynamics of both LinA variants with and without ligand bound and the computed free energy of activation of dehydrochlorination reactions, we have found that the conformational changes of the protein, C-terminus extensive movement in particular, should not be linked to the protein catalysis as it had been suggested previously.[23]

On the other hand, in-depth analysis of the binding pocket allowed us to demonstrate substantial differences between two protein variants reflected in a smaller volume and larger rigidity of the active site of the type II protein. These results, when combined with the computed energetic barrier determined for the H/Cl pair elimination from different ES complexes, led to a very interesting observation; namely, the more rigid is the protein, the less efficient is the catalysis. This conclusion is also related to the fact that LinA-type II seems to be less preorganized for the chemistry to take place, which is reflected in substantial conformational mobility of HCH molecule in its active site (non-existent in LinA-type I) and much weaker stabilization of dehydrochlorination products, the chloride leaving group in particular. Type I has been shown to be preorganized in a much better way by the presence of three instead of only one H-donating residues for stabilization of the eliminated chloride anion.

In the present study we have shown that not only factors such as active site plasticity, volume, and their positive effects on LinA preorganization are important. By calculating interaction energies between the HCH ligand and the first shell residues in each protein variant, we have pinpointed those enzyme–substrate interactions that may play a significant role in catalysis and that have not been mentioned in any earlier study comprising LinA and its substrates.

These findings have a potential to constitute a good starting point for studies aiming at improving the catalytic efficiency of LinA, an environmentally important enzyme.