Origin of the enantioselectivity of alcohol dehydrogenase

Alcohol dehydrogenases (ADH) are a family of enzymes that catalyse the interconversion between ketones/aldehydes and alcohols in the presence of NADPH cofactor. It is challenging to desymmetrise the substituted cyclopentane-1,3-dione by engineering an ADH, while the reaction mechanism of the metal independent ADH remains elusive. Here we measured the conversion of a model substrate 2-benzyl-2-methylcyclopentane-1,3-dione by Lb ADH and found it predominately gave the ( 2R, 3R ) product. Binding mode analysis of the substrate in Lb ADH from molecular dynamics simulations disclosed the origin of the enantioselectivity of the enzyme; the opening and closing of the loop 191-205 above the substrate are responsible for shaping the binding pocket to orientate the substrate, so as to give different stereoisomer products. Using QM/MM calculations, we elucidated the reaction mechanism of Lb ADH. Furthermore, we demonstrated the reaction profile corresponding to the production of different stereoisomers, which is in accordance with our experimental observations. This research here will shed a light on the rational engineering of ADH to achieve stereodivergent stereoisomer products.


Introduction
Alcohol dehydrogenases (ADH) can catalyse the reversible reduction of carbonyl compounds to yield alcohols in the presence of the cofactor NADPH.The desymmetrisation of substituted cyclopentane-1,3-dione was previously studied by chemical [1][2][3] and enzymatic catalysts 4,5 .Although high enantioselectivity was achieved by carbonyl reductase 4,5 , it is not known how the substrate would adapt itself in the catalytic pocket to give a certain stereoisomer product.Therefore, it is crucial to understand the preferred binding of the substrate in the enzyme in order to rationally engineer the enzyme to achieve stereodivergent products.
There are two types of ADHs, including Zinc containing and non-metal ADHs.The mechanism of Zinc-ADH catalysed reduction mechanism was reported previously 6,7 , however, to the best of our knowledge, the mechanism of the reduction of carbonyl compounds with metal-independent ADHs remains elusive.In this research, we expressed an ADH enzyme from Lactobacillus brevis (LbADH) and measured its conversion of the substrate.2-benzyl-2-methylcyclopentane-1,3-dione was chosen as a model substrate because the reduction would give four different enantiodivergent products and hence enable the desymmetrisation of the substrate.Our kinetic experiments showed that the (2R, 3R) and (2S, 3S) products were observed with a ratio of 9 to 1, while the other two chiral products (2S, 3R) or (2R, 3S) were not observed.
MD simulations disclosed that the loop 191-205 above the substrate could exhibit either open or closed conformation, which is responsible for the production of different stereoisomer products.Further hybrid quantum mechanics/molecular mechanics (QM/MM) calculations disclosed the mechanism of the reduction catalysed by LbADH.We elucidated that the reaction starts with the rate-limiting step, i.e., the concerted hydride transfer from NADPH to the substrate carbonyl carbon and proton transfer from Tyr155 to the substrate carbonyl oxygen.The reaction then proceeds with a concerted proton relay; deprotonated Tyr155 phenolate abstracts a proton from NADP+ ribose, which in turn abstracts a proton from the nearby catalytic acid Lys159.The second step is barrierless for the major product (2R, 3R), but it needs to overcome a notable barrier for the minor product (2S, 3S) and an even higher barrier for the inaccessible product (2S,3R) or (2R,3S), implicating the origin of the enantioselectivity of the enzyme.Additionally, we found the reaction coordinate in the reactant (i.e., the distance between the NADPH hydride to the carbonyl carbon) is correlated to the barrier for the rate-limiting step.

