Computational mutagenesis reveals the role of active-site tyrosine in stabilising a boat conformation for the substrate: QM/MM molecular dynamics studies of wild-type and mutant xylanases

Molecular dynamics simulations have been performed for non-covalent complexes of phenyl  xylobioside with the retaining endo-  -1, 4-xylanase from B. circulans (BCX) and its Tyr69Phe 10 mutant using a hybrid QM/MM methodology. A trajectory initiated for the wild-type enzyme-substrate complex with the proximal xylose ring bound at the  1 subsite (adjacent to the scissile glycosidic bond) in the 4 C 1 chair conformation shows spontaneous transformation to the 2,5 B boat conformation, and potential of mean force calculations indicate that the boat is ~30 kJ mol  1 lower in free energy than the chair. Analogous simulations for the mutant lacking one oxygen atom 15 confirm the key role of Tyr69 in stabilizing the boat in preference to the 4 C 1 chair conformation, with a relative free energy difference of about 20 kJ mol  1 , by donating a hydrogen bond to the endocyclic oxygen of the proximal xylose ring. QM/MM MD simulations for phenyl  -xyloside in water, with and without a propionate/propionic acid pair to mimic the catalytic glutamate/glutamic acid pair of the enzyme, show the


Introduction
Glycosidases 1 (glycoside hydrolases, GHs) are involved in various physiological as well as pathological processes and 25 are thus important therapeutic targets. 2 Much attention has been focussed over the past decade on the design and synthesis of glycosidase inhibitors. Successful examples include the influenza neuraminidase inhibitors, Tamiflu and Relenza, 3 and the anti-diabetes -glucosidase inhibitors, 30 Acarbose and Miglitol. 4,5 Other targets under developments are anti-cancer inhibitors for mannosidase 6 and for nucleoside hydrolases/phosphorylases. 7 The family G/11 endo-1,4--xylanase (BCX) from Bacillus circulans catalyses the hydrolysis of xylan 8 by means of a 35 double-displacement mechanism (Fig. 1), via oxacarbeniumlike transition states, with net retention of configuration at the anomeric carbon. 9 Being well characterized experimentally and relatively small, it is an ideal subject for computational modelling. Its 3D fold, its substrate conformation, and its 40 catalytic residues have been identified through a combination of sequence analysis, inhibition experiments and structural determination. [10][11][12][13][14][15] However, many questions still remain concerning details of the reaction pathway and catalytic mechanism. In this paper we attempt to answer the following 45 question: why does the proximal sugar-ring of the substrate distort from 4 C1 chair to 2,5 B boat conformation?
Substrate distortion appears to be a general feature of glycoside hydrolases, [16][17][18][19] as seen by X-ray crystallography of enzyme-inhibitor complexes or of complexes with mutated enzymes, but less easily observed with reactive substrates, for which computational modeling provides a powerful 55 investigative tool. Many modelling studies have confirmed that substrate ring distortion is a common feature among glycosidases. [16][17][18][19] Molecular dynamics (MD) simulations have shown that the boat conformation at the 1 subsite is critical in the mechanism of family 18 chitinases, 16 and other studies have demonstrated that the 1 sugar moiety in cellulase 5 Ce16A from Trichoderma reesi adopts a skew-boat conformation. 17 Similarly, modelling studies of galactosidases have provided evidence of substrate distortion. 18,19 However, most studies have relied exclusively on classical approximations to describe the interatomic 10 interactions: electronic factors, which could be crucial for substrate conformational analysis and subsequently for the mechanistic behaviour, have been neglected. To model accurately reactions occurring in enzymes (and other condensed phase systems) it is necessary to consider a very 15 large number of interactions which may influence the mechanism. It is not feasible to use quantum mechanical (QM) methods to treat the dynamics of entire megadimensional systems, whilst molecular mechanics (MM) methods cannot describe the electronic polarisability that is 20 probably a key component of substrate distortion, as exemplified by the studies of Rovira and co-workers on glucanase 20 and glucose. 21 In this paper, we apply a hybrid quantum mechanics/molecular mechanics (QM/MM) 22 approach in conjunction with MD simulations. 25

