Analysis of electrostatic coupling throughout the laboratory evolution of a designed retroaldolase.

The roles of local interactions in the laboratory evolution of a highly active, computationally designed retroaldolase (RA) are examined. Partial Order Optimum Likelihood (POOL) is used to identify catalytically important amino acid interactions in several RA95 enzyme variants. The series RA95.5, RA95.5-5, RA95.5-8, and RA95.5-8F, representing progress along an evolutionary trajectory with increasing activity, is examined. Computed measures of coupling between charged states of residues show that, as evolution proceeds and higher activities are achieved, electrostatic coupling between the biochemically active amino acids and other residues is increased. In silico residue scanning suggests multiple coupling partners for the catalytic lysine K83. The effects of two predicted partners, Y51 and E85, are tested using site-directed mutagenesis and kinetic analysis of the variants Y51F and E85Q. The Y51F variants show decreases in kcat relative to wild type, with the greatest losses observed for the more evolved constructs; they also exhibit significant decreases in kcat /KM across the series. Only modest decreases in kcat /KM are observed for the E85Q variants with little effect on kcat . Computed metrics of the degree of coupling between protonation states rise significantly as evolution proceeds and catalytic turnover rate increases. Specifically, the charge state of the catalytic lysine K83 becomes more strongly coupled to those of other amino acids as the enzyme evolves to a better catalyst.

partial order optimum likelihood

designed retroaldolase, series 95

theoretical microscopic titration curve shapes

nth central moment

INTRODUCTION

Enzymes, with their exquisite selectivity and catalytic power, have tremendous potential for use in the chemical industry and for diverse applications, such as bioremediation and biofuel production. Indeed enzymes have gained some use in industrial chemical processes; for example nitrile hydratase, which catalyzes the conversion of nitriles to amides, is used in the production of acrylamide on the kiloton scale yearly, at lower cost than conventionally catalyzed processes. , However, for many reactions involved in industrial chemical processes, a natural enzyme may not exist. The ability to design enzymes that can perform specific, unnatural reactions can open the door to greener chemical industry, with less energy consumption and fewer unwanted by‐products than thermal processes or processes that use more conventional catalysts. This requires expanded knowledge of the fundamental mechanisms by which both natural and designed enzymes work.

Very recently, Mazmanian, Sargsyan and Lim reviewed how the environment of enzyme active sites, namely local interactions, the presence of solvent molecules, and conformational flexibility, control enzymatic function. Here we focus on the role of local interactions , between catalytic residues and nearby residues on the activity of a series of designed retroaldolases of the RA95 family, which has been optimized by laboratory evolution. , , ,

Design of novel enzymes typically involves a new, computationally designed active site built into a natural protein scaffold. While this can be significantly more difficult than engineering natural enzymes for altered specificities, it opens up the functional space of enzymes to catalyze reactions that do not have a natural counterpart. Unfortunately, low activity of initial designs plagues most engineering attempts.

To increase the activity of the designed enzymes, directed evolution is employed; it relies on the screening of many enzyme variants, often with mutations introduced randomly across the sequence. Beneficial mutations are slowly accumulated and combined to produce better enzymes and the process can require many iterations. While this strategy can produce enzymes with activities orders of magnitude greater than the starting design, it is time‐ and resource‐intensive. To date such experiments have taken years to produce enzymes with useful levels of activity. Advances have been made to help directed evolution, such as the use of microfluidics to miniaturize the process and increase throughput, , , but improved initial designs could reduce substantially the time required to achieve effective activity.

A previously designed retroaldolase, RA95, has been the focus of much research as it has been evolved over several years and its activity increased over five orders of magnitude from its initial design. , , Kinetics data have been reported for the designs as they evolve from low activity to high activity, including data for multiple steps in between. Designed enzymes from this series were chosen for the present study because measured rate constants span a wide range of values and because crystal structures are available for multiple steps along the evolutionary path. Our computational tools allow us to understand better the sources of increased activity in the RA95 series and ultimately to help to establish the basic concepts for increased activity of designed enzymes. We have identified computed features that correlate with the increase in activity and have been able to attribute activity gains to specific residues in some instances. Using our computational tools such as THEMATICS , and POOL, , we have identified key mutations during evolution that contribute to increased activity along the evolutionary trajectory. This has also enabled the discovery of key metrics that mirror the activity differences between the enzymes, suggesting these values might be used to guide enzyme designs or to screen candidates computationally. These metrics are based on computed electrostatic and chemical properties of the individual amino acid residues in the protein structure.