Protein expression
100 μL stored bacteria containing LbADH gene in pET-28a was first inoculated in 5 mL LB medium (containing 50 μg/mL kanamycin) and was shaken at 200 rpm overnight as preculture.The preculture was used to inoculate large culture (250 mL LB + 50 μg/mL kanamycin in 1 L shake flasks) at 37 °C for about 4 hours until OD600 at 0.6-0.8.After cooling at 4 °C for 0.5 h, 1 mM isopropyl β-thiogalactopyranoside (IPTG) was added to induce LbADH expression.The culture was allowed to express at 18 °C for 18 h with shaking at 200 rpm.Then, cells were harvested by centrifugation at 9000 rpm and 4 °C for 10 min, and the supernatants were discarded.The cells were resuspended in 50 mM Tris-HCl Buffer (pH 7.5) containing 100 mM NaCl and 2 mM MgCl2 (1 g wet cell in 5 mL 50 mM Tris-HCl Buffer) and stored in -80 °C.The cells were repeated freezing and thawing for 3 times, and then released the target proteins by sonication.The cell debris was removed by centrifugation at 10,000 rpm for 20 min at 4 °C, and the crude enzyme solutions were stored at -80 °C.

Protein purification
The crude enzyme solution was filtered and loaded on a GE Healthcare HisTrap FF Crude column (5 mL) preequilibrated with 50 mM Tris-HCl buffer (pH 7.5) containing 100 mM NaCl, 2 mM MgCl2 and 25 mM imidazole.The enzyme was eluted by 50 mM Tris-HCl buffer containing 100 mM NaCl, 2 mM MgCl2 and 250 mM imidazole.The proteins were dialyzed by 50 mM Tris-HCl buffer (pH 7.5) containing 100 mM NaCl and 2 mM MgCl2 for 12 h at 4 °C.

Determination of kinetic parameters
Kinetic parameters of LbADH were measured by monitoring the decrease of NADPH following the absorbance at 340 nm using the spectrophotometry method (SHIMADZU-UV-2550).The activity assay was performed in a mixture (3 mL) containing 50 mM Tris-HCl buffer (pH 7.5), an appropriate amount of the enzyme and varying concentration of 2-benzyl-2-methylcyclopentane-1,3-dione (0-8 mM) with 5% (final) dimethyl sulfoxide (DMSO) as a cosolvent.NADPH was added to a final concentration of 100 μM as the last component to start the reaction.

Molecular Docking
In the docking studies, the crystal structure of wild-type Lactobacillus brevis alcohol dehydrogenase (LbADH) (PDB: 1ZK4) 8 was employed, which contains the NADPH cofactor and the acetophenone substrate.The substrate acetophenone was removed while the NADPH cofactor was retained to construct the apo-protein.A 100-ns MD simulation was conducted for the apo-protein using AMBER 20 9 , including energy minimization and equilibration steps.Representative structures were selected by cluster analysis of MD trajectories using CPPTRAJ 10 module of the Amber 20 for the docking study.The model substrate in this study 2-benzyl-2-methylcyclopentane-1,3-dione was constructed using GaussView 6 11 and optimized at the B3LYP/6-31G(d) level.The optimized substrate was docked in the simulated apo protein using Autodock 12 and then the top docked poses with favourable binding energies and appropriate orientation (that would lead to the respective stereoisomeric ketols) were selected for the subsequent MD simulations.