Isolated substrates
Unconstrained geometry optimization was performed for the substrate, starting from each of the two conformations ( 2,5 B and 4 C1), using the B3LYP density functional method 23 and 6-30 31+G(d) basis set, as implemented in the Gaussian03 package. 24

Enzyme-substrate initial structure
The X-ray crystallographic coordinates (Fig. 2a) of the covalent enzyme-inhibitor complex of wild-type BCX (PDB 35 accession code 1BVV) 25 were modified as follows: (1) pKa values for amino-acid residues at the pH optimum 5.7 were estimated using the H ++ programme. 26 (2) Hydrogen atoms were added in accordance with the predicted protonation state. 40 (3) The whole system (16475 atoms) was energy minimized to a residual gradient of less than 0.001 kJ mol 1 Å 1 with a QM/MM method using the DYNAMO package. 27 The QM region (68 atoms) comprised the covalently-attached substrate -a disaccharide of xylose (XYL) and 2-deoxyfluoroxylose 45 (DFX) -and the two catalytic residues, Glu78 and Glu172, and was described by AM1 semi-empirical Hamiltonian. 28 The MM region contained the rest of the enzyme, and was described by the OPLS-AA potential. 29  along the CC bonds of Glu78 and Glu172 (Fig. 2c). At this stage no additional water molecules were present beyond the crystallographic waters.
(4) A gas-phase QM/MM MD simulation (NVT, 300 K, 5 ps) was performed to pre-equilibrate the protein and to allow the 60 substrate to accommodate in the binding cavity. (5) The F atom at position 2 of the proximal sugar moiety of the covalently-attached substrate was changed to OH to simulate the natural substrate, xylose; the system was QM/MM energy-minimized and equilibrated by MD once 65 again. (6) The covalent bond between the anomeric carbon of the (anomer of the) attached substrate and the oxygen of Glu78 was broken, and a phenoxy (OPh) leaving group was inserted manually to form the -anomer. This change required care to 70 retain the integrity of the enzyme configuration prior to the whole system being freely optimized again and equilibrated by MD, still without solvating waters. (7) The whole system was enveloped in a cubic box of TIP3P water 30 of side-length 55.5 Å. First, all water molecules were 75 relaxed with a gradient minimizer, while keeping the protein structure frozen; then the system was subjected to a short MD pre-equilibration using mild constraining forces to maintain the desired interactions between the substrate and the catalytic residues. Next,  20 ps at 150 K, using the NVT ensemble, still with the frozen protein. Finally, the whole system was freely minimized without any constraints and subsequently equilibrated for 20 5 ps at 300 K.