RESULTS

Predicting chemical properties of protein active sites with POOL

Predictions by POOL , of biochemically active residues are based on metrics obtained from computed proton occupation functions C(pH) for each site with proton transfer capability, which are in turn based on the computed electrical potential function for the protein structure. For a residue that follows Henderson‐Hasselbalch behavior, the function C(pH), which always equals 1 at low pH and 0 at high pH, is a sigmoidal function with a sharp fall‐off at pH values near the pK . However, for residues active in catalysis and/or ligand recognition, C has an anomalous shape, arising from strong coupling to proton transfer equilibria on nearby residues, much like the behavior observed for polyprotic acids. To quantify the deviation from typical behavior, the C(pH) functions are expressed in terms of a first derivative function f, defined as:

The negative first derivative of C, f, is essentially a proton binding capacity, a measure of the change in concentration of a bound species per unit change in its chemical potential. For a typical titratable residue, the f function may be well approximated by a single Gaussian distribution function, with a peak at or near the pH equal to the pK . In contrast, for titratable residues involved in catalysis or ligand binding, the f function deviates significantly from Gaussian form and may be represented as a superposition of Gaussian functions. The area under the f function is always unity because of the 1 to 0 range of the C functions and f is nonnegative; therefore f(pH) resembles a probability density function, the shape of which can be characterized by summary statistics, such as the nth central moments μn, defined as: where m1 is the first moment, defined for n = 1 by the expression for the nth moment mn as: where the integrals in Equations (2) and (3) are over all real numbers (−∞ to +∞). We have established that the third and fourth central moments, μ3 and μ4 of Equation (2), are excellent predictors of the residues active in catalysis and/or ligand recognition for natural enzymes. , The third and fourth central moments may be calculated for each titratable residue, using the electrical potential function calculated from the three‐dimensional structure of the protein. These features are indicators of the degree of coupling between the proton binding reaction of a given titratable residue and those of other residues; for biochemically active residues in natural enzymes, these couplings are strong. , , , ,

There is considerable evidence that electrostatic effects, including contributions from distal residues, are important in enzyme catalysis. A 1991 density functional theory study was reported by Bajorath et al. of free dihydrofolate (DHF) and of DHF bound to dihydrofolate reductase. They showed that when DHF is placed in the electrical potential of the enzyme, the π orbitals of the reactive C=N bond of DHF are polarized, such that the π bond begins to split into two lobes that resemble atomic 2p orbitals on the C and N atoms. This facilitates reduction wherein σ bonds are formed with hydrogen 1 s orbitals; the polarizing potential was reported to arise significantly from both nearby and remote residues. Electric fields in protein binding sites can be substantial. Boxer et al. used vibrational Stark spectroscopy to measure local electric fields in protein binding sites and reported electric fields on the order of 100 MV/cm!

Experimental evidence of significant coupling by distal residues has been reported by us. , , , , , For all of these cases, the important couplings of the specific distal residues were predicted computationally. Although there is clear evidence of the importance of electrostatic effects and coupling of distal residues in enzyme catalysis, these effects are often not considered, including in the context of enzyme design.

Directed evolution often utilizes random mutagenesis to identify variants with a new function or higher activity, as was done in the RA95 family. It is often however not known exactly what causes the increased activity, specifically what residues are responsible or how they contribute to increased activity. The answers are generally not obvious. For instance, the introduction of a second lysine (K83) in the evolved variants of the RA95 series was originally believed to function by helping to modulate the pK a of the originally designed nucleophilic lysine (K210) as its function. However, x‐ray crystallography showed that K83 had become an active nucleophile itself. Understanding how and which mutations are beneficial is important to the development of guiding principles for enzyme engineering.

