Dynamic Preference for NADP/H Cofactor Binding/Release in E. coli YqhD Oxidoreductase

YqhD, an E. coli alcohol/aldehyde oxidoreductase, is an enzyme able to produce valuable bio-renewable fuels and fine chemicals from a broad range of starting materials. Herein, we report the first computational solution-phase structure-dynamics analysis of YqhD, shedding light on the effect of oxidized and reduced NADP/H cofactor binding on the conformational dynamics of the biocatalyst using molecular dynamics (MD) simulations. The cofactor oxidation states mainly influence the interdomain cleft region conformations of the YqhD monomers, involved in intricate cofactor binding and release. The ensemble of NADPH-bound monomers has a narrower average interdomain space resulting in more hydrogen bonds and rigid cofactor binding. NADP-bound YqhD fluctuates between open and closed conformations, while it was observed that NADPH-bound YqhD had slower opening/closing dynamics of the cofactor-binding cleft. In the light of enzyme kinetics and structural data, simulation findings have led us to postulate that the frequently sampled open conformation of the cofactor binding cleft with NADP leads to the more facile release of NADP while increased closed conformation sampling during NADPH binding enhances cofactor binding affinity and the aldehyde reductase activity of the enzyme.


Introduction
YqhD is a homo-dimeric, NADP(H)-dependent E. coli oxidoreductase identified in 2003 [1]. Since its discovery, it has received considerable attention for bioengineering efforts [2], being targeted for its utility in biomass conversion to biofuels [3] and chemical feedstocks [4]. The enzyme was initially identified during the E. coli structural genomics program and was shown to be widely distributed throughout the bacteria kingdom [5]. YqhD enzyme works as a homodimer with Zn 2+ in the active site, catalyzing the interconversion of alcohols and aldehydes with NADP and NADPH, respectively, as a cofactor. Although YqhD binds Zn 2+ , it has structural similarities to the iron-dependent (group III) alcohol dehydrogenase enzymes [6]. Initially, YqhD was proposed as NADP-dependent alcohol dehydrogenase with a weak affinity toward short-and medium-chain alcohols [7]. Later studies characterized this enzyme as a NADPH-dependent aldehyde reductase with a broad range of substrates, such as 3-hydroxypropionaldehyde [8], propanaldehyde, isobutylaldehyde [3], acrolein, and malondialdehyde [2,9]. The enzyme's biological role was also evident in the reduction of various reactive aldehydes derived from membrane lipid peroxidation [9] and mitigation of furfural toxicity [10][11][12]. Other biosynthetic efforts with YqhD have involved its utilization in the production of aromatic alcohols [13], ethylene glycol, 1-butanol [14], 1,4-butanediol, acetol [15,16], 1,2-propanediol [17] and, in a pathway using CO 2 , 1-butanol [18].

Figure 1. (a) YqhD dimer (comprised of monomers A and D) is in cartoon representation with the α-helical and
Rossmann-type domains colored in blue and pink, respectively. The shaded region shows the cofactor-binding site of monomer D. Zn 2+ is located at the interdomain face within the cofactor-binding site and is shown in orange colored VDW representation. The NADP cofactor is in licorice representation colored by elements (nitrogen, oxygen, hydrogen, carbon and, phosphorus in blue, red, white, carbon, and tan color, respectively). (b) NADP cofactor is present in the crystal structure as NADPH(OH)2 (with modification at C5 and C6 position in the nicotinamide moiety colored by elements and N and O atomic labeling) and in the MD simulations as NADPH and NADP in blue and red color, respectively. The black ovals highlight the nicotinamide moiety.