Sugar ring distortion
In order to understand the nature of ring distortion at the 1 subsite within the environment of the enzyme active site, it is necessary also to explore the conformational behavior of the 10 substrate in solution. To this end, we consider not only the final full system, as described in step 7 above, as model 3, but also two simpler models (Fig. 3): model 1 includes in the QM region only the proximal sugar ring of substrate in MM water, whereas model 2 additionally includes propionic acid and 15 propionate groups in the QM region, to represent the functionality of Glu78 and Glu172; models 1 and 2 both involve cubic boxes of TIP3P water with side-length 31.4 Å. Subscripts B and C below denote the 2,5 B boat and 4 C1 chair conformations, respectively. The QM region of models 1 and 20 2 does not involve any link atoms. Models 1 and 2 were allowed to equilibrate for 10 ps, and model 3 for 20 ps (same as above) of MD at 300 K in the NVT ensemble, prior to 30 ps of production MD for each model (1 fs timestep, switched cutoff radius of 16 Å applied 25 to all interactions). No constraints were applied at any stage of the dynamic runs. Atomic charges were monitored along the dynamic simulation. A free-energy pathway for interconversion of 2,5 B and 4 C1 conformers was generated using the potential of mean force (PMF) approach and the 30 weighted histogram analysis method (WHAM) as implemented in DYNAMO. 27 Before running PMF calculations, it was necessary to determine the best internal coordinate to describe this conformational change. Ideally the conformations of hexapyranose rings may be described by 35 combinations of distances and angles, such as the puckering coordinates of Cremer and Pople 31 used in a recent study of glucose, 21 but for simplicity we use a single coordinate. Five different dihedral angles were considered and the C3-C4-C5-O5 dihedral was found to be the best for the purpose. An 40 umbrella constraint force of 0.5 kJ mol 1 Å 2 was applied in each of the 31 windows; 1 ps of equilibration was performed prior to 10 ps of production MD in each window, starting with a structure perturbed in the distinguished coordinate from the final structure in the previous window. Block averaging 45 within each window suggests that the PMF is satisfactorily converged for its purpose in this study, although this does not mean that equilibrium has been achieved with respect to all hydroxyl group conformations.
Hexapyranose ring conformations were determined from 50 trajectories for each model and were described following the approach of Bérces et al. 32 Briefly, a conformation was specified in terms of a polar angle  (0° ≤  ≤ 180°), an equatorial angle  (0° ≤  < 360°), and an amplitude r.    local minimum, demonstrating that the boat is not a stable 10 conformer in the gas phase. The QM/MM MD trajectory (Fig. 6a) initiated from the 2,5 B conformation of phenyl -xyloside in water follows an itinerary through B3,0 to 2 S0 and 2,5 B over the course of 30ps, as shown by the plot of  vs. time; these conformations are 15 clustered in the region bounded by 70° <  < 110° and 120° <  < 45° (Fig. 6b). The corresponding trajectory initiated from the 4 C1 conformation in water remains stable: the apparent fluctuations in Fig. 6c and the spread of points in Fig. 6d reflect the ill-defined character of the angle  when  20  180°.

Effect of catalytic residues: model 2
The effects of the acid Glu172 and nucleophilic Glu78 5 residues may be mimicked simply by propionic acid and propionate moieties, respectively. The 2,5 B conformation with the additional groups located in positions taken from the structure of the enzyme-substrate complex optimizes to a local 2,5 B boat minimum in both gas-phase B3LYP/6-31+G(d) and 10 QM/MM aqueous-phase calculations. The QM/MM MD trajectory initiated from the 4 C1 conformation in water remains stable and behaves very similarly to that for model 1. The corresponding trajectory initiated from the 2,5 B boat in water fluctuates between this 15 conformation and the 2 S0 skew-boat during 30ps, as shown by the plot of  vs. time (Fig. 7a); these conformations are clustered in the region bounded by 75° <  < 105° and 135° <  < 90° (Fig. 7b). In comparison with model 1, the observed relative stability of 2,5 B boat may refer to the 20 stabilization by hydrogen bond between the OH group at C2 and the nucleophilic propionate residue (Fig. 8).