POOL PREDICTIONS AND DIRECTED EVOLUTION RETROSPECTIVE

For three constructs, RA95.5, RA95.5–5, and RA95.5–8, five residues have POOL scores higher than the average score; for RA95.5–8F, seven residues fall above this threshold. The amino acids ranked in the top 12 are shown in Table 1; residues with POOL scores below the average score are shown in italics. For RA95.5, RA95.5–5, and RA95.5–8 both active site lysines, K83 and K210, are predicted to be important for biochemical activity. K210 was the original catalytic residue introduced by design. Over the course of evolution, K83 was introduced and activity gradually shifted from K210 to K83. Thus, the original design only had K210; RA95.5, RA95.5–5, and RA95.5–8 have both. In RA95.5–8F, evolution replaced K210 with leucine.

TABLE 1

The top 12 POOL‐predicted biochemically active residues in rank order for members of the RA95 series. Residues in italics have POOL scores below the average score for each construct

POOL rank	RA95.5	RA95.5–5	RA95.5–8	RA95.5–8F
PDB	4A2S	4A2R	Model	6TFA
1st	K83	K83	K83	Y51
2nd	K210	Y51	Y51	K83
3rd	Y51	K210	K210	Y180
4th	D111	D111	D111	D212
5th	E85	E85	E85	D111
6th	S233	N110	R54	D61
7th	R19	L231	Y52	E85
8th	S110	R19	N110	E135
9th	F112	K135	S233	L210
10th	D165	L131	L131	R19
11th	L131	F112	F89	N110
12th	I232	Y52	R19	S233

The catalytic mechanism involves attack of the reactive lysine (K210 and/or K83) on the carbonyl group of the substrate to form a tetrahedral intermediate, which loses water to give a Schiff base adduct (Figure 1). The tyrosine residue Y51 is in contact with K83, located slightly behind K83 with respect to the substrate in RA95.5 but moves closer to the substrate in later constructs. It is ranked highly by POOL for all four constructs. Y51 has been proposed as a likely Brønsted acid/base, although a computational study has suggested that Y180 may assume this role in the highly evolved RA95.5–8F. Note that Y180, which emerged late in the evolutionary trajectory and is part of a hydrogen‐bonded network consisting of Y51, K83, and N110, is ranked third in RA95.5–8F; the earlier constructs all have phenylalanine in this position.

FIGURE 1

Reaction scheme for the retroaldolase

Located near the entrance to the binding pocket, along the rim of the beta barrel, E85 is predicted to be active in all four variants; no role in catalysis has been assigned previously to this residue. The acid side chain of E85 forms hydrogen‐bonds with the backbone N‐H groups of K87 and Y88, contributing to fold stability. D111, which is situated behind E85, is also predicted for all four constructs. D111 is conserved throughout the series and is critical to structural integrity; it forms a salt bridge with the conserved R19, helping to position the N‐terminal helix–loop, and has a charge‐dipole interaction with the alcohol group on the side chain of T84. Two additional aspartate residues, D212 and D61, are predicted for RA95.5–8F. D61 and D212 are located on the rim of the beta barrel at the entrance to the binding pocket, on the side opposite E85 and D111; these acidic residues at the entrance to the binding pocket may help to orient the substrate as it enters the pocket.

Residue N110, which provided a five‐fold boost in activity when first introduced by directed evolution, ranks in the top 12 for RA95.5–5, RA95.5–8, and RA95.5–8F, but has POOL scores below the average score, a somewhat arbitrary cutoff. Note that the original serine at this position ranks eighth in the early construct RA95.5. N110 is located within hydrogen‐bonding distance of K83 and its polar side chain may contribute to the catalytic power of this catalytic lysine.

Y51 and E85 were chosen for site‐directed mutagenesis studies here because they are consistently among the top POOL‐predicted residues across the series. The conservative mutation Y51F maintains the steric bulk of the aromatic side chain but loses the ionizable phenol group. The amide side chain in an E85Q variant can maintain the hydrogen‐bonding interactions with the N‐H groups of K87 and Y88. D111 is also consistently highly ranked but its structural importance is presumed to be too great to separate from a biochemical role in catalysis; it forms a key salt bridge with R19 and a charge‐dipole interaction with T84. The relative positions of residues Y51, K83, E85, and D111 in RA95.5–8F are shown in Figure 2.

