Cooperative Electrostatic Interactions Drive Functional Evolution in the Alkaline Phosphatase Superfamily.

It is becoming widely accepted that catalytic promiscuity, i.e., the ability of a single enzyme to catalyze the turnover of multiple, chemically distinct substrates, plays a key role in the evolution of new enzyme functions. In this context, the members of the alkaline phosphatase superfamily have been extensively studied as model systems in order to understand the phenomenon of enzyme multifunctionality. In the present work, we model the selectivity of two multiply promiscuous members of this superfamily, namely the phosphonate monoester hydrolases from Burkholderia caryophylli and Rhizobium leguminosarum. We have performed extensive simulations of the enzymatic reaction of both wild-type enzymes and several experimentally characterized mutants. Our computational models are in agreement with key experimental observables, such as the observed activities of the wild-type enzymes, qualitative interpretations of experimental pH-rate profiles, and activity trends among several active site mutants. In all cases the substrates of interest bind to the enzyme in similar conformations, with largely unperturbed transition states from their corresponding analogues in aqueous solution. Examination of transition-state geometries and the contribution of individual residues to the calculated activation barriers suggest that the broad promiscuity of these enzymes arises from cooperative electrostatic interactions in the active site, allowing each enzyme to adapt to the electrostatic needs of different substrates. By comparing the structural and electrostatic features of several alkaline phosphatases, we suggest that this phenomenon is a generalized feature driving selectivity and promiscuity within this superfamily and can be in turn used for artificial enzyme design.

Introduction

In recent years,[1] it has become widely accepted that catalytic promiscuity, i.e., the ability of many enzymes to catalyze the turnover of multiple chemically distinct substrates, plays a key role in the evolution of new functions, allowing for rapid responses to environmental changes.[2,3] Furthermore, interest in this phenomenon has exploded as it has been increasingly shown to be a powerful tool for gaining knowledge not just into the process of natural functional evolution,[2] but also as a factor that can be exploited in effective artificial enzyme design.[1,3,4] Such promiscuity appears to be highly pronounced among many phosphotransferases, such as the recently evolved bacterial phosphotriesterase (PTE),[5] serum paraoxonase 1 (PON1),[6] and members of the alkaline phosphatase (AP) superfamily,[7−9] to name a few examples. This latter superfamily has additionally played a central role as a model system for understanding enzyme catalytic promiscuity,[7,9−13] i.e., the ability of a given enzyme to catalyze more than one distinct chemical reaction.

The characterized members of the AP superfamily are highly promiscuous hydrolytic enzymes capable of interchangeable cleavage of P–O, S–O, and P–C bonds.[14,15] That is, they have been shown to catalyze the hydrolysis of a range of substrates that differ in the nature of their TS solvation and protonation patterns, and thus in their requirements for efficient catalysis (see discussion in refs (16−18)). Furthermore, all known AP superfamily members are metallohydrolases that employ similar catalytic scaffolds, which are comprised of at least one divalent metal ion in their respective active sites (Figure 1). This metal ion plays an important role in activating the nucleophile, which is generally thought to be an alcohol or alkoxide depending on the particular superfamily member,[15] by increasing the concentration of its active deprotonated form. Additionally, while there are a number of similarities between different known members of the superfamily, there are also broad differences in their metal requirements, overall structures, and specific choice of nucleophile, which can in turn be linked to changes in specificity patterns.[11] Despite these differences, a particular hallmark of this superfamily is crosswise-promiscuity, in that the native substrate for one member of the superfamily is often a promiscuous substrate for another,[11,15] in some cases with high (and almost comparable) proficiencies toward both the native and promiscuous substrates.[8,9,12] As a result, these enzymes provide a perfect showcase to generate a systematic roadmap of the process of functional evolution within an enzyme superfamily, as well as a broader model system for understanding the evolution of phosphohydrolase activity.[15]

Figure 1

An active site comparison of selected members of the AP superfamily. Shown here are the active sites of (A) Rhizobium leguminosarum phosphonate monoesterase (PDB code 2VQR),[8] (B) Escherichia coli alkaline phosphatase (1ALK),[20] (C) Pseudomonas aeruginosa arylsulfatase (1HDH)[19] and (D) Xanthomonas axonopodis nucleotide pyrophosphatase/phosphodiesterase (2GSN).[21] The figure highlights the presence of divalent metal ions as well as the conservation of some of the residues surrounding them.

Among the different superfamily members, the name giving member AP[20,22−25] as well as the very closely related nucleotide pyrophosphatase/phosphodiesterase (NPP)[21,25−27] have been the subject of extensive scrutiny. A lesser-studied subset of enzymes that stand out in this superfamily are those classified as phosphonate monoester hydrolases (PMHs), such as the enzymes from Rhizobium leguminosarum (RlPMH)[8] and Burkholderia caryophylli (BcPMH).[12] These highly promiscuous enzymes efficiently promote the hydrolysis of at least five different substrate classes (Figure 2) and stand out in particular as their promiscuous phosphodiesterase activity is almost as efficient as their native phosphonate monoesterase activity[8,12] (Table S1); note that the PTE activity reported in this work is ambiguous, as discussed in the Results and Discussion. Moreover, these PMHs provide the first example of biological PMH activity and are the only currently known enzyme capable of catalyzing the hydrolysis of xenobiotic sulfonate esters by direct S–O cleavage.[11,12] Note also that both enzymes are large homo tetramers with ∼56 kDa subunits and have extremely large binding pockets (≈10 × 20 Å2 wide and 15 Å deep).[8,12] Therefore, one would assume that such enzymes could easily accommodate a range of substrates of different shapes and sizes. Perhaps unsurprisingly, therefore, both PMHs are moderately efficient catalysts for the hydrolysis of the compounds shown in Figure 2 (kcat/KM values in the range of 103–104, see Table S1) and, in the case of RlPMH, apparently only marginally affected by mutations of the key active site residues with presumably multiple catalytic backups present in the active site (Table S2) that can take over the role of the mutated residues.

Figure 2

Structure of RlPMH (PDB ID: 2VQR) and the corresponding substrates studied in this work.[8] Both RlPMH and BcPMH are dimers of dimers, in which the monomeric units of each dimer communicate with its corresponding oligomeric unit through the C-terminal loop highlighted in this figure (which is in turn an adaptation from an analogous figure presented in ref (12)). This loop reaches into the adjacent active site, helping position key catalytic residues.[8] The substrates studied in this work are phenyl p-nitrophenyl phosphonate (PPP), ethyl p-nitrophenyl phosphate (PET), p-nitrophenyl sulfate (PNS), phenyl p-nitrophenyl sulfonate (PPS), and the p-nitrophenyl phosphate monoanion (PNPH).

A closely evolutionarily related enzyme in the AP superfamily is the arylsulfatase from Pseudomonas aeruginosa (PAS) (Figure 1C).[19] This enzyme only shares about 27% sequence similarity to RlPMH but has high structural similarity, in that 64% of the residues between the two enzymes structurally align with an RMSD of 2.54 Å.[8] This enzyme also has recently been the subject of extensive experimental[9,28] and computational[13,29] studies. Both PAS and the PMHs contain a mononuclear metal center with distorted octahedral conformation, which is most likely Mn2+ in the PMHs[8,12] and Ca2+ in PAS.[19] In addition, all three enzymes use an unusual geminal diol nucleophile[8,12,28] (Figure 1), a feature they share with all known sulfatases.[30] This noncanonical residue is a post-translationally modified cysteine or serine, which is first converted to an aldehyde and then hydrated to give rise to its active form.[31] Despite these apparent similarities, the two PMHs and PAS have very different specificity patterns. That is, PMHs are phosphonate monoesterases, while PAS is primarily a sulfatase, although all three enzymes have relatively low discrimination between native and promiscuous substrates.[8,9,12,28] In contrast, other superfamily members such as AP and NPP have dinuclear zinc centers in their catalytic sites (AP also possesses a third metal ion that appears to play an important role in determining the activity)[23] and utilize ionized serine or threonine residues as nucleophiles, respectively[10] (Figure 1). Thus, a direct atomic-level comparison between individual AP superfamily members and also related promiscuous phosphatases can provide better and broader understanding of the features that drive selectivity and promiscuity in these highly multifunctional systems.