Molecular dynamics (MD) simulations
The parameters for the substrate and the NADPH co-factor were calculated based on the optimized geometry at the B3LYP/6-31G(d) level of theory using Gaussian 16 13 with the PCM solvation model.Electrostatic potential (ESP) was derived at the same level according to the Merz-Singh-Kollman scheme and was fit using RESP in the antechamber module of AMBER 20 9 coupled with the general AMBER force-field GAFF2 14 .The H++ server 15 was used to predict the protonation states of titratable amino acids and the FF19SB 16 force field in AMBER 20 was used to obtain the parameters for the enzyme.
The enzyme-substrate complexes were solvated in a pre-equilibrated cuboid box of TIP3P 17 water molecules, and any protein atom was at least 10 Å to from the edge of the cuboid box.
The system was then neutralised by adding 8 Na + counterions using the tLEaP program implemented in AMBER 20.A harmonic restraint force constant of 100 kcal −1 Å −2 was applied to the solute molecules and ions to minimise the solvent molecules, followed by 1,000-step steepest descend and 1,000-step conjugate gradient unrestrained minimisation.A cut-off of 10 Å was applied for non-bonded Lennard-Jones potential and electrostatic interactions.Hydrogen bonds were constrained using the SHAKE algorithm during all MD simulations.Under constant volume and periodic boundary conditions, a progressive heating was performed from 0 K to 300 K in 100 ps (5000 steps with a step size of 0.02 ps), followed by 1 ns equilibration using NPT ensemble at 300 K.A harmonic restraint of 5 kcal mol -1 was applied to the solute at both equilibration stages.After the equilibration stages, a 200 ns production MD simulation was run using an NPT ensemble at 300 K and 1 bar.Three replicas of MD simulations were run for each system.

Hybrid quantum mechanics/molecular mechanics (QM/MM) calculations
The snapshots used for the QM/MM calculations were obtained from the CPPTRAJ cluster analysis 10 of the MD trajectory, and then three cluster structures of each system were optimized using 2,500-step steepest descend and 2,500-step conjugate gradient algorithms, with a harmonic restraint force constant of 100 kcal -1 Å -2 applied on the heavy atoms.
The substrate, truncated NADPH, Tyr155, Lys159, and Ser142 are included in the QM region.During the geometry optimization, residues within 6 Å of the QM region were allowed to move freely, while the remaining residues were kept frozen (Figure S1).The hydrogen link atoms were used to saturate the dangling bond at the QM/MM boundary.The QM/MM calculations were performed using the electronic embedding model within the Gaussian 13 ONIOM.Geometries of all intermediates and transition states were optimized at the B3LYP/6-31G(d) level.Frequency calculations within the harmonic approximation were used to verify the nature of all intermediates and transition states.Free energy and enthalpic corrections were carried out by computing harmonic frequencies analytically at 298.15 K.For each step on the reaction profile, thermochemical correction terms δEG were carried out as a difference of the reaction energy (ΔEB3LYP/6-31G(d)) and the corresponding free energy (ΔGB3LYP/6-31G(d)): B3LYP, M06-2X and the ωB97X-D functionals have been extensively used in previous studies on dehydrogenase with NADPH as the cofactor [18][19][20] .These functionals were also used in this study.To account for the dispersion effect, calculations were also performed using B3LYP-D3(BJ) and M06-2X-D3.
Potential energy surface scans were conducted to locate the transition states at the same level, i.e., B3LYP-D3(BJ)/6-31G(d).Transition states were confirmed through visual inspection of the imaginary frequency modes as well as intrinsic reaction coordinate (IRC) calculations.
The corrected free energies (∆G) were calculated as follows: where the electronic energy ∆ESP was computed at the B3LYP-D3(BJ)/6-31G(d)) level.
The single point energy of the optimized intermediates and transition states were also compared using the large basis set, i.e., 6-311+G(2d,2p) and def-2TZVPP.