FIGURE 2

Residues Y51, K83, E85, D111, and Y180 in the structure of RA95.5–8F, based on PDB 5AN7. The ligand is shown in orange. Image rendered in YASARA

Increasing values for the fourth central moment μ4 [see Equation (2) above] for the active site residues are observed as the designed enzymes evolve and their activity increases. The μ4 value is an indicator of the strength of electrostatic coupling between the subject residue and other residues. These calculated μ4 values are shown in Table 2 for K83 and for the residues most strongly coupled to it, together with the k cat values for each construct, reported previously by Obexer et al. and measured for the current study. More couplings between the catalytically active residues and nearby residues are built in as the enzyme evolves and thus the μ4 values rise. For comparison, Table 2 shows the μ4 values for the highly ranked residues E85 and D111. These values increase in the step from RA95.5–8 to RA95.5–8F as more acidic residues are coupled to them.

TABLE 2

Turnover , rates and μ4 values for K83 and strongly coupled residues Y51 and Y180 for wild‐type variants in the RA95 series; μ4 values are also shown for the highly ranked acidic residues E85 and D111

Variant	k cat (s−1) a	K83	Y51	Y180	E85	D111
RA95.5	0.00286 ± 0.00016	33 b	48 b	–	18	31
RA95.5‐5	0.172 ± 0.008	37	55	–	19	36
RA95.5–8	0.383 ± 0.024	41	60	–	20	35
RA95.5–8F	11.9 ± 0.60	104	135	70	32	47

These values were measured for the current study and are in agreement with the previously published k cat values. ,

Values averaged for A and B rotamers of K83.

The combination of residue scanning, where all possible single‐site variants are made in silico, and THEMATICS , , shows how different mutations affect the electrostatic features of nearby residues. This provides information to help to identify the new residues or mutations that are responsible for changes in the central moments of K83. Residue scanning results for multiple evolved RA95 family members are displayed in Figure 3. Contributions of specific ionizable residues to the catalytic power of K83 are evident when replacement of a given ionizable residue with any non‐ionizable residue type results in a significant change in the calculated μ4 value of the catalytic K83. Note the substantial reduction in the computed μ4 value of K83 for all four variants when the tyrosine at position 51 is mutated in silico to a non‐ionizable residue. (Sequence position 51 is marked with a red triangle in Figure 3.) The importance of the newly evolved Y180 in RA95.5–8F is also evident, with substantial decreases in μ4 calculated when this tyrosine is replaced with a non‐ionizable residue; sequence position 180 is marked with a purple triangle. Residues R19, Y69, and D165, all conserved through the series, plus a glutamate E135 newly introduced in RA95.5–8F, are all also coupled to K83. The overall trend of increasing μ4 for Lys83 with increasing activity is also apparent in Figure 3.

FIGURE 3

Mutational effects on the μ4 of Lys83 in the evolved variants of RA95. The position of the mutation in the sequence is shown on the x‐axis and the calculated μ4 for Lys83 is on the y‐axis. Mutations showing a change from an ionizable residue to a non‐ionizable residue are shown. The different colors represent the different structures or variants used: RA95.5 A conformation (dark blue); RA95.5 B conformation (yellow); RA95.5–5 (light blue); RA95.5–8 (orange); and RA95.5–8F (green). Key interacting residues can be identified by the changes they cause the catalytic Lys83. Triangles show sequence positions 51 (red), 85 (yellow), 111 (green), 180 (purple), and 210 (blue)

EXPERIMENTAL TESTING OF PREDICTED RESIDUES Y51 AND E85

