|Home | About | Journals | Submit | Contact Us | Français|
Inosine monophosphate dehydrogenase (IMPDH) catalyzes an essential step in the biosynthesis of guanine nucleotides. This reaction involves two different chemical transformations, an NAD-linked redox reaction and a hydrolase reaction, that utilize mutually exclusive protein conformations with distinct catalytic residues. How did Nature construct such a complicated catalyst? Here we employ a “Wang-Landau” metadynamics algorithm in hybrid quantum mechanical/molecular mechanical (QM/MM) simulations to investigate the mechanism of the hydrolase reaction. These simulations show that the lowest energy pathway utilizes Arg418 as the base that activates water, in remarkable agreement with previous experiments. Surprisingly, the simulations also reveal a second pathway for water activation involving a proton relay from Thr321 to Glu431. The energy barrier for the Thr321 pathway is similar to the barrier observed experimentally when Arg418 is removed by mutation. The Thr321 pathway dominates at low pH when Arg418 is protonated, which predicts that the substitution of Glu431 with Gln will shift the pH-rate profile to the right. This prediction is confirmed in subsequent experiments. Phylogenetic analysis suggests that the Thr321 pathway was present in the ancestral enzyme, but was lost when the eukaryotic lineage diverged. We propose that the primordial IMPDH utilized the Thr321 pathway exclusively, and that this mechanism became obsolete when the more sophisticated catalytic machinery of the Arg418 pathway was installed. Thus, our simulations provide an unanticipated window into the evolution of a complex enzyme.
Many enzymes have the remarkable ability to catalyze several different chemical transformations. For example, IMP dehydrogenase catalyzes both an NAD-linked redox reaction and a hydrolase reaction. These reactions utilize distinct catalytic residues and protein conformations. How did Nature construct such a complicated catalyst? While using computational methods to investigate the mechanism of the hydrolase reaction, we have discovered that IMP dehydrogenase contains two sets of catalytic residues to activate water. Importantly, the simulations are in good agreement with previous experimental observations and are further validated by subsequent experiments. Phylogenetic analysis suggests that the simpler, less efficient catalytic machinery was present in the ancestral enzyme, but was lost when the eukaryotic lineage diverged. We propose that the primordial IMP dehydrogenase utilized the less efficient machinery exclusively, and that this mechanism became obsolete when the more sophisticated catalytic machinery evolved. The presence of the less efficient machinery could facilitate adaptation, making the evolutionary challenge of the IMPDH reaction much less formidable. Thus our simulations provide an unanticipated window into the evolution of a complex enzyme.
Textbooks extol the extraordinary catalytic power and specificity of enzymes, yet the ability of many enzymes to promote several different chemical transformations is even more remarkable. In examples such as the polyketide synthases, the substrate is tethered to a flexible linker and swings gymnastically between separate active sites . The evolutionary path to the assembly of such enzymes seems reasonably straightforward: gene duplication and recombination, followed by optimization of a promiscuous activity [2–6]. In contrast, enzymes such as IMP dehydrogenase (IMPDH) move around a stationary substrate, restructuring the active site to accommodate different transition states . Such enzymes pose an evolutionary conundrum: it seems unlikely that Nature could simultaneously install multiple sets of catalytic machinery into the ancestral protein. IMPDH controls the entry of purines into the guanine nucleotide pool, which suggests that the origins of IMPDH are primordial, so the ancestral IMPDH probably utilized a simpler catalytic strategy.
IMPDH catalyzes two very different chemical transformations: (1) a dehydrogenase reaction between IMP and NAD+ that produces a Cys319-linked intermediate E-XMP* and NADH, and (2) a hydrolysis reaction that releases XMP (Figure 1A) [7,8]. A mobile flap is open during the hydride transfer reaction, permitting the association of NAD+. After NADH departs, this flap occupies the dinucleotide site, carrying Arg418 and Tyr419 into the active site and converting the enzyme into a hydrolase (Figure 1B). Thus, the dehydrogenase and hydrolase reactions utilize mutually exclusive conformations of the active site.
All enzymes that catalyze hydrolysis reactions have some strategy to activate water. This strategy has been difficult to recognize in IMPDH because the hydrolytic water interacts with three residues that are usually protonated at physiological pH: Thr321, Arg418, and Tyr419 (Figure 1C) . The rate of the hydrolysis step decreases by a factor of 103 when Arg418 is substituted with Ala or Gln, whereas a decrease of approximately 20 is observed when Tyr419 is substituted with Phe [10,11]. Neither Arg418 nor Tyr419 is involved in the dehydrogenase reaction, as expected, given their position on the mobile flap. In contrast, Thr321 is found on the same loop as the catalytic Cys319, and both the dehydrogenase and hydrolysis reactions are decreased by a factor of 20 when this residue is substituted . These observations suggest that Arg418 is the most likely candidate for the role of general base in the IMPDH reaction [11,12].
We performed a series of hybrid quantum mechanical/molecular mechanical (QM/MM) simulations to further investigate the mechanism of the hydrolysis reaction of IMPDH. Surprisingly, these simulations find that IMPDH possesses two mechanisms to activate water: the Arg418 pathway as previously proposed, and a second pathway utilizing Thr321. Phylogenetic analysis indicates that the Thr321 pathway was present in the ancestral enzyme. These observations suggest that the primordial IMPDH used the Thr321 pathway exclusively, and elimination of the Arg418 pathway by mutation of modern IMPDH creates an enzymatic atavist.
We have applied computational methods to further investigate the mechanism of water activation in IMPDH, employing a “Wang-Landau” metadynamics algorithm  in hybrid quantum mechanical/molecular mechanical (QM/MM) simulations [14–21]. The Wang-Landau recursion procedure adaptively updates the height of the basis Gaussian, which allows the metadynamics algorithm to be realized in a more robust and efficient fashion. The simulation models were derived from the crystal structure of Tritrichomonas foetus IMPDH in complex with mizoribine monophosphate (MZP) (Protein Data Bank [PDB] accession code 1PVN), which describes the closed conformation of the hydrolysis reaction . E-XMP* was modeled based on the guanine residue topology and parameters in CHARMM 22 . The atoms in the reaction centers (colored in red in Figures 2–4) were treated quantum mechanically using the self-consistent charge density-functional tight-binding (SCCDFTB) potential , and the rest of the system (colored in blue in Figures 2–4) was treated classically using the CHARMM 22 force field . The molecular dynamics simulations were performed with generalized solvent boundary conditions (GSBP) [24,25]. Further details on the simulations are provided in Figure S1.
When Arg418 is deprotonated in the starting condition, the lowest energy pathway for the hydrolysis reaction involves the transfer of a proton to the neutral Arg418 (the Arg418 pathway, Figure 2A–2C). Proton transfer is virtually complete at the transition state, and the developing hydroxide is stabilized by interactions with Tyr419, Thr321, and another water molecule. Importantly, a stable hydroxide intermediate is not observed; the developing hydroxide instantaneously reacts with E-XMP*. These results are in remarkable agreement with experimental observations: solvent isotope effects (SIE) demonstrate that proton transfer is rate limiting (SIE = ~2 ), and Bronsted analysis indicates that proton transfer is virtually complete in the transition state (β = ~1 ). However, the calculated energy barrier is only 8.0 kcal/mol, much less than the experimentally observed barrier of 16 kcal/mol . This difference may reflect uncertainties in the calculation, but we believe this is unlikely. A more intriguing source of discrepancy arises from the starting condition of neutral Arg418; if only a small fraction of the enzyme exists in this state, the energy barrier will be correspondingly increased. Indeed, if the pKa of Arg418 is 12.5, as for a typical Arg residue, the barrier would be increased by approximately 6 kcal/mol.
The pKa of a Tyr residue is usually two units lower than an Arg, which suggests that a deprotonated Tyr419 might activate water while Arg418 remains protonated. Further simulations argue against such a mechanism; instead, the deprotonated, negatively charged Tyr419 interacts strongly with positively charged Arg418 and cannot interact with water. Therefore, Tyr419 is unlikely to play the role of general base in the wild-type enzyme. However, the situation changes when Arg418 is substituted with Gln: now the Tyr419 phenolate can accept a proton from water. The barrier is approximately 17 kcal/mol (Figure 3). Assuming a pKa of 10, as is usual for a Tyr residue, then deprotonation of Tyr419 will further increase the barrier to 21–22 kcal/mol, which is very similar to the barrier observed in the reactions of the Arg418Gln and Arg418Ala variants (~20 kcal/mol [10,11]). As above, the simulations suggest that proton transfer is rate limiting and essentially complete in the transition state. Whereas the landscape contains a shallow valley suggesting the presence of a hydroxide intermediate, the barriers are less than a kT, so the intermediate would not have a finite lifetime. This simulation is generally consistent with experiments, where SIEs of 3–5 are observed when Arg418 is substituted . However, the magnitude of these SIEs is greater than expected if the transition state is indeed late as suggested by the simulations. Interestingly, no activity is observed in the Arg418Ala/Tyr419Phe double mutant, though this fact may equally well be attributed to the inability to form the closed conformation required for the hydrolysis reaction as to the loss of the general base catalyst . Together, the simulations and experiments suggest Tyr419 may act as a surrogate general base in the absence of Arg418.
Similar surrogate residues have been invoked to explain residual activity in other enzyme systems . In RNase T1, His40 residue assumes the role the general base when Glu58 is substituted with Ala . Similarly, in ketosteroid isomerase, Asp99 may catalyze proton transfers in the Asp38Ala variant . Water or buffer molecules can also replace the function of missing catalytic residues [29,30]. These examples illustrate the resilience and plasticity of enzyme catalysis.
Surprisingly, the simulations suggest a second pathway for water activation when the starting condition is protonated Arg418: Thr321 abstracts a proton from water while simultaneously transferring its own proton to Glu431 (Figure 4). As in the Arg418 pathway, the developing hydroxide is stabilized by Tyr419 and another water molecule, and the hydroxide attack occurs instantaneously; the protonated Arg418 also stabilizes the developing hydroxide by 1–2 kcal/mol. The calculated free energy barrier for the Thr321 pathway is approximately 20 kcal/mol. The proton transfers are simultaneous and rate limiting. When a simulation was performed with Glu431 treated as a molecular mechanical (MM) residue, which eliminates the possibility of proton transfer while maintaining electrostatic interactions, the energy barrier increases to at least 35 kcal/mol (Figure S2). Likewise, when Glu431 is substituted with Gln, the barrier increases to at least 38 kcal/mol. Therefore, the presence of Glu431 is essential for the operation of the Thr321 pathway.
The simulations suggest that the Thr321 pathway is favored at low pH, whereas the Arg418 pathway becomes dominant at high pH, which predicts that the pH-rate profile will shift to the right when the Thr321 pathway is disrupted by the Glu431Gln mutation. This prediction was confirmed experimentally (Figure 5): the Glu431Gln mutation shifts the pKa from 7.2 ± 0.1 to 7.6 ± 0.1, but has only a small effect on the pH-independent value of kcat (kcat = 2.2 and 1.4 s−1 for wild type and Glu431Gln, respectively; these values are in good agreement with previous reports [31,32]). Assuming that the pKa shift is entirely attributable to the loss of the Thr321 pathway, the barrier for the Thr321 pathway is approximately 19 kcal/mol, as predicted by the simulations (see Figure S3).
When Arg418 is substituted with Gln, the barrier for the Thr321 pathway is approximately 21 kcal/mol, which is similar to the barrier observed experimentally in the Arg418Ala and Arg418Gln variants [10,11]. Therefore, both the Thr321 pathway and the Tyr419 pathway can account for the residual activity of the Arg418Ala and Arg418Gln variants. However, since the Thr321 pathway involves the simultaneous transfer of two protons, this pathway can account for the large solvent isotope effects observed in the Arg418 variants (SIE = 3–5 [10,11]). Therefore, we constructed the Arg418Gln/Glu431Gln variant, which should disrupt the Thr321 pathway but leave the Tyr419 pathway intact. The simulations predict that the activity of this variant should be approximately the same as the Arg418Gln, but that the solvent isotope effect should be reduced. These predictions were confirmed in subsequent experiments: (1) the value of kcat for Arg418Gln/Glu431Gln is decreased by 50% relative to that of Arg418Gln, as expected if the Thr321 pathway was lost (0.0020 ± 0.0002 s−1 and 0.0040 ± 0.0004 s−1, respectively); and (2) though the errors on the SIE are larger than ideal, nonetheless, a smaller SIE is observed in the reaction of Arg418Gln/Glu431Gln, consistent with the loss of the Thr321 pathway (SIE = 2.1 ± 0.3 and 2.3 ± 0.4 for Arg418Gln/Glu431Gln in two independent determinations versus 2.9 ± 0.5 for Arg418Gln and 5 ± 2 for Arg418Ala [10,11]). These experiments confirm the operation of the Thr321 pathway in IMPDH.
To the best of our knowledge, the presence of dual mechanisms for water activation in an enzyme active site is unprecedented. Why would an enzyme have two pathways to accomplish the same task? We believe the Thr321 pathway may be vestige of evolution, and phylogenetic analysis is consistent with this hypothesis (Figure 6; see Figure S4 for the complete phylogenetic tree). The closest relative of IMPDH is GMP reductase (GMPR), which catalyzes the conversion of GMP to IMP and ammonia with concomitant oxidation of NADPH (Figure 6) . Cys319, Thr321, and Glu341 are also conserved in GMPR, which suggests that these residues were present in the IMPDH/GMPR ancestor. X-ray crystal structures show that the conserved Cys, Thr, and Glu display similar interactions in both GMPR and IMPDH (Figure 6), suggesting that these residues may have similar functions in both enzymes.
To confirm that GMPR activity depends on the presence of Cys186, Thr188, and Glu289, we tested the effect of mutations of these residues on the activity of Escherichia coli GMPR in a complementation assay (Figure 7). E. coli H1173 requires both adenosine and guanosine for growth due to mutations in both purA and guaC  (Figure 7). Growth on guanosine alone is restored with plasmid pGS682, which carries the wild-type E. coli guaC gene . However, mutations in Cys186, Thr188, and Glu289 clearly compromise the ability of pGS682 to restore GMPR activity, demonstrating that selective pressure exists to conserve these residues.
Although the mechanism of the GMPR reaction has not been characterized, some clear parallels can be drawn with the IMPDH reaction, and E-XMP* may well be an intermediate. Importantly, if E-XMP* forms as proposed, the active site must be constructed to prevent the hydrolysis reaction. Kinetic and structural experiments clearly indicate that the reaction only proceeds when NADPH is bound in the active site and can block the access of water [33,36,37]. Moreover, GMPR does not contain the Arg418-Tyr419 dyad, and the flap is truncated relative to the corresponding region of IMPDH, as expected, given that the hydrolysis of E-XMP* must be avoided during the GMPR reaction. Therefore, the Arg418-Tyr419 dyad could have been installed as IMPDH optimized. Alternatively, the dyad may have been present in the ancestral IMPDH/GMPR, but was subsequently remodeled in the GMPR lineage; since the flap binds in the same site as NAD+, this scenario suggests that the ancestral IMPDH/GMPR was a hydrolase. While we cannot rule out the latter scenario, we note that IMPDH is a member of the FMN oxidoreductase superfamily of (β/α)8 barrel proteins (unfortunately, none of these proteins is sufficiently similar to permit rooting of the tree) [38–40]. Therefore, it seems more likely that the ancestral enzyme was a promiscuous dehydrogenase, and the flap carrying the hydrolase activity was the later addition.
In contrast, the Thr321 pathway was likely present in the ancestral IMPDH/GMPR. All IMPDHs and GMPRs contain Thr321 (Figures 6 and S4, and Text S1). As noted above, Thr321 also plays a role in the dehydrogenase reaction of IMPDH , which suggests that Thr321, like Cys319, was inherited from the ancestral redox enzyme. Glu431 is conserved among GMPRs, suggesting that the Thr321 pathway has a crucial function in this reaction, perhaps operating in the reverse to protonate the ammonia leaving group. Curiously, although Glu431 is highly conserved among IMPDHs, it is substituted with Gln in the eukaryotic branch as well as in a few prokaryotic IMPDHs. We suggest that the ancestral IMPDH/GMPR utilized the Thr321 pathway exclusively, but this pathway became expendable once the Arg418 pathway was established. Phylogenetic analysis is consistent with this view: maximum likelihood analysis indicates that the ancestral enzyme almost certainly contained Glu at position 431 (probability = 0.87) .
Why then is Glu431 conserved in the majority of prokaryotic IMPDHs? The presence of the Thr pathway increases turnover, which may be important in maintaining the high concentration of guanine nucleotides required to support the rapid proliferation of prokaryotes. More intriguingly, Glu431 provides 5–10-fold resistance to mycophenolic acid, a natural product that specifically inhibits IMPDH . Approximately 5% of microorganisms contain some mechanism to modify mycophenolic acid, which suggests that this compound is reasonably prevalent in the environment . Indeed, the extraordinary divergence of the adenosine subsite of IMPDH may be a response to the assault of natural product inhibitors such as mycophenolic acid and mizoribine . This divergence occurs despite the multiple functional constraints imposed by interactions with both NAD+/NADH and the flap. The presence of the Thr pathway could facilitate this adaptation, making the evolutionary challenge of the IMPDH reaction much less formidable.
Plasmid pGS682, a pUC plasmid carrying the 1.4-kb guaC insert from pGS89 , was a generous gift from Simon Andrews (University of Sheffield). E. coli strain H1173 was obtained from the E. coli Genetic Stock Center (Yale University).
Atoms within a radius of 22 Å around the reaction center were treated as the dynamic region; this region was propagated with regular Newtonian dynamics by applying leapfrog integrator and 1-fs time step. The atoms in the layer between the radii of 22 Å and 25 Å were treated as the buffer region; the heavy atoms in this region were harmonically restrained with the force constants scaled linearly with the distance from the sphere center. The force constants around the boundary of the 25 Å sphere were set as implied by the B factors of the crystal structure. In the buffer region, Langevin dynamics were applied with the friction coefficients also linearly scaled with the distance from the sphere center. The friction coefficients around the boundary 25 Å sphere were set as 60. CHARMM 22 force fields  were utilized as the molecular mechanical potentials in these simulations (colored in blue in Figures 2–4) and SCCDFTB (self-consistent charge density-functional tight-binding) method was applied as the quantum mechanical potential on the atoms involved in the chemical reactions (colored in red in Figures 2–4). For the nonbonded interactions, an extended electrostatic treatment was applied with the electrostatic interactions within 12 Å described by group-based coulombic interactions.
IMP, acetylpyridine adenine dinucleotide (APAD+), Tris, and MES were purchased from Sigma. DTT was purchased from Research Organics. Wild-type and Glu431Gln IMPDH from T. foetus were expressed in E. coli and purified as described previously [10,32]. All assays were performed at 25 °C. The release of NADH is partially rate limiting [11,31]. Therefore, to ensure that hydrolysis is completely rate limiting, these experiments used APAD+ . Pre-steady-state experiments were performed to demonstrate that hydride transfer and APADH are rapid over the entire pH range ( and unpublished data). Standard IMPDH assays contained saturating concentrations of IMP (2 mM) and varying concentrations of APAD+ in 100 mM KCl, 1 mM DTT, and 50 mM of the appropriate buffer (MES for pH 5.0–7.0, and Tris-HCl for pH 7.3–9.3). Activity was measured by monitoring the absorbance of APADH at 363 nm on a Hitachi U-2000 UV-visible spectrophotometer. Steady-state parameters with respect to APAD+ were derived at saturating IMP concentrations by plotting the initial velocity against APAD+ concentration and fitting to an equation describing uncompetitive substrate inhibition using SigmaPlot (SPSS):
where (kcat)app are the apparent values at each pH, (kcat)indep are the pH-independent values, and Ka is the acid dissociation constant for the most acidic ionization.
IMPDH/GMPR amino acid sequences (IMPDH IPR005990, GMPR1 IPR005993, and GMPR2 IPR005994) were retrieved from the InterPro database (http://www.ebi.ac.uk/interpro/). Additionally, BLAST  searches with the T. foetus IMPDH (P50097) and human GMPR1 (P36959) amino acid sequences were performed. Sequences from the BLAST search that were already part of the InterPro dataset were removed, and an initial multiple sequence alignment was performed with MUSCLE . A neighbor joining tree (unpublished data) was constructed in PAUP* 4.0b10 , and 95 sequences were selected for a Bayesian phylogenetic analysis. The sequences of this subset were realigned with Espresso [47,48]. A Bayesian phylogenetic analysis was performed with the parallel version of MrBayes 3.1.2 [49,50]. Amino acid substitution rates and state frequencies were fixed to the WAG parameters . A uniform (0.0, 200.0) prior was assumed for the shape parameter of the gamma distribution of substitution rates , an unconstrained exponential prior with rate 10.0 for branch lengths, and all labeled topologies were a priori equally probable. Two independent MCMC analyses were run, each with one cold chain and three heated chains, with the incremental heating schema implemented in MrBayes (λ=0.2). Convergence was assumed after the topology samples from the two cold chains had reached an average standard deviation of split frequencies of less than 0.01 (after 1,610,000 generations). Accession numbers, detailed results, and the full tree are found in Text S1.
E. coli strain H1173 (F-, guaC23, tonA2, proA35, lacY1, tsx-70, supE44?, gal-6, l-, trp-45, tyrA2, rpsL125, malA1 (lR), xyl-7, mtl-2, thi-1, purH57) contains mutations in purH and guaC, and therefore requires both adenosine and guanosine for growth. Bacteria were transformed with pGS682 carrying either the wild-type guaC gene or variants containing C186A, T188A, and E289Q mutations. Cultures were grown overnight in LB or LB/ampicillin and 5 μl of 1/20 serial dilutions were plated on M9 minimal media containing 0.5% casamino acids, 100 μg/ml l-tryptophan, 0.1% thiamin, 50 μg/ml guanosine, and/or 50 μg/ml adenosine.
In order to optimize the Wang-Landau metadynamics conditions, three setups with the final Gaussian heights 1.0 kcal/mol, 0.06 kcal/mol, and 0.01 kcal/mol were executed. The 1.0 kcal/mol simulation yielded a result with large uncertainties and gave the free energy barrier of 14 kcal/mol at the end of a 1-ns simulation. The 0.01 kcal/mol simulation yielded a nicely converged free energy diagram with a barrier of 8 kcal/mol, but required more than 20 ns. The 0.06 kcal/mol also yielded a free energy barrier of 8 kcal/mol with acceptable fluctuations, but required only 5 ns. Based on these benchmark results, 0.06 kcal/mol was utilized as the final Gaussian height throughout all the simulations.
(3.5 MB TIF)
Methods as described in Figure S1
(1.75 MB TIF)
Assuming that E431Q mutation disables the Thr pathway, but has no effect on the pH dependence of the Arg pathway, then the pH-rate profile of the wild-type enzyme is described by:
where Ka = 10−7.6 as determined from the pH dependence of E431Q. The pH-rate profile of the wild-type enzyme could be reasonably described with the above equation when the value of kThr is 0.15 s−1, which corresponds to an energy barrier of approximately 19 kcal/mol. For easier visualization, the pH-rate profiles of Figure 3A were normalized so that the pH-independent values of kcat = 1 for both wild type and E431Q, where x becomes the normalized value of kThr.
(2.05 MB TIF)
The unrooted tree was inferred with MrBayes (including posterior probabilities) . Organism names are followed by their sequence accession codes. IMPDH, GMPR1, and GMPR2 refer to members of InterPro accession codes IPR005990, IPR005993, and IPR005994, respectively.
(4.39 MB TIF)
We thank Karen Allen, Margaret Phillips, and Gregory Petsko for comments on the manuscript. Jeffrey Boucher constructed the GMPR mutants.
Funding. This work was supported by the National Institutes of Health GM54403 (LH).
Competing interests. The authors have declared that no competing interests exist.