In the present work, we have performed an extensive number of empirical valence bond (EVB)[32,33] simulations (total simulation time of ∼4 μs) of both the native and several characterized promiscuous activities of BcPMH and RlPMH, reproducing key experimental observables such as activation barriers, qualitative predictions from the pH-rate profiles for BcPMH activity, and the effect of mutations on RlPMH activity.[8,12] We demonstrate that, despite their broad promiscuity, the PMHs studied in this work hydrolyze all substrates through a unified mechanism with similar substrate binding positions, transition states, and electrostatic contributions to transition-state stabilization. Additionally, we showcase the importance of compensatory and cooperative electrostatic effects, which allow for an electrostatically flexible active site environment that can accommodate a range of substrates with different charge distributions, transition-state geometries, and requirements for efficient catalysis.

Finally, in order to test whether these observations are general to other members of the superfamily, we provide a detailed comparison of a range of AP superfamily members, in terms of active site shape, volume, and polarity. From this analysis we find a strong correlation between these properties and both substrate charge preference and number of known promiscuous activities, once again emphasizing the central role of the electrostatic environment of the active site in determining enzyme specificity and promiscuity. It is commonly accepted that enzymes achieve their tremendous catalytic proficiencies through an exquisite network of interactions that preferentially stabilize their transition states over their ground states,[34] and it has been argued that this is achieved through preorganization of the catalytic residues into an optimal conformation for transition-state stabilization.[35,36] This has been demonstrated for a wide range of systems through both experimental and computational work.[37−39] We illustrate here that while having an optimal electrostatic environment is clearly important to the catalysis of these enzymes toward individual substrates, one should also take into account the cooperativity between these residues, where the effect of the combined electrostatic environment from all relevant residues on the transition state stabilization is greater than the effect of each residue determined individually. We demonstrate that this cooperativity renders the active site electrostatic environment sufficiently flexible to accommodate a broader range of substrates with different electrostatic needs for efficient catalysis (without necessarily altering either substrate binding position or enzyme conformation). Additionally, our comparative analysis of different alkaline phosphatases shows that the higher the number of polar active side residues, the greater the propensity toward catalytic promiscuity. This highlights the importance of such cooperative electrostatic interactions as a common feature to functional evolution among members of the AP superfamily,[13] illustrating the power of subtle amino acid substitutions to drive very different solutions for the same chemical problem.

Methodology

Initial System Setup

Initial structures for both RlPMH[8] and BcPMH[12] (1.42 and 2.40 Å resolution, respectively) were obtained from the Protein Data Bank[40] (accession codes 2VQR and 2W8S). As the deposited structure for RlPMH contains only the monomeric unit without the transformation matrix, the structure of the full tetramer was obtained directly from the authors.[8] Potential flips of histidine, asparagine, and glutamine side chains were evaluated using the MolProbity server,[41] and those suggested by the software were applied to the structure. In all cases, the substrates were placed manually in the active site in such a way to optimize nonbonded interactions between the substrate and nearby amino acid side chains, including charge–charge, hydrogen-bonding and hydrophobic interactions. Structures for the corresponding Q13A, N78A, Y105A, T107A, H218A, and K337A variants presented in ref (8) were generated by manual truncation of the relevant side chains starting from the wild-type crystal structure, and structures were equilibrated using the same protocol as for the wild-type enzymes in order to allow the active site to adapt to the perturbation introduced.

Both PMHs are metalloenzymes with a single metal per active site of the tetramer (4 total), and the most likely candidate for this role has been identified as being a divalent manganese ion.[8,12] We recently presented a set of force-field-independent parameters to describe a range of alkali earth and transition-metal centers[42] based on Åqvist and Warshel’s original cationic dummy model,[43] which describes the metal as a delocalized charge spread over a number of dummy atoms placed around the metal center (in this case six particles in octahedral coordination, as shown in Figure S1). These particles are bonded to the central atom and to each other, and the frame is allowed to freely rotate in its coordination sphere without the need for external constraints or artificial bonds. We demonstrated[42] that this model also allows one to capture subtle structural effects upon metal substitution without the need for the artificial restraints that need to be imposed in a fully bonded model, while simultaneously capturing key electrostatic properties of the metal center. We have successfully used our Ca2+ model in simulations of the selectivity of PON1,[44] and the Mn2+ model presented in our original paper[42] has been used in the present work to describe the catalytic metal center in the PMHs.

All relevant reactions were simulated in the active sites of both enzyme species and in a 24 Å water droplet, in order to quantify the catalytic effect of the enzyme compared to background reaction in aqueous solution. To model the reaction in solution, we used truncated residues to model the nucleophile (acetaldehyde hydrate as a model for the formylglycine) and the relevant general acids (ethylamine and ethylimidazole for Lys and His, respectively). In the enzyme simulations, one of the main computational difficulties encountered comes from the fact that truncating the 16 C-terminal residues of RlPMH causes the enzyme to lose its tetrameric structure, with a corresponding loss of activity. This strongly suggests that interactions at the subunit interfaces can be important to catalysis, as can also be observed from the protrusion of the interfacial loop almost into the active site of the adjacent subunit (see Figure 2). Thus, it is necessary to include the entire (2056 amino acid) tetramer, which creates substantial computational cost.

To simplify this problem and reduce computational cost, the system was divided into three layers: the EVB (reacting) atoms, an active region encompassing all residues within a 24 Å sphere of the reacting atoms centered at the metal center, and an external layer in which the remainder of the system was present, but the atoms were constrained to their crystallographic positions (as is commonly done in similar studies, see e.g., refs (24 and 45)). The simulation sphere encompassing the active region was centered on the catalytic Mn2+ ion, and all crystallographic water molecules within 18 Å of this center were retained in our simulations, with the exception of any crystal waters clashing directly with the substrate once it was placed in the active site. The solvation sphere was then completed and extended to 24 Å using TIP3P[46] water molecules subjected to the surface constraint all atom solvent (SCAAS) spherical boundary conditions.[47] A 10 kcal·mol–1·Å–2 harmonic restraint was applied to the outer layer of the active region and associated solvent molecules (15%, 3.6 Å), in order to ease the transition between the active and constrained regions, which is why only an 18 Å of crystallographic water molecules were retained for the simulation. All forces on the constrained atoms were set to zero, in order to prevent them from distorting the dynamics of the active region. Ionizable residues within ∼18 Å of the center were ionized during the course of the simulation, leading to a total system charge of −1 (without including the substrate). The protonation states of histidine side chains were investigated using the MolProbity server,[41] PROPKA 3.1,[48,49] and by visual inspection. All other residues, in particular those outside the active region, were set to their neutral form for system stability. The Mn2+ ions in the adjacent monomeric units (the positions of which were kept constrained) were removed in order to avoid the presence of residual charge outside the simulation sphere (we note that the adjacent active sites all fall within the constrained external layer and all surrounding residues are therefore not allowed to move). In contrast, the catalytic metal center in the active region was described using a 7-pointed dummy model with distributed charges as described above.

Molecular Dynamics Equilibration of the Systems of Interest

All molecular dynamics (MD) and EVB simulations in this work were conducted using the OPLS-AA force field[50] implemented in the Q simulation package (Version 5.0.6).[51] For the substrate and nucleophile, OPLS-AA compatible force field parameters were generated with Macromodel 9.1 (OPLS-AA force field, 2001, Schrödinger LLC).[52] The only exceptions were the force field parameters for the carbon and oxygen of the deprotonated geminal diol, which were available in the literature and obtained from ref (53). Partial charges for the reacting atoms were generated at the HF/6-31G* level of theory using the Gaussian 09 simulation package,[54] followed by the standard RESP procedure.[55]