Previously, Y51 was shown to be important in RA95.5–8F and is predicted to be important in all members of the series studied here. It is believed to function as a proton shuttle in the catalytic mechanism. Y51 and E85 were chosen for experimental investigation here because both scored consistently highly in the POOL rankings throughout the series. To test the effects of the predicted residues Y51 and E85 on the members of the RA95 series, the catalytic activity was measured for variants with mutations of the POOL‐predicted residues Y51 and E85. Conservative mutations, that is, Y51F and E85Q, were made. Results of kinetics measurements for the wild type constructs are shown in Table 3, for the Y51F variants in Table 4, and for E85Q variants in Table 5. Results for the wild type RA95.5–8F and its Y51F variant are consistent with values reported previously. Michaelis–Menten plots are shown in Supplementary Material.

TABLE 3

Kinetics constants for wild type members of the RA95 series

	WT
	k cat (s−1)	K M (μM)	k cat/K M (M−1 s−1)
RA95.5	0.00286 ± 0.00016	279 ± 35	10.3 ± 1.9
RA95.5–5	0.172 ± 0.008	295 ± 32	(5.83 ± 0.90) 102
RA95.5–8	0.383 ± 0.024	275 ± 35	(1.39 ± 0.26) 103
RA95.5–8F	11.9 ± 0.60	308 ± 34	(3.86 ± 0.62) 104

The loss of Y51 in the Y51F variants results in significant reduction in the μ4 values for K83 for all members of the series. Table 6 gives the μ4 values of key residues for homology model structures of the Y51F variants. The E85Q variants have a small effect on the μ4 values for K83. A graph of log10(k cat) as a function of μ4 of K83 for the wild type, Y51F, and E85Q variants is shown in Figure S12 of Supplementary Material.

TABLE 6

Calculated μ4 values for key residues in homology models of the Y51F variants

Variant	K83	Y51	Y180	E85	D111	K210
RA95.5 Y51F	17	–	–	22	41	16
RA95.5–5 Y51F	19	–	–	15	24	14
RA95.5–8 Y51F	21	–	–	20	36	19
RA95.5–8F Y51F	36	–	40	33	52	–

Thermal shift assays for the wild type and Y51F and E85Q variants are shown in Supplementary Material. In all four constructs, the Y51F variant shows a small reduction in melting temperature of 1–3°C. The effect is larger for the E85Q variants, with reductions in melting temperature of 5–7°C. Melting temperatures are in the range 70–79°C for Y51F and 65–74°C for E85Q.

DISCUSSION

Evolution builds stronger coupling between protonation states of active residues

The rising μ4 values of key residues shown in Table 2 indicate increasing coupling of the protonation equilibria of biochemically active residues with other ionizable residues along the evolutionary trajectory. Strong coupling between protonation states gives rise to expanded buffer ranges, , , , as in polyprotic acids, enabling the active residues to have significant population of both the protonated and deprotonated states across a wide pH range. Therefore, the side chains can be in the correct protonation state for activity and residues that serve as Brønsted‐Lowry acids and bases are able to regenerate for the next turnover cycle. ,

Y51

Table 4 shows that Y51 plays an important role in catalysis for all four constructs, as significant activity is lost upon mutation to phenylalanine. In RA95.5–5 it is located behind K83 with respect to the substrate, with the OH atom of Y51 situated 4.5 Å away from the NZ atom of K83. In RA95.5–8F, it is in direct contact with the ligand and the OH‐NZ distance is shorter in the crystal structure (3.2 and 4.0 Å for the two rotameric states). Y51 has been identified by molecular dynamics simulation as a residue that contributes to shifting the conformational equilibria toward catalytically competent conformations ; for some states along the dynamic trajectory, partial hydrogen bonding may be possible between K83 and Y51 in the earlier constructs. Even in the absence of direct hydrogen bonding, in the highly evolved RA95.5–8F, K83 and Y51 are connected via a hydrogen bonding network with N110 and Y180. For RA95.5–8F, it is reported to be a part of a catalytic tetrad. In the proposed mechanism, K83 forms a protonated Schiff base, attacking the carbonyl moiety of the substrate to form a carbinolamine intermediate; Y51 is in position to deprotonate the hydroxyl group of the substrate in the Schiff base intermediate. In addition, the high, similar intrinsic pK as of the lysine and tyrosine side chains, with the spatial proximity of Y51 and K83, result in strong coupling of the proton transfer equilibria of the two residues, giving the side chains of both residues a wide buffer range. The wide buffer range enables the primary amine group of K83 to have significant population of its deprotonated state at neutral pH, a feature necessary for attack on the substrate carbonyl group. Calculated energies of interaction between Y51 and K83 are given in Table S1 of Supplementary Material.