Results and discussion
Enantioselectivity of Alcohol Dehydrogenase 2-benzyl-2-methylcyclopentane-1,3-dione was chosen as a model substrate for LbADH, because the reduction would introduce desymmetrisation of the substrate.The wild-type LbADH was expressed and purified and then kinetic parameters of the enzyme were measured by monitoring the consumption of NADPH.
The stereoisomeric distribution of the products for WT-LbADH was measured as following.
Sodium phosphate buffer (1 ml, 100 mM, pH 7.0) containing 10% DMSO (v/v), 10 mM substrate, 30 mM glucose, 0.5 mM NADP+, 0.5-2 mg of LbADH and 4 U GDH were shaken at 30 ℃ for overnight.Then the reaction was stopped by the addition of an equal amount of ethyl acetate.The organic phase was separated and the solvent was removed; the resulting sample was analysed by chiral HPLC to determine the enantiomeric excess value of the alcohol products.
The reduction of the substrate yielded 90% and 10% stereoisomeric (2R,3R) and (2S,3S) products, respectively, whereas no (2S, 3R) nor (2R,3S) product was obtained (Figure 1, Scheme 1).A catalytic tetrad comprised of Asn113, Ser142, Tyr155 and Lys159 was proposed in the previous crystal structure study of LbADH 8 .However, the relative distances between Ser142, Lys159 and Tyr155 are not ideal for the reaction to take place, indicating the crystal structure may not be catalytically active conformation.The binding mode of the substrate 2benzyl-2-methylcyclopentane-1,3-dione in the catalytic site of the enzyme was decided by combining molecular docking and molecular dynamics simulations.Molecular docking was performed and the substrate poses that would lead to different stereoisomer products were retained, including the most dominant (2R,3R) and the less dominant (2S,3S) product, as well as the (2S, 3R) which was not obtained by the WT enzyme.These enzyme-substrate complexes were then subjected to MD simulations.RMSD analysis showed the simulations reached equilibrium (Figure S2).
In the proximity of catalytic sites of all the structures, Ser142 stabilizes the substrate carbonyl group by a hydrogen bond (Figure 2, Figure S3); Tyr155 is positioned toward the substrate carbonyl to serve as an H-bond donor and also forms H-bond with the NADPH ribose hydroxyl group, which is, in turn, H-bonded to Lys159.The binding of the substrate in the enzyme is mediated by Glu144, which results in different substrate orientations.In the substrate pose that leads to the main (2R, 3R) product, Glu144 forms a hydrogen bond with Tyr189 (Figure 2A, Figure S3), which is stabilized by the methyl group of Met205.The benzene ring of the substrate is nested between the surrounding Ala93 and Leu194, forming hydrophobic interactions with them.Thr192 forms H-bond with the amide NH2 of the NADPH nicotinamide.The loop 191-205 adjacent to Tyr189 displays a closed conformation, with Asp196 interacting with Lys191 by ionic interactions, which in turn tethered by Glu202 (Figure 2).
In the substrate pose that leads to the minor (2S, 3S) product (Figure 2B), the sidechain of Glu144 is turned away from the catalytic pocket.The loss of the H-bond between Glu144 and Tyr189 caused the loop K191-M205 to display a significant conformational change to adopt an open conformation.The opening of the loop makes Leu194 become far away from the substrate, such that the substrate was able to move freely in the binding pocket and eventually stabilized by the favourable hydrophobic interactions with the surrounding Leu152 and Ala193 with its benzene ring.The opening of the loop also caused the loss of the H-bond between Thr192 and NADPH.Meanwhile, Asp196 on the loop becomes exposed and its interaction with Lys191 is lost, leaving the latter interacting only with Glu202.A water flux moved in the catalytic site with the opening of the loop (Figure 2B).
Thus we demonstrated that the loop191-205 region is largely stabilized in the substrate pose leading to the dominant (2R, 3R) product (Figure S2B), compared to other substrate poses that lead to minor product (2S, 3S) or do not give the corresponding chiral product such as (2S, 3R).