All the simulations performed herein used time steps of 1 fs, while the temperatures of the system were regulated using the Berendsen thermostat[56] (with a 100 fs bath relaxation time). The systems were initially heated from 1 to 300 K over a short 80 ps simulation, applying a 200 kcal·mol–1·Å–2 harmonic force constant on the solute atoms to restrain them to their crystallographic positions. This allowed for the solvent molecules to equilibrate around the protein and the removal of initial contacts due to substrate placement. The system was then cooled down to 5 K for another 10 ps and then gradually heated to 300 K for 90 ps of simulation time, while the force constants of the harmonic restraint were gradually decreased from 200 to 0.5 kcal·mol–1·Å–2. Subsequently, a 5 ns equilibration was performed at 300 K for both wild-type and mutant enzyme simulations (300 ps for the reference reaction in solution) using a 0.5 kcal·mol–1·Å–2 position restraint on the substrate atoms, the side chain of the nucleophile, the catalytic metal center, and the side chain of the general acid (H218 or K337, depending on the mechanism being considered) to keep the reacting atoms in place. An RMSD plot of the active monomer for the wild-type enzyme and each enzyme variant is shown in Figures S2 and S3. As shown in this figure, due to the fixed excluded region, these systems equilibrated rapidly, with RMSD of < 0.5 Å from the reference crystal structure. After the final equilibration step, we ran an additional 500 ps of molecular dynamics, during which 10 snapshots of the whole system were taken every 50 ps to be used as starting points for subsequent EVB simulations. Finally, although we remained as faithful as possible to a fully nonbonded model for the catalytic metal center, we introduced an angle parameter (50 kcal·mol–1·rad–2, equilibrium angle 180°) between the center of the Mn2+ dummy model and the residue D324 (Mn2+–Ometal–Ofree, where Ometal corresponds to the oxygen atom closest to the Mn2+, and Ofree the one not coordinated to it), which would otherwise become bidentately coordinate to the metal center and make the active site unstable. As can be seen from the Results and Discussion, despite the inclusion of this extra parameter, we are able to systematically reproduce the activation energies of both wild-type and mutant forms of these PMHs with different substrates with good agreement to experimental data.

Empirical Valence Bond Calculations

Our methodology of choice in this work to model chemical reactivity was the EVB approach of Warshel and co-workers.[32,33] This is an empirically-based multiscale valence-bond/molecular mechanics approach that is fast enough to allow for the extensive sampling required to obtain convergent free energies for complex biochemical processes, while having a proven track record as a powerful tool for quantifying and rationalizing the catalytic power of native and mutant enzymes.[36,44] All EVB calculations were performed using the standard EVB free energy perturbation/umbrella sampling (EVB-FEP/US) procedure outlined in refs (33 and 57), as implemented in the Q simulation package.[51]

The reaction under study was described in terms of two valence bond structures, as illustrated in Section S4 of the Supporting Information. It should be pointed out that as all atoms in the two valence bond states are treated using the same force field, the only differences between the reacting (EVB) and nonreacting atoms are the use of Morse rather than harmonic potentials to describe bonds that are being broken and formed during the reaction (see the Supporting Information) and the fact that unlike the rest of the protein, the EVB atoms do not have a cutoff for calculation of the nonbonded interactions. All EVB-FEP/US simulations were performed at 300 K, using 51 mapping windows of 200 ps per window, resulting in 10.2 ns of simulation time for each individual trajectory, sampling over 10 starting conformations per system (102 ns per system) and ∼4 μs cumulative simulation time for all systems studied in this work. All MD and EVB simulations were performed using a 1 fs time step, and long-range effects were treated using the local reaction field (LRF) approach.[58]

Finally, as outlined above, we also modeled the corresponding uncatalyzed reaction for each substrate of interest in this work, as we needed these calculations not only for the calibration of the EVB parameters, but also to compare the reactions for different environments. Coordinates for each system were based on the equilibrated enzyme system, but in the absence of the enzyme itself, and the reactions were modeled using the relevant substrate as well as acetaldehyde and protonated ethylamine as models for the nucleophile and the general acid, as described in the Initial System Setup section. All equilibration and EVB protocols were the same as for the full enzyme system, with the exception that the background reaction in the absence of the enzyme was only equilibrated for 300 ps rather than 5 ns. As with the corresponding enzymatic reactions, a weak position restraint of 0.5 kcal·mol–1·Å–2 was placed on all solute atoms (nucleophile, substrate, and model for general acid) in order to keep the reacting fragments in the center of the simulation sphere. This weak restraint is sufficient to keep the nitrogen atom of the general acid within 3.5 Å of the leaving group oxygen throughout the 300 ps equilibration of the background reaction, in part due to electrostatic interactions between the charge on the general acid and the substrate. The relevant background reaction was calibrated based on estimations using experimental data, as described in detail in the Supporting Information. Furthermore, the parameters describing the relative positions of the VB parabolas and the coupling between them were then transferred unchanged to the enzyme in order to be able to quantify and correctly predict the catalytic effect of wild-type and mutant enzymes.

Results and Discussion

Modeling the Catalytic Mechanism for the Wild-Type Enzymes

A unique feature of the PMHs being considered in the present work is that they are the only nonsulfatase enzymes known to date to possess a post-translational modification from a cysteine to an aldehyde (formylglycine, fGly).[8,12,19] The current proposed mechanism[8] for both native and promiscuous PMH activities is shown in Figure 3. In a first step (I → II), hydration of the post-translationally modified aldehyde yields a reactive geminal diol, which can act as a nucleophile. This geminal diol is activated by the catalytic metal center and exists in its alkoxide form. Following substrate binding (II → III), this geminal diol then attacks the phosphorus/sulfur center (III → IV) of the relevant substrate (Figure 2) to give rise to a hemiacetal intermediate (IV). In the final step (IV → I), this intermediate is hydrolyzed by hemiacetal cleavage to regenerate the aldehyde and yield the final product.

Figure 3

(A) Overlay of the active sites in RlPMH and PAS, illustrating conservation of active site structure between the two enzymes. (B) A simplified version of the proposed catalytic mechanism for both PMHs considered in this work, based on refs (8 and 12). Since the catalytic metal is suggested to be Mn2+, which is a hard Lewis acid, the nucleophile could be stable as an alkoxide, in agreement with the pH-rate profiles shown in Figure S4.

The two PMHs considered in this work show an absolute dependence on divalent metal ions, with the most likely candidate for fulfilling this role being Mn2+, based on experimental data presented in refs (8 and 12). A transition metal would be expected to substantially decrease the pKa of the metal-bound nucleophile to yield an alkoxide, as also suggested by the acid limb of the pH-rate profiles shown in Figure 2 of ref (12) (reproduced as Figure S4), which most likely corresponds to the deprotonated nucleophile (see discussion in ref (12)). The pH-rate profiles, which are coincidental for all substrates except the phosphate triester, also suggest the involvement of an acid catalyst, most likely either H218 or K337[12] (see also Figures 3 and S4).

It has been argued that steps IV → I of Figure 3 can play an important role in facilitating promiscuity.[12] Specifically, harnessing hemiacetal cleavage allows for a common mechanism irrespective of the functional group used in the intermediate,[12] while simultaneously providing a thermodynamically less challenging route to facilitate C–O cleavage, compared to the repeated cleavage of an extremely stable P(S)–O bond.[59] However, kinetic data on base-catalyzed hemiacetal cleavage in aqueous solution demonstrate that this reaction is extremely fast.[60] Additionally, as all substrates will be broken down by a common mechanism through a common intermediate, the selectivity will be already determined in steps III → IV, which is also therefore the focus of the present work (Figure 3).