Structural and Dynamical Properties of YqhD Enzyme
As mentioned in the Introduction section, the crystal structure has a modified NADPH(OH)2 cofactor, but in this work, simulations were performed for the YqhD homodimer with oxidized NADP and reduced NADPH cofactors (separately) for comparison (see Figure 1 for details). First, we performed an equilibration step involving 20 ns MD simulation in isothermal-isobaric (NPT) ensemble with each state. During the equilibration step, the YqhD protein adopts a conformation suitable for the binding of NADP Figure 1. (a) YqhD dimer (comprised of monomers A and D) is in cartoon representation with the α-helical and Rossmanntype domains colored in blue and pink, respectively. The shaded region shows the cofactor-binding site of monomer D. Zn 2+ is located at the interdomain face within the cofactor-binding site and is shown in orange colored VDW representation. The NADP cofactor is in licorice representation colored by elements (nitrogen, oxygen, hydrogen, carbon and, phosphorus in blue, red, white, carbon, and tan color, respectively). (b) NADP cofactor is present in the crystal structure as NADPH(OH) 2 (with modification at C5 and C6 position in the nicotinamide moiety colored by elements and N and O atomic labeling) and in the MD simulations as NADPH and NADP in blue and red color, respectively. The black ovals highlight the nicotinamide moiety.
It has been observed previously that differences in cofactor oxidation states induce protein conformational changes [25,28] and also influence the cofactor binding affinities [29]. Along these lines, structural differences in cofactor binding for FucO oxidoreductase, a NADH-dependent enzyme with similarities to YqhD, were recently reported [30], while conformational changes in cofactor were shown to be critical in the pathway toward the transition state for horse liver alcohol dehydrogenase [31]. Recent work on product release, a rate-limiting step in dihydrofolate reductase (DHFR), indicates that NADPH induces product release through steric repulsion during conformational sampling of the closed excited state [32]. There is some variation in conformational dynamics (rates between exchange of states) when the DHFR-product complex is bound to oxidized (NADP) vs. reduced (NADPH) cofactor, but conformational exchange rates measured by nmR are on the same order of magnitude (1890 ± 80 s −1 bound to NADPH vs. 1420 ± 70 s −1 when bound to NADP + ).
The present work has gone beyond structural effects to computationally study how dynamics are affected by cofactor oxidation/reduction. To gain insight into dynamical preference of the YqhD protein towards NADPH over NADP cofactor, the article is organized as follows: the details of MD simulations are reported in the Methods section; the results section has the analysis of MD trajectories, with the main focus on structural and dynamical properties of YqhD enzyme and NADP/H cofactor binding; next, we discuss our results, focusing on the effect of cofactor oxidation/reduction on structure-dynamics and activity of the enzyme.

Structural and Dynamical Properties of YqhD Enzyme
As mentioned in the Introduction section, the crystal structure has a modified NADPH(OH) 2 cofactor, but in this work, simulations were performed for the YqhD homodimer with oxidized NADP and reduced NADPH cofactors (separately) for comparison (see Figure 1 for details). First, we performed an equilibration step involving 20 ns MD simulation in isothermal-isobaric (NPT) ensemble with each state. During the equilibration step, the YqhD protein adopts a conformation suitable for the binding of NADP and NADPH cofactors. Starting with equilibrated structures, we performed five 200 ns long independent MD simulations for each state, starting with different initial velocities (generated through a random number seed) to check the structural convergence. First, analysis of structure and dynamics was done using the crystal structure as reference for three different types of structural data subsets within each state of the NADP/H-cofactor: the entire dimer, the monomers comprising dimer, and the domains (α-helical and Rossmann-type domains) within each monomer, in addition to dynamics during the simulations (see Figure 1 for cofactor-bound crystal structure of YqhD).

YqhD Dimer
During the simulations, the populated dimer conformations have an average root mean square deviation (RMSD) of 0.28 ± 0.05 nm and 0.24 ± 0.04 nm when bound with NADP and NADPH cofactors, respectively. Radius of gyration (Rg) of dimer for the crystal structure is 2.71 nm. During the simulations, the YqhD dimer has conformations with an average 2.75 ± 0.01 nm Rg, which is slightly less compact than the crystal structure and is consistent with solution phase dynamical sampling. The cofactor binding and protonation states affect the dimer compactness during the simulations resulting in a less compact structure with NADPH and NADP cofactors than with NADPH(OH) 2 .

YqhD Domain
A Rossmann-type domain is involved in NADP cofactor binding. With NADP-bound, the domain shows an average RMSD of 0.17 ± 0.03 nm with slightly higher variations than NADPH-bound (average RMSD of 0.15 ± 0.02 nm). The α-helical domain is a metalbinding domain and has an average RMSD of 0.11 ± 0.02 nm for NADP-bound and 0.13 ± 0.03 nm for NADPH-bound YqhD. Monomers bound to reduced cofactor show slightly higher deviations than monomers bound to oxidized cofactor. Both domains within each monomer become slightly more compact in the simulations with Rg values of 1.57 ± 0.01 for Rossmann-type and 1.63 ± 0.01 for α-helical domains during the simula- tions, compared to 1.60 nm for Rossmann-type and 1.65 nm for α-helical domains in the crystal structure.