Effect of protein environment: model 3
The QM/MM MD trajectory (Fig. 9a and 9b) initiated from the 2,5 B conformer shows this to be a stable species which 25 does not evolve towards a skew conformation. The C2-C1-O5-C5 dihedral angle of the 2,5 B-substrate in the enzyme achieves a better coplanarity (0-30 o ) than in model 2.
Interestingly, after initial energy minimization, the 4 C1 chair conformation of the proximal sugar ring within the 30 enzyme active was not found as a local minimum but converted to an E1 envelope and then transformed quickly to   5 four nearby residues that interact directly with the xylose ring at the 1 subsite, namely Glu172 and Glu78, Arg112 and Tyr69 (Fig. 10). The acidic group Glu172 donates a hydrogen bond to the glycosidic O1atom of the substrate, whereas the nucleophilic group Glu78 accepts a hydrogen bond from the 10 OH group at C2 of the sugar ring, and Arg112 forms two hydrogen bonds with OH group at C3 of the substrate. We have monitored these hydrogen-bond distances along the 30ps MD trajectories for the chair and the boat conformations of the proximal sugar ring of the substrate, but have found 15 essentially no difference between the two conformers. On the other hand, Tyr69 may form a hydrogen bond with either the endocyclic ring oxygen O5 or with Glu78 (Fig. 11a): when one hydrogen-bond distance is short, the other is long (Fig.  11b, red and purple lines). In the QM/MM MD simulation 20 started from the 4 C1 chair, the H1 … O5 distance between Tyr69 and the substrate is initially > 3.5 Å and distinct from the short distance for the 2,5 B boat (Fig. 11c), but after 5 ps the two trajectories behave in the same way (Fig. 11b, blue line).
The remarkable change in the hydrogen bond between 25 Tyr69 and the proximal sugar ring suggests that this interaction is required to assist the conformational rearrangement from the 4 C1 chair to the 2,5 B boat in the enzyme active site. To test this hypothesis we made the Y69F mutant, replacing tyrosine by phenylalanine, in order to 30 eliminate this hydrogen bonding interaction with the substrate.
Using the same computational procedure as before, QM/MM optimization showed the 4 C1 conformation to be a local minimum, and subsequent MD simulation confirmed that it remained as a stable species during the 30 ps trajectory (  To quantify the extent of this stabilization, we computed free energy profiles (Fig. 13) for the interconversion of the 2,5 B boat and the 4 C1 chair for both the wild-type enzyme and Y69F mutant enzyme-substrate complexes, using the C3-C4- 10 C5-O5 dihedral as the reaction coordinate. The (Helmholtz) free changes along this coordinate are evaluated relative to the well-defined boat conformer with C3-C4-C5-O5  45° for both profiles. Consequently, the energetic stabilization of the boat conformer in the wild-type due to the hydrogen bond 15 between the substrate and Tyr69 appears in Fig. 13 as an apparent stabilization of the chair conformer in the mutant by  20 kJ mol 1 . Note that, for the wild-type, the dihedral angle C3-C4-C5-O5  60° does not correspond purely to a 4 C1 chair: the population of conformers sampled in the simulation 20 for values of the reaction coordinate in this range includes also envelope and half-chair conformations. It is therefore probably safer to estimate the degree of stabilization by comparison between the local energy minima at  60° for the mutant and  40° for the wild-type. Clearly, the barrier for 25 conversion of the 4 C1 chair to the more-stable 2,5 B boat in the wild-type enzyme-substrate complex is significantly lower than it is for the Y69F mutant.
The sequence of conformations populated by the wild-type BCX-substrate complex as it is driven along the C3-C4-C5- 30 O5 dihedral angle reaction coordinate in the PMF calculations is presented in Fig. 14 H5 and E3 on the way towards 4 C1. Note that this sequence of clusters does not define a dynamical trajectory since it is driven by the dihedral angle constraint and is not a 40 function of time. Also note that the "polar" 4 C1 region is distorted by the Mercator projection in regard to area and distance relative to the other sugar-ring conformations.

Development of oxacarbenium-ion character
Mulliken atomic charges for selected atoms (Table 1) averaged over the 30 ps MD trajectory for each of the three 55 models suggest a gradual increase in the oxacarbenium-ion character in the proximal xylose ring as the substrate interacts with the carboxyl/carboxylate pair in solution and within the wild-type enzyme active site. This is shown by increases in positive charge on C1 and H1 and a decrease of negative 60 charge on O5; these changes are accompanied by increased polarization of the C1-O1 bond as shown by an increase of negative charge on O1. Note that these charge changes are more pronounced for the boat than for the chair conformer. Moreover this development of oxacarbenium-ion character is 65 demonstrated by shortening and lengthening, respectively, of the average C1-O5 and C1-O1 bonds (Table 2); curiously, these changes are larger for the chair conformer. The effect of the Y69F mutation is to reduce the degree of oxacarbenium-ion character in the proximal xylose ring of the 70 enzyme-substrate complex, as evidenced by the charges and bond lengths shown in Tables 1 and 2. The non-bonded distances On .... C1 and Ha .... O1 in Table 2 indicate that greater oxacarbenium-ion character is accompanied by a shorter distance to the nucleophilic carboxylate group but a longer 75 distance to the acidic carboxyl group.

Discussion
The suggestion that a pyranoside substrate might preferentially adopt a 2,5 B conformation during enzymecatalysed glycoside hydrolysis was first made by Hosie and Sinnott 33 on the basis of kinetic isotope effects for yeast glucosidase catalysed hydrolysis of aryl glucosides and 5 glucosyl pyridinium ions. Their data required that breaking of the bond to the aglycone moiety of their substrates was preceded by a kinetically discrete non-covalent transformation of the initial enzyme-substrate complex, which they identified with a change to the 2,5 B conformation in which the C5-O5-10 C1-C2 atoms are approximately co-planar. The latter arrangement would facilitate the formation of a transition state with significant oxacarbenium-ion character.
In this paper we have taken one step toward investigating the reason for, as well as the structural implications of, 15 substrate ring distortion in -1,4-xylanase, by applying hybrid QM/MM molecular dynamics and free energy calculations. How reliable are our results, based as they are upon use of the semiempirical AM1 hamiltonian within the QM/MM method? Momany and co-workers 34 have shown that the B3LYP/6-20 311++G** level of theory gives consistently reliable geometries, conformations, and energies for carbohydrates in vacuo, and have noted the importance of hydroxyl group orientations and interactions. Intramolecular interactions between these groups are less important in aqueous solution 25 because of hydrogen bonding with solvent molecules, and in an enzyme active site interactions with the protein environment may drastically reduce the number of significant hydroxyl group rotameric forms. In regard to the relative energies of the 4 C1 and 2,5 B conformations of the xylose ring 30 in each of our condensed phase models 1, 2 and 3, one may well suppose that there is an underlying systematic error due to AM1 as compared to a large-basis DFT, or other high-level theoretical method, but it is not a straightforward task to estimate its magnitude. Moreover, our concern in this study is 35 not with the energy difference between the boat and chair conformers, but rather with the influence of enzyme environment upon that energy difference. We expect that the change from tyrosine to phenylalanine (both in the MM region) will polarize higher-level QM wavefunctions in a 40 qualitatively similar (though undoubtedly quantitatvely different) fashion to AM1. We therefore consider our AM1/OPLS estimate of the free energy of stabilization of the boat conformer due to the hydrogen bond with Tyr69 in BCX to be at least qualitatively reliable. A quantitative assessment 45 of error, as compared to an appropriate high-level QM method, would be meaningless unless it were evaluated for completely converged populations of conformations accessible under condensed-phase conditions, since there are large fluctuations between energy differences determined at 50 arbitrary "snapshot" structures taken from the MD trajectories.
The first structural evidence for sugar-ring distortion to the boat conformer was obtained by Brayer and co-workers 25 who obtained the 1.8 Å resolution structure of the covalent 55 glycosyl-enzyme intermediate formed between BCX and a 2deoxy-2-fluoro--xylobioside. (This structure was, of course, the starting point for the modelling studies described in the present paper.) The proximal xylose ring was found to be heavily distorted from 4 C1 to 2,5 B, thus achieving near 60 coplanarity of the C5-O5-C1-C2 atoms, whereas the distal xylose ring remained in the chair conformation. It was noted that an itinerary for the conformational rearrangement that were to proceed via 4 C1  2 H3  2 SO  2,5 B would involve motion (relatively) of only C5 and O5, and that these atoms 65 do not bear substituents which might add to the energetic cost of the conformational change. These workers also obtained the X-ray crystal structure for the catalytically inactive Y69F mutant BCX, with no substrate in the active site, and reported that it has a very similar structure to that of the covalent 70 intermediate. In particular, the position of Glu78 was almost identical, thus refuting an earlier proposal 35 that Tyr69 might have a crucial catalytic role in correctly positioning the nucleophile. On the other hand, it was noted 25 that the noncovalent enzyme-substrate complex observed by Wakarchuk 75 et al. 35 between the catalytically incompetent E172C mutant BCX and a xylotetraose substrate has the proximal xylose ring (occupying the 1 subsite) bound in the 4 C1 conformation.
What none of these prior experimental studies has been able to ascertain is the nature of the conformation of a xylose 80 substrate bound non-covalently in the 1 subsite of a catalytically competent BCX, but that is precisely what we have been able to achieve by means of the computational modelling carried out in the present study. Our QM/MM model for the enzyme-substrate complex between wild-type 85 BCX and phenyl -xylobioside clearly prefers the proximal ring to adopt the 2,5 B conformation. An MD trajectory initiated from a structure with this ring in the 4 C1 conformation spontaneously transformed into the 2,5 B boat. The PMF for interconversion of these conformers along the 90 dihedral angle C3-C4-C5-O5 coordinate suggests that the boat is lower in free energy than the chair site by ~30 kJ mol 1 within the active site of wild-type BCX.
We have demonstrated the key role of hydrogen bond donation from Tyr69 to O5 of the proximal xylose ring in 95 stabilizing the boat over the chair by means of a comparative QM/MM MD study for the Y69F mutant BCX which clearly shows the 4 C1 conformation to be a stable local energy minimum in the absence of that hydrogen bond. The PMF for interconversion of boat and chair along the dihedral angle C3-100 C4-C5-O5 coordinate suggests that the boat is still lower in free energy than the chair site by ~10 kJ mol 1 within the active site of the Y69F mutant, implying a preferential stabilization of the boat by wild-type BCX by about 20 kJ mol 1 .
The Y69F mutant has been reported 35 to have less than 0.01% of the activity of wild-type BCX at 40 C, which would imply a free energy difference between the rate-determining transition states of > 24 kJ mol 1 . It has been hypothesized 25,36 that the catalytic role of Tyr69 is provide a stabilizing dipolar 5 interaction with the partial positive charge on O5 of the oxacarbenium-ion like transition state for the rate-determining glycosylation of BCX by the substrate. We will present the results of QM/MM PMF calculations for the free energies of activation for the glycosylation step in wild-type BCX and the 10 Y69F mutant in a forthcoming paper, 37 and so will not speculate further on this point here. However, it is worth noting an alternative possibility that arises by analogy with the experimental 38 and theoretical 39 studies by Schramm, Schwartz and co-workers on human purine nucleoside 15 phosphorylase (hPNP). The glycosyl transfer reaction catalysed by this enzyme has a transition state with oxacarbenium-ion character and the ribofuranoside substrate possesses a hydroxyl group C5. It is suggested that the neighbouring His257 provides a mechanical push upon O5 20 towards the endocyclic O4 in a compressive motion with the phosphate nucleophile such that the build-up of electron density stabilizes the oxacarbenium-like transition state and facilitates the reaction. 38,39 Regardless of whether this particular idea is correct, it is striking to note the importance 25 of a suitably located hydroxyl group which exerts on effect upon the endocyclic oxygen of a glycoside by means of its oxygen rather than its proton. Of course, the xylose ring in the 1 subsite of BCX has no hydroxymethyl substituent at C5; indeed, it has been noted that there is no space around C5 to 30 accommodate any substituent. 25 Instead, however, the active site of wild-type BCX presents Tyr69 in close proximity to O5. We suggest that the OH groups of Tyr69 in BCX and of the hydroxylmethyl substituent of the ribofuranoside substrate of hPNP may have similar roles, although we do not wish to 35 speculate here on what exactly those roles may be.

Conclusion
Molecular dynamics simulations using a hybrid QM/MM potential have demonstrated that wild-type BCX preferentially binds a phenyl -xylobioside substrate in the 2,5 B 40 conformation at the 1 subsite adjacent to the scissile glycosidic bond. This result complements and extends the earlier experimental report that the proximal xylose ring of a covalently-bound 2-deoxy-2-fluoroxylobioside complex with BCX adopts the same 2,5 B boat conformation. Analogous 45 simulations for the Y69F mutant confirm the key role of Tyr69 in stabilizing the boat in preference to the 4 C1 chair conformation with a relative free energy difference of about 20 kJ mol 1 . 5 By donating a hydrogen bond to the endocyclic oxygen of the sugar ring, Tyr69 lowers the free energy of the 2,5B conformer.