Our mechanistic model assumes an anionic nucleophile and general acid catalysis from either K337 or H218 to protonate the departing leaving group. As discussed below, in this work K337 was chosen as the general acid based on empirical pKa calculations and experimental data. That is, the experimental pH-rate profiles[12] suggest a two-pKE model, with a pKE1 of 7.0–7.2 (5.8 for the sulfate monoester) and a pKE2 of 7.5–8.1 (see Figure S4). The first pKa is likely to correspond to the nucleophile, contributing to catalysis in its deprotonated form as discussed above, and the second pKa to the general acid. The pKE2, which is very close to the pKa of around 8 suggested for K337 by PROPKA, led to the choice of this residue as the putative general acid, as also suggested by refs (8 and 12). The only exception to this model is the p-nitrophenyl phosphate monoester, of which the dianionic form will be extremely resistant to attack by an anionic nucleophile (see ref (62)). However, the pKa of the already basic nonbridging oxygens of this substrate (pKa ∼ 5.0[63]) is likely to be substantially elevated due to the close proximity of the anionic nucleophile. This, in turn raises the possibility that this substrate binds as a monoanion, as has already been demonstrated by simulations, for example, for protein tyrosine phosphatase 1B.[64] Note also that, as shown in Figure S4 and ref (12), the monoanionic sulfate and dianionic phosphate monoesters give rise to very different pH-rate profiles that not only have different slopes but also are shifted by 2 pH units. Thus, we have herein considered a mechanism involving an anionic nucleophile attacking a monoanion phosphate, which yields excellent agreement with experiment as discussed below.

In the present work we have not examined phosphate triester hydrolysis, which shows a very different kcat/KM pH-rate profile to all other substrates studied (Figure S4). The inverted pH-rate profile observed for this substrate suggests either the involvement of a completely different set of residues or a completely different mechanism of catalysis. Additionally, BcPMH shows extremely poor activity toward this substrate (kcat/KM of 1.6 × 10–2 M–1·s–1), which is actually slower than the corresponding uncatalyzed alkaline hydrolysis of the model substrate paraoxon.[65] Taken together, this suggests that the hydrolytic mechanism of this substrate, if it at all binds in the same active site, is impossible to prove conclusively through calculations due to lack of concrete experimental data.

As mentioned before, the reactions examined in this work correspond to the third step (III → IV) of the catalytic cycle shown in Figure 3. Since experiments show that the hemiacetal cleavage is a fast step, we focused only on this second step, which is the most chemically challenging step of the cycle, being thus the one related to the measured kinetic parameters. Figure 4 shows a comparison between our calculated and, where available, experimental activation free energies (derived from kcat, which provides an upper limit for the reaction rate).[8,12] The corresponding tabulated values can be found in Table 1. From our results, it can be seen that the model used in the present work reproduces the experimental activation free energies within an accuracy of 1.7 kcal·mol–1 for all substrates. It has additionally been argued[8] that the PMHs considered in this work can accept both diesters and phosphonates with such high proficiency in the same active site due to similar geometrical and steric demands for the respective substrates and transition states. To probe this further, we have examined transition-state geometries for all uncatalyzed and enzyme-catalyzed reactions considered in this work. Table 2 shows a comparison of P(S)–O distances to the oxygen atoms of the incoming nucleophile and departing leaving group for all substrates and reactions. Representative transition-state structures in the BcPMH active site are also illustrated in Figure 5. From these results, it can be seen that the PMHs hydrolyze all substrates through a unified mechanism with similar substrate binding positions and transition states. With the exception of the phosphonate, little change is seen in transition-state geometry upon moving from aqueous solution to the enzyme active sites, in agreement with related experimental work by Herschlag and co-workers,[22] as well as theoretical analysis by Hou and Cui,[25] on alkaline phosphatase. Even in the case of the phosphonate, the overall transition-state size (considering the distance between Onuc–Olg) stays very similar, and the main change is that the symmetry of the transition states changes, with P–Onuc becoming slightly elongated and P–Olg slightly compressed compared to the corresponding uncatalyzed reaction. Hence, as suggested in previous works[22,25] for alkaline phosphatase, we find very little effect on the transition-state geometries of moving to the enzyme active when compared with those obtained through modeling the corresponding uncatalyzed reaction in aqueous solution.

Figure 4

Calculated and experimentally derived activation energies for the enzyme-catalyzed reactions of the five substrates studied here by the wild-type forms of (A) RlPMH and (B) BcPMH.[61]

Table 1

Calculateda and Observedb Activation (ΔG⧧) and Reaction Free Energies (ΔG°) for the Hydrolysis of the Five Substrates by the Wild-Type Forms of RlPMH and BcPMHc

	ΔG⧧		ΔG°	ΔG⧧		ΔG°
substrate	calc.	expt.	calc.	calc.	expt.	calc.
	RlPMH			BcPMH
PPP	16.3 ± 0.7	15.9	–16.9 ± 1.1	17.2 ± 1.6	17.0	–15.7 ± 2.8
PET	15.1 ± 1.8	16.4	–6.3 ± 2.0	15.9 ± 0.7	16.7	–6.3 ± 1.1
PNS	19.5 ± 1.2	n.d.	–2.0 ± 1.5	21.3 ± 1.3	19.7	–6.2 ± 1.9
PPS	19.4 ± 1.7	n.d.	1.3 ± 2.1	19.7 ± 2.5	20.2	3.2 ± 3.6
PNPH	19.1 ± 0.3	19.1	–9.2 ± 1.0	20.0 ± 0.8	20.7	–7.8 ± 1.2

“Expt.” and “calc.” denote experimental and calculated values respectively, and “n.d.” refers to values that have not been experimentally determined. All energies are given in kcal·mol–1 and are averages and standard deviations based on 10 individual EVB simulations generated from different starting structures, as outlined in the Methodology section.

ΔG⧧(expt.) corresponds to experimental values of the enzyme-catalyzed reaction, based on the kinetic data presented in refs (8 and 12).

The corresponding EVB parameters are presented in the Supporting Information.

Table 2

Average P(S)–O Distances, in Å, at the Transition State for the Relevant Group Transfer Reaction in Water and in the Wild-Type forms of RlPMH and BcPMHa

	water		RlPMH		difference RlPMH – water		BcPMH		difference BcPMH – water
	P/S–Onuc	P/S–Olg	P/S–Onuc	P/S–Olg	P/S–Onuc	P/S–Olg	P/S–Onuc	P/S–Olg	P/S–Onuc	P/S–Olg
PPP	1.89 ± 0.08	2.14 ± 0.12	2.09 ± 0.11	1.96 ± 0.09	0.20	–0.18	2.07 ± 0.11	1.95 ± 0.08	0.18	–0.19
PET	2.09 ± 0.29	2.04 ± 0.12	2.06 ± 0.11	1.94 ± 0.08	–0.03	–0.10	2.09 ± 0.13	1.97 ± 0.09	0.00	–0.10
PNS	2.13 ± 0.17	2.06 ± 0.13	2.11 ± 0.13	2.08 ± 0.12	–0.02	0.02	2.12 ± 0.13	2.08 ± 0.12	–0.01	0.02
PPS	1.85 ± 0.07	2.03 ± 0.09	1.88 ± 0.07	1.99 ± 0.08	0.03	–0.01	1.84 ± 0.06	1.99 ± 0.08	0.03	–0.01
PNPH	1.99 ± 0.25	2.14 ± 0.18	2.03 ± 0.11	2.0 ± 0.11	0.04	–0.14	2.03 ± 0.10	2.03 ± 0.10	0.04	–0.11

All data are averages and standard deviations over 10 individual simulations as outlined in the Methodology section. For an extended version of this table, including both P(S)–O and H–N/Olg distances, we refer the reader to the Table S4.

Figure 5