YqhD Monomer
When comparing monomers in the 2 µs aggregated trajectory, it can be seen that the oxidation state of the NADP/H cofactor affects the overall monomer conformation significantly. NADP-bound monomer shows an average backbone RMSD of 0.27 ± 0.06 nm with slightly higher variations than NADPH-bound monomer, which has an average RMSD of 0.23 ± 0.05 nm. The distribution of backbone RMSD and Rg values for YqhD monomers is reported in Figure 2. YqhD monomers bound to NADP show higher deviation from the reference crystal structure. The RMSD distribution ranges from 0.13-0.45 nm, with a main peak at 0.31 nm and two additional peaks at 0.23 and 0.17 nm. NADPH-bound YqhD has a higher probability of remaining closer to the reference structure, showing a main peak at 0.19 nm and two additional sharp peaks at 0.25 and 0.15 nm; however, the RMSD distribution is broader, ranging from 0.09-0.47 nm. The Rg value of monomer is 2.07 nm in the crystal structure. NADP-bound monomers show an average Rg of 2.14 ± 0.02 nm, while NADPH-bound monomers have an average Rg of 2.11 ± 0.03 nm. For NADPHbound monomers, Rg values show a broad distribution ranges from 2.04-2.27 nm with three peaks at 2.07, 2.10 and 2.13 nm. NADP-bound monomer exhibits a narrower Rg distribution, ranging from 2.04-2.25 nm, with a main peak at 2.15 nm and a wider shoulder at 2.09 nm. Overall, YqhD monomer bound with reduced NADPH cofactor populates conformations (two peaks) with lower RMSD (<0.2 nm) and Rg (<2.1 nm) values, i.e., close to the crystal structure that has the modified NADPH(OH) 2 cofactor. Thus, some of the structural differences in YqhD monomer relative to the crystal structure can be attributed to the changes in the oxidation states of the cofactor (reflecting relevant, functional states of NADP/H bound to the enzyme); the molecular details are discussed later.   Figure 3a shows the per-residue backbone RMSD for the YqhD domains illustrated in Figure 3b. High deviations are restricted mainly to loop regions and the N-and Cterminus. Both the monomers show significant deviations in the loops that are present at the dimer interface (between monomers) and the interface between domains within each monomer (marked by the cyan colored horizontal bar in Figure 3a). Monomers bound to NADPH show a higher deviation in loop regions α6/α7, α7/α8, α9/α10 and α12/α13 of the metal-binding helical domain then the NADP-bound monomer. It is worthwhile to keep in mind that the RMSD and Rg values are calculated using the crystal structure as a reference, which may not reflect the functional enzyme structure due to its modified cofactor (NADPH(OH) 2 ). Still, the crystal structure gives a good starting point to witness how the enzyme samples different conformations during the binding of oxidized and reduced cofactors. In the simulations, YqhD shows slightly higher structural flexibility, populating diverse conformations with a more compact yet broader, distinctive distribution of Rg values when bound to NADPH (required for reductase function) compared to NADP (released for enzyme turnover). The population of conformations towards lower RMSD and Rg values in the simulations with reduced cofactor indicates the slight change in protein conformation in NADPH binding relative to NADP binding. Figure 3a shows the per-residue backbone RMSD for the YqhD domains illustrated in Figure 3b. High deviations are restricted mainly to loop regions and the N-and Cterminus. Both the monomers show significant deviations in the loops that are present at the dimer interface (between monomers) and the interface between domains within each monomer (marked by the cyan colored horizontal bar in Figure 3a). Monomers bound to NADPH show a higher deviation in loop regions α6/α7, α7/α8, α9/α10 and α12/α13 of the metal-binding helical domain then the NADP-bound monomer.  Figure 3a shows the per-residue backbone RMSD for the YqhD domains illustrated in Figure 3b. High deviations are restricted mainly to loop regions and the N-and Cterminus. Both the monomers show significant deviations in the loops that are present at the dimer interface (between monomers) and the interface between domains within each monomer (marked by the cyan colored horizontal bar in Figure 3a). Monomers bound to NADPH show a higher deviation in loop regions α6/α7, α7/α8, α9/α10 and α12/α13 of the metal-binding helical domain then the NADP-bound monomer.

Hydrogen Bonding in YqhD Enzyme
Inter-monomer and interdomain hydrogen bonds were calculated in the crystal structure and during the simulations of YqhD bound to NADP/H. Hydrogen bonds were calculated using a distance between acceptor and hydrogen donor of 3.5 Å and an angle of <30° among acceptor, donor, and hydrogen as criteria. The dimer has ~8 inter-monomer hydrogen bonds, and monomers A and D have ~7 interdomain hydrogen

Hydrogen Bonding in YqhD Enzyme
Inter-monomer and interdomain hydrogen bonds were calculated in the crystal structure and during the simulations of YqhD bound to NADP/H. Hydrogen bonds were calculated using a distance between acceptor and hydrogen donor of 3.5 Å and an angle of <30 • among acceptor, donor, and hydrogen as criteria. The dimer has~8 inter-monomer hydrogen bonds, and monomers A and D have~7 interdomain hydrogen bonds during the simulations (averaging over data collected every 10 ps). The hydrogen bonds discussed below are present for a cumulative total of >30% of the time points analyzed in the trajectories.