TABLE 4

Kinetics constants for Y51F variant members of the RA95 series, with fold change in k cat and in k cat/K M compared to wild type

	Y51F
	k cat(s−1)	K M(μM)	k cat/K M (M−1 s−1)	Fold change k cat	Fold change k cat/K M
RA95.5	0.00154 ± 0.00028	1700 ± 422	9.06 ∙ 10−7 ± 3.9 ∙ 10−7	1.9	11
RA95.5–5	0.00379 ± 0.00017	301 ± 30	1.26 ∙ 10−5 ± 0.18 ∙ 10−5	45	46
RA95.5–8	0.00522 ± 0.00011	133 ± 8.5	3.92 ∙ 10−5 ± 0.19 ∙ 10−5	73	35
RA95.5–8F	0.162 ± 0.0021	36.7 ± 2.1	4.41 ∙ 10−3 ± 0.31 ∙ 10−3	73	8.8

The Y51F mutation leads to reductions of 1.9‐, 45‐, 73‐, and 73‐ fold in k cat compared to wild type for RA95.5, RA95.5–5, RA 95.5–8, and RA95.5–8F respectively, with 11‐, 46‐, 35‐, and 8.8‐ fold reductions, respectively, in k cat/K M. The level effect on k cat between RA95.5–8 and RA95.5–8F could be because of a change in rate‐determining step along the evolutionary path; the rate‐determining step is reported to be early in the mechanism (the first step in Figure 1) for the early constructs, but shifts to the protonation of the enamine intermediate just prior to product release for the more highly evolved constructs.

The introduction of Y180 in proximity to K83 (about 4.4 Å away) in RA95.5–8F, resulting from an F180Y substitution on evolving from RA95.5–8 to RA95.5–8F, adds redundancy to the role of Y51. This could be one contributing factor to the observation that the reduction in k cat is level between RA95.5–8 and RA95.5–8F. Indeed Y180 has been identified as possibly serving as the catalytic base in the highly evolved RA95.5–8F. , A previous site‐directed mutagenesis study (Table S1) showed that mutation of either of these two tyrosine residues individually in RA95.5–8F had a relatively modest effect on k cat/K M, while simultaneous mutation of both tyrosine residues resulted in substantially reduced activity (200‐fold loss in catalytic efficiency for the double mutant Y51F/Y180F). Note that Y51, K83, and Y180 are the three top residues in the POOL rankings for RA95.5–8F.

E85

E85 is located at the entrance to the active site pocket, on the rim of the beta barrel. Tables 5, 6 show that the conservative mutation E85Q leads to small but significant reduction in catalytic efficiency for RA95.5, RA95.5–5 and RA95.5–8 but has little effect on k cat. It plays some role in thermal stability, as shown in the thermal shift data (Supplementary Material); in RA85.5–8F its side chain hydrogen‐bonds to the backbone amide groups of K87 and Y88, helping to stabilize the 85–91 loop. It is one of the acidic residues, D61, E85, D111, E188, located along the base of the beta barrel and conserved across the RA95.5 – RA95.5–8F series. These residues are most likely involved in orienting the substrate in the active site pocket. Additional acidic residues in this location along the rim of the beta barrel appear in the evolutionary trajectory, D212 in RA95.5–8 and E135 in RA95.5–8F. The presence of additional acidic residues adds redundancy and diminishes the role of E85.

TABLE 5

Kinetics constants for E85Q variant members of the RA95 series, with fold change in k cat and in k cat/K M compared to wild type

	E85Q
	k cat(s−1)	K M(μM)	k cat/K M (M−1 s−1)	Fold change k cat	Fold change k cat/K M