Representative transition-state structures for the BcPMH catalyzed hydrolysis of (A) PPP, (B) PET, (C) PNS, (D) PPS, and (E) PNPH. Values presented here correspond to the averages over 10 trajectories.

To locate the origin for the differences in the observed catalytic activity, we have performed a comparison of the electrostatic contributions of individual residues to the hydrolysis of each substrate, calculated using the linear response approximation following previous works (e.g., refs (66 and 67); see Figure 6). This analysis shows that, strikingly, despite the calculations of residue interactions being completely independent of each other, with different charge distributions and different transition states, all residues that make substantial electrostatic contributions to the calculated activation barrier are conserved among different substrate. This is similar to the observations from our previous computational work on the arylsulfatase from PAS.[13] However, although qualitatively similar, there are some key quantitative differences between the different substrates, most notably in the case of D12, D324, H325, and E327. While this itself is hardly surprising, considering these are electrostatically quite different substrates, it highlights that the active site pre-organization is not “perfect”, but rather cooperative electrostatic interactions render this preorganization flexible enough to readily adapt to the electrostatic needs of different substrates (see also the related discussion of catalytic backups in serum paraoxonase 1).[68]

Figure 6

Electrostatic contribution of key residues to the calculated activation barrier for the hydrolysis of the five substrates for wild-type form of BcPMH.[72]

To further explore this observation, we have also examined the charge change on the central P/S atom and all atoms bonded to it upon moving from reactant to transition state for the wild-type BcPMH catalyzed hydrolysis of the different substrates studied in the present work. These atoms were split into three fragments: a central fragment comprising the P/S atom, the nonbridging oxygens of the substrate, and the C atom of the phenyl group of phenyl p-nitrophenyl sulfonate (PPS) and phenyl p-nitrophenyl phosphonate (PPP) connected to the central P/S atom as well as the oxygen atoms of the departing leaving group (Olg) and attacking nucleophile (Onuc) as individual fragments (see Table S3). As the transition states for the reactions studied involve partial bond formation to the nucleophile and leaving group, we have summed up the charges on the central atoms (P/S, nonbridging oxygens, and carbon) and treat this as one unit, which we will henceforth refer to as “central fragment” for simplicity. The schematic for this division is shown in Table S3 which also provides absolute charges for each fragment at the reactant and transition states as well as the charge shift upon moving from reactant to transition state. These have then in turn been ranked against the measured kcat/KM for BcPMH for each substrate.[12] From this table, it can be seen that while there is little trend in the charge shift on the leaving group oxygen (which is partially protonated by the general acid), there are subtle but clear trends in the partial charges of the nucleophile oxygen and the central fragment. That is, for the native substrate, PPP, there is a substantial charge shift corresponding to a loss of +0.2728 au on Onuc. This charge shift gradually decreases across the series, correlated with a reduction in kcat/KM. In parallel to this, for the native substrate, there is a small buildup of negative charge (−0.0468 au) on the central fragment, which increases to −0.0884 au for the promiscuous activity with the lowest observed kcat/KM (the sulfate monoester). This suggests a subtle preference for minimizing negative charge at the transition state. Therefore, there appears to be a correlation between the calculated charge shift at the transition state and the subsequent catalytic efficiency of the enzyme. This also ties in with experimental observations that other alkaline phosphatases such as AP and NPP clearly discriminate on the basis of substrate charge.[69−71] The most radical example of such charge discrimination in this superfamily, in fact, is in the arylsulfatase from Pseudomonas aeruginosa (PAS), where the hydrolysis of large bulky monanionic diesters such as bis-p-nitrophenyl phosphate shows only 100-fold lower values than the monanionic substrate p-nitrophenyl sulfate (with kcat/KM = 4.9 × 107 M–1·s–1 for the sulfate monoester and 2.5 × 105 M–1·s–1 for the phosphate diester).[9] In contrast, the dianionic analogue of the sulfate monoester, p-nitrophenyl phosphate, is a much poorer substrate than the sulfate monoester (kcat/KM = 790 M–1·s–1),[9] despite having the same ground-state geometry and similar predicted transition-state geometries to the sulfate monoester. Even further examples of such charge discrimination have been seen by Baxter and co-workers[73] in studies of aluminum and magnesium fluoride transition-state analogues (TSA) of phosphoryl transfer enzymes, where they showed clear preference for preserving anionic charge at the expense of TSA geometry over a broad range of pH.

Clearly, these enzymes have evolved to provide the key active site interactions that optimally stabilize the transition state for the native reaction, which in turn leads to the observed preference for the native substrate. However, in the electrostatic cooperativity model we present in this work, these interactions are sufficiently flexible to accommodate the electrostatic needs of other substrates, although the same residues can make quantitatively different contributions as shown in Figure 6. These differences in ability to stabilize different transition states would in turn lead to the selectivity displayed by these enzymes for different substrates.

Probing Key Active Site Mutations

Hollfelder and co-workers[8] have performed a detailed alanine scan of the RlPMH active site residues, testing against both the phosphonatase (PPP) and phosphodiesterase (PET) activities of the enzyme. Both substrates appear to be highly insensitive to active site single mutations, with the individual substitution of each key active site residue in RlPMH leading to, at worst, a ∼20-fold reduction in kcat. An exception to this is modifying the nucleophile to alanine, but even in this extreme case, these enzymes still show some activity.[8,12] To probe the origin of the seeming resilience of these enzymes to substitution of individual active site residues, we have performed EVB calculations to obtain free energy profiles for the chemical step for the hydrolysis of substrates PPP and PET of Figure 2 by a range of Ala-substituted forms of RlPMH presented in ref (8). Critically, when moving from wild-type to mutants, we used exactly the same parameter set, unchanged, allowing us to in parallel rigorously validate our valence bond model for the reaction mechanism catalyzed by these enzymes and its predictive power (see the Supporting Information for theoretical background). The only exception to this is the K337A mutant that has the general acid functionality of K337 removed, and for which we model the reaction as proceeding with H218 as the general acid instead (for all other calculations, H218 is kept in its neutral form as close proximity to the K337 side chain and the catalytic metal ion will depress its pKa). Note again that although we have focused on K337 as a putative general acid in this work, due to the agreement between the predicted pKa of K337 and the experimental pH-rate profiles as outlined in the previous section, in practice either of these two residues could fulfill the role of general acid.

Table 3

Calculateda and Observedb Activation (ΔG⧧) and Reaction Free Energies (ΔG°) for the Hydrolysis of PPP and PET in Both the Wild-Type and Different Mutant Forms of RlPMHc

	ΔG⧧		ΔG°	ΔG⧧		ΔG°
system	expt.	calc.	calc.	expt.	calc.	calc.
	PPP			PET
WT	15.9	16.3 ± 0.7	–16.9 ± 1.1	16.4	15.1 ± 1.8	–6.3 ± 2.0
Q13A	15.6	15.6 ± 1.5	–18.0 ± 1.4	16.0	14.8 ± 1.0	–10.5 ± 1.4
N78A	17.3	17.6 ± 1.2	–18.4 ± 1.7	17.7	17.0 ± 0.9	–5.8 ± 2.0
Y105A	17.0	17.2 ± 3.8	–13.8 ± 3.6	18.8	19.7 ± 1.7	–5.6 ± 1.9
T107A	17.1	17.3 ± 0.9	–9.9 ± 1.8	18.1	17.5 ± 0.7	–3.3 ± 1.8
H218A	17.6	17.3 ± 1.1	–13.4 ± 1.2	18.5	17.5 ± 0.7	–5.7 ± 1.2
K337A	16.1	16.2 ± 0.9	–7.1 ± 1.3	17.1	17.5 ± 0.8	–2.2 ± 1.9

“Expt.” and “calc.” denote experimental and calculated values, respectively. All energies are given in kcal·mol–1 and are averages and standard deviations over 10 individual trajectories using different starting conformations, as outlined in the Methodology section.