Cluster Analysis of YqhD Enzyme
Conformational differences observed due to the cofactor oxidation/reduction state were further quantified using cluster analysis on monomers (see Section 4.4 for details). In the combined trajectory, a total of 54 (NADP) and 43 (NADPH) clusters were obtained for the monomer using an RMSD cutoff of 0.13 nm. Figure 4 represents the first three clusters of NADP-bound and NADPH-bound monomers in red, blue and green colored ribbons, respectively. These first three clusters comprise 62.7% (30.8%, 19.3%, and 12.6%) and 72.9% (46.0%, 13.8%, and 13.1%) of the total monomer conformations occupied by YqhD for NADP-and NADPH-bound enzyme, respectively.  By looking at the superimposed mean cluster structures and crystal structure illustrated in Figure 4, it can be seen that the primary differences in monomer conformations arise from the position of the β6/β7 loop region (circled) of the Rossmann-type domain and the α8/α12 helix of α-helical domain. Figure 5a shows the superimposed crystallographic structures of apo-and holo-proteins. The β6/β7 loop region is in open conformation in the apo form of the YqhD crystal structure [7] (and closed in the holoenzyme), indicating that they comprise an interdomain cleft which opens and closes for cofactor binding and release. In this side-by-side comparison of the monomer, it can be seen that, By looking at the superimposed mean cluster structures and crystal structure illustrated in Figure 4, it can be seen that the primary differences in monomer conformations arise from the position of the β6/β7 loop region (circled) of the Rossmann-type domain and the α8/α12 helix of α-helical domain. Figure 5a shows the superimposed crystallographic structures of apo-and holo-proteins. The β6/β7 loop region is in open conformation in the apo form of the YqhD crystal structure [7] (and closed in the holoenzyme), indicating that they comprise an interdomain cleft which opens and closes for cofactor binding and release. In this side-by-side comparison of the monomer, it can be seen that, with both cofactors, these loops are significantly dynamic. The highly populated conformations of NADPH-bound involve both monomers within the dimer remaining with a partially closed cleft (cluster 1), closed cleft ( Figure 5b by vertical cyan and blue colored lines for the holo-and apo-enzymes, respectively. The NADPH-bound monomer remains in more closed and partially closed conformation, indicated by two main peaks for the distances between domains at 2.68 nm and 2.73 nm, respectively. A third peak is also observed at 2.83 nm, which is close to the one observed in the apoprotein, at 2.86 nm. In contrast, the NADP-bound monomer has the main peak at 2.90 nm for distances between domains, indicating a population of more open conformations, and a second small peak at 2.72 nm for sampling partially-closed structures; the molecular details are discussed later.

Interdomain Cleft
The cofactor binding site in each YqhD monomer spans the interdomain face. The residues of β6/β7 loop region act as the clamp at the mouth of the cofactor binding region, which remains in open conformation in the apo form [7] and closed in the holoenzyme crystal structures. The closed conformation of holoenzyme has hydrogen bond Asp159-Lys359/His363 involving residues D159 of β6/β7 loop region in the Rossmann-type domain and residues K359/H363 of α12 helix in the α-helical domain. For a better understanding of cleft opening/closing and its relationship with cofactor binding/release, we defined interdomain cleft using residue pair D159-H363 at the mouth of the cofactor binding site, on the interface of the Rossmann-type and α-helical domains (see Figure 5a). In the crystal structure, the holoenzyme has a distance of 0.16 nm for cleft, while the apo-enzyme has distances of 1.25 nm for cleft [7], which indicates the involvement of cleft opening/closing in cofactor binding and release. Figure 5a shows the residue pair D159-H363 which is involved in interdomain cleft formation and clamps the Figure 5. (a) Crystallographic structure of the YqhD monomer as apoprotein and holoprotein is in tan and sky-blue colored cartoon representation, respectively, with bound NADP (green licorice representation) and Zn 2+ (violent sphere) in the active site [7]. Cleft residues and Zn 2+ binding residues are labeled. (b) The distribution of distances between Rossmann-type and α-helical domains is shown using the center of mass. (c) The distribution for cleft minimum distances using residue pair D159-H363 is shown for NADP-and NADPH-bound monomers during the simulations. The cleft distances observed in the apo-and holo-YqhD protein crystal structures are marked with blue (apo) and cyan (holo) bars.

Interdomain Cleft
The cofactor binding site in each YqhD monomer spans the interdomain face. The residues of β6/β7 loop region act as the clamp at the mouth of the cofactor binding region, which remains in open conformation in the apo form [7] and closed in the holoenzyme crystal structures. The closed conformation of holoenzyme has hydrogen bond Asp159-Lys359/His363 involving residues D159 of β6/β7 loop region in the Rossmanntype domain and residues K359/H363 of α12 helix in the α-helical domain. For a better understanding of cleft opening/closing and its relationship with cofactor binding/release, we defined interdomain cleft using residue pair D159-H363 at the mouth of the cofactor binding site, on the interface of the Rossmann-type and α-helical domains (see Figure 5a). In the crystal structure, the holoenzyme has a distance of 0.16 nm for cleft, while the apo-enzyme has distances of 1.25 nm for cleft [7], which indicates the involvement of cleft opening/closing in cofactor binding and release. Figure 5a shows the residue pair D159-H363 which is involved in interdomain cleft formation and clamps the NADP/H cofactor in the binding site. The interdomain cleft shows hydrogen bonding between residue pairs D159-K359/H363 and K160-H271/E272 in the crystal structure and simulations. YqhD has~5% D159-H363 hydrogen bond existence in the NADPH-bound monomer trajectory; however, its occurrence is negligible in NADP-bound monomers. Figure 5c shows the distribution of cleft distances using the residue pair D159-H363. Cleft distances observed in the crystal structures are shown in Figure 5c by vertical cyan and blue colored lines for the holo-and apo-enzymes, respectively. NADP-bound monomer shows a wider distribution range of cleft distances, from 0.14-2.0 nm, compared to those in NADPH: 0.14-1.7 nm. It can be seen in Figure 5c Figure 5a), but it is a cooperative motion (involving the cleft forming regions of both domains). Cleft distance data indicate that residue D159 of cleft works as the clamp to monitor cleft opening and closing for cofactor binding and release, and its movement is affected by cofactor oxidation state.