RA95.5	0.00281 ± 0.00021	1,320 ± 144	2.13 ∙ 10−6 ± 0.39 ∙ 10−6	1.0	4.8
RA95.5–5	0.0555 ± 0.0034	578 ± 69	9.60 ∙ 10−5 ± 1.73 ∙ 10−5	3.1	6.1
RA95.5–8	0.164 ± 0.0071	509 ± 44	3.22 ∙ 10−4 ± 0.42 ∙ 10−4	2.3	4.3
RA95.5–8F	10.75 ± 0.50	493 ± 45	2.18 ∙ 10−2 ± 0.30 ∙ 10−2	1.1	1.8

CONCLUSIONS

While POOL was created and optimized to predict the biochemically important residues for natural enzymes, it shows predictive capability for this designed enzyme. Note that K83 ranks first in Table 1 for the early construct RA95.5 and ranks first or second for the later constructs. Y51 ranks among the top three for all four constructs. K210 ranks highly for all constructs where it is present. The top three residues for the highly evolved RA95.5–8F are Y51, K83, and Y180, all shown to contribute to catalysis. N110, previously shown to be important in catalysis, ranks among the top 12 for RA95.5–5, RA95.5–8 and RA95.5–8F; S110 also ranks eighth in RA95.5.

The calculated values for μ4 of the catalytic K83, a measure of the degree of coupling between its protonation equilibrium and those of other residues, are 33, 37, 41, and 104, for wild type RA95.5, RA95.5–5, RA95.5–8, and RA95.5–8F, respectively. For the original construct RA95.0, with a measured k cat of 1.0 • 10−4 s−1, its catalytic K210 has a μ4 of only 15. For contrast, the active site lysine residues of the natural enzyme human fructose 1,6‐bisphosphate aldolase have μ4 values of 54, 272, and 190 for K107, K146, and K229, respectively. The larger scaffold of this natural enzyme provides more opportunities for favorable coupling partners to be built in by natural evolution.

Strong coupling between the catalytic K83 and nearby residues, including Y51 and Y180, and also weaker coupling with more distant residues (5–10 Å away), including R19, Y69, and D135, improves the catalytic power of K83. Thus K83 in the highly evolved enzyme more closely resembles the catalytic lysine residues of natural enzymes than its counterpart in the earlier designs. These couplings expand the buffer range of the side chain and enable the deprotonated state to have significant population at the pH at which the enzyme functions. These couplings are necessary features to build into designed enzymes to produce activity that rivals natural enzymes. Mutations to introduce such couplings can be determined computationally; the difficulty is finding the mutations that do not diminish other vital properties.

MATERIALS AND METHODS

Computational

Input structures for POOL predictions were obtained as follows: The RA95.5 (PDB ID 4A2S) and RA95.5–5 (PDB ID 4A2R) structures were downloaded from the PDB. A homology model was built for RA95.5–8, as no structures have been reported. The comparative model was built using the homology model feature in the YASARA suite of programs. A total of five models were built using five different templates. The highest‐scoring model came from the RA95.5 structure (PDB ID 4A2S). To evaluate the model, a Ramachandran diagram was drawn using the MolProbity server; the model has zero Ramachandran outliers and 98.4% of residues are in favored regions. The model passed the 3D‐1D profile assessment by VERIFY3D, with 94.4% of the residues at or above the threshold value of 0.2 (>80% is required for the structure to pass). For RA95.5–8F, the structure (PDB ID 5AN7) was downloaded from the PDB; the missing loop for residues 58–61 was built in using the homology model module in YASARA. Structures for all four Y51F variants were built using the homology model module in YASARA.

POOL calculations were performed using electrostatic properties from THEMATICS , , and surface topology properties from the structure‐only version of ConCavity as input features, using the method of Somarowthu. The average POOL score for all residues in a given structure was used to define the cutoff, similar to the method of Brodkin used for natural enzymes; residues with POOL scores higher than the average score are designated, somewhat arbitrarily, as predicted residues for present purposes.