ΔG⧧(expt.) corresponds to experimental values of the enzyme-catalyzed reaction, based on the kinetic data presented in refs (8 and 12).

The corresponding EVB parameters are presented in the Supporting Information.

A comparison between calculated and experimental activation barriers for the hydrolysis of these substrates by each key RlPMH variant is shown in Figure 7 with the corresponding energetics presented in Table 3. It can be seen that we are able to reproduce the experimentally observed effect of all active site mutants (ΔΔG⧧WT→mut) to within an average error margin of ∼1 kcal·mol–1. The most challenging of these mutations to model is the Y105A mutation (as can also be seen from the large error bar shown in Figure 7), as this mutation leads to a larger perturbation in the active site, for example repositioning residues such as T107 and Y215. This results in a larger standard deviation for these calculations in the case of the phenyl phosphonate than in other simulations. However, the average over 10 trajectories is still in good agreement with experimental results. Taken together with our other simulations, this provides support for the quality of our calculations, the suggested mechanism, and our assumption that the group transfer is the key step in determining the specificity. Additionally, an examination of the corresponding P–O distances at the transition state for each variant and substrate shows that, as we move from the background reaction to the enzyme (Table S5), the single mutants in the active site have little effect on the transition-state geometry.

Figure 7

Calculated and experimentally derived activation energies for the RlPMH catalyzed reactions of (A) PPP and (B) PET. Shown here is data both for the reaction catalyzed by the wild-type enzyme as well as several mutant forms of the enzyme. The data plotted in this figure are presented in Table 3, and the error bars represent standard deviations over 10 independent trajectories.

Here, we will provide a brief discussion of our simulations of each mutant below and refer the reader to Figure 5 for an overview of how each residue interacts with the substrates of interest in the equilibrated wild-type enzyme.

Q13A

As suggested in ref (8), this residue seems to play a key role by holding K337 in place. For the simulations of the Q13A variant, we observe that both K337 and H218 are perturbed, and the RMS displacement of K337 after equilibration compared to the wild-type enzyme is 0.80 and 0.78 Å for PPP and PET, respectively. However, the effect of losing this interaction translates to only a 0.4 kcal·mol–1 reduction in activation barrier for both substrates both experimentally and from our simulations (see Table 3). Note that, as shown in Table S2, this mutation primarily affects KM rather than kcat for both substrates considered here.

N78A

N78 is a key active site residue, as it provides a hydrogen-bonding interaction to the nonbridging oxygen of both PPP and PET. This interaction stabilizes the substrate and also helps to optimally position the substrate in the active site. Loss of this interaction results in a 1.4 kcal·mol–1 increase in barrier for both substrates. As shown in Table 3 and Figure 7 we are able to reproduce the detrimental effect of the N78A mutant, which is slightly larger for PET than for PPP. This could be in part due to the presence of the phenyl ring in PPP, which can help to position the substrate in the active site even in the absence of the hydrogen bond from N78.

Y105A

In the wild-type enzyme this residue does not directly interact with the nucleophile or the substrate. However, it is part of a hydrogen-bond network that keeps D324 and R61 correctly positioned for catalysis (Figure 1). It has been experimentally shown that this mutation drastically reduced the kinetic efficiency of the enzyme. As can be seen from Table 2, while we are able to reproduce the experimental activation barrier within 1 kcal·mol–1, for both the phosphonate monoester and the phosphate diester, we obtain larger standard deviations for this variant in the case of the phosphonate monoester. In the case of the phosphate diester, upon truncating Y105 to alanine and equilibrating the system, we see an increase in the number of water molecules around the nucleophile as well as repositioning of other residues, such as T107 and Y215. Specifically, the interaction between the nucleophile and T107 is broken, and Y215 occupies the space left by mutation. In the case of the phosphonate monoester, however, the large hydrophobic phenyl ring of the phosphonate (compared to the smaller ethyl group in the diester) blocks water access to the nucleophile and restricts the movement of other residues around it.

T107A

T107 is another key active site residue, as it directly interacts with the nucleophile, helping optimally position it for catalysis (Figure 1). Unsurprisingly, truncating this residue to alanine leads to an increase in activation barrier of 1.2 kcal·mol–1 for the phosphonate and 1.7 kcal·mol–1 for the diester (a trend we reproduce computationally, see Table 3 and Figure 7). This is primarily due to loss of the hydrogen bonding interaction with T107 as well as the resulting subtle repositioning of the nucleophile in the active site.

H218A

Interestingly, this mutation appears not only to impact the catalytic activity, but also to increase KM from 2.8 mM in the wild-type enzyme to 15 and 57 mM in the mutant form for both phosphonate and diester substrates.[8] It has been argued that this increase in KM is due to substrate binding, suggesting that this residue as well as K337 are directly involved in this step. H218 could also play a role as a general acid, due to its close proximity with the leaving group. However, in its unprotonated form, H218 also plays a key role in positioning K337 for optimal leaving group stabilization (the distance between Nε of H218 and K337 is 3.44 Å in the RlPMH crystal structure).

K337A

The role of K337 is two-fold: it helps to position the substrate in the active site in the Michaelis complex (through a hydrogen bonding interaction to one of the nonbridging oxygens of the substrate) as well as to stabilize the leaving group upon departure by acting as a general acid and protonating it. One would expect, then, that mutation of this key residue to alanine would result in a substantial increase in activation barrier. However, the experimentally observed increase is only 0.2 kcal·mol–1 for the phosphonate monoester and 0.7 kcal·mol–1 for the phosphate diester, which is lower than for example either the N78A or T107A mutations. This suggests that another positively charged residue is taking up the role of K337 in leaving group stabilization. As discussed above, the role of general acid could be fulfilled by either K337 or H218, and in absence of the close proximity of the K337 positive charge upon mutation (the distance between the Nε of H218 and the nitrogen of K337 is 3.44 Å in the wild-type crystal structure),[8] one would presume that H218 is more likely to be protonated than in the wild-type. Therefore, we tested modeling H218 as a general acid, demonstrating that this in fact provides activation barriers in very good agreement with experiment (Table 3), suggesting that in the absence of K337, H218 takes up the role of this residue in leaving group stabilization (either through hydrogen-bonding/charge–charge interactions with the anionic leaving group or as a general acid).

So far, we have not yet discussed the details of the proton transfer to the leaving group from either K337 or H218. In the present work, we have used a two-state valence bond model to describe this process, as outlined in the Methodology section. In our model, for both the lysine- and histidine-catalyzed mechanisms (wild-type and K337A mutants, respectively), the group transfer and proton transfer reactions take place in a single, concerted but slightly asynchronous reaction step (Figures S5 and S6 and Table S4–S6). As seen from these figures and associated table, when the general acid is modeled as being K337, the transition state is dominated by the group transfer reaction, with the proton transfer reaction taking place very slightly after the group transfer. This would tie in with the fact that the p-nitrophenol leaving group is sufficiently basic to not a priori need protonation to depart. In the case of H218 as the general acid (in the K337A mutant), the proton transfer becomes more concerted with the group transfer reaction (Figure S6). This would be in agreement with our previous DFT study of the hydrolysis of phosphate and sulfate monoesters,[18] in which we carefully examined all proton transfer steps involved and demonstrated that they either immediately precede or succeed the group transfer step, but along the same reaction coordinate (without the need for discrete intermediates). Such a model also agrees with high-level QM/MM calculations of the phosphoryl transfer reaction catalyzed by dUTPase,[74] which show a similar coupling between proton transfer and group transfer reactions, suggesting that a two-state VB model is adequate for capturing the key features of the relevant reaction mechanisms, and this is also borne out by the agreement with the experimental data.