Reaction profile of enantioselective reduction by ADH
So far, the reaction mechanism of NADPH-dependant ADH remains elusive.To elucidate the reaction mechanism of the LbADH and the reaction profiles of ADH for desymmetrisation of the substituted cyclopentane-1,3-dione to produce different stereoisomers, the representative structures from MD simulations of the enzyme-substrate complexes were minimized and the minimized structures are similar to the initial representative structures from cluster analysis based on their RMSD values (Supplementary information Table 1).These structures were selected as the starting structures for the subsequent mechanism study by QM/MM calculations.
LbADH is homologous to RasADH 3 and there is 5 AA difference between the two homologous enzymes.Previously, Chen et al. studied RasADH and suggested that Ser137, Tyr150, and Lys154 (corresponding to Ser142, Tyr155, Lys159 in LbADH) may participate in the catalytic process, along with the nicotinamide ring of NADPH 3 .However, the exact catalytic process of ADH is not known.
From our MD simulations, Tyr155 and Lys159 are located in the proximity of the substrate and NADPH; therefore, these two residues were included in the QM region, along with the substrate, the NADPH cofactor.Tyr155 and Lys159 were truncated to keep the sidechain and NADPH was truncated to keep the nicotinamide riboside part in the QM/MM calculations.
Potential energy surface scans were conducted for the enzyme-substrate complex structures that would lead to different stereoisomers, i. e. the most dominant (2R, 3R) (90% distribution among all products), the less dominant (2S,3S) (10% distribution), and also a chiral product (2S, 3R) that was not observed from the kinetic experiments.
Interestingly, for the substrate pose that leads to the minor product (2S,3S) or the unachievable product (2S, 3R), the deprotonated Tyr155 in I1 needs to overcome a notable barrier to be reprotonated (Figure S4b & Figure S5).
The dispersion in the DFT method may affect the energy barrier [21][22][23] .For example, QM/MM calculations for cytochrome P450 catalysed reactions showed that dispersion correction may reduce the barrier of hydrogen abstraction significantly by around 5 kcal/mol 21 .To examine if dispersion is necessary to consider the polarisation effect for ADH, we conducted bench marking calculations with dispersion effect using B3YLP-D3(BJ) and M06-2X-D3 with the 6/31G* basis set.Since M06-2X implicitly includes dispersion energy, as expected, including Grimme's correction did not cause significant change, whereas notable differences were observed with the B3LYP functional (Figure 3 & S4), highlighting the importance of including dispersion for B3LYP.

Reaction mechanism of ADH catalyzed reduction
Our QM/MM calculations disclosed that the ADH catalyzed reduction reaction proceeds with a hydride transfer from the C4 methylene of NADPH nicotinamide to the carbonyl carbon atom of the substrate that is accompanied by the proton transfer from Tyr155 to the substrate carbonyl oxygen, to give the reduced product and deprotonated Tyr155 (I1) (Scheme 2).The reaction proceeds with an abstraction of a proton from the nicotinamide riboside hydroxyl by Tyr155 to resume its protonated state while the deprotonated sugar abstracts a proton from the nearby Lys159, which functions as a catalytic acid (Scheme 2).Notably, the previous hypothesis for RasADH 4 proposed that the reprotonation of Tyr150 is fulfilled by the proton transfer from Lys154 (Lys159 in LbADH) to the phenolate oxygen of Tyr150 (Tyr155 in LbADH).Our calculations disclosed that the reprotonation of Tyr155 is facilitated by the nicotinamide riboside OH with the aid of the nearby lysine residue Lys159, thus presenting a complete reaction profile of the ADH-catalyzed reduction.
and observed a linear correlation (Figure 6).This correlation was consistently observed across all the optimization methods employed including B3LYP, B3LYP-D3(BJ), M06-2X, M06-2X-D3, and ωB97X-D.Thus, the reaction coordinate for the rate-determining step in the substrate can be directly used to estimate the reaction energy barrier and predict the reaction kinetics without the need to conduct a complete PES scan using computingdemanding QM/MM calculations.

Figure 1 .Scheme 1
Figure 1.The stereoisomeric distribution of the products from WT-LbADH catalyzed

Figure 2 .
Figure 2. MD simulated LbADH in complex with the substrate (A) Substrate in closed loop

Figure 4 .
Figure 4. Free energy profile for the enzyme-substrate complexes with the substrate poses

Figure 6
Figure 6 Energy barrier for the rate-limiting step is linear correlated with the reaction