The program FoldEx , was used for residue scanning (Figure 3) to generate all possible single‐site variants in silico. POOL was performed on the in silico variants as above.

Materials

Chemicals, solvents, buffer components, media, and antibiotics were purchased from commercial sources (Sigma‐Aldrich, Merck, Acros, Fluka) and used without further purification. Phusion DNA polymerase and the DpnI restriction enzyme were purchased from New England Biolabs (USA). A plasmid purification kit from ZYMO Research (USA) and Ni‐NTA agarose from Qiagen (Germany) were used.

Site‐directed mutagenesis

pET29b plasmids that contain the genes coding for RA95.5, RA95.5–5, RA95.5–8 and RA95.5–8F with a C‐terminal His6‐tag were available from previous work. , , The mutations coding for Y51F or E85Q were introduced in all of the templates according to the Quikchange® standard protocol (Agilent). DNA oligo synthesis and DNA sequencing were performed by Microsynth AG (Switzerland).

Protein production and purification

The 12 different RA95 variants were produced in E. coli BL21 Gold (DE3) and purified by nickel affinity chromatography as previously described. , , The enzymes were stored in buffer A (25 mM HEPES pH 7.5 and 150 mM NaCl), which was also used for all kinetic assays.

Synthesis of enantiopure (R)‐methodol

(R)‐Methodol was synthesized enzymatically from acetone and 6‐methoxy‐naphthaldehyde (6‐MNA) as reported previously, using RA95.5–8F as the biocatalyst. , The reaction mixture containing 2 mM 6‐MNA, 2 M acetone and 0.1 μM RA95.5–8F in 50 ml buffer A was gently shaken for 3 hr at 29°C and subsequently saturated with 15 g solid sodium chloride and extracted with 3 × 50 ml ethyl acetate. The organic phase was dried over sodium sulfate before removing the solvent under vacuum. The crude material was purified by flash chromatography. (R)‐Methodol was obtained in 58% yield. H‐NMR, HR‐MS and chiral HPLC were performed for quality control as described previously. ,

Steady‐state enzyme activity assays

Kinetic measurements were carried out using a Perkin Elmer UV/Vis spectrometer (Lambda35 with Peltier system and cell changer). (R)‐Methodol cleavage was measured at 29°C in buffer A supplemented with 2.7% acetonitrile as described previously. , , The formation of 6‐MNA was monitored over time by absorption at 350 nm (ε350 = 5,970 M−1 cm−1). The enzyme concentrations used in the assay ranged from 30 nM for the most active variant to 12 μM for the least active variant. The (R)‐methodol concentration in the assay buffer (25 to 900 μM) was determined by absorption at 330 nm (ε330 = 1,390 M−1 cm−1) prior to initiating the reaction by addition of enzyme. All measurements were performed in triplicate. The steady‐state parameters k cat and KM were determined by fitting the data to the Michaelis–Menten equation.

Protein stability measurements using a thermal shift assay

Melting curves were measured according to the standard protocol of the Protein Thermal Shift Kit (Thermo Fisher Scientific) in 96‐well PCR plates (V = 20 μL per sample) using a StepOnePlus Real‐Time PCR thermocycler (Thermo Fisher Scientific) and SYPRO Orange dye (Thermo Fisher). The melting curves were recorded in triplicate from 25–99°C. The protein concentration was 20 μM in buffer A. The fluorescence signal derivative was plotted against temperature to extract the midpoints of the melting transitions.

AUTHOR CONTRIBUTIONS

Timothy A. Coulther: Conceptualization; formal analysis; investigation; writing‐review & editing. Moritz Pott: Conceptualization; investigation. Cathleen Zeymer: Conceptualization; formal analysis; investigation; writing‐original draft; writing‐review & editing. Don Hilvert: Conceptualization; formal analysis; funding acquisition; resources; supervision; writing‐review & editing. Mary Jo Ondrechen: Conceptualization; formal analysis; funding acquisition; methodology; project administration; resources; software; supervision; visualization; writing‐review & editing.

CONFLICT OF INTEREST

The authors declares no conflict of interest.

Appendix S1: Supplementary Information

Click here for additional data file.