Finally, a comparison of the electrostatic contributions of individual residues to catalysis for both substrates and all variants (Figure 8) shows that, as with the wild-type reactions, changes in activity correspond to cooperative electrostatic effects, where the active site residues are able to compensate the absence of key active site residues and stabilize the transition state of different substrates. This effect was also seen in our previous computational studies of the evolutionarily related PAS[13] and has also been alluded to in other recent works.[38,75] A similar phenomenon has been observed in experimental studies of serum paraoxonase 1,[68] suggesting that such electrostatic flexibility is a feature of multiple enzymes that catalyze phosphoryl transfer.

Figure 8

Electrostatic contribution of key residues to the calculated activation barrier for group transfer reactions of (A) PPP and (B) PET for different RlPMH variants.

Examining Other Plausible Contributions to the Observed Selectivity and Promiscuity

Although our data strongly point toward electrostatic cooperativity as the origin for the observed selectivity and promiscuity among members of this superfamily, it is important to also examine other possible origins of this effect. Before proceeding further in this discussion, it is worth mentioning that there are several different ways to define this concept,[76] that range from an enzyme performing distinct chemistry using a similar set of residues and the same mechanism to cases where different sites of an enzyme are used to perform different chemical reactions (a form of “protein moonlighting”). Common to all these definitions, however, is the fact that promiscuity can be regarded simply as a converse of specificity, in that a highly specific enzyme would only be able to perform a single chemical reaction, whereas a catalytically promiscuous enzyme would be able to perform multiple distinct chemical reactions.[3,77] In addition, while the enzymes studied in the present work are multifunctional, with very high proficiencies for both phosphonate monoester and phosphodiester hydrolysis, they nevertheless show high selectivity and an order of preference between these and other promiscuous reactions that they catalyze (see Table S1). Therefore, in the present discussion, we will use “specificity” as a converse to promiscuity (i.e., referring to the number of reactions the enzyme catalyzes) and “selectivity” to indicate the discrimination between different reactions catalyzed by the same enzyme.

An important point to take into account in the present work is that, with the exception of the monoesters, the reactivity of the substrates examined herein is substantially lower by several orders of magnitude under neutral conditions than at high pH (see Table S1). Therefore, it is plausible to consider that a part of the broad substrate specificity might result because the enhanced reactivity of the active site nucleophile by the catalytic metal center already provides substantial rate accelerations for a broad range of substrates. This would be consistent with the small effects of < 10-fold (on kcat) on the enzymatic activity of the mutation of they key residues that interact with the substrate oxygens, specifically N78 and K337.[8] However, while having an activated nucleophile is clearly important for the overall activity toward different substrate, this in itself is not fully sufficient to describe the observed promiscuity, as there are clear variations in rate acceleration between the different substrates, even when considering the alkaline reaction as the relevant reference state for the uncatalyzed reaction (see values presented in Table S1).

Additionally, as discussed and demonstrated in several of our previous works,[17,18,78] despite the superficial similarities between the different transition states involved (substitution of P for S, adding or removing functional groups), they have very different charge distributions and thus solvation patterns, leading to very different requirements for efficient catalysis. This can also be seen in both the quantitative differences in the residue contributions shown in Figure 6, and the fact that although multifunctional, BcPMH (for which more kinetic data is available with different substrates, see Table S1 and ref (12)) shows up to ∼25,000-fold differences in kcat/KM values toward different substrates. These range from 0.59 M–1·s–1 for the sulfate monoester to 1.5 × 104 M–1·s–1 for the phosphonate monoester. Following from this, there is also large sensitivity among alkaline phosphatases to the nature of the leaving group. For example, in the case of BcPMH, simply changing the leaving group to phenol substantially reduces the catalytic activity for all substrates by up to 350-fold (in terms of kcat/KM). In the case of AP, for which linear free energy relationships do exist, these also show moderate-to-strong leaving group dependence (see, e.g., refs (22, 79, and 80)). In particular, reported literature values for the AP-facilitated hydrolysis of sulfate monoesters phosphorothioates, phosphate monoesters and phosphate diesters all show steep leaving group dependence, with βlg values in the range from −0.76 to −0.95.[22,79,80] Therefore, although the transition states are superficially very similar, the enzyme is actually highly selective between the different substrates and leaving groups.

As these enzymes are all metalloenzymes, the catalytic metal center plays a major role in substrate positioning in all cases, guiding the ultimate orientation of the electrophile relative to the nucleophile. This is further facilitated by the involvement of a number of key residues that are strategically positioned to assist in overall substrate positioning, which for the PMHs studied here are N78, H218, and K337. These residues primarily interact with the leaving group or nonbridging oxygens of the relevant substrates and keep the electrophile in similar positions relative to the nucleophile for different substrates and electrophiles. However, due to the large binding pocket (≈10 × 20 Å2 wide and 15 Å deep),[12] there is extensive space to accommodate variations in binding conformations of spectator and leaving groups, which in turn would facilitate the accommodation of a broader range of substrates shapes. Such substrate repositioning through diversity in placement of leaving and spectator groups would therefore also play a role in determining the resulting overall catalytic efficiency and promiscuity. This is in line with our computational evidence, which highlights the importance of cooperative electrostatic interactions, brought about by active site plasticity, in accommodating a range of different substrates.

Tying in with this, if plasticity is important for facilitating promiscuity in these enzymes, one can ask whether flexibility would be reduced for catalysis by enzymes that show high specificity for their physiological substrates. A classical set of proteins where rigidity is very important in ligand binding specificity vs promiscuity are the periplasmic binding proteins, as illustrated by the structure of the cellulose-binding protein from Thermotoga maritima.[81] This protein has a bipartite active site, comprised of a solvent excluded region involved in highly specific ligand binding and is adjacent to a second and more flexible solvent-filled cavity in which semi-specific ligand binding occurs. Detailed studies of these systems as well as exploration of the role of water molecules have provided important insight into how the interplay between flexibility and rigidity allows both specificity and promiscuity to be encoded into a single binding site, moving beyond a single highly specific and fixed protein scaffold.

Other examples of highly specific enzymes that show comparably “rigid” active sites (in the substrate-bound conformation) are orotidine 5′-monophosphate decarboxylase[82] and β-phosphoglucomutase,[83] among others. These enzymes can exist in more than one conformational form and undergo ligand-gated conformational changes, engulfing their substrates upon ligand binding. However, once the substrate is bound, crucial tight binding hydrogen-bonding networks are involved in keeping the key catalytic residues in place, and, for example, truncation of various functional groups on the respective substrates can lead to tremendous reductions in catalytic activity due to the loss of key stabilizing interactions (see, for example, Richards’ “substrate-in-pieces” studies of these systems that quantify the contribution of different parts of the substrate to binding and catalysis). Yet another example is a recent study of the evolution of β-lactamases from their promiscuous ancestral variants to their modern specific counterparts (such as TEM-1 β-lactamase).[84] This work demonstrates that, within these enzymes, evolution from a generalist to a specialist enzyme is coupled to a loss of conformational flexibility.

Finally, as pointed out by a reviewer, it should be noted out that the experimental work on which the present study is based[8,12] was performed using nonphysiological substrates, where the leaving group is the weakly basic nitrophenoxide anion (as the physiological substrate has not been identified).[12] The relatively small requirement for stabilizing negative charge at weakly basic anions will in turn increase the contribution of the enhanced reactivity of the active site nucleophile to the total rate accelerations outlined in Table S1, facilitating the turnover of a broader range of substrates.

Implications for Evolution in the Alkaline Phosphatase Superfamily