Interdomain Opening-Closing Cleft Dynamics
To characterize the cleft dynamics, we evaluated the minimum distance data of cleft opening/closing residues D159 and H363. For Yqhd monomer, we set a distance of <0.55 nm for closed states and >1.0 nm for open states. By assigning the open, partiallyclosed, and closed states a value of 0, 0.5 and 1, respectively, we were able to generate a graph of state vs time, shown in Figure 6a,b. The state-time analysis is performed on the 2 µs long aggregated trajectory involving 10 sets of NADP/H-bound monomer data. Within the aggregated data, the NADP-bound monomer is in the open state for 41% of the time, partially closed 49%, and closed for 10%. In contrast, the NADPH-bound monomer highly populates the closed state (44% of time points) and exists in the transitional conformation 39% and open state for 17% of the time points in the aggregated trajectories. The state population data indicate an efficient sampling of protein conformations during the 10 independent sets of simulations. From the simulation data, the effect of cofactor oxidation state is clearly evident on the occurrence of one state over another (see Figure 6a,b). For example, the NADP-bound monomer is prone to access open states, while NADPH-bound monomer remains in closed states more frequently. We also evaluated cleft minimum distance data to calculate waiting times for the enzyme to remain in a conformation before switching cleft state (i.e., switching conformation from closed to open, and vice versa). , which is consistent with the experimentally observed properties of the biocatalysts such as the higher affinity of YqhD enzyme for NADPH cofactor indicated by its lower K M value of 0.008 mM [11] over NADP cofactor, K M = 0.150 mM [21].

Cofactor Binding and Release in YqhD Enzyme
The cofactor was found to be quite flexible during simulations, which can be quantified with variations in cofactor RMSD and Rg values.