As discussed in the Introduction, the phenomenon of catalytic promiscuity appears to be common among members of the AP superfamily.[15] Although the members of this family (which include AP, NPP, PMH, and PAS) have diverged considerably, they still share considerable similarities in active site architecture and key catalytic residues, as highlighted in Figure 1 and discussed in ref (15). In the present work, we have performed a detailed computational study of the hydrolysis of a range of substrates by two PMH, demonstrating the key role of cooperative electrostatic interactions at the PMH active site. Moreover, there is no significant conformational change in the active site alongside the reaction coordinate when comparing the different reactions studied. This suggests an important role for electrostatic rather than conformational[86] active site plasticity in facilitating the observed promiscuity in these enzymes. In order to examine this effect, we have studied the geometry of the active site cavity for different members of the AP superfamily, using the Fpocket 2 software package.[87] We have inspected several members of this superfamily, as well as related promiscuous phosphatases, in the search for a possible correlation between the physical properties of the different active sites and their corresponding catalytic promiscuity. The corresponding results are summarized in Figure 9 and Tables S7 and S8. From this data, it can be seen that when comparing different AP superfamily members, clear trends emerge with respect to the active site volumes and the number of activities that have been reported for each enzyme to date.[7−10,12,21,28,71,80,88,89] Specifically, according to this analysis, pockets with a larger polar surface allow the enzyme to exploit distinct residue conformations to create an optimal electrostatic environment and accommodate different transition states. In addition to this, the large volume of the different pockets would allow for the accommodation of a more diverse range of substrates, which can then be hydrolyzed through cooperative enzyme–substrate electrostatic interactions in the corresponding active sites.

Figure 9

Correlation between the number of known catalytically activities (in parentheses) and the total volume and SASA for several members of the AP superfamily.[85]

Through this comparison, we find that AP, NPP, and the PMHs have the largest active site volumes and polar solvent accessible surface areas (SASAs), mostly due to the width of their active sites. Tying in with this, the arylsulfatases and other members of the superfamily have comparably smaller and narrower pockets. Interestingly, this can be directly correlated to the number of reported activities in the literature (Figure 9 and Table S7) where AP has both the largest and most accessible active site as well as the highest number of reported activities.[7,10,80,88,89] This is closely followed by NPP and the PMHs, with five clear activities each (not including the anomalous PTE activity of BcPMH for reasons outlined above). Finally, a BLAST search with RlPMH against known protein structures yielded the related alkaline phosphatase, PAS, as the protein with the highest sequence identity.[8] As can be seen from Figure 9, upon moving from the PMHs to PAS, the active site starts to reduce in volume to give a much narrower pocket, in line with the lower number of reported activities for PAS.[9,28] The remaining enzymes all have smaller active site volumes compared to AP, NPP, and PMH, and all of them, according to Braunschweig Enzyme Database,[90] have so far been reported to have only one activity each. When this observation is combined with the similar mechanisms and highly polar active site residues that members of this superfamily possess (Figure 10) as well as with the experimental evidence for the existence of catalytic backups in PON1,[68] this strongly suggests that the cooperative electrostatic flexibility observed in PMH and related enzymes is a common feature for evolution in the AP superfamily as well as related phosphotransferases.

Figure 10

Surface representation of the active sites of (A) Burkholderia caryophylli PMH, (B) Escherichia coli AP, (C) PAS, and (D) Xanthomonas axonopodis NPP (PDB IDs: 2W8S, 1ALK, 1HDH, and 2GSN, respectively), displaying the polar character inside the pocket. (B) and (D) show a strong presence of negatively charged residues (red) near the metal site. For both (A) and (C), there are more apolar residues present (white), although BcPMH still has a larger number of polar residues in its active site, mostly due to the very large size of the binding pocket.

One limitation of this analysis, however, is that while we consider the total number of characterized activities, it neither takes into account the possibility of further as-yet uncharacterized activities in these enzymes nor the relative proficiency of these enzymes toward their promiscuous substrates. For example, even though AP has the highest number of known activities among the enzymes we examine, only two of those six activities (phosphate and phosphothioate monoester hydrolysis) are particularly proficient with kcat/KM values of 3.3 × 107 and 2.0 × 104 M–1·s–1, while the other activities can have kcat/KM values as low as 10–3 M–1·s–1.[8] However, clearly, a very high number of polar residues in the active site, as shown in Figure 10 and Table S7, as well as very large active sites, would allow for the presence of multiple distinct catalytic backups and a shifting electrostatic field upon substrate binding that can accommodate substrates of other shapes and charge distributions. Additionally, it is of course useful to consider not only the total binding interactions available for transition-state stabilization but also the binding interactions that are actually needed to account for a given observed enzymatic rate acceleration, because nonspecificity will be favored whenever these total possible interactions greatly exceed the number of required interactions. For example, in the case of phosphate monoester hydrolysis, strong interactions between the enzyme and heavily charged reacting phosphate may be more than sufficient to account for the enzymatic rate acceleration. This would favor a lack of specificity for the leaving group and, perhaps, also floppiness in transition-state binding, provided there is no strict requirement of the precise placement of enzymatic side chains around the phosphate. In contrast, there may be a greater requirement for precision in the binding of the less highly charged transition state for sulfate monoester hydrolysis, which would also tie in with the experimental observation that PAS, a native sulfatase, is a more proficient phosphatase than the corresponding phosphatases in the AP superfamily are sulfatases.[15]

This would be supported by our observed correlation between larger active site volume/polar SASA and a great number of activities and is in sharp contrast to enzymes such as orotidine 5′-monophosphate decarboxylase, which is highly selective for the decarboxylation of orotidine monophosphate, because the enzyme makes use of every possible interaction with OMP in the stabilization of the decarboxylation state.[91] Interestingly, there appears to be also some sort of correlation between tertiary structure and the number of reported activities, in that the enzymes with the highest number of characterized activities shown in Figure 9, namely AP, NPP, and PMH, are dimers and a tetramer, respectively, whereas all other enzymes are monomeric. One would assume that having a large number of polar residues in the active site would result in electrostatic strain that would have to be compensated elsewhere in the structure, which could potentially be correlated to the oligomeric states of these proteins. Overall, however, the clear correlation between increased promiscuity and a larger active site volume and SASA highlights the crucial importance of substrate charge and active site electrostatics in facilitated selectivity and evolution among these highly promiscuous enzymes.

Overview and Conclusions

In the present work, we have performed a detailed EVB study of both the native and several promiscuous activities of two PMH, BcPMH and RlPMH, as well as RlPMH variants with mutations in key active site residues. Our calculations can reproduce key experimental observables such as experimentally observed activation barriers for the wild-type reactions of all substrates and qualitative mechanistic predictions based on examining pH-rate profiles as well as energetic trends upon mutation of key active-site residues in RlPMH. We demonstrate that despite their broad promiscuity, both PMHs studied in this work hydrolyze all five chemically distinct substrates through a unified mechanism, binding substrates in similar positions and without the need for any significant local or global conformational changes. Additionally, we demonstrate that the apparent resilience of these enzymes to active site mutations as well as the overall promiscuity is due to compensatory electrostatic effects from different residues, allowing enough flexibility in the electrostatic environment of the active site to accommodate multiple substrates with distinct transition states and charge distributions. Finally, we provide a detailed structural and physical comparison of a range of highly promiscuous members of the AP superfamily.

These results demonstrate the strong correlation between the structural and electrostatic features of these enzyme’s active sites and the corresponding variations in both substrate charge preference and the number of known promiscuous activities. This further supports our hypothesis by strongly suggesting that active site shape, size, and more critically number of polar residues available can be directly correlated the ability to accommodate increasing numbers of promiscuous activities. Our simulations and comparative analysis therefore highlight the importance of cooperative electrostatic interactions and an electrostatically flexible active site as a common feature in the evolution of promiscuous side reactions among members of the AP superfamily. In the present work we demonstrate that, in addition to the electrostatic preorganization originally suggested by Warshel in 1978,[35] the active site can also electrostatically reorganize to accommodate the needs of different substrates. This provides a classical example of protein flexibility allowing the reaction to occur, as these enzymes do not know in advance what substrate is going to bind. Rather, they adjust their active site environment to a given substrate after the binding step. These insights, in turn, helps us not only to understand protein evolution within a superfamily at the molecular level but also highlights a concrete feature that can be manipulated in targeted artificial enzyme design.