Cofactor Binding and Release in YqhD Enzyme
The cofactor was found to be quite flexible during simulations, which can be quantified with variations in cofactor RMSD and Rg values. Subsequently, we observed an average minimum distance in the monomer ries between protein and cofactor, and cofactor and Zn 2+ of 0.16 and 0.45 nm, resp to the crystal structure (see Figure 8a,b). The distances between cofactor and pro significant differences during the simulations and are related to the binding of co each monomer. Zn 2+ has an average distance of 0.25 nm from residues Asp194, His267 and His281 during the simulations as observed in the crystal structure (se 5a). A minimum of one water and Ala141 remain within 0.19 nm of the Zn 2+ course of both NADP and NADPH simulations. Within 0.35 nm distances from average of 12 ± 3 and 10 ± 3 water molecules were present in NADP and NADPH YqhD monomer, respectively. A higher number of water molecules were present Zn 2+ in open conformations than the closed one: 15 ± 2 water molecules for NAD 15 ± 7 for NADP-bound YqhD during the time-interval 800-1000 ns, which sh population of open conformations indicated by higher distances between cofac protein in Figure 8a.
NADP-bound YqhD showed cofactor release in three out of five sets of sim trajectories, while NADPH-bound monomer exhibited cofactor release only onc 3. The NADPH-bound monomer has a more rigid cofactor binding, which coincid a conformational change in the cofactor binding pocket. The changes in the cofa tances from the protein are related to the dynamic nature of cofactor in the bind Subsequently, we observed an average minimum distance in the monomer trajectories between protein and cofactor, and cofactor and Zn 2+ of 0.16 and 0.45 nm, respectively, to the crystal structure (see Figure 8a,b). The distances between cofactor and protein had significant differences during the simulations and are related to the binding of cofactor in each monomer. Zn 2+ has an average distance of 0.25 nm from residues Asp194, His198, His267 and His281 during the simulations as observed in the crystal structure (see Figure 5a). A minimum of one water and Ala141 remain within 0.19 nm of the Zn 2+ over the course of both NADP and NADPH simulations. Within 0.35 nm distances from Zn 2+ , an average of 12 ± 3 and 10 ± 3 water molecules were present in NADP and NADPH-bound YqhD monomer, respectively. A higher number of water molecules were present close to Zn 2+ in open conformations than the closed one: 15 ± 2 water molecules for NADPH and 15 ± 7 for NADP-bound YqhD during the time-interval 800-1000 ns, which shows the population of open conformations indicated by higher distances between cofactor and protein in Figure 8a.
lower Rg and indicates the population of a bent conformation of NADP/H cofactor in open cleft states before release shown in Figure 8d. Overall, these results indicate that cofactor release is associated with an increase in the interdomain cleft (>0.55 nm), followed by loosened cofactor binding at cleft distances of 0.55-1.0 nm and loss of hydrogen bonds involving nicotinamide moiety, a change in cofactor conformation to a more compact structure (lower Rg), and, finally, release. Figure 8. The minimum distance is in nm between (a) protein and cofactor, and (b) cofactor and zinc ion as a function of time in the aggregated trajectory of monomer bound to NADP and NADPH. Large distances correspond to cofactor release events. (c) Cofactor binding site shows residues forming hydrogen bonds in the starting structure of YqhD monomer. Cofactor binding residues are colored based on their hydrogen bonding partner, being adenine moiety in orange, nicotinamide in green color, and nucleotide or phosphate group in purple colors. (d) Snapshots of NADP-bound YqhD monomer showing cofactor binding site starting from the equilibrated structure in Set 2 trajectory followed by every 25 ns interval up to cofactor release event in the second set of YqhD simulation.

Discussions and Conclusions
Molecular dynamics simulations were performed on the YqhD dimer, with oxidized and reduced NADP/H cofactors bound in aqueous solution. The starting point of the simulations came from the only solved crystal structure of holo YqhD [7], which has a modified NADPH(OH)2 cofactor, and Zn 2+ in the active site of only one of the monomers within the functional dimer unit. Simulations were run on the functional structure of YqhD, with Zn 2+ present in each monomer, and with NADP + and NADPH cofactors, using careful methodology to prepare these structures (i.e., repairing the crystal structure artifacts). Good conformational sampling of YqhD bound to NADP/H was obtained by NADP-bound YqhD showed cofactor release in three out of five sets of simulation trajectories, while NADPH-bound monomer exhibited cofactor release only once, in Set 3. The NADPH-bound monomer has a more rigid cofactor binding, which coincides with a conformational change in the cofactor binding pocket. The changes in the cofactor distances from the protein are related to the dynamic nature of cofactor in the binding site, as observed earlier with deviations in RMSD and Rg values of monomer and cofactor during the simulations. A total of 15 hydrogen bonds were observed in the starting structure of the cofactor binding site to the cofactor. Figure 8c shows the residues involved in hydrogen bonding with NADP/H cofactor. Hydrogen bonds involve the adenine moiety with residues Thr138, Tyr179, and Thr182; the dinucleotide moiety and phosphate groups with residues Gly38, Ser40, Gly95, Ser96, His281, and Lys160; and the nicotinamide moiety with residues Asp99 and Gly149. During the simulations, Lys160 has no hydrogen bond with NADP/H cofactor, while His281 shows hydrogen bond existence in 5% of NADP-bound monomer trajectory frames. Excluding these latter two hydrogen bonds present in the starting structure, NADP/H-bound structures have conserved hydrogen bonds involving the adenine moiety, nucleotide and phosphate group in 50% of the aggregated trajectory frames. The NADP-bound monomer has only one hydrogen bond involving Leu279 and nicotinamide moiety for 6% existence in the aggregated trajectory. However, NADPH bound monomer forms three potential hydrogen bonds involving Asp99, Ser144, and Gly149 in~30% of the aggregated trajectory. In both NADP/H-bound structures, hydrogen bonds involving the adenine moiety were observed more frequently (50% time points) than with nicotinamide moiety (6% in NADP and 30% in NADPH-bound monomers), which also supports our earlier observations of lower Rg and indicates the population of a bent conformation of NADP/H cofactor in open cleft states before release shown in Figure 8d. Overall, these results indicate that cofactor release is associated with an increase in the interdomain cleft (>0.55 nm), followed by loosened cofactor binding at cleft distances of 0.55-1.0 nm and loss of hydrogen bonds involving nicotinamide moiety, a change in cofactor conformation to a more compact structure (lower Rg), and, finally, release.

Discussion and Conclusions
Molecular dynamics simulations were performed on the YqhD dimer, with oxidized and reduced NADP/H cofactors bound in aqueous solution. The starting point of the simulations came from the only solved crystal structure of holo YqhD [7], which has a modified NADPH(OH) 2 cofactor, and Zn 2+ in the active site of only one of the monomers within the functional dimer unit. Simulations were run on the functional structure of YqhD, with Zn 2+ present in each monomer, and with NADP + and NADPH cofactors, using careful methodology to prepare these structures (i.e., repairing the crystal structure artifacts). Good conformational sampling of YqhD bound to NADP/H was obtained by performing five sets of simulations that were assigned different initial velocities from a Maxwell-Boltzmann velocity distribution at 300 K. The structures remain conserved throughout the simulations, regardless of the oxidation state of NADP/H cofactor. So, we compiled a combined trajectory of the 200 ns from the set of five simulations with both NADP/H cofactors and obtained sampling of various conformations representing opening and closing of the cofactor-binding interdomain cleft. Each monomer of the YqhD dimer showed coordinated cleft opening/closing, with both oxidized and reduced NADP/H, via the movement of β6/β7 and α12 regions in the interdomain cleft. The cleft remains open in the apo form of the YqhD crystal structure [7], indicating that it opens and closes for cofactor binding. Observations of cofactor release during MD simulations confirm the movement to an open-cleft conformation prior to the release of both NADP and NADPH.
The sampling of open-cleft conformations depends on the cofactor oxidation state. NADPbound monomers tend to sample more open conformations, with cleft distances ranging from~0.7-1.8 nm (cf.~0.7-1.5 for NADPH, 1.25 nm for apoprotein crystal structure [7]).
The dynamics of cleft opening/closing were also found to depend on the oxidation/reduction of the cofactor. With both cofactors, YqhD undergoes periods of rapid cleft opening and closing, sampling open conformations similar to the apoenzyme [7]. NADPbound monomer scarcely populates conformations with a closed cleft in 10% of time points and remains in a partially closed state for 49% of time points or adopts open cleft conformation in 41% of time points. NADPH-bound monomers remain in closed conformation for 44% and spend only 17% in open state conformations. The maximum waiting time of transition from closed to open state was observed to be~8 ns in NADP-bound monomer and 68 ns in NADPH-bound monomer. However, NADP-bound monomer underwent a total of 750 transitions between open and closed states, compared to 318 total cleft state transitions observed for NADPH-bound monomer, indicating a higher propensity for cleft opening dynamics when the enzyme is bound to the oxidized cofactor. Out of five independent simulations, cofactor release was observed in three sets of NADP-bound and only one set of NADPH-bound structures. This suggests an entropic, dynamics-based preference for cofactor release when NADP is bound.
Significant differences were observed even in the conformation of NADP/H cofactor, concomitant with conformational differences in the monomers. The differences in cofactor-enzyme interactions and YqhD conformational dynamics may rationalize differences in the Michaelis constant, K M , which depends on cofactor binding (k f ) and release (k r ) substrate-binding rates as K M = (k r + k cat )/k f . The higher number of cofactor-monomer hydrogen bonds between NADPH-bound vs. NADP-bound YqhD may lead to higher binding rates and/or slower release rates for NADPH, resulting in the higher affinity for NADPH indicated by its lower K M value (0.008 mM [11] vs. 0.150 mM [21] for NADP). Furthermore, the less frequent sampling of open conformations with NADPH-bound YqhD may hinder NADPH release, dropping k r and subsequently K M for NADPH. Meanwhile, conformational dynamics appear to promote the more facile release of NADP, following hydride transfer from NADPH to aldehydes. Thus, the conformational differences induced by cofactor oxidation state, dynamical effects of more-frequent cofactor cleft opening with NADP, and differences in hydrogen bond motifs may lead to preferential kinetics for the aldehyde reductase activity of YqhD. These findings raise questions about whether enzymes with higher alcohol dehydrogenase activity show a dynamic preference for the release of reduced cofactor. Results of this study enhance our basic understanding toward this class of enzyme, with the possible application of guiding the rational design of YqhD to enhance substrate affinity and biocatalyst efficiency.

Starting Coordinates
The crystal structure of YqhD (PDB ID: 1OJ7, 2.0 Å resolution) [7] was used to extract the starting coordinates for setting up MD simulations. The starting coordinates include a dimer (crystallographic monomers A and D, shown in Figure 1a) with bound Zn 2+ , modified NADP cofactor as NADPH(OH) 2 due to oxidative stress, and crystallographic waters. In the crystal structure, only one monomer (D) of the YqhD homodimer has Zn 2+ in the active site. Herein, monomer D is used to model the starting structure of YqhD homodimer with bound Zn 2+ and NADPH(OH) 2 cofactor using PyMOL molecular graphics system [33]. Monomer D was duplicated and superimposed onto monomer A using the root mean square deviation (RMSD) minimization criteria of PyMOL to generate a homodimer structure with two bound Zn 2+ ions for MD simulations.

Modeling of NADP/H Cofactor
In the crystal structure of YqhD homodimer, [7] NADP cofactor is present as NADPH(OH) 2 (see Figure 1b), with modification at the C5 and C6 positions of the nicotinamide moiety. The native oxidized cofactor (NADP) was modeled by removing the hydroxyl moieties at the fifth and sixth positions in the nicotinamide ring. The reduced cofactor (NADPH) was prepared by modifying the oxidation state of the C4 atom of the nicotinamide ring. Forcefield parameters come from CHARMM 36 force field [34] for Zn 2+ , NADP, and NADPH cofactors. Throughout this paper, NADP/H is used to signify the NADPH cofactor generally, whether in its oxidized (NADP + ) or reduced (NADPH) state. When a specific oxidation state is indicated, NADP and NADPH are specified.

Molecular Dynamics Simulations
The CHARMM 36 force field [34] was used for the simulations summarized in Table 1 using the GROMACS software package version 5.1.3 [35]. The YqhD dimer bound to Zn 2+ and NADP/H cofactors was centered in a cubic periodic box (~11 nm 3 ) and set to have a distance larger than 1 nm from any side of the box. Solvent molecules having any atom within 0.15 nm from the protein were removed. TIP3P model [36,37] was used for explicit water solvent. Sodium counter ions were added by replacing the solvent molecules at the sites of most negative electrostatic potential to provide the box with a total charge of zero. The protonation state of residues was assumed to be the same as that of the isolated amino acids in solution at pH 7. The LINCS [38] algorithm was used to constrain all bond lengths and the SETTLE [39] algorithm was used for the water molecules. Electrostatic interactions were calculated using the Particle Mesh Ewald method [40]. For the calculation of long-range interactions, a grid spacing of 0.12 nm combined with a fourth-order B-spline interpolation was used to compute the potential and forces between grid points. A nonbonded pair-list cutoff of 1.4 nm was used and updated at every five time-steps. V-rescale thermostat [41] was used to keep the temperature at 300 K through a weak coupling of the system to an external thermal bath with a relaxation time constant τ = 0.1 ps. The pressure of the system was kept at 1 bar using Berendsen's barostat [42] with a time constant of 1 ps. A time step of 2 fs was used to integrate the equations of motion. First, the simulated systems were energy minimized, using the steepest descent algorithm, for at least 5000 steps to remove clashes between atoms that were too close. After energy minimization, all atoms were given an initial velocity obtained from a Maxwell-Boltzmann velocity distribution at 300 K to start the MD simulations. The system was initially equilibrated by 30 ps with position restraints on the heavy atoms of the dimer to allow relaxation of the solvent molecules. After the equilibration procedure, position restraints were removed, and the system was gradually heated from 50 K to 300 K during 200 ps of simulation. The equilibrated structure was used to perform MD simulations in the NPT ensemble for 20 ns at 300 K. The final conformation of the 20 ns isothermalisobaric NPT ensemble was used to set up five sets of 200 ns independent simulations in the canonical NVT ensemble initiating with five different velocities.
The starting coordinates were adopted from the crystal structure which has modified cofactor (NADPH(OH) 2 ); hence, the equilibrated conformations of YqhD bound to NADP and NADPH cofactors were obtained after a 20 ns NPT simulation, which was used as the starting point for the final production run of 200 ns long NVT simulations. The starting crystallographic coordinates of the YqhD homodimer were used as the reference structure for the analysis of the trajectories to characterize the conformational changes induced by a change in the oxidation state of NADP/H cofactor. The analysis was focused on three sets/sub-structures within the simulation data, considering dimer, monomer, and domain in a 1 µs aggregated trajectory for each cofactor oxidation state (NADP and NADPH). Structural analysis includes the root mean square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyration (Rg), and cofactor binding interactions (Zn +2 and NADP/H), comparing these with respect to the crystal structure.

Cluster Analysis
The conformational diversity of the structures generated during the MD simulations was characterized using the Gromos [43] clustering algorithm. In this method, an RMSD cutoff criterion is used to assign a structure in a cluster based on the root-mean-square differences of selected atoms among the conformations obtained from the simulations. An RMSD cutoff of 0.13 nm was used to determine neighboring backbone atom conformations of YqhD monomers to discriminate among the less varying conformations. For the cluster analysis, a total of 52,815 structures were evaluated for each monomer within the dimer over the 2 µs aggregated trajectory, using time intervals of 100 ps. [33]